Evolutionary Approaches to the Identification of Dynamic Processes in the Form of Differential Equations and Their Systems ^†

Abstract

Evolutionary approaches are widely applied in solving various types of problems. The paper considers the application of EvolODE and EvolODES approaches to the identification of dynamic systems. EvolODE helps to obtain a model in the form of an ordinary differential equation without restrictions on the type of the equation. EvolODES searches for a model in the form of an ordinary differential equation system. The algorithmic basis of these approaches is a modified genetic programming algorithm for finding the structure of ordinary differential equations and differential evolution to optimize the values of numerical constants used in the equation. Testing for these approaches on problems in the form of ordinary differential equations and their systems was conducted. The influence of noise present in the data and the sample size on the model error was considered for each of the approaches. The symbolic accuracy of the resulting equations was studied. The proposed approaches make it possible to obtain models in symbolic form. They will provide opportunities for further interpretation and application.

1. Introduction

Nowadays, evolutionary algorithms are applied in solving problems of optimization, classification, clustering, forecasting, and regression in various practical fields [1,2,3,4].

In addition, evolutionary algorithms are an efficient tool for solving identification problems [5,6]. System identification is the process of constructing mathematical models of dynamic systems based on measurements that describe input and output variables of the process under consideration. One of the most informative ways to represent models of dynamical systems are differential equations and their systems [7,8]. This way of describing the processes makes it possible to present the relationship between variables in an interpretable form as well as to reflect the dynamics.

The problem of parametric identification is mainly solved with the help of evolutionary algorithms. Thus, evolutionary approaches to parametric identification are considered in the following papers [9,10].

Parametric identification determines the parameter values of a model with an already known structure. This class of models makes it possible to describe a lot of real processes. For example, there exist processes described by a homogeneous differential equation in economics, biology and chemistry [6].

The presented papers consider a genetic programming algorithm applied to construct a model in the form of an ordinary differential equation (ODE). However, they are all from a class of homogeneous equations. Consequently, the method is not universal due to its applicability to only one class of ODEs [11,12].

However, there is a need to identify processes with a previously unknown structure with the increasing complexity of technological processes. In this case, there is a need for automated methods to construct models of dynamic systems of arbitrary structure. The need to search for such models is faced by specialists in such fields as metallurgy, chemical industry and thermal power plants.

Various methods have been developed due to the current significance of the problem being solved. Nonparametric methods and artificial neural networks are of the greatest importance [13,14,15,16,17].

However, these methods are often unable to construct a model that reflects the true nature of the process in an understandable and interpretable form. The paper presents approaches to the identification of dynamical systems in the form of ODEs and their systems.

2. Materials and Methods

The paper solves an inverse problem of mathematical modelling. The inverse problem involves constructing a system model based on a sample of measurements of inputs and outputs. The approaches proposed in the paper fulfil the search for models in the form of ODEs and their systems. It is possible to apply a genetic programming algorithm to solve the inverse problem of mathematical modelling in order to construct a model in symbolic form. The application of the genetic programming algorithm becomes possible when considering the identification problem as a problem of symbolic regression [18].

2.1. The Evolutionary Approach EvolODE for the Identification of Dynamical Systems in the Form of an ODE

The authors’ approach based on a self-configuring genetic programming algorithm to identify a system in ODE form of arbitrary structure is presented in detail in [19]. The paper denotes this approach as EvolODE (evolutionary approach, ordinary differential equation).

Let us consider a statement of the problem in more detail and highlight the key features of the approach under consideration.

It is necessary to construct a model in ODE form of arbitrary order (1) based on the measurement data of the input x_ti and output y_i characteristics of the process (where t is the number of inputs, i = 1, …, n, where n is sample size):

$${\widehat{y}}^{\left(k\right)}=F\left({\widehat{y}}^{\left(k-1\right)},\dots ,{\widehat{y}}^{\prime },\widehat{y}, x\right).$$

(1)

The authors propose to apply the genetic programming algorithm to construct models in ODE form. Consequently, the terminal set has been changed to encode an ODE in the form of a tree. Changes to the terminal set were made relative to the classical GP algorithm for solving the symbolic regression problem. The terminal set in the proposed approach was represented by a set of all input (x) and output (y) variables, a set of constants, and derivatives ($y^{\prime }, \dots , y^{\left(k-1\right)}$). The root vertex is supplemented with information about the maximum order of the derivative for the current individual at the current iteration. The functional set is standard for the problem of symbolic regression. It contains a set of functions for constructing dependencies. Thus, Figure 1 shows an example of the ODE encoded as a tree.

