Evolutionary modeling-based approach for model errors correction

The inverse problem of using the information of historical data to estimate model errors is one of the science frontier research topics. In this study, we investigate such a problem using the classic Lorenz (1963) equation as a prediction model and the Lorenz equation with a periodic evolutionary function as an accurate representation of reality to generate “observational data.” On the basis of the intelligent features of evolutionary modeling (EM), including self-organization, self-adaptive and self-learning, the dynamic information contained in the historical data can be identified and extracted by computer automatically. Thereby, a new approach is proposed to estimate model errors based on EM in the present paper. Numerical tests demonstrate the ability of the new approach to correct model structural errors. In fact, it can actualize the combination of the statistics and dynamics to certain extent.


Introduction
Monthly, seasonal, annual and inter-annual climatic predictions became the next targets of the frontier research of atmospheric sciences since the successful implementation of short-term weather forecasts.However, with the discovery of chaos, it is well known that the climate system is a complex nonlinear system, and climate prediction is a great challenge to researchers.Moreover, inherent defects of statistical methods and dynamical methods result in significant uncertainties of their applications in climate prediction.Statistical methods are mainly used to search prediction clues from historical observational data, based on the historical behavior and statistical regularities of the climate system.The results of statistical model are sometimes satisfactory, but such statistical regularities and the prediction models derived from such statistical methods are often unstable due to the effects of nonlinearity of the climate system (Gu, 1958).
The prediction problem based on dynamical methods can be regarded as a problem of initial value of differential equations.It implies that it is necessary to ensure sufficiently precise initial values and boundary conditions as well as the accuracy of the prediction model in order to obtain a reliable prediction.In fact, these conditions cannot be satisfied completely.Despite the continuous improvement and optimization of numerical model and data assimilation systems, it is difficult to improve the time limit of weather forecasts beyond two weeks (Lorenz, 1965(Lorenz, , 1969) ) due to the complexity of the atmospheric motion.Therefore, it is crucial to develop alternative ways to improve the capability of the climate prediction model (Schubert, 1985;Vannitsem and Toth, 2002;Chou, 2003a, b;Li and Ding, 2011).
The dynamics-statistics approach is an important step for the improvement of climate prediction models.Extensive research has been carried out using a combination of statistics and dynamics.Qiu and Chou (1988) noted that the observational data could be regarded as sufficiently precise solutions for atmospheric models, and could be used to correct the errors of model parameters by solving an inverse problem.Cao (1993) and Feng (2004) suggested that the atmospheric motion is an irreversible process, and they introduced a memory function that can make use of historical observational data to deduce the self-memorization equation for the atmospheric motion including multi-time observational data.Chou (2007) argued that past daily weather changes, especially the recent evolutionary status of the atmosphere, contain the information of numerical model errors.It should be used to correct the model errors.Some related mathematical and numerical issues in the geophysical fluid dynamics and climate dynamics have been discussed (Li and Wang, 2008).These methods have shown to be useful techniques for climate prediction in numerical experiments and applications.It is an evidence of the effectiveness of the dynamics-statistics approach.
In fact, the idea of the correction of the model errors is old.Many studies provided some objective methods of estimating model errors, but their use for improving model performance, is comparatively small.Examples are given by Schubert (1985), Klinker and Sardeshmukh (1992).D' Andrea and Vautard (2000) proposed a methodology for the correction of systematic errors in a simplified atmospheric general circulation model, and confirmed that this improvement actually stems from the flow dependence of the model error.Such flow dependence was found in the Euro-Atlantic sector, while similar attempts to establish this relation in other sectors of the globe failed.Vannitsem and Toth (2002) investigated the short-term dynamics of model errors by means of numerical analysis of the Lorenz (1984) low-order atmospheric system, and found that the short-term mean square error evolution is mainly characterized by an initial quadratic or linear behavior, depending on the dynamical properties of the model error source terms.Also, there are studies of the impact of model errors associated with parameter errors.Duan and Zhang (2010) investigated the effect of initial errors and model parameter errors on a significant "spring predictability barrier" (SPB) for El Niño events.They inferred that initial errors, rather than model parameter errors, may be the dominant source of uncertainties that cause a significant SPB for El Niño events.Mu el al. (2010) suggested that conditional nonlinear optimal perturbation approach (CNOP) can be applied to estimate model parameter errors for El Niño events, but did not consider other kinds of model errors at present.Despite this, it is expected that CNOP will play an important role in the studies of atmospheric and oceanic sciences.
However, most of the works mentioned above remain in the theoretical stage, and many issues need to be addressed for their application in weather forecasts (Chou et al., 2007).For example, traditional methods for solving inverse problems of differential equations face one essential problem of ill-posed characteristic, such as the instability of approximate solutions.The models describing the complex climate system can only approximately describe the major dynamic processes of atmospheric motion.Therefore, not all the details are described in these prediction models.In fact, climate change can be viewed as a long-range evolutionary process with self-adaptation.So, it is unnecessary to describe all the characteristics of the problem in detail, and all it requires is a solution in accordance with the evolutionary laws of the nature.
In the case that only limited information is known on a dynamic system, a possibility is to replace human intelligence with computational intelligence in some steps of the traditional modeling, including the development of assumptions, the construction and the calculation of the model.This is the idea of so-called evolutionary modeling (EM; Cao et al., 2000).Evolutionary computation is a mathematical modeling approach based on the evolutionary laws of nature (Back et al., 1997;Cao et al., 2000).In EM, simple coding techniques can be used to represent a variety of complex mathematical structures.The approximate solutions of inverse problem can be automatically searched by computer programs.This step is iterated for some number of generations until the termination criterion of the run has been satisfied.The characteristics of evolutionary algorithms are selforganizing, self-adaptive and self-learning.Natural selection, namely survival of the fittest, and evolution strategies provide many conveniences and advantages for solving inverse problems of complex differential equations.
In the present work, the problem of how to use combination of the dynamics and statistics to correct model errors has been studied.On the basis of the EM method, a new approach is proposed to correct model errors.Using the new approach, dynamic information of model errors can be extracted from historical observational data.The results of numerical experiments based on the Lorenz (1963) model have been preliminarily validated in terms of their ability to correct model errors.It must be noted that the model error studied here mainly refers to the dynamic structural error in a prediction model.

