Combining Grammatical Evolution with Modal Interval Analysis: An Application to Solve Problems with Uncertainty

: Complex systems are usually affected by various sources of uncertainty, and it is essential to account for mechanisms that ensure the proper management of such disturbances. This paper introduces a novel approach to solve symbolic regression problems, which combines the potential of Grammatical Evolution to obtain solutions by describing the search space with context-free grammars, and the ability of Modal Interval Analysis (MIA) to handle quantiﬁed uncertainty. The presented methodology uses an MIA solver to evaluate the ﬁtness function, which represents a novel method to manage uncertainty by means of interval-based prediction models. This paper ﬁrst introduces the theory that establishes the basis of the proposed methodology, and follows with a description of the system architecture and implementation details. Then, we present an illustrative application example which consists of determining the outer and inner approximations of the mean velocity of the water current of a river stretch. Finally, the interpretation of the obtained results and the limitations of the proposed methodology are discussed.


Introduction
Continuous improvements in algorithmic problem-solving, together with an increase in the availability of high-performance computing, have resulted in a new generation of precise and highly detailed mathematical models [1]. The modeling of the physical properties within these systems is complex, because they are generally unknown or indeterminate. Explicitly incorporating uncertainty in the predictions of a model from the beginning is vital. The identification and management of uncertainty is essential to manage the complexity that occurs in numerous engineering applications [2]. Methods that model these uncertainties seek to develop robust mechanisms designed to remain flexible and resilient to appropriately react to new situations.
In general, the uncertainty of a system can arise from many sources, such as disturbances from the physical environment (e.g., variability in the flow controlling a water level), noise from devices that collect data (e.g., blood glucose readings from a continuous blood glucose monitor [3]), the representation of errors in physical quantities (e.g., the error in the gravitational acceleration constant g = 9.807 ± 0.027 m/s 2 ), truncation errors in floating point operations, etc. These uncertainties will lead to undesired system behaviors if not addressed properly. Therefore, methodologies able to manage such uncertainty in a robust manner are required.
The purpose of this paper is to present a computational framework for solving symbolic regression problems, which is able to cope with the uncertainty encountered in real-life problems, and to make predictions accounting for such uncertainty. The proposed approach is based on interval analysis and data-driven optimization techniques. In particular, it combines the advantages of a grammar-based search algorithm [4] with the arithmetic of Modal Interval Analysis (MIA) [5]. The excellent flexibility of Grammatical Evolution (GE) for generating models, together with the inherent ability of MIA to manage uncertainty, led us to propose the development of a combined method to deal with the generation of predictive models from imprecise, or uncertain, data, without losing interpretability in the results. This study contributes to the literature by providing an evolutionary algorithm capable of working directly with MIA, which provides a mathematical framework to deal with quantified uncertainty. While the GE algorithm delimits the space of the solutions and guides the search, MIA manages uncertainty naturally and provides semantics to transform logical formulas into inclusion relationships, which are used to build the algorithm loss function. In order to illustrate the utilisation of the proposed approach, the resolution of a real-life problem, consisting of determining the meanstream velocity of the waters of a river stretch, is presented. This paper is organized as follows. Uncertainty handling in predictive modeling is briefly reviewed in the next Subsection, while Sections 2 and 3 summarize the foundations of GE and MIA theory, respectively. Section 4 introduces our proposed interval-based grammatical evolution approach. Section 5 presents an illustrative example, the required implementation steps, and the results. Finally, in Section 6, the obtained results are discussed and some conclusions provided.