The main stages of the evolutionary process of the EvolODE algorithm were also modified. Its scheme is presented in Figure 2:

The authors modified such steps of the genetic programming algorithm as the initialization, fitness calculation, and mutation to implement the proposed approach.

The population is generated at the initialization stage. Here, each individual encodes an ODE of the k-order not higher than K, where K is the maximum order of the derivative specified by a researcher. The algorithm determines an order for each individual automatically, but the value of the order does not exceed the value of K. Thus, the order of the equation is not fixed initially, but restricted.

The fitness calculation is carried out in accordance with Equations (2) and (3):

$$fitness=\frac{1}{1+error},$$

(2)

$$error=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\frac{\sqrt{{\left({\widehat{y}}_i-y_i\right)}^2}}{\max \left(y_i\right)-\min \left(y_i\right)},$$

(3)

where n is the sample size, ${\widehat{y}}_i$ is the value of the individual at the i-th point, and y_i represents values from the original sample.

The peculiarity is that mutation could be at the root node. The node contains the maximum order of the derivative for the given individual. In this case, the order of all derivatives included in the tree is randomly changed so that their order does not exceed a new maximum order of the derivative.

The authors of the paper chose the achievement by the algorithm of a given number of generations as a stop criterion.

The authors apply the method of differential evolution to optimize the numerical constants included in the equation and initial conditions [20]. The study presents the efficiency of a hybrid evolutionary approach that combines genetic programming and differential evolution [21].

The presented approach makes it possible to obtain a model of a dynamic system in ODE form of arbitrary structure.

2.2. The Evolutionary Approach EvolODES for the Identification of Dynamical Systems in the Form of an ODE System

An output of a dynamic system is often characterized by several variables in real processes. It is necessary to use the representation in the form of a system of ODEs to develop a model of such processes.

As provided by this problem statement, it is obligatory to construct a model in the form of an ODE system of arbitrary order (4) according to the measurement data of the input x_ti and output y_si characteristics of the process (where t is the number of input and output variables, i = 1, …, n, where n is sample size):

