Synchronization of non-smooth chaotic systems via an improved reservoir computing

The reservoir computing (RC) is increasingly used to learn the synchronization behavior of chaotic systems as well as the dynamical behavior of complex systems, but it is scarcely applied in studying synchronization of non-smooth chaotic systems likely due to its complexity leading to the unimpressive effect. Here proposes a simulated annealing-based differential evolution (SADE) algorithm for the optimal parameter selection in the reservoir, and constructs an improved RC model for synchronization, which can work well not only for non-smooth chaotic systems but for smooth ones. Extensive simulations show that the trained RC model with optimal parameters has far longer prediction time than those with empirical and random parameters. More importantly, the well-trained RC system can be well synchronized to its original chaotic system as well as its replicate RC system via one shared signal, whereas the traditional RC system with empirical or random parameters fails for some chaotic systems, particularly for some non-smooth chaotic systems.

www.nature.com/scientificreports/At the aspect of synchronization, Ibáñez-Soria et al. 36 make use of RC to successfully detect generalized synchronization between two coupled Rössler systems.Guo et al. 33 investigate the transfer learning of chaotic systems via RC method from the perspective of synchronization-based state inference.Weng et al. 34 adopt RC model to reconstruct chaotic attractors and achieve synchronization and cascade synchronization between well-trained chaotic systems by sharing one common signal.Hu et al. 35 train chaotic attractors via RC method, and realize synchronization between linearly coupled reservoir computers by control the coupling strength.However, all above mentioned works aim at smooth chaotic systems, few investigations involve in the prediction and learning of non-smooth chaotic systems whose dynamical behaviors may be more complicated than the smooth ones since the non-differential or discontinuous vector field in non-smooth systems easily lead to the singularity.No coincidentally, our extensive simulations show that the traditional RC method has poor prediction effect for some non-smooth chaotic time series, especially for discrete-time chaotic systems, please see Table 1 and Fig. 1.On the other hand, non-smooth chaotic systems have many applications in engineering fields, and they are frequently used to model circuit systems, such as Chua circuit system.Motivated by the above aspects, the paper proposes a hybrid RC method with optimal parameter selection algorithm for improving the prediction ability for non-smooth as well as smooth chaotic systems, and based on the improved RC model we carry out synchronization between well-trained RC systems via the variable substitution control.

Results
In this section, based on the improved Reservoir Computing (RC) machine learning model, we will use chaotic data to check that the proposed RC model with SADE algorithm can work well for the non-smooth chaotic systems as well as the smooth chaotic systems, and to verify the availability of PC synchronization between the constructed and trained RC systems.To begin with, the fourth-order Runge-Kutta method is used to calculate numerical solutions of smooth chaotic systems (such as, Lorenz and Rössler systems shown in 37,38 ) and piecewise smooth chaotic systems (such as Chua and PLUC systems), the observation time step t = 0.02 , and the numerical solution is transformed into the data at the range of [−1, 1] for convenience.For all the used chaotic systems, the first 3000 observations are used to train the RC system after discarding the leading appropriate amount of observations to eliminate transient states, and the first 100 observations of the trained RC system are discarded for the steady running.
Furthermore, the basic parameters of RC model are set as follows.In Eq. 9, σ = 1 for which the weight W in is randomly selected from the uniform distribution of (−σ , σ ) , the leakage rate of the reservoir is α = 0.2 , and = 1 × 10 −8 in Eq. 10.