Algorithm of EM
EM is a calculation strategy and method on the basis of the principles of biological natural selection and genetic inheritance (Back et al., 1997;Cao et al., 2000).Complex problems can be solved by means of computer simulations of the natural evolutionary process in EM.In accordance with rules of natural evolution, such as the "survival of the fittest", the complexity of the computing and solution can be controlled in EM.Evolutionary computation provides an effective approach to deal with those problems (Holland, 1975), which are complex and intractable problems by traditional mathematics.Generally, differential equations can be expressed as binary trees in EM (Fig. 1).In these binary trees, each tree structure uniquely corresponds to a complex or simple function within the solution space of the problem.These functions generally consist of arithmetic operations, logic operations, variables, constants, and elementary functions, such as sine, cosine, and exponential function.The tree structure can be changed dynamically by using genetic operators, such as crossover and mutation.The evolution of these tree structures continues from one generation to another generation until the termination criterion of the run has been satisfied.The complexity of a function can be controlled by setting the depths of the layers of a binary tree.The detailed algorithm of EM can be found in Cao et al. (2000).As an example, the classic three-dimensional Lorenz (1963) model can be written as the form: (1) Here, the three parameters, a, r and b, are positive, and are called the Prandtl number, the Rayleigh number and a physical proportion related to convective scale, respectively.In the present study, the three parameter values are 10.0, 28.0 and 8/3, respectively.When r > 24.74, the Lorenz system exhibits chaotic behavior.The binary tree form of the Lorenz model is shown in Fig. 1.
In evolutionary algorithms, an individual mathematical model structure corresponds to a unique binary tree (Fig. 1).Three popular methods can be used to implement evolutionary computation, namely, genetic algorithm (GA), gene expression programming (GEP; Ferreira, 2001), and genetic programming (GP; Saiedian, 1997).The GP and GEP methods will be selected to implement EM program in the present study.