$$\{\begin{array}{c}{y}_1^{\left(k\right)}=f_1\left(t, x_1, \dots ,x_t,y_1,\dots ,y_s\right)\\ y_2^{\left(k\right)}=f_2\left(t, x_1, \dots ,x_t,y_1,\dots ,y_s\right)\\ \dots \\ y_s^{\left(k\right)}=f_s\left(t, x_1, \dots ,x_t,y_1,\dots ,y_s\right)\end{array}.$$

(4)

The order of ODE k will be considered as a limited K-value set by a user.

The paper presents a new approach based on the EvolODE evolutionary approach described above. Thus, according to the new proposed approach, a system of ODEs should be encoded by a set of trees (Figure 3).

The evolutionary approach to the identification of dynamical systems in the form of a system of ODEs will be denoted as EvolODES (Evolutionary, Ordinary Differential Equations System).

The peculiarity of the proposed EvolODES approach is the implementation of a combination of several genetic programming algorithms (k genetic programming algorithms are required for k equations in the system) [22]. Each of them searches for one equation of the system. The number of equations n corresponds to the number of output variables in the data (Figure 4).

Next, a population of individuals is generated by each genetic programming algorithm. Individuals are potential representations of each of the k ODEs in the system.

The fitness is calculated for all individuals of each population (one-time full enumeration of all possible combinations of individuals without regard to order). It is based on the obtained values in each population; the best individual is determined. The fitness calculation for each individual of each population is conducted by substitution into each of the equations of the system. The best individuals are selected as the remaining equations. The calculation is performed in accordance with Equations (5) and (6).

$$fitness=\frac{1}{1+error},$$

(5)

$$error=\frac{{{\displaystyle \sum }}_{i=1}^k{{\displaystyle \sum }}_{j=1}^n{\left(y_{ij}-{\widehat{y}}_{ij}\right)}^2}{sn},$$

(6)

where n is the sample size, s is the number of equations, $y_{ij}$ are values from the original sample, and ${\widehat{y}}_{ij}$ is the value of the model.

According to the new approach, there are no restrictions on the number of ODEs included in the system, their order or their structure. The models obtained in symbolic form can be interpreted in the context of the subject area under consideration.

Both proposed approaches apply a self-configuring type of evolutionary algorithms (genetic programming and differential evolution) [19,21]. The success history adaptation algorithm is applied for genetic programming and differential evolution parameters containing numerical values [23]. The PDP method is applied for parameters with the selection of the type [24]. For the experiment, the authors determined the number of generations for differential evolution as 30. The population contains 10 individuals.

3. Results

Let us consider the proposed approaches in detail. Firstly, consider the efficiency on test tasks represented by ODEs and their systems.

The tasks presented in

Table 1. Test tasks for the evolutionary approach to the identification of dynamical systems in the form of an ODE.

No	Differential Equation	Initial Sample Point
1	$y^{\prime }=\frac{y+x^2\cos x}{x}$	y(x0) = 0.06
2	$y^{\prime }=-\frac{y\ln y}{x}$	y(x0) = e
3	$y^{″}=\frac{2{\left(y^{\prime }\right)}^2}{y-1}$	y(x0) = 2
4	$y^{″}=6y^{\prime }-9y+6x{e}^{3x}$	y(x0) = 3
5	$y^{‴}=4y^{\prime }+24e^{2x}-4\cos \left(2x\right)+8\sin \left(2x\right)$	y(x0) = 1
6	$y^{IV}=y+2\cos \left(x\right)$	y(x0) = 0
7	$y^V=y^{IV}+2x+3$	y(x0) = 2
8	$y^{\prime }=\frac{y+x_1^2\cos x_2}{x_3}$	y($\overline{x}$0) = 0.06
9	$y^{″}=\frac{y^{\prime }+x_1^{2}y{y}^{\prime }}{x_2^2}$	y($\overline{x}$0) = 0
10	$y^{″}=-2y^{\prime }-2y+2x_1^2+8x_2+6$	y($\overline{x}$0) = 3

were used to test the evolutionary approach of EvolDE. ODEs represent tasks of various types and orders. The exhaustive testing and study of the efficiency of this algorithm on tasks with some input variables and with input variables represented by functions are reflected in [25].

Table 2. Testing results of the evolutionary approach to the identification of dynamical systems in the form of an ODE.

No	Sample Size
	150			100			50
	Noise, %
	0	5	10	0	5	10	0	5	10
1	0	0	0.001	0	0	0	0	0	0.0001
2	0	0	0	0	0	0	0	0.0001	0.0001
3	0.0001	0.0008	0.0008	0.0001	0.0008	0.0008	0.0002	0.0009	0.0012
4	0.0002	0.0003	0.0003	0.0002	0.0003	0.0003	0.0003	0.0004	0.0005
5	0	0.0005	0.0007	0	0.0006	0.0007	0.0001	0.0008	0.0008
6	0	0.0001	0.0002	0	0.0002	0.0003	0.0001	0.0004	0.0005
7	0	0.0008	0.0009	0.0001	0.0009	0.0009	0.0001	0.0013	0.0014
8	0.0046	0.0055	0.0059	0.0074	0.0092	0.0104	0.0113	0.0201	0.0287
9	0.0017	0.0047	0.0053	0.002	0.0049	0.0051	0.0033	0.0058	0.0075
10	0.0011	0.0033	0.0041	0.0023	0.0037	0.0042	0.0024	0.0037	0.0039

presents the results of identification according to the data generated by the tasks from Table 1. Data of different volumes were generated without noise and with noise for each task. Noise is distributed normally. The authors conducted 50 runs for each task. Table 2 shows the results by error values (averaged over 50 runs) depending on the sample size and the presence of noise in the data.

The authors consider the dependence of the rate of a chemical reaction on temperature to study the possibility of using this approach in solving applied tasks. A set of values describing the process for a chemical reaction is known. It is necessary to construct a model that describes a change in the rate of a chemical reaction depending on temperature.

It is necessary to estimate the dependence of the reaction rate constant on temperature for the reaction H₂ + I₂ → 2HI. A sample of temperature (T) and rate constant (k) values was applied to construct a model. The sample size was 12 points. As a result, the authors obtained a symbolic model. It corresponds to the well-known Arrhenius Equation (7):

$$\frac{dlnk}{dT}=\frac{E_a}{RT},$$

(7)

where E_a is the activation energy, and R is the universal gas constant [26].

Figure 5 shows graphs corresponding to the original data and points obtained from the model. The mean squared error of the fit was 0.

The tasks presented in

Table 3. Test tasks for the evolutionary approach to the identification of dynamical systems in the form of an ODE system.

No	System of Differential Equations	Initial Sample Point
1	$\{\begin{array}{c}\frac{d{y}_1}{dt}=2y_1-5y_2+3\\ \frac{d{y}_2}{dt}=5y_1-6y_2-1\end{array}$	y1(0) = 6 y2(0) = 5
2	$\{\begin{array}{c}\frac{d{y}_1}{dt}=2y_1+y_2\\ \frac{d{y}_2}{dt}=3y_1+t{e}^t\end{array}$	y1(1) = 1 y2(1) = 2
3	$\{\begin{array}{c}\frac{d{y}_1}{dt}=y_1+2y_2+e^{-2t}\\ \frac{d{y}_2}{dt}=4y_1-y_2\end{array}$	y1(0.1) = 5 y2(0.1) = 8
4	$\{\begin{array}{c}\frac{d{y}_1}{dt}=y_2-y_1^2-y_1\\ \frac{d{y}_2}{dt}=3y_1-y_2^2-y_1\end{array}$	y1(0) = 0 y2(0) = 1
5	$\{\begin{array}{c}\begin{array}{c}\frac{d{y}_1}{dt}=-y_2\\ \frac{d{y}_2}{dt}=y_1\end{array}\\ \begin{array}{c}\frac{d{y}_3}{dt}=y_1-y_4\\ \frac{d{y}_4}{dt}=y_2+y_3\end{array}\end{array}$	y1(0) = 1 y2(0) = 0 y3(0) = 0 y4(0) = 0
6	$\{\begin{array}{c}\begin{array}{c}\frac{d{y}_1}{dt}=2y_1-y_2-y_3\\ \frac{d{y}_2}{dt}=3y_1-2y_2-3y_3+2t\end{array}\\ \frac{d{y}_3}{dt}=2y_3-y_1-y_2-t^2\end{array}$	y1(0) = 2 y2(0) = 3 y3(0) = 2
7	$\{\begin{array}{c}\begin{array}{c}\frac{d{y}_1}{dt}=-10y_1+10y_2\\ \frac{d{y}_2}{dt}=25y_1-y_2-y_1{y}_2\end{array}\\ \frac{d{y}_3}{dt}=y_1{y}_2-10y_3\end{array}$	y1(0) = 1 y2(0) = 0 y3(0) = −1
8	$\{\begin{array}{c}\frac{d{y}_1}{dt}=3y_1+2y_2\\ \frac{d^2{y}_2}{d{t}^2}=-10y_1-y_2\end{array}$	y1(0) = −1 y2(0) = 7
9	$\{\begin{array}{c}\frac{d^2{y}_1}{d{t}^2}=-t{y}_2\\ \frac{d^2{y}_2}{d{t}^2}=y_2-2\frac{d{y}_1}{dt}\end{array}$	y1(1) = 4 y2(1) = −4
10	$\{\begin{array}{c}\frac{d^2{y}_1}{d{t}^2}=\frac{d{y}_1}{dt}-\frac{d{y}_2}{dt}+e^{-t}+\cos t\\ \frac{d^2{y}_2}{d{t}^2}=\frac{d{y}_1}{dt}-\frac{d{y}_2}{dt}-2e^t-\sin t\end{array}$	y1(0) = 2 y2(0) = 0

are applied to test the evolutionary approach of EvolDES. The systems of ODEs of various types and orders represented these tasks.

The analysis of the error of the obtained models was conducted in a similar way to the previous testing scheme. This depended on the amount of data and the presence of noise (

Table 4. Testing results of the evolutionary approach to the identification of dynamical systems in the form of an ODE system.

No	Sample Size
	150			100			50
	Noise, %
	0	5	10	0	5	10	0	5	10
1	0.034	0.036	0.039	0.032	0.037	0.039	0.03	0.038	0.04
2	0.018	0.023	0.024	0.018	0.024	0.024	0.019	0.025	0.024
3	0.018	0.051	0.057	0.018	0.052	0.057	0.019	0.055	0.058
4	0	0	0.001	0	0	0.001	0.002	0.002	0.003
5	0	0	0.002	0	0	0.002	0	0.001	0.003
6	0.022	0.024	0.025	0.031	0.035	0.39	0.038	0.039	0.038
7	0.027	0.031	0.032	0.028	0.029	0.033	0.031	0.032	0.039
8	0	0	0.002	0	0.001	0.002	0.002	0.003	0.003
9	0	0.002	0.002	0	0.002	0.002	0.001	0.002	0.003
10	0.041	0.042	0.053	0.042	0.045	0.054	0.039	0.045	0.053

).

Let us consider the application of this approach in solving a practical task from the field of biology. Let values describing the change in the population size of predators and prey over time be given in some ecosystem. The sample size is 50 points. It is necessary to calculate the change in the number of prey and predator populations over time.

As a result, an expression was obtained represented by a system of differential equations.

$$\begin{gathered} \{\begin{array}{c}\frac{dx}{dt}=0.998x-\frac{xy}{1.023+2x-x}-0.124x^2+0.009x^2\\ \frac{dy}{dt}=\frac{\left(0.255-0.101\frac{y}{x}+0.08\right)y^2}{y}\end{array}= \\ =\{\begin{array}{c}\frac{dx}{dt}=\left(0.998-\frac{y}{1.023+x}-0.115x\right)x\\ \frac{dy}{dt}=\left(0.175-0.101\frac{y}{x}\right)y\end{array} \end{gathered}$$

(8)

The presented system corresponds to the best solution found among 50 runs. It corresponds to an error of 0.0019. The best solution was found 24 times from 50 runs. The averaged error over 50 runs was 0.0031.

The presented system of differential equations corresponds to the Holding–Tanner model. Generally, it is represented as follows [27]:

$$\{\begin{array}{c}\frac{dx}{dt}=\left(a-\frac{cy}{g+x}-ex\right)x\\ \frac{dy}{dt}=\left(b-f\frac{y}{x}\right)y\end{array}$$

(9)

where x is the size of the predator population, y is the size of the prey population, $\frac{cy}{g+x}$ is the predation coefficient, and a, b, f are positive parameters.

Figure 6 shows graphs corresponding to the original data and points obtained from the model represented by the system of differential Equation (8). The mean squared error of the fit is 0.

4. Discussion

The authors conducted the study of two evolutionary approaches for constructing models of dynamical systems in the form of differential equations (EvolDE) and their systems (EvolDES). According to the results, we can observe high accuracy of the models obtained as a result of testing. Changing the sample size within presented limits does not have a significant effect on the error.

An important peculiarity of these approaches is the possibility of obtaining models in symbolic form. However, such representation is expected to be consistent with the true nature of the process. Therefore, it was determined that the models obtained as the identification of test tasks belong to three types: symbolically accurate (full compliance of the structure and parameters of the resulting equation or system with known true ones), symbolically conditionally accurate (the resulting structure corresponds to the known true one after performing the simplest mathematical operations), and approximate ones (solutions that do not belong to two previous types). Figure 7 and Figure 8 show the number of symbolically accurate, symbolically conditionally accurate, and approximate solutions for each of the 50 runs for each task. Figure 7 contains information about the results for the evolutionary approach to the identification of dynamical systems in the form of an ODE. Figure 8 contains information about the results for the evolutionary approach to the identification of dynamical systems in the form of ODE systems.

The successful search for a structure that corresponds to the true nature of the process is also proven by the above models of processes from chemistry and biology.

The paper presents two evolutionary approaches that combine a genetic programming algorithm, differential evolution, and methods for the self-configuring parameters of evolutionary algorithms. The proposed approaches make it possible to identify dynamical systems in the form of ODEs and their systems. In the process of identification, both approaches automatically select the order, structure, and parameters of an ODE, both single and included in the system. The developed approaches make it possible to operate with systems in the form of ODEs of arbitrary order, a number of input and output variables.