Parameter identification for discrete memristive chaotic map using adaptive differential evolution algorithm

Since the concept of discrete memristor was proposed, more and more scholars began to study this topic. At present, most of works on the discrete memristor are devoted to the mathematical modeling and digital circuit implementation, but the research on its synchronization control has not received much attention. This paper focuses on the parameter identification for the discrete memristive chaotic map, and a modified intelligent optimization algorithm named adaptive differential evolution algorithm is proposed. To deal with the complex behaviors of hyperchaos and coexisting attractors of the considered discrete memristive chaotic maps, the identification objective function adopts two special parts: time sequences and return maps. Numerical simulations demonstrate that the proposed algorithm has the best performance among six existing algorithms, and it can still accurately identify the parameters of the original system under noise interference.


Introduction
Memristor is a kind of nonlinear electronic component which connects the relationship between magnetic flux and electric charge in the circuit [1]. Since Strukov et al. [2] successfully designed the real memristor in 2008, the research on memristor becomes a hot topic. And it has received extensive attentions in the fields of materials [3], artificial intelligence [4,5], circuits [6,7], etc.
The memristor is completely different from the traditional nonlinear electronic devices. Due to its special nonlinearity, memristance and switching mechanism, it can produce a variety of complexity behaviors when building the oscillation circuits [8]. Consequently, the memristor is very suitable for constructing continuous chaotic circuits to generate numerous special behaviors and apply to the engineering applications. So far, many interesting phenomena are reported in the continuous memristor-based chaotic systems, like hyperchaos [9], extremely multi-stability [10,11], hidden attractors [12,13] or offset-boosting [14,15]. In addition, some secure communication applications based on memristor-based chaotic systems are constantly proposed [16,17]. These studies are basically based on the continuous-time domain, but there are not enough discussions in the discrete-time domain.
As a matter of fact, the structure of discrete nonlinear system is easier for the logic design of digital hardware circuit. And it avoids the disadvantage of param-eter sensitivity in continuous nonlinear system, so the discrete system is more concerned in the practical engineering. Inspired by the mathematical model of continuous memristor, the discussion of discrete memristor becomes a topical issue in recently [18][19][20][21][22][23][24][25][26]. He and Peng et al. proposed the integer-order and fractionalorder discrete memristor mathematical models and introduced these discrete models into the sine map [18] and two high-dimensional chaotic maps [19,20]. Simulations show that these systems can not only enhance the complexity of the original system, but also have abundant dynamic behaviors such as hyperchaos and coexistence attractors. Bao's group presented several hyperchaotic maps based on the discrete memristor [21][22][23][24][25] and verified the numerical results in digital circuits, which laid the foundation for the practical applications of discrete memristor. Kong et al. [26] also reported a new discrete memristive chaotic map, and the offset-boosting and amplitude control are discussed in detail by linear transformation. Although there are more and more studies about the discrete memristive chaotic map, most of them focus on the system modeling and the circuit implementation. The application of the discrete memristive chaotic map has not attracted much attention, especially in its control and synchronization.
Chaos synchronization [27] has great significance in the control of chaotic systems, and it plays an important role in chaotic secure communication. In many real-world applications, the parameters cannot be directly obtained, so parameter identification becomes the key step for the successful synchronization control of chaotic systems. With the rise of artificial intelligence, it is very popular for utilizing intelligent optimization algorithm to identify the parameters of a chaotic system. It converts the parameter identification into a typical multi-dimensional optimization problem [28], and the prerequisite is just some time sequences generated by the original chaotic system. Compared with the traditional identification method, it is simple and not sensitive to the considered system. For the Lorenz, Chen, Rössler and other classical continuous chaotic systems, it emerges the differential evolution (DE) algorithm [29], particle swarm optimization (PSO) algorithm [30], JAYA algorithm [31], artificial bee colony (ABC) algorithm [32], bird swarm algorithm (BSA) [33] to identify the parameters of these systems. In contrast, only a few reports focus on the parameter identification of discrete chaotic maps [34][35][36]. Due to the stronger sensitivity of the parameters in discrete nonlinear systems, it is a challenge to identify the parameters of discrete chaotic maps, especially for the discrete memristive chaotic maps with coexisting attractors (dynamics guided by the initial state) [37]. This motivates us to develop an intelligent method to accurately identify the parameters of the discrete memristive chaotic maps with coexisting attractors.
The main contributions of this paper are summarized: (1) a modified differential evolution algorithm called adaptive differential evolution (ADE) algorithm is proposed. (2) The ADE algorithm is used for the parameter identification of two discrete memristive chaotic maps. As far as the authors know, this is the first time to study the parameter identification of the discrete memristive chaotic map. (3) In order to solve this complex parameter identification problem, the identification process is divided into two special steps to deal with the global search and the local search, respectively. The remainder of this paper is organized as follows. The preliminaries for parameter identification are introduced in Sect. 2. Section 3 introduces the ADE algorithm. Numerical simulations are presented in Sect. 4. Finally, we summarize the conclusion and gives the prospects for the future researches in Sect. 5.