Correction scheme for model errors
Numerical prediction models are the dynamic equations constructed according to the basic laws of atmospheric and climate motion as currently understood.Undoubtedly, these models use simplified and approximate equations.It can only describe main features of atmospheric motion.In fact, the problem is that the disturbance terms have to be overlooked, or certain important dynamic elements are not included due to the limitation of the existing knowledge.One of the current measures for correcting errors is a post-processing strategy for model prediction results, i.e. statistically correct on the prediction results.However, the omission of dynamic characteristics of the model error affects the quality of the results of the statistical correction.

Description of the inverse problem
An inverse problem is a general framework that is used to convert observed measurements into information about a physical object or system that we are interested in.For example, if we have measurements of the Earth's gravity field, then we might ask the question: "given the data that we have available, what can we say about the density distribution of the Earth in that area?"The solution to this problem (i.e. the density distribution that best matches the data) is useful because it generally tells us something about a physical parameter that we cannot directly observe (Keller, 1976; http://en.wikipedia.org/wiki/Inverseproblem).
In contrast, with a forward problem, an inverse problem of differential equations is based on existing results to deduce the cause of the results.In other words, on the basis of the solution or partial solution of one differential equation, the unknown components of the equation can be deduced.In terms of differential equations, Leverentiev and his colleagues argued that the inverse problem can be defined as follows: an inverse problem for a partial differential equation is any problem involving the determination of the coefficients or right-hand side of a partial differential equation on the basis of certain functionals of the solution of the equation (Leverentiev et al., 2003).Without loss of generality, the general form of differential equations can be written as (2) Here, L is the differential operator, and u(x, y, t) is the solution of the differential equation.The function f (x, y, t) is the right-hand side source term of the equation.When L is unknown, this is called an inverse problem of operator identification; furthermore, when the right-hand source term f (x, y, t) is unknown, this is called a source-term inverse problem.The ill-posedness and nonlinearity of an inverse problem make the related theories and solutions more complex and difficult than those of a forward problem.
In order to investigate the inverse problem based on the use of the information of historical data to estimate model errors, Eq. ( 1) is regarded as an approximate model derived from practice, while the corresponding true model is represented by (3) Here, E(x, y, z) is the potential unknown error term, which is a composite function, such as trigonometric, exponential and power functions.The symbol " " represents connecting nodes.Generally, the symbol can only be basic arithmetic operators, such as "+, −, ×, ÷".If L(x, y, z) is the main dynamic structure of the function, " " can be determined.In this case, the symbol " " usually is "+".Even for a simple dynamic equation, such as Eq. ( 1), it is difficult to determine the specific mathematical expression of E(x, y, z) using traditional methods.On the basis of the definition of the inverse problem, determining the error term E(x, y, z) in Eq. ( 3) is a classic inverse problem, namely, solving the source term on the right-hand side.
The limitation of those traditional methods is that model structure must be determined in advance, and then model parameters can be estimated.However, it generally depends on personal experience to construct model structure.It is a complicated problem, especially when the data space takes on the form of hyper-surfaces with multi-parameters and multivariables.