Related Work
Two approaches are typically used to handle uncertainty in predictive modeling: explicit, or knowledge-driven, models and data-driven models. Techniques such as Bayesian networks, fuzzy logic, interval analysis, and rough sets have been successfully integrated into explicit, or knowledge-driven, models to address systems with uncertainty [6][7][8][9]. The advantage of these models is that they provide an explicit representation of the internal operation of the system to predict, simulate, and explain its behavior from the structure, causality, functional, and behavior of its components. However, these techniques are often computationally prohibitive or inaccurate for modelling complex systems that depend on a large number of variables. One of the more important techniques in this group is the well-know classical interval analysis. Classical intervals have shown a wide range of abilities, such as addressing problems using optimization to avoid overestimation of the optimal range of values [10][11][12]. In addition, the study of the solutions of a linear system when the coefficients and the independent term are considered as intervals has received much attention in the interval community, since it appears in the control design of different physical systems [13][14][15][16][17]. Pursuing the same line, in 2014 Modal Interval Analysis (MIA) [5], was presented as an extension of classical interval analysis. MIA allows uncertainty to be introduced, as classical intervals, in the generation of mathematical models by representing it using an interval (i.e., a pair of real numbers bounding a real domain). In addition to a range, model intervals are defined by a logical quantifier (∀ or ∃), that affects the represented uncertainty. Hence, MIA allows for a richer semantic interpretation of the results.
On the other hand, data-driven methods such as deep neural network techniques have focused on the extraction of predictive models from uncertain data [18]. Such models are especially effective when knowledge-driven models are difficult to build (e.g., owing to insufficient understanding of the underlying processes). Data-driven models cannot catch the physics of the modeled process; they just capture the relationships between the relevant input and output variables. However, such data-driven methods could potentially be more accurate than knowledge-driven models because they are based on objective information (i.e., data). The main drawback of all these black-box approaches is that they lead to a loss of result interpretability, and causality cannot be inferred. Machinelearning methods have been prompted by a recent wave of interest due to improvements in the algorithms (e.g., activation function) that have made such techniques more efficient. Finally, the combinations of more traditional approaches, such as control engineering theory, with artificial intelligence, have resulted in a number of hybrid approaches [19,20].
This paper introduces an approach to solve symbolic regression problems, combining the potential of MIA to handle uncertainties with the ability of GE to obtain solutions. MIA is a mature mathematical technique that has been widely applied in control engineering [21], modeling and simulation [22], civil engineering for structural control [23], biomedical applications [24], computer graphics [25], and finance for insurance problems [26]. On the other hand, GE is a relatively new machine-learning methodology with a modular architecture that is applied in the design and optimization of predictive models [27]. GE shows good performance in complex environments, such as financial forecasting [28] or blood glucose predictions [4].
To the best of our knowledge, evolutionary algorithms have not been previously used in combination with MIA and we have not found any work using interval methods combined with GE. However, different studies investigating hybrid methods applying interval analysis techniques to deal with uncertainties and evolutionary algorithms to guide the search of a solution have been proposed. For instance, hybrid approaches based on multi-objective evolutionary algorithms, either based on the conversion of an interval multi-objective evolutionary algorithm to a deterministic single-or multi-objective optimization problem [29][30][31], or based on the interval dominance relation [32,33]. Furthermore, hybrid approaches on multiple types of evolutionary algorithms have also been presented. For example, Sun et al. [34] implemented a particle swarm optimization method to power lithium-ion batteries, Femia et al. [35] applied a combination of affine arithmetic (an extension of interval arithmetic) and genetic algorithms for worst-case analysis in circuit tolerance analysis, while Claudio et al. presented an hybrid approach applying interval analysis in a cellular evolutionary strategy [36]. Finally, one of the closer studies we have found is the work performed by Keijer, which proposed to enhance genetic programming using interval arithmetic [37]. The author presented an approach to produce solutions which avoided undefined behaviours and which performed better than standard approaches.

Predictive Modeling by Grammatical Evolution
One of the most compelling forms of evolutionary algorithms can be found in the GE methodology. GE is a search-based heuristic method with a modular design. It is widely used to solve optimization problems and is well-suited to the design of predictive models [27], which are, essentially, multivariate optimization problems. The goal of this search algorithm is to exploit historical information to guide it toward the regions of better performance within the search space. The general process of functioning in GE algorithms remains similar to those applied in other evolutionary algorithm approaches. Essentially, the algorithm initializes a population of solutions and initiates an iterative process that applies selection and transformation operators until the algorithm reaches an optimal solution or a predefined number of iterations.
The candidate solutions are presented as a sequence of codons. Typically, a codon comprises 8 bits representing an integer. Codons determine which specific grammar rules will be used.
The array of codons is called a genotype, while the code derived from the values is the phenotype. GE uses context-free grammar to enable a many-to-one mapping process, from genotypes to phenotypes. The genetic algorithm handles the generation of populations, whereas, at the phenotype level, solutions are evaluated in terms of a fitness function. Therefore, the genotype-phenotype mapping process involves the decoding of variablelength integer array (chromosomes) into expressions or programs in an arbitrary language. This feature allows the phenotype to be as complex as necessary, since all genetic operators are applied to the genotype. The general operation of the algorithm and the mapping process are outlined in Figure 1.