Preliminaries of parameter identification
Preliminary knowledge of parameter identification is introduced in this section. The purpose of parameter identification is to make the parameters of the identification system match the parameters of the original system by minimizing the error of dynamic behavior (i.e., objective function value) between the original system and the identification system. In this paper, two different errors of dynamic behavior are utilized for the identification process, so this section is divided into two subsections: the parameter identification by time sequences and the parameter identification by return maps.

Parameter identification by time sequences
Parameter identification by time series is the traditional identification method [38]. Its objective function is set to the error of generated time sequences between orig- Fig. 1 The main principle of parameter identification by time sequences inal system and identification system. The main principle is illustrated in Fig. 1.
As presented in Fig. 1, considering there is an original chaotic map where and M is the number of time sequences for identification (1 ≤ M ≤ t). t is the iteration number, and x 0 is the initial state. θ = (θ 1 , θ 2 , ..., θ D ) T ∈ R D denotes the original system parameter and D is the number of identified system parameter.
The identification system is The objective function is J , which is calculated by where X k andX k (k=1,2,...,M) are the state at the k-th time of the original system and the identification system, respectively. Thus, the parameter identification by time sequences can be regarded as a multi-dimensional optimization problem which is solved by the intelligent optimization algorithm.

Parameter identification by return maps
Return maps is new approach for parameter identification problem [39]. Although it also collects a period of time sequences as the identification sample, its objective function is based on the similarity of chaotic attractors between the original and the identification system. Compared with the traditional parameter identification by time sequences, the return maps approach is less disturbed by the sensitivity of system parameters, but the disadvantage is high time-consuming [40]. Here, the return maps is only used for the local search to deal with the dilemma of local optimum. The main principle of parameter identification by return maps is shown in Fig. 2. Considering an original system X = F(X, t, θ) and a identification system X = F(X , t,θ). L 1 and L 2 are the number of time sequences in which the original system and the identification system generate attractor structures, respectively. The meaning of other symbols is consistent with those in Sect. 2.1. The specific process of Step (a) and Step (b) in the figure is described as follows.
Step (a): For each point in the original system, find their nearest neighbors in the identification system, and calculate the Euclidean distance between them.
Step (b): For each point in the identification system, find their nearest neighbors in the original system, and calculate their Euclidean distance like Step (a).
Add all the Euclidean distances obtained in Step (a) and Step (b) to get the sum of distance D, and the objective function is calculated by The above process is repeated until all the ranges of parameters are searched done. Finally, the identified parameters with objective function of value 0 are the parameters of the original system.

Adaptive differential evolution algorithm
Here, the differential evolution (DE) algorithm is employed as the intelligent optimization algorithm for parameter identification. DE algorithm is a famous nature-inspired metaheuristics algorithm for solving multi-dimensional optimization problem [41]. It has four main steps: initialization, mutation, crossover and selection.

Process of classical DE algorithm
Initialization: The first stage is to initialize the population in the multi-dimensional variable space.
., n m i,g ), and g, G is the current number and the maximum number of generation, respectively. The lower search bound and the upper search bound are defined by N min = (n 1 min , n 2 min , ..., n m min ) and where j = 1, 2, 3, ..., m. rand(0, 1) is a random number generated between 0 and 1. Mutation: In the nature-inspired metaheuristics algorithm, mutation is taken as a sudden change or perturbation with the random element. Each individual generates a new individual by the mutation step, and the new individual is called the mutant vector V i,g . The most common mutation definition is where r 1 , r 2 and r 3 are randomly selected from [1, N P]. F is the scaling factor, which is a real number between [0, 1]. Crossover: The crossover is employed to generate the trial vector . In this step, mutant vectors are combined with the original individuals of the current population to form the trial vectors. In general, the binomial recombination is performed on each variable where j rand is randomly selected to ensure that the trial vector could get at least one component from the mutant vector ( j rand ∈ [1, m]). C R is called the crossover probability. Selection: The selection is employed to generate the next generation population. Just like the survival of the fittest in biology, individuals with smaller (or bigger) value are preserved. Obviously, the parameter identification is a minimum optimization problem, so individuals with smaller objective function value are retained, and the selection operation is where f (·) is the objective function.