Error correction scheme of evolutionary algorithm-based prediction model
The main steps for correcting prediction model error terms through historical data are as follows: an approximate prediction model and an accurate prediction model are constructed, respectively.The solutions of the approximate model are regarded as the prediction results of a specific system, and the solutions of the accurate model are regarded as "observational data".Subsequently, the EM algorithm is used to automatically search for the error term.During modeling, a large ensemble of mathematical model structures is randomly formed by means of the search in the basic function library using the evolution program.Simultaneously, the parameters of the individual model structures will be continuously estimated and optimized by the parameter optimization module.According to certain criteria, poor individual model structures will be eliminated, and good model structures will be retained, which can be continuously used to EM of the next generation.Finally, several mathematical models with relatively minimal errors will be obtained.Firstly, in Eq. (3), for example, the function L(x, y, z) must be unchanged in the EM.Then, the function will be expressed as the binary tree hierarchy (Fig. 1), which can be represented in a certain program code, and then can be placed into the evolutionary process.Obviously, obtaining the expression of the error function E(x, y, z) is an inverse problem that solves the right-hand side of the differential equation.Here, L(x, y, z) is the approximate model, and the correction scheme of the prediction equation is as follows.
1. Constructing Gene(L), the binary tree of the main function L(x, y, z).
2. Specifying randomly the node symbols and the constants contained in the error terms E(x, y, z), and the nodes include {+, −, ×, ÷, sin, cos, ln, exp} in the present study.The populations Pop with a certain size will be generated randomly.
3. Forming several binary tree structures Gene i (E) of the error functions using the evolutionary algorithm GP and GEP.Meanwhile, Gene i (E) will concurrently evolve with the main binary tree Gene(L) by applying evolutionary operations, including selection, hybridization and mutation.By controlling the depths of the binary tree layers of Gene i (E), several complex-controllable composite binary trees (i.e. the revised function) can be obtained.
4. Computing and evaluating the adaptive values of these composite binary trees in combination with "observational data".Retaining the superior individuals Gene i (E) in the generation, and then these superior individuals will be used as EM of the next generation.
5. Repeating the above procedures until a new function body meets the predetermined conditions (i.e. the termination condition).