Context-Free Grammar
The key feature of predictive modeling via GE is the use of a context-free grammar (shortened to "grammar" hereinafter). This structure should incorporate knowledge related to the addressed problem. It comprises a set of rules that defines the structure of the generated candidate models. A grammar is defined by a 4-tuple {N, T, P, S}, where N and T are the set of non-terminal symbols and terminal symbols, respectively, P represents the production rules by which the symbols in N can be derived, and S the starting non-terminal symbol that must appear in N. The core of a grammar is typically defined in the Backus normal form (BNF) [Symbol] → production 1 | ... | production k Thus, given a grammar G = {N, T, P, S}, non-terminals N must be defined by a derivation rule, whereas terminals T are the lexicons of the predictive model language. The derivation rules in P comprise a non-terminal on the left-hand side [Symbol] and its possible derivations on the right-hand side. Each rule definition comprises one or more alternatives separated by the symbol "|". A single production is composed of a sequence σ = (σ 1 , · · · , σ n ) where σ i ∈ N ∪ T and #σ ≥ 0. Therefore, if we consider a production rule {[Symbol 1 ] → γ Symbol n β} ∈ P where [Symbol n ] ∈ N and β, γ ∈ N ∪ T, we define that [Symbol n ] is a direct derivation of [Symbol 1 ], denoted as [Symbol 1 ] ⇒ [Symbol n ], if a sequence of direct derivations exists where n ≥ 0.
A complete derivation tree is constructed by expanding all the non-terminals of the candidate solutions. Particularly, the process derives the leftmost non-terminal, as it is selected according the corresponding codon of the genotype. We repeat the process of derivation until the candidate solutions have no symbols remaining in N. During genotype mapping, we may reach the end of the genotype and be depleted of values before all non-terminals are transformed into terminals. At this point, the algorithm assesses the candidate model as an invalid solution. Therefore, candidate models are penalized if their genotypes lead to large phenotypes. Invalid solutions are assigned with poor fitness values, to discard in subsequent iterations.