Adaptive strategy
According to the "No free lunch" theorems [42], the classical DE algorithm cannot handle all kinds of optimization problems well. In the face of some highdimensional complex problems, the classical algorithm would meet the problems of slow convergence rate or low precision [43], so it is still necessary to improve the DE algorithm. There are two main control parameters of the DE algorithm: the crossover probability C R and the mutation scale factor F. Especially for the value of C R, it is very important to the search process of DE algorithm. Large C R value is more conducive to the global search, but it is easy to get premature convergence. On the contrary, small C R is more conducive to local search, but the search rate is relatively slow. According to the aggregation degree of individuals in the search process, an adaptive strategy for C R is proposed as where Gb and W b represent the current global optimal vector and the global worst vector, respectively. In addition, according to the situation that the algorithm is easy to fall into premature in the search process, a random mutation factor F is also proposed as Based on the two improvements, the new DE algorithm is named the adaptive differential evolution (ADE) algorithm. Finally, the execution process of ADE algorithm is presented in Algorithm 1. Update C R by Eq. (9) and F by Eq. (10) 6: Update the individual's mutant vector and trial vector by Eq. (6) and Eq. (7), respectively 7: Calculate the individual's objective function value 8: if f (U i,g ) ≤ f (N i,g ) then 9: N i,g+1 is replaced by U i,g 10: else 11: N i,g+1 is replaced by N i,g 12: g = g + 1

Numerical simulations
To verify the effectiveness and advantage of the proposed ADE algorithm, simulations were carried out in two discrete memristive chaotic maps. Five intelligent optimization algorithms, including DE algorithm [43], PSO algorithm [44], JAYA algorithm [45], ABC algorithm [46] and BSA [47] were chosen for the comparison. Simulations were done based on MATLAB 2020a in Intel(R) Core(TM) i7-7700HQ CPU @2.80-GHz with 8 GB RAM.

Identified systems
The second-order discrete memristive chaotic map is [21] x n+1 = x n − kx n cos(y n ), y n+1 = y n + εx n , where k and ε are two positive system parameters. In general, ε is fixed to ε = 1.
Another second-order discrete memristive chaotic map is [26] x n+1 = a sin(x n ) sin(y n ) + (1 − 0.1 | y n |)x n , y n+1 = 0.4x n + y n , where a is the system parameter and a = 1. When the parameters of system (11) are set to k = 2.75, x(0) = 0.5, y(0) = −0.5 (Lyapunov exponents: 0.1804, 0.1123) and the parameters of System (12) are set to a = 3.5, x(0) = 1, y(0) = −2 (Lyapunov exponents: 0.2593, 0.0361), both of the two systems exhibit hyperchaotic behavior, and the phase diagrams of hyperchaotic attractor are shown in Fig. 3. Furthermore, the two discrete memristive chaotic maps have been proved to have the special characteristic of coexisting attractors. So, the dynamic behaviors of the two systems under this set of parameters are very complex. In the following simulations, these parameters k = 2.75, x 0 = 0.5, y 0 = −0.5 and a = 3.5, x 0 = 1, y 0 = −2 were selected as the original system parameters.

