=Paper=
{{Paper
|id=Vol-2964/article_166
|storemode=property
|title=Data-based Discovery of Governing Equations
|pdfUrl=https://ceur-ws.org/Vol-2964/article_166.pdf
|volume=Vol-2964
|authors=Waad Subber,Piyush Pandita,Sayan Ghosh,Genghis Khan,Liping Wang,Roger Ghanem
|dblpUrl=https://dblp.org/rec/conf/aaaiss/SubberPGKWG21
}}
==Data-based Discovery of Governing Equations==
https://ceur-ws.org/Vol-2964/article_166.pdf
Data-based Discovery of Governing Equations
Waad Subber1∗ , Piyush Pandita1 , Sayan Ghosh1 , Genghis Khan1 , Liping Wang1 , Roger Ghanem2
1
GE Research, Niskayuna, NY 12309
2
University of Southern California, Los Angeles, CA 90007
Abstract scribes causal mechanisms by simplified mathematical for-
mulations, while the ML approach seeks to establish a sta-
Most common mechanistic models are traditionally presented
in mathematical forms to explain a given physical phe- tistical relationship between inputs and outputs. These two
nomenon. Machine learning algorithms, on the other hand, approaches should not be seen as direct competitors (Baker
provide a mechanism to map the input data to output with- et al. 2018). The advantages of one approach should be used
out explicitly describing the underlying physical process that to complement its counterpart, which suggests that mod-
generated the data. We propose a Data-based Physics Dis- ern research efforts in the scientific machine learning field
covery (DPD) framework for automatic discovery of gov- should be directed towards enabling a symbiotic relation-
erning equations from observed data. Without a prior defi- ship between both approaches (Baker et al. 2019; Jain and
nition of the model structure, first a free-form of the equa- Singh 2003; Baker et al. 2018). A synergy framework be-
tion is discovered, and then calibrated and validated against tween machine learning and mechanistic approach can be
the available data. In addition to the observed data, the DPD
built based on integrating multiple sources of information
framework can utilize available prior physical models, and
domain expert feedback. When prior models are available, (such as field and lab data), prior domain knowledge, phys-
the DPD framework can discover an additive or multiplica- ical constraints and expert feedback in one unified frame-
tive correction term represented symbolically. The correction work. The main feature of such an approach is manifested
term can be a function of the existing input variable to the in the symbolic representation of the predictive model. The
prior model, or a newly introduced variable. In case a prior symbolic description of the prediction mechanism by math-
model is not available, the DPD framework discovers a new ematical expressions can provide explainability for its pre-
data-based standalone model governing the observations. We dictions, facilitate integrating expert feedback, and fuse the
demonstrate the performance of the proposed framework on state-of-art domain knowledge in an explicit manner.
a real-world application in the aerospace industry.
To this end, we propose a Data-based Physics Discov-
ery (DPD) framework for automatic discovery of govern-
Introduction ing equations from observed data. Symbolic regression and
Modern machine learning (ML) methods are aimed at pro- Bayesian calibration are utilized in discovering the physi-
viding a statistical mechanism to predict the outcome of a cal laws governing the data. The approach is based on inte-
system under new conditions. This statistical mechanism is grating multiple sources of data, domain knowledge, phys-
constructed based on exploring the correlation between in- ical constraints and expert feedback in one unified frame-
puts and outputs that is embedded in data (Jain and Singh work for model discovery. In our previous work (Atkin-
2003). However, in many engineering applications, the in- son et al. 2019), we introduced a method to infer a sym-
puts, outputs, and ML model structure are not selected such bolic representations of differential operators from data, in
that learning elucidates insight into the underlying physical a free-form manner as opposed to methods such as SINDy
process that generated the data beyond a black-box function (Brunton, Proctor, and Kutz 2016) which require a user to
approximation. Thus, knowledge is discovered in a data- postulate a library of terms from which only linear com-
driven manner without fully explaining the physics of the binations might be considered. In this work, we propose a
problem. The mechanistic modeling approach, on the other framework to discover governing equation from data uti-
hand, proceeds from a starting point of scientific laws and lizing prior physical domain knowledge and constraints. In
axioms, producing by way of logical deduction formal mod- case a prior physical model is available, the DPD frame-
els of the physics underlying a phenomenon and measure- work discovers a multiplicative or additive correction term;
ments of it. Typically, mechanistic modeling approach de- otherwise a standalone model based on data only is pro-
∗
Corresponding Author: Waad Subber, email:
posed. The discovered correction term or standalone model
Waad.Subber@ge.com can be a function of existing input variable or a newly intro-
Copyright c 2021, for this paper by its authors. Use permitted un- duced variable. The symbolic representation of the discov-
der Creative Commons License Attribution 4.0 International (CC ered model helps in explaining the effect of a new variable
BY 4.0). on the current physical process. The discovered model is
calibrated using Bayesian Hybrid Modeling (GEBHM) ap- Field & lab data
discovering
proach (Ghosh et al. 2020; Zhang et al. 2020). The GEBHM calibration
Input testing
is a probabilistic ML method that enables calibration, valida- Observations,
tion, multi-fidelity modeling and uncertainty quantification. Prior physics models and Constraints
In addition, the framework equipped with an optimal experi- Physics discovery
BHM calibration
ment design tool to propose a new experiment for enhancing Data
Hypothesis
+
the model accuracy. Here, the Intelligent Design and Anal- 𝑢 𝑢 !
Fitness
Model
× 𝑓 𝑢, 𝑡 = 0
ysis of Computer Experiments (IDACE) (Kristensen et al. Prior
exp
2019) methodology is utilized. The IDACE is a technique 𝑡
that adaptively generates new experiment setting for improv- Discover a novel model
that better describes the
ing the accuracy of the model based on the estimated uncer- data and prior knowledge IDACE updating
tainty and output desirability.
We provide a technical description of the approach and
demonstrate its use on a real-world problem with relevance Output
Updated or new model
to the aerospace industry. The application is focused on dis-
covering a predictive model describing corrosion, given the
Figure 1: The workflow of the Data-based Physics Discov-
availability of a prior physical model. We demonstrate the
ery (DPD) framework. The inputs to the framework are ob-
ability of the proposed framework to fuse prior knowledge
servation, and optional prior model and constraints. First,
with a discovered correction term, improving the model’s fi-
depending on the setting, a standalone model or a correction
delity to measurements.
term is discovered in the physical discovery part. Second,
the discovered model is calibrated. Third, V&V and UQ are
Data-based Physics Discovery (DPD) performed to the discovered model and, if needed, a new set
framework of experiments are proposed to improve the accuracy. The
The proposed Data-based Physics Discovery (DPD) ap- output is a calibrated and validated model.
proach for integrating domain knowledge, physical con-
straints and observational data in one framework is devel-
oped to predict the response (with confidence bounds) of Here θ = {θp , θd , ρ}, where θp denote the prior model
complex engineering systems under novel conditions. This parameters, θd represent the ephemeral random constants
framework enables a relationship between data-based ma- generated during the evolution of the symbolic regres-
chine learning and causality-based mechanistic modeling sion (Fortin et al. 2012), and ρ is a weighting parameter.
approaches. The DPD framework serves as a human’s re- This setting can be readily generalized to the case of multi-
search assistant with two goals in mind: enhance the cred- plicative correction as:
ibility of prediction and enable learning new physics. This g(x, θ) = f (x, θp ) × ρδ(x, θd ), (4)
constitutes a new model-building paradigm whereby ma-
chine learning algorithms and scientists work together to and to the case where a prior model is not available:
formulate a symbolic mathematical form explaining a novel g(x, θ) = δ(x, θd ). (5)
physical phenomenon. Similar to the physical experiments
Note that we do not a priori define a model structure of the
performed by scientists to discover the interaction relation-
correction term δ(x, θd ). In this step, the scaling parame-
ship between the input and output variables, we propose a set
ter ρ is set to unity and the prior model parameters θp are
of variables to our symbolic regression framework and ask
set to their maximum likelihood estimation values, and θd is
for a relationship govern them. An overview for the work-
created and set during the evolutionary computation at run
flow of the proposed approach is shown in Fig. (1).
time (Koza and Koza 1992), (ii) After discovering the model
Given set of observations D = {xi , yi }ni=1 , (and if
structure, the model parameters θ = {θp , θd ρ} in the newly
available, a prior physical model f (x, θp ) and a constraint
discovered model g(x, θ) are calibrated next using Bayesian
φ(x) = 0), we seek a new or updated model that bet-
approach, (iii) In case the calibrated model does not sat-
ter describes the observed data. The process as shown in
isfy the V&V requirements (e.g., test data are within 95%
Fig. (1) can be summarized as follows: (i) Let a prior model
confidence interval of the model prediction), a new optimal
f (x, θp ), where θp are the calibration parameters, be avail-
experiment for enhancing the model accuracy is proposed
able and an additive correction term is required, then the
using an intelligent design of experiments. Once the discov-
problem is posed as: Find δ(x, θd ) such that
ered model satisfies the V&V and UQ requirements it will
J = kg(x, θ) − yk2 → min, (1) be proposed as an explainable model to the domain experts
for their feedback.
subject to: In essence, the proposed workflow embeds domain
φ(x) = 0, (2) knowledge and physical laws into machine learning algo-
rithms to provide an interpretable predictive model. Salient
where features of the DPD framework are: 1) enhancing the fore-
casting ability of the machine learning algorithms by embed-
g(x, θ) = f (x, θp ) + ρδ(x, θd ). (3) ding physics principles enabling extrapolation for regions
were data is not available, 2) addressing the limitations in the 1.0 train
traditional mechanistic modeling approach in incorporating test
multiple sources of information, 3) minimizing the size of 0.8
corrosion current
the training data required for machine learning approach by 0.6
incorporating domain knowledge.
0.4
Real-World Application 0.2
An aerospace industrial application is considered to demon-
strate the potential practicality of the proposed framework. 0.0
The problem is focused on developing a predictive model 0.0 0.2 0.4 0.6 0.8 1.0
time index
for the corrosion process of an aircraft structural material.
In this application, we demonstrate: 1) the ability of DPD
in discovering a correction term to an existing physical Figure 2: The data is halved equally into training and test-
model for better forecasting ability, 2) the improved perfor- ing sets. The decomposition is based on the time index (i.e.,
mance obtained by the DPD compared to an existing phys- training data t ≤ 0.5 and testing data t > 0.5).
ical model for the case of small training data set, 3) the
performance of the DPD in discovering a standalone model
is used as a prior model for the temperature T . The correc-
based on only the training data set without incorporating a
tion term is set as a function of the relative humidity H and
prior physical model.
is denoted by δ(H, θd ). The symbolic regression in the DPD
Corrosion problem framework will discover δ(H, θd ) such that:
Extending the service lifetime of a corrosion protection coat- J = kg(T, H, θ) − yk2 → min, (8)
ing system necessitates a fundamental understanding of the and
materials under harsh operation conditions; the sensitivity
of the corrosion process to material features and service g(T, H, θ) = f (T, θ̂p ) × ρ̂δ(H, θd ), (9)
conditions constitutes a challenge to the development of a
where y is the measured corrosion current. Here θ̂p are set
novel corrosion-resistant materials (Olajire 2018). The pre-
a priori to the maximum likelihood estimation values, and
diction of the corrosion under long-term operation condi-
ρ̂ = 1. Once the most suitable structure of the correction
tions and unforeseen environmental events is a crucial pro-
term is discovered, a Bayesian calibration is performed to
cess for maintenance scheduling of the aircraft structural. To
estimate the model parameters θ = {θp , θd , ρ}. The equa-
this end, our proposed framework will provide a calibrated,
tion discovery and the model calibration are performed on
validated and uncertainty quantified predictive model for the
the first 50% portion of the data, whereas the prediction and
corrosion as a function of the environmental and operations
testing of the model are performed on the forthcoming 50%
conditions. The symbolic representation of the discovered
of the data as shown Fig. (2).
corrosion model enables interpretability of the prediction
model, which is required for understanding the corrosion Baseline model: First, using the training data set t ≤ 0.5,
process under novel conditions. a gradient based optimization is performed to obtain the
maximum likelihood estimation (MLE) of the parameters in
Problem setting: The experimental data is provided by
the Butler-Volmer equation. The model, which is a function
the AFRL as part of the DARPA AIRA challenge prob-
of temperature only, is validated on the remaining test data.
lems (AIRA 2018). The measurements are represented by
The validation plot for the test data is shown in Fig. (3). The
time-series data of the corrosion current, temperature and
performance metrics R-squared R2 = 0.691 and root mean
relative humidity are reported every 5 minutes. The objec-
square error RM SE = 0.109 are calculated using the test
tive here is to develop a forecasting model for the corrosion
data.
current, given the environmental conditions.
In order to improve the performance of the Butler-Volmer
For this application, we pose the problem as finding a cor-
equation without requiring extra training data, a new correc-
rection term to an available prior physical model as:
tion term to the existing model is introduced. The goal of
g(T, H, θ) = f (T, θp ) × ρδ(H, θd ), (6) the correction term is not only to enhance the model perfor-
mance, but also to study the effect of newly introduced vari-
where T is the temperature, H is relative humidity, and the able on the corrosion current. Commonly, the former goal is
θ = {θp , θd , ρ} are the calibration parameters. Note that the achieved by increasing the size of the training data set, and
available prior model is a function of temperature T only, the latter goal is carried out by a domain expert in mecha-
while the discovered term is set to be as a function of humid- nistic way. Next, we will show how DPD framework can im-
ity H only. This setting is to show how the DPD framework prove the performance of an existing model by an automated
can introduce a new variable to an existing physical model in discovery of a new governing equation using the available
order to improve the forecasting accuracy. Specifically, the training data set.
Butler-Volmer equation of the form:
DPD with a prior model: Utilizing the existing model
f (T, θp ) = θ0 exp(θ1 /T ) + θ2 exp(θ3 /T ) (7) f (T, θ̂p ) with the a priori estimated parameters θ̂p and the
1.0 1.0
0.8 0.8
0.6
predicted
predicted
0.6
0.4 0.4
0.2 0.2
R2(test) = 0.887
R2(test) = 0.691 RMSE(test) = 0.066
RMSE(test) = 0.109 0.0
0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
observed observed
Figure 3: Validation of the baseline model f (T, θp ) on the Figure 4: Validation of the discovered model g(T, H, θ) on
testing data set, t > 0.5. the testing data set, t > 0.5.
Table 1: Discovered correction terms to the Butler-Volmer 1.0
equation for the training data set t ≤ 0.5. Performance met- observed
rics are based on the testing data set t > 0.5. 0.8 forecasted
corrosion current
0.6
BIC R2 δ(H, θd ) V[%]
-2078.2 0.787 0.481H exp(H 3 ) 14 0.4
-2073.9 0.774 H(2H − 0.73) 12
0.2
–2064.3 0.763 2H − 0.821 10
0.0
0.5 0.6 0.7 0.8 0.9 1.0
time index
training data set, Table (1) lists some of the discovered mul-
tiplicative correction terms. In addition to the R2 listed for Figure 5: The forecasting of the discovered model
each term, the Bayesian information criterion (BIC) and g(T, H, θ) for t > 0.5. The confidence bounds are defined
the gained improvement are presented. The percentage im- as (µ ± σ).
provement is defined based on the R2 value as V [%] =
(Ri2 /Rb2 − 1) × 100, where Rb2 = 0.691 is value of the
baseline model and Ri2 is the value of the new discov- size (t ≤ 0.08), the performance metrics of the baseline
ered model, before Bayesian calibration. The parameters model are much lower than that obtained when using half
θd of the correction term δ(H, θd ) are not calibrated with of the data (t ≤ 0.5) as shown in Fig. (3). Precisely a reduc-
Bayesian approach yet. A forecasting enhancement ranging tion of (89%) in the R2 value is observed. This observation
between 10%-14%, depending on the complexity of the ob- highlights the effect of the training data size on the perfor-
tained term, can be achieved by introducing a new variable mance of the model. Next, using the available training set
to the existing model. (t ≤ 0.08) we show an improvement can be obtained us-
ing the DPD framework. The improvement is achieved by
Model calibration: In order to capture any missing
introducing a correction term as a function of a new input
physics, the Kennedy and O’Hagan Bayesian calibra-
variable.
tion (Kennedy and O’Hagan 2001) is carried out next using
For the given training data set (t ≤ 0.08), Table (2) lists
the GEBHM framework. The validation results of the cali-
the discovered correction terms δ(H, θd ) with the corre-
brated model is reported in Fig. (4). An improvement of 28%
sponding BIC, R2 and the percentage value of improvement
can be achieved by the new discovered model with parame-
obtained, before calibrating the model. By introducing a new
ters being calibrated using GEBHM. Fig. (5) shows the fore-
variable to the existing model, the DPD discovered a correc-
casting of the discovered model as well as the observed data
tion term (e.g., δ(H) = 0.391 exp(H)) that can lead to an
for t > 0.5. The gray shaded area is the confidence bounds
improvement of 78%. The percentage improvement V [%] in
defined as µ ± σ. All data are normalized between 0 and 1.
Table (2) is calculated based on the baseline model shown in
The model prediction agrees quite well with the measured
Fig. (7).
data.
Next, the new discovered model g(T, H, θ) is calibrated
Reduction of the training data size: To study the ef- using using GEBHM. The forecasting and validation results
fect of the training data size on the model performance, are shown in Fig. (9) and Fig. (8), respectively. The predic-
we use only the first 0.08% of the available data for train- tion of the current discovered model follows the trend of
ing and the remaining portion are set for testing as shown the observed data; however, shows an overestimation to the
in Fig. (6). Fig. (7) shows the validation plot of the cur- corrosion current between the peaks. Nevertheless, an im-
rent baseline model (the Butler-Volmer with maximum like- provement of 239% as shown in Fig. (8) can be achieved
lihood estimation of the parameters). For the training data compared to the baseline model shown in Fig. (7). The large
1.0 Table 2: Discovered correction terms to the Butler-Volmer
train equation for the training data set t ≤ 0.08. Performance
0.8 test metrics are based on the testing data set t > 0.08..
corrosion current
0.6
BIC R2 δ(H) V[%]
0.4 -5620.2 0.132 0.391 exp(H) 78
0.2 -5609.6 0.098 0.99H + 0.083 32
-5608.8 0.095 0.274H 2 + 0.509H + 0.293 28
0.0
0.0 0.2 0.4 0.6 0.8 1.0
time index
1.0
Figure 6: The data is partitioned based on the time index into 0.8
training (t ≤ 0.08) and testing (t > 0.08) sets.
0.6
predicted
1.0 0.4
0.8 0.2
R2(test) = 0.251
0.0 RMSE(test) = 0.172
predicted
0.6
0.0 0.2 0.4 0.6 0.8 1.0
0.4 observed
0.2
R2(test) = 0.074 Figure 8: Validation of the discovered model g(T, H, θ) on
0.0 RMSE(test) = 0.192 the testing data set, t > 0.08.
0.0 0.2 0.4 0.6 0.8 1.0
observed
standalone model based on only the observed data, in case a
Figure 7: Validation of the baseline model f (T, θp ) on the prior model is not available. A Bayesian approach is utilized
testing data set, t > 0.08. for the model calibration and validation under uncertainty.
A real-world application in the aerospace industry is consid-
ered to show the practicality of the proposed framework.
improvement obtained here indicates that although the dis-
covered model has improved the performance by 78% (as
shown previously), it is still missing some physics which Acknowledgments
can be captured by GEBHM based calibration. The authors thank the United States Air Force Research
DPD without a prior model: Utilizing the training data Laboratory for providing the corrosion data for public re-
set t ≤ 0.5 only, Table (3) shows some of the discov- lease under Distribution A as defined by the United States
ered standalone models without any prior physical model. Department of Defense Instruction 5230.24. This material
Here the discovered model takes the form g(T, H, θ) = is based upon work supported by the Defense Advanced
δ(T, H, θd ). Large improvement can be achieved with more Research Projects Agency (DARPA) under Agreement No.
complicated model structure. These models explain the data HR00111990032. Approved for public release; distribution
in mathematical expressions, and they are ready for inter- is unlimited.
pretation by the domain experts. Bayesian calibration of the
discovered model References
g(T, H, θ) = H 2 (θ0 H 2 − θ1 T ), (10) AIRA. 2018. Artificial Intelligence Research Asso-
ciate AIRA. https://www.darpa.mil/program/artificial-
gives the validation results shown in Fig. (10), with an im- intelligence-research-associate.
provement of 28%. The forecasting performance of the dis-
covered standalone model is shown in Fig. (11). Atkinson, S.; Subber, W.; Wang, L.; Khan, G.; Hawi, P.;
and Ghanem, R. 2019. Data-driven discovery of free-
Conclusions form governing differential equations. arXiv preprint
arXiv:1910.05117 .
A Data-based Physics Discovery (DPD) framework for an
automatic discovery of governing equations from observed Baker, N.; Alexander, F.; Bremer, T.; Hagberg, A.;
data is presented. The framework can discover a correction Kevrekidis, Y.; Najm, H.; Parashar, M.; Patra, A.; Sethian,
term to an existing physical model in order to improve the J.; Wild, S.; et al. 2019. Workshop report on basic research
prediction performance. The correction term can be a multi- needs for scientific machine learning: Core technologies for
plicative or additive, and a function of the current input vari- artificial intelligence. Technical report, USDOE Office of
able or new variable. In addition, the DPD can discover a Science (SC), Washington, DC (United States).
1.0 observed 1.00
0.8 forecasted
0.75
corrosion current
predicted
0.6 0.50
0.4 0.25
R2(test) = 0.884
0.2 0.00 RMSE(test) = 0.067
0.0 0.0 0.2 0.4 0.6 0.8 1.0
0.2 0.4 0.6 0.8 1.0
time index observed
Figure 9: The forecasting of the discovered model Figure 10: Validation of the discovered standalone model
g(T, H, θ) for t > 0.08. The confidence bounds are defined g(T, H, θ) on the testing data set, t > 0.5.
as (µ ± σ).
1.0 observed
Table 3: Discovered standalone models for the training data forecasted
set t ≤ 0.5. Performance metrics are based on the testing 0.8
corrosion current
data set t > 0.5.
0.6
BIC R2 δ(T, H) V[%] 0.4
0.537HT (H 2 T + H − T + 0.2
-2235.6 0.862 25
0.457(H − T ) exp(H(H + T )))
0.0
-2220.5 0.833 H 2 (H 2 − 0.456T ) 21 0.5 0.6 0.7 0.8 0.9 1.0
H 3 (−0.479T + time index
-2186.2 0.825 19
0.479 exp(H) − 0.304) Figure 11: The forecasting of the discovered standalone
-2195.5 0.802 0.527H 4 exp(H − T ) 16 model g(T, H, θ) for t > 0.5. The confidence bounds are
defined as (µ ± σ).
Baker, R. E.; Peña, J.-M.; Jayamohan, J.; and Jérusalem, A.
2018. Mechanistic models versus machine learning, a fight Kristensen, J.; Subber, W.; Zhang, Y.; Ghosh, S.; Kumar,
worth fighting for the biological community? Biology letters N. C.; Khan, G.; and Wang, L. 2019. Industrial applications
14(5): 20170660. of intelligent adaptive sampling methods for multi-objective
optimization. In Design Engineering and Manufacturing.
Brunton, S. L.; Proctor, J. L.; and Kutz, J. N. 2016. Discov- IntechOpen.
ering governing equations from data by sparse identification
of nonlinear dynamical systems. Proceedings of the national Olajire, A. A. 2018. Recent advances on organic coating sys-
academy of sciences 113(15): 3932–3937. tem technologies for corrosion protection of offshore metal-
lic structures. Journal of Molecular Liquids 269: 572–606.
Fortin, F.-A.; De Rainville, F.-M.; Gardner, M.-A.; Parizeau,
Zhang, Y.; Ghosh, S.; Pandita, P.; Subber, W.; Khan, G.; and
M.; and Gagné, C. 2012. DEAP: Evolutionary Algorithms
Wang, L. 2020. Remarks for scaling up a general gaussian
Made Easy. Journal of Machine Learning Research 13:
process to model large dataset with sub-models. In AIAA
2171–2175.
Scitech 2020 Forum, 0678.
Ghosh, S.; Pandita, P.; Atkinson, S.; Subber, W.; Zhang, Y.;
Kumar, N. C.; Chakrabarti, S.; and Wang, L. 2020. Ad-
vances in Bayesian Probabilistic Modeling for Industrial
Applications. ASCE-ASME J Risk and Uncert in Engrg Sys
Part B Mech Engrg 6(3).
Jain, S. K.; and Singh, V. P. 2003. Water resources systems
planning and management. Elsevier.
Kennedy, M. C.; and O’Hagan, A. 2001. Bayesian calibra-
tion of computer models. Journal of the Royal Statistical
Society: Series B (Statistical Methodology) 63(3): 425–464.
Koza, J. R.; and Koza, J. R. 1992. Genetic programming:
on the programming of computers by means of natural se-
lection, volume 1. MIT press.