Main parameters selection and model-free prediction
This part will employ SADE algorithm mentioned in the previous section to select the appropriate parameters (N, p, r) for sparse Erdös-Rényi (ER) random network of RC systems, make model-free prediction of the trained RC system with the selected parameters, and compare prediction effect with those with empirical parameters as well as randomly-selected one.
Based on the presented SADE algorithm, the selected parameters in ER random networks and the prediction duration of RC for different chaotic systems are listed in Table 1, where each in 50 experiments has the same prediction duration for the improved RC at the given precision threshold.It can be seen that for PLUC chaotic system, the corresponding trained RC system can well predict up to 20 seconds (about 31 units of Lyapunov time) at the square mean prediction error of least 10 −3 order.But for smooth chaotic systems, Lorenz and Rössler systems, surprisingly, the trained RC for Rössler system can predict up to 110 seconds (about 9 units of Lyapunov time) at the prediction error of least 10 −8 order.Similarly, the trained RC also performed well for Lorenz system, and the prediction time is up to 20 units of Lyapunov time at the error of 10 −4 order, more 7 and 11 units than the empirical and randomly-selected parameter approaches, respectively.For the discrete Hénon and Ikeda chaotic systems, the trained RC can predict up to 60 steps (about 25 and 8 units of Lyapunov time, respectively) at the prediction error of least 10 −3 order and 10 −5 order, respectively.
More importantly, for all the considered chaotic systems including discrete, piecewise smooth, and smooth chaotic systems, the prediction time duration of RC with optimally-selected parameters is far longer at the same of prediction precision threshold than those with empirical parameters ( N = 500 , p = 0.25 and r = 0.95 ) as Table 1.A list of the prediction time duration of RC systems for different chaotic systems at the given precision threshold.Here, symbols * and # represents the prediction time using the empirical and randomly- selected parameters, respectively.Lyapunov time is defined as the reciprocal of the largest Lyapunov exponent of a system.Each prediction duration of improved RC systems averages 50 realizations in the SADE algorithm, and the listed optimally-selected parameters are one of 50 experiments.well as those with randomly-selected parameters, see the last two columns in Table 1, and the corresponding state variable evolution shown in Fig. 1.Specifically, the Lorenz-trained RC system with optimal parameters can predict up to 15-seconds time span, about 20 units of Lyapunov time, at the error of 5.67 × 10 −4 , but that with empirical parameters (random parameters) just can predict up to about 10 (6.6)-seconds time span, about 13 (9) units of Lyapunov time, at the same error.Surprisingly, for another smooth chaotic system, Rössler system, the trained RC system with optimal parameters can reach up to 110-seconds prediction time span, about 9 units of Lyapunov time, with high accuracy ( 8.89 × 10 −8 order error), but those with empirical parameters and with randomly-selected parameters has much less prediction time at the same of error, 22.96 seconds and 31.88 seconds (about 2 and 3 units of Lyapunov time), respectively.Please see Fig. 1e,f for the state evolution.

Systems
Particularly for non-smooth chaotic systems, the trained RC systems with optimal parameters have the overwhelming prediction effect.For example, at the same error, the Chua-trained RC system with optimal parameters has 40-seconds prediction time, far bigger than those with empirical parameters (18.78-seconds prediction time) and with randomly-selected parameters (11.66-seconds prediction time), as shown in Fig. 1c.For PLUC system, there are huge difference of prediction effect between using optimal parameters and using empirical (or randomly-selected) parameters, the former reaches up to 20-seconds prediction time at the error listed in Table 1, and the latter just predicts less than 1-second time span at the same error, as shown in Fig. 1d.For discrete Hénon and Ikeda chaotic systems, their trained RC systems with optimal parameters can predict up to 60 steps, three time as much as those with empirical (or randomly-selected) parameters where both have the same prediction steps, exactly as cyan and purple vertical lines overlap each other in Fig. 1a,b.
In brief, compared with the traditional RC systems (using empirical or randomly-selected parameters), the improved RC systems based on SADE parameter selection algorithm can make a far better prediction not only for non-smooth chaotic systems, but for smooth chaotic systems (such as Lorenz system and Rössler system).In other words, the improved RC system can better learn the chaotic systems including smooth and non-smooth www.nature.com/scientificreports/ones.In next section, the advantage will be shown again for PC synchronization behaviors between chaotic system and its trained RC system as well as between two trained RC systems.

PC synchronization between chaotic systems and trained RC systems
Next, we use the proposed RC model based on SADE algorithm to train the RC systems for all the chaotic systems under consideration, and perform PC synchronization between the chaotic system and the trained RC system, between two trained RC systems, and between a chaotic system and two trained RC systems through the corresponding drive variables shown in Table 2. Consequently, synchronization behavior of four non-smooth chaotic systems are shown in Figs. 2, 3, and 4 respectively for the cases between the chaotic system and its trained RC system, between two trained systems, and between Ikeda system and its two trained RC systems where the x component of Ikeda system is used to drive the first RC system, and the y component of the first RC system is used to drive the second RC system.The results demonstrate the good realization of PC synchronization, implying the availability of the improved RC model with SADE algorithm for non-smooth chaotic systems, of course as well as for smooth chaotic systems (synchronization evolution plots are not provided here).
However, using empirical parameters and randomly-selected parameters, not all four non-smooth chaotic systems can be synchronized to their corresponding trained RC systems.As shown in Figs. 5 and 6, for PLUC system and Ikeda system, the synchronization error does not tend to zero as time goes, implying the failure of synchronization between the chaotic system and its trained RC system, as well as between two trained RC systems.Further, for Chua system, synchronization also fails between two trained RC systems with randomlyselected parameters (see Fig. 6a).These show the advantage of the proposed RC model with SADE algorithm and the importance of main parameters of ER network in RC model.
Interestingly, whatever using the empirical parameter or randomly-selected parameter, let alone the optimal parameter, the RC system trained from Hénon system works well, and it is able to be well synchronized with the original Hénon system as well as with its replicate RC system.A possible reason is the fact that the error w ≡ 0 in (8) as x is the drive variable.Just because of the perfect prediction ability for chaotic systems, the well-trained RC system is more closely approximated to the original chaotic system, and thus it can successfully demonstrate synchronization behavior occurring in the drive-response chaotic system described in the second section.Theoretically, if the RC is sufficiently approximated to the original system, the variable driving the original drive-response system into synchronization is able to drive the trained RC system and its replica (or the original system) into synchronization.www.nature.com/scientificreports/Conversely, the variable not driving the original one into synchronization should not drive two trained RC systems into synchronization.Although the simple grid search approach may yield a as at least good result as the SADE approach in this work does, it spends far more time to find the optimal parameters than the SADE method.On the other hand, compared with the empirical or randomly-selected parameter approaches, the proposed method has significant advantage and performs well even for non-smooth chaotic systems, as shown in the above section.As a trade-off method, the proposed RC model with SADE algorithm is a good choice for training chaotic data and achieving PC synchronization.It may be applicable to many other chaotic systems including some complicated piecewise smooth chaotic systems only if avoiding the overfitting problem.Besides PC synchronization in the driveresponse system where the coupling is unidirectional, the proposed RC model may be used to realize synchronization in the bi-directionally coupled non-smooth chaotic systems even the coupled networked chaotic systems.