Parameter identification by different algorithms
In this subsection, the parameter identification by time sequences was carried out under six intelligent optimization algorithms. The maximum number generation of each algorithm was set to G = 100, and other settings are presented in Table 1. The sample number of time sequences for parameter identification was M = 3. To eliminate the stochastic nature of the algorithms, we performed 100 consecutive algorithm runs. The search ranges of the original system parameters were set as: 0 ≤ k ≤ 4, −4 ≤ x 0 , y 0 ≤ 4 for system (11) and 0 ≤ a ≤ 4, −4 ≤ x 0 , y 0 ≤ 4 for system (12).
The simulation results for parameter identification by time sequences under different algorithms are listed in Table 2. The statistical results of the best J value, average J value and accurate identification rate are given in this table. Generally, the smaller the value of objective function J is, the higher the identification accuracy is. When the error between each identification parameters and the corresponding original parameters is less than 1e-04, it is considered as the accurate identification. From Table 2, it is shown that DE, PSO, BSA  and ADE algorithms can accurately identify the parameters of the original system. But the JAYA algorithm and ABC algorithm are invalid for the two discrete memristive chaotic maps. According to the accurate identifica-tion rate, the ADE algorithm has the best precision, and the accuracy of the modified algorithm is much higher than that of the classical DE algorithm. PSO algorithm and BSA algorithm are ranked in the second and the  ( third, respectively. In addition, it should be noticed that the identification difficulty of system (12) is higher than that of system (11). Although four algorithms have probability to accurately identify the parameters of the original system, they inevitably fall into the local optimum (i.e., the final identified parameters cannot converge to the original parameters). To check the detailed identification results, the final identified parameters of the ADE, PSO and BSA are given in Table 3. It shows that these algorithms are very accurate for the identification of system parameters, but there are many local optimal solutions for the initial state identification. The PSO algorithm and BSA algorithm have more than two local optimal solutions, while the ADE algorithm has only one local optimal solution. Therefore, the ADE algorithm has higher accuracy and better stability.
When the system parameters are fixed to k = 2.75, a = 3.5 and the initial states are set as variable, the J value of parameter identification by time sequences for system (11) and (12) is plotted in Figs. 4 and 5, respectively. In these figures, dark blue indicates that the value of the objective function is low (i.e., the smaller error). The figures illustrate that there is a large area of local optimal region near the global optimum, and the difference of J value between local optimal region and global optimum is only about 1e-02 to 1e-05. Thus, due to the special property of coexisting attractors in the two maps, the change of initial state has big influence on the dynamic behavior, resulting in the situation of multiple local optimal solutions. To solve this problem, the return maps approach is introduced in the following.

Hybrid method for parameter identification
Because there is the local optimal solution problem when using time sequences for parameter identification, the return maps method is employed for the local search. Due to the better performance of ADE algorithm, other algorithms are not considered in this subsection. Thus, a hybrid method is proposed (ADE-Return maps), which consists of two steps: the ADE algorithm of parameter identification by time sequences is utilized for global search, and the parameter identification by return maps is utilized for local When the system parameters are fixed and the initial states are set to the variable, the J value of parameter identification by return maps for system (11) and (12) is plotted in Figs. 6 and 7, respectively. Obviously, the change of J value in return maps method is more drastic than that of the parameter identification by time sequences (it is similar to the difference between the high-rise buildings in the city and the rolling hills). The reason is that the change of the initial states leads to the sudden change of the attractor structure, which is far greater than the error between time sequences. Furthermore, the difference of J value between local optimal region and global optimum is about 0.1; this is also much larger than that of the parameter identification by time sequences. Therefore, the global optimal solution can be easily identified from other solutions in this region.
As listed in Table 3, there is always only one local optimal solution by the ADE algorithm, so the local search by return maps can be carried out according to the local optimal solution (here the search range is set to the local optimal solution ±0.5, which is given Parameter identification by ADE-Return maps for system (11) with noise a SNR 20 dB b SNR 15 dB in Algorithm 2). Hence, such the mechanism ensures that the ADE-Return maps method can find the global optimal solution, as shown in the white star in Fig. 8.

Noise interference
Numerical simulation of parameter identification by ADE-Return maps under noise interference is investigated in this subsection. Noise is a troublesome factor that cannot be ignored in practical engineering, so the additive white Gaussian noise (AWGN) was considered in the time sequences and the attractor structure obtained by the original system.
The numerical simulation results for system (11) are presented in Fig. 9. It is illustrated that the ADE-Return maps successfully find the global optimum (namely x 0 = 0.5, y 0 = −0.5) with the signal-to-noise ratio (SNR) 20 dB. However, it cannot find the correct global optimum with SNR 15 dB. The results for system (12) are illustrated in Fig. 10. It also shows that the ADE-Return maps can identify the global optimum (namely x 0 = 1, y 0 = −2) until the SNR reaches 15 dB. The experiments show that the proposed approach has good anti-noise performance. For the two discrete memristive chaotic maps, it can resist noise larger than SNR 20 dB. In contrast, if the objective function of ADE-Return maps method is changed into the error of time sequences (called ADE-Time sequences), the identification results under AWGN interference are given in Fig. 11. The figure shows that the algorithm cannot find the global optimal solution for systems (11) and (12) when the SNR only reaches 45 and 40 dB, respectively. This is also a reason why the RM method is chosen for the local search.

Conclusion
In this paper, an improved intelligent optimization algorithm called ADE algorithm is proposed for the parameter identification. The target system is selected as two discrete memristive chaotic maps with complex hyperchaotic behaviors and coexisting attractors. In order to solve the problem of falling into the local optimal solutions, the ADE algorithm based on the two objective function of time sequences error and return maps is employed. Numerical simulation results lead to the following conclusions.
(a) For the parameter identification by time sequences, the proposed ADE algorithm outperforms the other five intelligent optimization algorithms, due to the highest identification accuracy and only one local optimal solution. (b) By utilizing the return maps as the objective function, the ADE algorithm successfully jumps out of the local optimal solution and find the global optimum. (c) After considering the interference of AWGN, the proposed approach still finds the global optimum until SNR reaches 15dB.
In the numerical simulations, the factors in realworld applications are studied, including unknown system parameters, unknown initial values and noise. These simulations were therefore theoretically constitute a possibility for the synchronization control of the discrete memristive chaotic maps. What is more important is that the introduced approach may be a useful tool for the future applications of discrete memristor. Our next work is to try to apply this method in the secure communication.