In interval analysis, a classical interval [a, b] is defined as the set of numerical values
x that lies between a and b, and operations are performed using intervals instead of real numbers. [ In the system of intervals the real arithmetic operators are introduced. An interval extension of a continuous function Because the interval united extension of a continuous function can not always be computed, an interval syntactic extension f R(X 1 , . . . , X n ) is defined as its corresponding real function f (x 1 , . . . , x n ) replacing their numerical arguments x 1 , . . . , x n by the interval arguments X 1 , . . . , X n and their "real" arithmetic operators by their corresponding interval operators, such that where f R(X 1 , . . . , X n ), is computed from the bounds of the intervals X 1 , . . . , X n . The interval syntactic extension is useful for computing the range of a function because it guarantees the result. However, the exact range is not obtained in the general case. This is because, in the case that a variable appears more than once in the expression of the function, each appearance is considered independent of the others, which leads to an overestimation of the real range.

Modal Interval Analysis
Modal interval analysis is a logical and algebraic completion of classical interval analysis, defined to overcome its semantic and algebraic deficiencies. The key concept is the inclusion of a logical quantifier in the definition of a modal interval. This provides semantic interpretation to the interval computations. For a more complete introduction to modal interval analysis, see [5,38].
A modal interval (see Figure 2), A, is a pair formed by a classic interval A = [a, b] , its domain, and a logic quantifier ∀ or ∃, which is its modality, that is, A = (A , Q A ), where Q A is one of the logic quantifiers. Modal intervals of the type A = (A , ∃) are called proper intervals and modal intervals of the type A = (A , ∀) are called improper intervals. A set of modal intervals is denoted by I * (R). A modal interval can be represented using its canonical coordinates in the form For example, the interval [1,3] is equal to ([1, 3] , ∃) and the interval [3,1] is equal to ([1, 3] , ∀).
A convenient symmetry between proper and improper intervals is established by the duality operator, which changes the modality of a modal interval. Therefore, The construction of modal intervals is completed with the concept of a modal quantifier Q defined by that allows us to define the set of real predicates accepted as a modal interval X = (X , Q X ). The identification of a modal interval with the set of those real predicates that introduce the definition of the modal interval inclusion (see Figure 3). If X, Y ∈ Pred(R), Using their canonical coordinates, the inclusion between two modal intervals A = [a 1 , The lattice operations meet and join (see Figure 4) on I * (R) for a bounded family of modal intervals A(I) := {A(i) = [a 1 (i), a 2 (i)] ∈ I * (R) | i ∈ I} (I is the index's domain) are defined as a function of the interval bounds as The inequality relations (see Figure 5) lead to the lattice operators min and max.

Consider a bounded family of modal intervals satisfying
Hence, min and max are defined by the following equations and computationally: Important binary operators involving two operators, A and B, are the arithmetical operators +, −, * , /, together with other operators, namely log, exp, sin, · · · , and some metric operators, e.g., the width of an interval A = [a 1 , a 2 ], defined by wid(A) = a 2 − a 1 , or distance, defined as which allows for the inclusion to be characterized by the 0 value.

Modal Interval Extensions
Two modal interval extensions exist for a continuous function f from R n to R in the variables x = (x 1 , · · · , x n ) to the n-dimensional interval X = (X 1 , · · · , X n ), and a split in its proper and improper interval components X = (X p , X i ), which is called the *-semantic extension of f , and the following, which is its **-semantic extension Both are related by the equality

Primary Theorems
The following two key results, named the * and ** semantic theorems, provide a logical interpretation to these semantic extensions.
*-semantic theorem: Given a continuous real function f : R n → R and a modal vector and **-semantic theorem: Given a continuous real function f : R n → R and a modal vector For a given n-dimensional interval A, computing f * (A) (or f * * (A)) can be a difficult problem, except for simple cases. When the function f is a one-or two-variable operator, the computations are easy to obtain and, as for the arithmetic operators, the results are the same as the ones obtained via Kaucher arithmetic [39].
The values of semantic extensions, f * and f * * , can not be obtained by direct computation, except for simple real functions. If f is a R n to R continuous with the syntactic tree function, its syntactic extension to the modal intervals X 1 , . . . , X n , is represented by f R(X 1 , . . . , X n ), is the function f R from I * (R n ) to I * (R) defined by the computational program indicated by the syntax tree of f . Modal syntactic functions are easy to compute but are generally not interpretable. To overcome this problem, semantic theorems provide inclusion relations to the corresponding *-and **-semantic extensions. Computations with f R(X) must be performed with the external truncation of each operator to obtain the inclusions f * (X) ⊆ f R(X), and with the inner truncation to obtain the inclusions f R(X) ⊆ f * * (X). In many cases, the rational extension f R(X) is optimal, i.e., Modal interval analysis provides results for these inclusions or equalities. Coercion theorems provide the conditions and the method to obtain the optimal extensions according to the monotony.

f * Algorithm
When these theorems and their results are not applicable, the so-called f * algorithm, based on branch-and-bound techniques, yields the inner and outer approximations to f * . Consider a continuous function f (u, v), and its associated proper and improper interval vectors. (U, V).
To obtain the inner (Inn) and the outer (Out) approximations of f * (U, V) or f * * (U, V), we divide the initial domain into cells such that, in each of them, the monotony conditions necessary to apply the optimality theorems of modal interval analysis are met. The f * algorithm applies a set of strategies to minimize the number of bisections, and to obtain better local approximations of the resulting partitions. A tolerance function, defined as the distance between inner and outer approximation, provides the criteria to stop the algorithm where Inf and Sup are the left and right bounds of the corresponding approximations, respectively.
If we assume a desired tolerance ϕ for the output, the algorithm stops when Tolerance (Inn, Out) ≤ ϕ. The algorithm also stops when the width of all cell dimensions is smaller than a fixed precision (wid(X) ≤ ). A standalone version of this algorithm can be found online, incorporated in the modal interval calculator package [40,41].

Interval-Based Grammatical Evolution
Evolutionary algorithms are are well-known because their high scalability and flexibility. They have been largely studied, and countless approaches can be found in the literature. There are other algorithms in this category, such as genetic algorithms and genetic programming, to name just two of the most popular members. The idea of an MIA-based extension of an evolutionary algorithm should be applicable in a more traditional setting than GE, such as in the aforementioned methodologies. GE has been chosen because of its versatility and adaptability to extract models from data and our experience working with it. Figure 6 shows a schematic representation of the general architecture and the operation of the proposed methodology. Initially, the system is powered by a database (A) that provides the required information for system training and the subsequent validation of the method. In addition to the above-mentioned data sources, there are two other primary system inputs that must be defined according to the addressed problem. First, the definition of a problem-specific objective function (B), which evaluates the solutions. Secondly, a customized grammar (C), which defines the structure of the generated solutions. The module implementing the arithmetic of modal intervals is responsible for the interval operations associated with the solutions generated by the system (E). Solutions are iteratively combined to create new and better solutions, aiming to incrementally improve the quality of solutions (F) and reach a final solution that minimizes the fitness function satisfactorily. Once a final solution (G) is generated, the prediction model is evaluated using data from the remaining database information. The following subsections describe the complete setup of this methodology in detail. Figure 6. General diagram of the general methodology, where A is the database; B, C and D are the inputs of the grammatical evolution methodology; E is the core of the method combining the interval arithmetic and the evolutionary operators; F is the full set of candidate models and G is the selected final model.

Including Uncertainty in Grammatical Evolution
Given a set of the experimental measures y, we shall consider a predictive model defined by a continuous function fŷ whereŷ is the state variable corresponding to y, y 0 the initial values, and p are the model parameters. If we assume uncertainty in the variables and the experimental measures of this model, we obtain the interval model defined bŷ where P, Y 0 and, consequently, the model output (Ŷ) are proper or improper intervals. The generation of models that satisfy the aforementioned equation would require the implementation of mechanisms to achieve the following • Manage intervals as inputs of the system; • Generate intervals as parameters; • Compute intervals as basic operations.
Therefore, the integration of uncertainty management into the method of GE directly affects its two essential features: the use of context-free grammar and the fitness function.

Uncertainty in Context-Free Grammars
First, before beginning the search-based process of GE (see Figure 1), the method incorporates a series of preprocessing mechanisms to enable intervals as a native data type in the computational framework. Interval inputs I associated with the explored methodology can be defined using the lower and upper bounds of the interval following Equation (3) or by defining an increment in the discrete values as follows Therefore, the method encompasses the typical cases and formats to gather data. Meanwhile, the data input may involve uncertain factors provided by the intervals; this is typical in cases where inaccuracy should be determined via repetitive trials, e.g., air-pollutant concentrations [42]. Nevertheless, the data entries may be defined as the imprecision of measuring devices; this is a typical and realistic scenario in which inaccuracy is associated to an error rate, e.g., in a continuous blood glucose monitor [3].
Additionally, the framework aimed to define the associated grammar that provides the possibility of generating interval parameters as an alternative to constant values. Thus, the predictive-modeling search engine can include intervals in the form [x ± ∆x] into the candidate models. This feature is intended to represent the uncertainty of a parameter without an a priori associated error or deviation values, e.g., the circadian variability of a patient with diabetes [43].

Intervalized Fitness Functions
The role of the inclusion of uncertainty and the integration of intervals in the assessment criteria has important implications. Because the outcome of the candidate model evaluations implies an interval representation, two generic approaches are available to guide the definition of a fitness function. These can be defined depending on the inclusion relationship between the experimental data and the output of the model. If the pursued model requires the inclusion of the experimental data in the model estimates, thenŶ = f * * (Y 0 , P) ⊇ Y with, e.g., Y 0 and P as proper intervals. Therefore, we obtain the following expression by applying the **-semantic theorem (∀y ∈ Y ) (∃y 0 ∈ Y 0 ) (∃p ∈ P ) y = f (y 0 , p).
According to (17), this inclusion is obtained by minimizing an assessment criterion Tolerance(Y, Y ∧Ŷ) or Tolerance(Ŷ, Y ∨Ŷ) (29) Otherwise, if the requirement is that the estimated solutions are to be contained in the experimental dataset, i.e.,Ŷ = f * (Y 0 , P) ⊆ Y, we apply the *-semantic theorem to obtain Therefore, the evaluation criteria to be minimized in this case are These generic criteria might be sufficient to handle most predictive modeling problems. However, according to the semantic theorems and the inclusion conditions considered, new semantics can be proposed by varying the modalities of the associated intervals.
Finally, the proposed methodology supports the complete set of interval arithmetic operations, including the required f * algorithm, provided that multi-incidence is considered in the generation of the models. Both of them are presented in Section 3 and are essential for computing the final outcomes and training fitness values.

Predictive Modeling Solutions
GE is an evolutionary algorithm that mimics the natural selection process and, in this mimicked process, randomness is crucial. The stochastic nature of GE makes it a method well suited for problems with high requirements, but it should be remembered that GE is a probabilistic stochastic search algorithm, and solution models may vary across different executions. This method can assess a large number of candidate prediction models that iteratively converge to the final model. However, GE algorithms have not been proven to converge to the global optima of the fitness function in a finite time. Thus, an efficient evolutionary algorithm, which typically implies a balance of exploration and exploitation of the solutions, is usually expected to converge to multiple sub-optimal solutions rather than to the global optima. Different strategies can be implemented to select the sub-optimal prediction model that can be used to provide new estimations from unseen data inputs; however, they will not be discussed herein, as they are out of the scope of this study.
The proposed combination of interval methods and GE aims to obtain a soft interval band of estimated intervals that contains, is contained, or is simply the closest to the input intervals that are to be modeled. However, a discrete approximation, e.g., an independent modeling of the upper and lower boundaries of a set of intervals, would aim to obtain two trajectories that coincide or are close to the upper and lower boundary points, regardless of whether the other extreme coincides or is close to the other extreme of the interval. Thus, a non-interval approach could not generate a soft-band solution if the trajectories intersect. In addition, operating directly with intervals allows us to enter conditions and constraints in the form of logical conditions. By defining the semantics, the modal intervals allow us to transform these logical formulas into their equivalent form as inclusion relationships between the modal intervals. Therefore, the interval-based methodology presented herein is indispensable to generate solutions that follow such semantics.

Illustrative Example
In this section, we present an illustrative example in which interval GE is applied to the problem of determining the average water velocity of a river stretch. We used the interval GE methodology to build a model based on the bed slope and flow rate of the river. The method aims to capture the dynamics of the water velocity while managing the uncertainty associated to the physical system. Thus, the expected output of the model is the velocity of the river, whereas the model parameters are the river flow and the bed slope intervals. The interval model that we present in this section is presented as a mere example of the ability of the interval GE methodology to generate models from intervalized data.
The experimental dataset was obtained from a study aimed to generate a water-quality model of the Bajo Ter river (Spain) [44]. The measurements include the bed slope (m/m), the river flow rate (m 3 /s), and the average surface velocity (m/s) of the river. To encompass a wide range of possible flow rates and a variety of bed slope values that could comprise the river's extreme values, the measurements of the experimental values were performed at different locations along the river and at different times, and their values were modeled as intervals.

Modeling River Velocity
The grammar structure used in this case study is simple but flexible, and is adequate for illustrating the mechanism of the method. The grammar defined to estimate river velocity values is presented in Algorithm 1.
The non-terminal [λ] is defined by a derivation rule, with no production involved. This non-terminal, aside from be a path to avoid the generation of any production of [OperatorComplex], allows for the termination of the recursion originated by the production [OperatorComplex][OperatorComplex]. The expressions generated by this grammar comprise a variable number of mathematical operations using the bed slope and flow-rate variables. In order to illustrate the modeling process, Equation (32) whereV t is the predicted velocity of the waters, Bs is the bed slope and Fr is the river flow rate for a given step t.
The goodness of the achieved fit relies on the previous grammar definition, and also on the definition of a criterion to evaluate each model of the population. The definition of the fitness function guides the methodology toward an objective that will be used to model the final solutions. In this case study, we analyze different fitness functions based on the error calculation between the experimental measures and the estimated intervals. The goal is to achieve the most accurate approximation according to the experimental measurements. Here, we propose adjusting our estimates using two general approaches for an interval problem: an external and internal approximation to experimental values. First, we consider an external approximation of the estimates, i.e., we seek solutions that contain the experimental measures, thus fulfilling the following condition Therefore, we define the first assessment criteria as the minimum distance between the experimental and estimated intervals constrained by condition (33) where K is a constant for weighting constraint (33). Secondly, we consider an internal approximation of the estimates, i.e., we seek solutions that are contained in the experimental measures, thus fulfilling the following condition Therefore, we define the first assessment criteria as the minimum distance between the experimental and estimated intervals constrained by condition (35) where K is, once again, a constant for weighing the constraint.
The methodology adjusts the model parameters using previously collected historical data, i.e., the experimental measurements of the river. Additionally, the method has to be adjusted using a set of hyperparameters that express a higher level of structural settings and are not explicitly learned from the data. Table 1 defines the values of these hyperparameters, which were selected empirically and are associated with the operators used in this implementation.

Results
The aim of this case study is to generate and evaluate the solutions for a modeling problem scenario that uses information contained in an experimental database. This section presents the estimated results of the two presented approaches and their assessment in terms of the typical performance metrics used to evaluate predictive accuracy.
The first step in analyzing this case study is to examine the convergence of the evolutionary algorithm. Figure 7 illustrates the average, maximum, and minimum fitness achieved from the best individuals of the proposed method using the assessment criteria C 1 . The fitness value is represented with a logarithmic scale on the y-axis and the number of generations on the horizontal axis. Results show the mean values of the best individual, the fittest 50 individuals and the complete population. Similar to all the experiments presented in this study, performance was assessed after 20 runs, because the method is a probabilistic stochastic search algorithm. Next, we evaluate the estimation results of the approach using the two previously defined assessment criteria. The results are presented in Figures 8 and 9, which are a graphical representation of the outcomes and illustrate the differences between the internal and external approximations. They show the estimated average velocity values of the waters in the analyzed river stretch and the standard deviation (σ) of the computed intervals. The interval results are the mean average of 20 executions, and they estimate a range of possible values, thus including the uncertainty of the system in the estimations. Meanwhile, the standard deviation values indicate the variability of the set of executions in each of the interval bounds.  To evaluate the impact of interval operations on runtime, we implemented a version of the algorithm that does not use interval arithmetic in the assessment evaluations. This version is limited to training models that minimize the root-mean-square error (RMSE) of the mean value between the upper and lower limits of the experimental intervals. Figure 10 represents the mean times' ± standard deviation of 20 trials. Considering that the variation in the computation time between C 1 and C 2 is negligible, we used condition C 1 for these trials.

Discussion
In this study, we presented an interval-based GE methodology, providing both its theoretical foundations and implementation details. We also introduced an innovative approach aiming to resolve regression problems affected by uncertainty. Finally, a series of experiments to test and illustrate the proposed method have been presented. Figure 7 shows how the functioning of the algorithm provides a population that is incrementally fitted through generations. The mean values of the best individuals achieve the best solution in the 80th generation, while the 50 fittest individuals give shape to a typical exponential function that converges at the mean minimum value at approximately the 100th generation. The curve representing all the individuals does not converge at the same values because the functioning of the crossover operator implements a measure to avoid premature convergence of the population. The operator avoids the generation of identical individuals in the same population by replacing these individuals with new randomized chromosomes. Therefore, there is a variable percentage of the population that is generated randomly each generation, and increments the mean value of fitness values.
As mentioned above, the objective is to analyze two complementary and general approaches for an interval problem, an external and an internal approximation with the targets C 1 and C 2 , respectively. Although the defined criteria were not fulfilled in all the experimental points, the obtained intervals suggest a fitted representation of the models for both of them. TheV C1 : [a, b] provides mean standard deviation values of 0.017 and 0.018 andV C2 : [c, d] mean standard deviation values of 0.002 and 0.001, which are satisfactorily low compared with our cost function, indicating a competitive design of the GE algorithm. Figure 7 also shows how the evaluation times of both implementations decrease as the population converges, and that new assessments are not required. After comparing this with the version that does not use intervals, we observed a significant increase in the evaluation times, with an average multiplier of 3.7 over the version that does not use intervals.
The method is strongly influenced by the nature of the addressed problem; this implies a trade-off between the complexity of the solution-seeking space and the convergence of the algorithm. Moreover, it is not exempt from the typical drawbacks of evolutionary computation methods, such as the reproducibility of the results, the risk of premature convergence, and the high computational effort required. The latter is of significant importance, mainly due to the f * function evaluation, which leads to a significant increase in the computational power required by the system. This is especially due to the branch-and-bound nature of such algorithms, which aims to reduce the so-called interval overestimation problem. The use of the f * function involves an inherent algorithmic complexity that grows with the number of non-monotonic variables, the individuals, and the generations.
Typically, modern machine-learning techniques have a high computational costs due to the nature of the involved algorithms. An important example is deep-learning methods, which have recently attracted significant attention from the scientific community and industry, which are willing to invest in the required high-performance hardware systems. A serious drawback of GE is its inefficiency when implemented sequentially, as is the case here. However, GEs have inherent parallel properties that allow them to be successfully parallellized and obtain considerable accelerations (e.g., calculate the fitness of all chromosomes in the current population in parallel). The methodology presented herein still offers ample room for algorithm optimization and a parallelizable architecture, which could lead to a more mature and efficient machine-learning method.
The experiments conducted were performed based on an instance of a real case, which was presented as a practical approach to show the proposal. However, the real utility of the proposed method would be oriented to model more complex real-world problems from imprecise or uncertain data. A good example could be the glycemic monitoring of type 1 diabetics patients with continuous glucose sensors. The monitoring of blood glucose permits subjects detect unusual behavior of their glucose dynamics or even take preventive actions before the occurrence of an adverse events, reducing the risk of complications. Since continuous glucose monitoring devices were launched, the gathered historical data are used by data-driven models to capture the dynamics of blood glucose levels [20]; however, the measurement accuracy of these devices is far from optimal [45]. In addition to this lack of accuracy, we must add the delays that are intrinsically linked with the measurements on these devices [46,47]. Here, machine-learning methods could be applied to acquire and maintain knowledge based on treatment information; however, models will be always biased by the aforementioned measurement errors. This is an open problem in the diabetes technology community, and we hope to contribute to the solution by using the methodology presented in this paper.

Conclusions
Completely eliminating uncertainty from physical systems is often impossible; hence, a proactive approach to manage uncertainty is required. The presented approach acknowledges the existence of uncertainty and uses intervals to power a GE system architecture aiming to resolve regression problems. The methodology involves an evolutionary engine based on a GE approach that implements a set of features to operate with MIA. The context-free grammar itself allows for the attribution of knowledge that delimits the search space of the solutions, whereas the intervals allowed for the uncertainty to be managed naturally and MIA provide semantics that allows to transform logical formulas into inclusion relationships.
This work has explored the first evolutionary algorithm that enables interval arithmetic, thus handling uncertainties by a search algorithm. Furthermore, it has studied the role of uncertainty and integration of intervals in the objective functions of the algorithms, introducing mathematical definitions of fitness functions in a modal interval context. The performed experiments targeted an instance of a real case, which was presented as a practical approach to showcase the proposal. The final utility of the proposed methodology would be oriented to modeling complex real-world problems from imprecise or uncertain data. A deeper analysis of the methodology results would require the modeling of a complex system in real-life scenarios. Although this illustrative example might fall a bit short in providing an indication of performance in more complex problems (usually involving a higher number of records and variables), the convergence of the population, the execution times, and the final estimations results make this approach promising for modelling problems that are influenced by diverse sources of uncertainty.