Non-smooth chaotic systems
Here briefly reviews the used non-smooth chaotic systems, including piecewise smooth chaotic systems (Chua system and a piecewise linear unified chaotic system) and discrete chaotic systems (Hénon map and Ikeda map).
Hénon Map 41 : Ikeda Map 42 : The phase diagrams of Hénon Map and Ikeda Map are shown in Fig. 8.

PC synchronization principles
In early years, Pecora and Carroll presented a variable substitution method 43,44 for achieving synchronization between chaotic systems.The basic idea is to split the original system into two subsystems, copy one of the subsystems, and then replace the corresponding variable of the replicated subsystem with the variable (called the drive (1)  variable) of the other subsystem, and eventually realize synchronization between the replicated subsystem and the original system.The diagram is clearly shown in Fig. 9, and the mathematical model is described as follows.
For an n-dimensional chaotic system, Divide arbitrarily it into two subsystems, for instance, and where v = (u 1 , u 2 , . . ., u m ) and w = (u m+1 , u m+2 , . . ., u n ).Now create a new system dw ′ dt = h(v ′ , w ′ ) identical to the subsystem of dw dt = h(v, w) , substitute the variable v ′ in the function h(v ′ , w ′ ) with the corresponding variable v (called the drive variable), and get the called response system dw ′ dt = h(v, w ′ ) .Together with the drive system (5), it follows that The subsystems of w and w ′ is synchronized each other, called PC synchronization here, if the error w = w − w ′ → 0 as t → ∞ , which can be theoretically checked through the negative conditional Lyapunov exponent of the drive-response system (8).Through our checking, the system can achieve synchronization when the appropriate variables listed in Table 2 are selected as the drive variable.www.nature.com/scientificreports/ In recent years, RC model is frequently applied and studied on account of its excellent predictive effect for nonlinear time series.As many researchers know, the reservoir consisting of ER random network structures plays a key role in this machine learning model, which can well characterize the connected neurons in the recurrent neural network, and is with strong dynamical characteristics.
In general, the parameters (network size N, connected probability p and network spectral radius r) of ER random networks in the reservoir are selected by experience, but the parameters have significant impact on the prediction effect, and therefore the inappropriate selection may lead to the poor prediction.In this section, we introduce a simulated annealing-based differential evolution (SADE) method for parameter selection, and present a reservoir computing system with SADE algorithm.Furtherly, we introduce how to realize PC synchronization via the presented reservoir computing model.The reservoir computing system is composed of the input layer, the reservoir part, and the output layer.Denote u(t) = (u 1 (t), u 2 (t), . . ., u n (t)) T the input vector in the input layer, r(t) = (r 1 (t), r 2 (t), . . ., r m (t)) T the state vector in the reservoir, and y(t) = (y 1 (t), y 2 (t), . . ., y l (t)) T the output vector in the output layer.Then, the run of reservoir computers can be divided into three parts: Initialization, Training and Prediction, which are described as follows:

Reservoir computing systems with SADE parameter selection algorithms
• Initialization: Generate the reservoir's adjacency matrix A and the input weight matrix W in .In general, A is a sparse ER random network matrix with the appropriate network size N, connected probability p and the spectral radius r, and its non-zero entries are randomly chosen from the uniform distribution of [−1, 1] .The elements of matrix W in are randomly chosen from the uniform distribution of (-σ , σ). • Training: Train the output weight matrix W out by using the training data set.Here we update the reservoir network state vector r(t) according to the following iterative law: where α is the "leakage" rate limited in the interval of [0, 1], tanh is the hyperbolic tangent function, b in = 1 is the offset, and A and W in are the mentioned above adjacency matrix and input weights, respectively.To eliminate the transient state of the reservoir computer with the initial r(0) , the first τ time-steps states are removed out.Next, the state vector and the training data denoted {s(t)|t = τ , τ + 1, . . ., τ + T} are used to train the output W out .Here the training data is taken as the input data {u(t)|t = τ , τ + 1, . . ., τ + T} where u(t) is used to predict next forward u(t + 1) .According to the reference 20 , the output weights W out can be calculated by Where I is an identity matrix and is a ridge regression parameter for avoiding overfitting.X is the matrix whose t-th column is (b out , s(t), r(t)) T , and Y is the one whose t-th column is s(t + 1).• Prediction: With the trained output weights, the reservoir computing system based on Eqs. 9, 10, 11 can run autonomously, and use the current output vector y(t) of ( 11) to predict u(t + 1) at the next moment, and so on.
As mentioned in the beginning of this section, the number of nodes N, the connected probability p and the spectral radius r in ER random network of the reservoir, have a great impact on the prediction effect particularly for some chaotic systems.Our experiments have shown that the parameters selected by experience lead to poor prediction for some discrete chaotic systems and non-smooth chaotic systems.
For better prediction, we will introduce the simulated annealing-based differential evolution algorithm (SADE) to select the optimal parameters via which we train the reservoir computing system and then realize synchronization between trained reservoir computing systems.   is the trained RC system where the input u ′ = (x ′ t , y ′ t , z ′ t ) T , and create the replicate RC system denoted as where the input u ′′ = (x ′ t , y ′′ t , z ′′ t ) T for which the initial value y ′′ 0 and z ′′ 0 are different from y ′ 0 and z ′ 0 , and x ′ t is the drive variable come from the system (12).
For the cascade synchronization between a chaotic system and two RC systems shown in Fig. 11c, take x and y components as the drive variables, and suppose the chaotic system and two trained RC systems are described as where u(t) = (x t , y t , z t ) T , u ′ (t) = (x t , y ′ t , z ′ t ) T and u ′′ (t) = (x ′′ t , y ′ t , z ′′ t ) T , and the initial values u 0 = (x 0 , y 0 , z 0 ) T , u ′ 0 = (x 0 , y ′ 0 , z ′ 0 ) T and u ′′ 0 = (x ′′ 0 , y ′ 0 , z ′′ 0 ) T .Here, the chaotic system (14a) is called the drive system, RC systems (14b) and (14c) are called the first and second response system, respectively.

Data availability
All data generated and all programs used during the current study are available from the corresponding author on reasonable request.

Figure 1 .
Figure 1.Comparison of model-free prediction effect between optimal parameters (red dashed line) based on SADE algorithm, empirical parameters (purple dash-dot line) and random parameters (cyan dotted line) for Hénon (a), Ikeda (b), Chua (c), PLUC (d), Lorenz (e) and Rössler (f) chaotic systems.Here the cyan, purple, red solid vertical lines represent the prediction time (or step) span respectively for random, empirical and optimal parameters at the same square mean error listed in Table1.

Figure 2 .
Figure 2. The phase diagram of synchronization state evolution between the chaotic drive system and its trained RC response system where x (y) component is the drive variable in a,c,d (b).(a) Chua system, (b) PULC system, (c) Hénon system, (d) Ikeda system.

Figure 3 .
Figure 3. State variable evolution plots of synchronization between two RC systems trained from (a) Chua system, (b) PLUC system, (c) Hénon system, and (d) Ikeda system, respectively.Here the drive variable is the same as in Fig. 2.

Figure 4 .Figure 5 .
Figure 4. State variable evolution plots of cascade synchronization between (a) Ikeda system and the second trained RC system and (b) Ikeda system and the two trained RC systems.

Figure 6 .
Figure 6.Synchronization error between two trained RC systems with different types of parameters where RC models are trained from (a) Chua system, (b) PLUC system, (c) Hénon system, and (d) Ikeda system, respectively.Here the drive variable is the same as in Fig. 2.

Figure 9 .
Figure 9.The schematic diagram for PC synchronization.

Figure 11 .
Figure 11.Diagrams for PC synchronization between different systems.(a) between chaotic system and its trained RC system, (b) between two trained RC systems, (c) cascade PC synchronization between chaotic system and its two trained RC systems.

Table 2 .
A list of driving and response variables for four non-smooth chaotic systems.