Prediction of dynamical systems by symbolic regression

We study the modeling and prediction of dynamical systems based on conventional models derived from measurements. Such algorithms are highly desirable in situations where the underlying dynamics are hard to model from physical principles or simplified models need to be found. We focus on symbolic regression methods as a part of machine learning. These algorithms are capable of learning an analytically tractable model from data, a highly valuable property. Symbolic regression methods can be considered as generalized regression methods. We investigate two particular algorithms, the so-called fast function extraction which is a generalized linear regression algorithm, and genetic programming which is a very general method. Both are able to combine functions in a certain way such that a good model for the prediction of the temporal evolution of a dynamical system can be identified. We illustrate the algorithms by finding a prediction for the evolution of a harmonic oscillator based on measurements, by detecting an arriving front in an excitable system, and as a real-world application, the prediction of solar power production based on energy production observations at a given site together with the weather forecast.

Figure 1

Illustration of the relation of the methods mentioned in the introduction. We only sketch a part of machine learning, which is relevant for our work.

Figure 2

Illustration of the genetic programming mutation and crossover. The upper left expression tree describes the function $f(x,y)=\sqrt{0.981+sin\left(x\right)}$. Mutation is conducted by picking a random subtree, here the single terminal node $0.981$, and replacing it with a new random expression tree. Similarly, the crossover operator (right) takes two expression trees and swaps two random subtrees.

Figure 3

Illustration of the ranking for Pareto optimization. For a given set of points in the “objectives plane” $({\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_2)$ one labels first all points with rank zero that possess at least one minimum objective, the Pareto front (filled circles, solid line). Then, these points are eliminated from the set and the remaining points are ranked the same way: the rank 1 model emerges (filled diamonds, dashed line). This procedure is repeated until all points are ranked.

Figure 4

Harmonic oscillator study, method FFX: NRMSE (6) versus noise level $\sigma$ for different training set lengths $N$ and fixed $\tau =10$. Sufficiently small noise does not worsen the predictability; i.e., the prediction algorithm stops at the target training NRMSE of 1%. After 0.3 the error does not increase further, since the noise covers the signal completely. Dashed line marks $\sigma =0.17$.

Figure 5

Harmonic oscillator study, method FFX: In solid blue: normalized root mean squared error vs training set length $N$ for $\sigma =0.17$. Dashed green: $e^{-2}/\sqrt{N}$. The error decreases slightly with $N$, but the scaling is much less rapid than $1/\sqrt{N}$.

Figure 6

Harmonic oscillator study, method FFX: Average complexity of the chosen model vs noise amplitude $\sigma$. The form-free structure allows for overfitting. For small noise, the true solution with complexity 2 is found; for higher noise levels, the algorithm starts to fit the noise and more terms are added, reflected by a higher complexity.

Figure 7

Space-time plot of the pulse evolution. $v_i$ is color coded. The front velocity is $v_f=1.28$. Pulse width (full width half maximum) ${\tau }_P=8.4$.

Figure 8

Coupled spiking oscillators, method FFX: Pareto fronts for the different spatiotemporal samplings of the network data. For this plot we use $\rho =0.95$. This leads to the nonconvex shape of the front based on the most information. The models are reevaluated on the testing set.

Figure 9

Coupled spiking oscillators, method FFX. For each feature set, the most accurate model is used as the predictor ${\widehat{v}}_0$. In (a) we show the difference $\delta v_0=v_0-{\widehat{v}}_0$. The upper two curves are shifted by 0.25 and 0.5, respectively. In (b) we compare the time derivative (approximated by the finite difference quotient) of the most accurate model overall (mixed) and the real data. To better visualize the curves we only plot every seventh point. For details see text.

Figure 10

Coupled spiking oscillators, method GP. Averaged Pareto fronts; for each spatiotemporal sampling option, 10 runs are conducted and the resulting complexities and errors are averaged. Error bars represent the standard deviation. For the spatially extended and mixed data sets the errors are smaller than the circle size. The models are reevaluated on the testing set.

Figure 11

Coupled spiking oscillators, method GP. For each feature set, the most accurate model is used as the predictor ${\widehat{v}}_0$. In (a) we show the difference $\delta v_0=v_0-{\widehat{v}}_0$. The upper two curves are shifted by 0.25 and 0.5, respectively. In (b) we compare the time derivative (approximated by the finite difference quotient) of the most accurate model overall (mixed) and the real data. Prediction and true data cannot be distinguished by eye. To better visualize the curves we only plot every seventh point.

Figure 12

Solar power study, Pareto front obtained using FFX. The results for FFX are as accurate as the ones obtained with GP. Test and training set are, however, nicely aligned. This demonstrates not only consistency of the models, but less variability of the models found.

Figure 13

Solar power study, method FFX. (a) Time series of the predicted $\left(\widehat{P}\right)$, and observed $\left(P\right)$ data. We display the results of the first week of August 2015. Similarly to the GP prediction extrema are not particularly well predicted. For the linear model, even the zero values are not well hit. The reason for this is the regression to mean values and the inability of powers to stay at zero for a sufficient time. (b) Histogram of the residuals $\varepsilon =P-\widehat{P}$. Despite different formulas, the histogram of the residuals is asymmetric around zero with a trend to underpredict as well.

Figure 14

Solar power study, average Pareto front obtained using GP: with increasing complexity, training and testing data first behave similarly, then testing deviates strongly indicating overfitting or too small testing data set, respectively. The peaks around complexity 20 are due to two reasons: there are only a few models on the fronts (1–3), and one of them is an extreme outlier.

Figure 15

Solar power study, method GP. (a) Real and predicted time series. We display the results of the first week of August 2015. Prediction used model of complexity 4 which had lowest error on the test set. (b) Histogram of the residuals $\varepsilon =P-\widehat{P}$. The distribution is asymmetric around zero. The model tends to underpredict.

Figure 16

Harmonic oscillator study, method GP: NRMSE (6) versus noise level $\sigma$. We fixed $N=50$ and $\tau =10$. The results are averaged 30 times. The absolute values for the NRMSE are higher compared to the results obtained with FFX (Fig. 16).