Numerical experiments
In order to test the performance of the error correction scheme presented in this study, it is necessary to generate "observational data".To deal with this issue, a classic threedimensional Lorenz model is regarded as an inaccurate prediction equation L(x, y, z).The Lorenz model with an error term E x = 5 sin(sin(x)) is regarded as an accurate prediction equation, and connecting node = "+ , as follows: Thus, Eq. ( 3) can be written as 700 integral steps of Eq. ( 4) are regarded as "observed data" (solid line in Fig. 2).The corresponding prediction data can be generated by the classic Lorenz model (Eq.1), and the solution of the approximation model has been shown in Fig. 2 (dash line).Comparing the prediction data with the observational data, it is clear that there is some small difference between the prediction data and the observational data when the integral time is less than about 180 steps.Thereafter, the differences gradually become larger, and even the phases of evolution of the prediction data and the observational data are sometimes in opposition.
In accordance with the correction scheme, the GP and GEP are used in the modeling approach.The details of the related parameters are already listed in Table 1, and a relative error function is used as the fitness function, Equation ( 5) is called the evaluation function, which is a quantitative indicator of the quality of model, and also provides the termination condition of the EM program.Here, x i is the "observational data", and x i is the solutions of the revised model by EM.It is the most crucial element of the evolutionary algorithm for the problem to be solved.The selection of the fitness function has a significant impact on the convergence speed and results of the evolutionary algorithm.The evaluation function (Eq.5) is comparatively consistent with the law of evolution, the "survival of the fittest".It retains a suitable adaptive range and can create a state space that can be addressed by the evolutionary algorithm.A suitable evaluation function can improve the convergence x , E (2) x and E x , E (5) x , E x , respectively.velocity of the evolutionary algorithm and avoid excessively rapid or slow evolution of populations.Thus, it can ensure to obtain several reliable solutions with high quality.A few prediction equations have been listed in Table 2 obtained by the error correction scheme.
Generally, it is difficult to have the knowledge on the forms and positions of the error terms.The error corrections have been implemented for the three functions L x (x, y), L y (x, y, z), L z (x, y, z), respectively.Based on the results of the EM, seven representative error functions are shown in Table 2 for the error corrections on the first function L x (x, y) in the Lorenz model (Eq.1).The error function E (0) x obtained by EM is identical to the real error term because the function E (0) x can be transformed into the original error function.Indeed, it is no doubt that the probability for obtaining an identical error function is very low in climate models based on the existing technology.The reason is that the variables and physical mechanisms of the atmospheric motion are more complicated than those of the current example.The other six error functions can be divided into two categories.The first three error functions only contain time variable t, and the style of these functions is similar to that of the real error function.In contrast, there are more than one variable contained in the other three error functions.It is easy to find that most of these error functions obtained by EM are similar to the given error term E x , namely, trigonometric functions are the main form in these error functions.
The fitting curves and prediction curves are shown in Fig. 3.The sample size of the training data is 500, and that of the prediction data is 200.Obviously, they are all good x, z ln(cos(e x ) + cos(cos(z + 0.5)) + 5) 0.03 0.06 fitting results for different error functions by EM, and the prediction errors are largely reduced by the correction before reaching the first 150 steps.Even though the prediction error becomes slightly larger in the final 50 steps, the overall varying trend remains unaffected.It is also indicated that there exists a time limit in the correction efficiency of these error functions.In Fig. 3a and b, the adaptive value calculated by Eq. ( 5) is less than 2, and the maximum of the predication error is only 8.41.Obviously, it is a successful error correction.Figure 3c and d show the results of the revised equation corrected by three relatively more complex error items: E x , E (5) x and E (6) x .They contain more than one variable, including components (x, y, z) and t.Even though they are more complicated than the original error function, the effects of the error correction are still satisfactory.Except for the prediction result of the revised equation by E (5) x , the average prediction errors are very small.The corrected prediction results are similar for the prediction of the variables y and z by error correction on the first function L x (x, y) in Eq. ( 4).It shows the ability of EM to reduce model error using historical observational data.
The results of the error corrections on the second function L y (x, y, z) of the Lorenz model are listed in Table 3.The original error term E x is in the first function L x (x, y) in Eq. ( 4), but the error corrections on the second function L y (x, y, z) can achieve the effect of reducing prediction error by the automatic searching process of the EM and the correction effect is almost the same as with the function L x (x, y).Apart from the error functions E (1) y and E (2) y , the fitting values of other four revised functions are less than 1.0.The fitting value of E (1) y is the largest in the six cases; however, it can be seen that the prediction error also can be largely reduced (Fig. 4).For the other three error functions E  for the improvement of prediction precision.In future work, the idea presented in this paper will be gradually applied to higher dimensional ordinary differential and partial differential equations.The ultimate goal is to apply the present method to a variety of numerical prediction models.Currently, EM can be optimized by means of parallel computing or cloud computing.It can be realized by implementing multithread parallel computing not only for the correction of model errors of single dynamic equation but also for that of multi-equations.Based on this, the computational efficiency of EM can be improved greatly.In fact, the numerical calculations used in this paper all are based on multicore multi-threaded algorithms, which obviously improves the computational efficiency of EM.So we believe that the present study could be applied to more complex systems for the correction of model errors including the model errors of large-scale numerical prediction models.

Fig. 1 .
Fig. 1.Hierarchical structure diagram of binary tree of the Lorenz model.

Fig. 2 .
Fig. 2. The "observed value" (x 0 , z 0 ) and the corresponding prediction data (x, z) of the prediction model, where the sample size is 700.(a) Variable x in the Lorenz model; and (b) variable z in the Lorenz model.

Fig. 3 .
Fig. 3. Error correction based on the first component L x (x, y) using EM.(a) and (b) are the fitting and prediction data of the corrected models using the error function E (1) and (d) are the same as (a) and (b), but for E (4) average prediction errors of the revised function are very small.The prediction errors of the revised model by the error function E y become relatively larger when the integration time is greater than about 150 units.The results further indicate that there exist time limits in the correction efficiency by these error functions.It Nonlin.Processes Geophys., 19, 439-447, 2012 www.nonlin-processes-geophys.net/19/439/2012/

Table 1 .
Parameters in the EM procedure.

Table 2 .
Error functions and errors of the equations obtained by EM (for variable x).

Table 3 .
Error functions and errors of the equations obtained by EM (for variable y).

Table 4 .
Error functions and errors of the equations obtained by EM (for variable z).