Vol 17 No 4, April 2008 1674-1056/2008/17(04)/1196-06

Chinese Physics B

c 2008 Chin. Phys. Soc. ° and IOP Publishing Ltd

Parameters estimation online for Lorenz system by a novel quantum-behaved particle swarm optimization∗ Gao Fei(高 飞)a)† , Li Zhuo-Qiu(李卓球)b) , and Tong Heng-Qing(童恒庆)a) a) Department

of Science, School of Science, Wuhan University of Technology, Wuhan 430070, China b) Department of Engineering and Mechanics, School of Science, Wuhan University of Technology, Wuhan 430070, China (Received 16 April 2007; revised manuscript received 8 October 2007)

This paper proposes a novel quantum-behaved particle swarm optimization (NQPSO) for the estimation of chaos’ unknown parameters by transforming them into nonlinear functions’ optimization. By means of the techniques in the following three aspects: contracting the searching space self-adaptively; boundaries restriction strategy; substituting the particles’ convex combination for their centre of mass, this paper achieves a quite effective search mechanism with fine equilibrium between exploitation and exploration. Details of applying the proposed method and other methods into Lorenz systems are given, and experiments done show that NQPSO has better adaptability, dependability and robustness. It is a successful approach in unknown parameter estimation online especially in the cases with white noises.

Keywords: parameter estimation online, chaos system, quantum particle swarm optimization PACC: 0545

1. Introduction Recently, the traditional trend of analysing and understanding chaos has evolved a new phase in investigation: controlling and utilizing chaos[1] where detecting the unstable periodic orbits and estimating the unknown parameters of the chaos are of vital importance. However, many methods in chaos control and synchronization are based on the cases with accurate deterministic parameters. And they are invalid in the cases with parameters unknown for some special reasons, for example, in secret communications. Then we have to estimate the unknown parameters in chaotic system.[2−4] Actually, parameter estimation of chaotic systems is the chief task to be resolved and of vital significance. These are inverse problems in chaotic dynamic systems, using the data in theory and observations, with deterministic dynamic equations but some specific parameters unknown to be estimated generally. They belong to the classic black-box problems. Many methods have been proposed in estimating the un∗ Project

known parameters of systems.[2−6] The concept of observers for identification of unknown parameters is proposed in Ref.2 and an effective identification for Lorenz system is also given. Through an appropriate fitness function to translate the problems of parameter estimation into that of a parameter optimization[3,4,6] proposed a parameter estimation method based on genetic algorithms (GA), chaotic ant swarm algorithms and particle swarm optimizations (PSO). PSO is a relatively new computational intelligence tool relation to artificial neural nets, Fuzzy Logic, and Evolutionary Computation,[7] developed by Eberhart and Kennedy in 1995,[8,9] inspired by social behaviour of bird flocking or fish schooling. However, unlike GA, PSO does not have genetic operators, such as crossover and mutation. The dynamics of population in PSO resembles the collective behaviour and self-organization of socially intelligent organisms. With the idea leading the population into the most potential district by the best individual, inspired by the PSO and quantum mechanics,[10,11] a quantumbehaved PSO (QPSO)[12] has been put forward recently.

supported by the National Natural Science Foundation of China (Grant No 10647141). [email protected] http://www.iop.org/journals/cpb　http://cpb.iphy.ac.cn

† E-mail:

No. 4

Parameters estimation online for Lorenz system by a novel quantum-behaved · · ·

1197

In this paper, a novel quantum PSO (NQPSO) is proposed to estimate the unknown parameters of the chaotic systems. By means of three techniques to QPSO: contracting the searching space selfadaptively, boundaries restriction (BR) and substituting the particles’ convex combination for their centre of mass, NQPSO is expected to maintain an effective search mechanism with fine equilibrium between exploitation and exploration. As the experiments done with other methods show, the put way is a robust method. The rest of this paper is organized as follows. Section 2 describes the features of classical PSO and QPSO for continuous optimization respectively, while Section 3 explains the proposed NQPSO with the two techniques, contracting the searching space selfadaptively and BR. In Section 4 different methods for estimation of the unknown parameters of Lorenz system are given in detail and the experimental results of different methods are analysed. Section 5 summarizes the paper.

of each particle toward its pbest and gbest locations (PSO without neighbourhood model). Acceleration is weighted by a random term, with separate random numbers being generated for acceleration toward pbest and gbest locations.

2. Classical and behaved PSO

2.2. QPSO

quantum-

2.1. PSO PSO belongs to the category of Swarm Intelligence methods closely related to the methods of Evolutionary Computation.[7−9] The PSO system is initialized with a population of random solutions and searches for optima by updating generations. The procedure for implementing the global version of PSO is given as following. At step k, Each particle Xi (k) = (xi,1 (k), . . . , xi,D (k)) (i = 1, 2, . . . , M ) keeps track of its coordinates in the problem space which are associated with the best solution (fitness) it has achieved (self best) so far. (The fitness value is also stored.) This value is called pbest Pi (k) = (pi,1 (k), . . . , pi,D (k)) . Another ‘best’ value that is tracked by the particle swarm optimizer is the best value, obtained so far by any particle in the neighbours of the particle, called lbest Li (t) = (li,1 (t), . . . , li,D (t)). when a particle takes all the population as its topological neighbours, the best value is a global best and is called gbest G(k) = ((g1 (k), . . . , gD (k)). The PSO concept consists of, at each time step, changing the velocity Vi (k) = (vi,1 (k), . . . , vi,D (k))

Ai,d (k) = rand(0, c1 ) · [pi,d (k) − xi,d (k)], Bi,d (k) = rand(0, c2 ) · [gd (k) − xi,d (k)], vi, d (k + 1) = χ[w · vi, d (k) + Ai,d (k) + Bi,d (k)], xi, d (k + 1) = xi, d (k) + vi, d (k + 1),

(1)

where w is called inertia weight, c1 , c2 are cognition and social acceleration constant respectively, χ is constriction factor, normally χ = 0.9. c1 determines the effect of the distance between the current position of the particle and its best previous position Pi on its velocity. c2 plays a similar role but it concerns the global best previous position G(k), attained by any particle in the neighbourhood. rand(a, b) denotes random in [a, b], in this way, the randomicity is introduced to PSO.

Though PSO converges fast, sometimes it relapses into local optimum easily. Inspired by the classical PSO method and quantum mechanics, an improved edition[12] of PSO (QPSO) with the ideas from quantum computing[11] is also proposed to ameliorate convergence through maintaining more attractors. For each particle Xi (k), let ϕ1 = rand (0,1) , ϕ2 = rand (0,1) , Ti (k) = ϕ1 · Pi (k) + (1 − ϕ1 ) · G(k).

(2)

Using the idea from the centre-of-mass position, we define a new particle S(k) for all Xi (k) as mbest S(k) =

M X Pi . M i=1

(3)

And the j-th dimension xi j (k + 1) of the particle Xi (k + 1) is updated as following Wi j (k) = −β · |Si j (k) − Xi j (k)| · ln (1/u) , xi j (k + 1) = Ti j (k) + Wi j (k)(−1)round(t) .

(4)

where β is the extended coefficient which decreases linearly in interval [1.3, 0.6], u, t are random in (0, 1), round(t) rounds the elements of t to the nearest integers.

1198

Gao Fei et al

3. NQPSO In this section, firstly two techniques, contracting the searching space self-adaptively and BR strategy, are introduced. Then a novel QPSO is proposed and experiments are done for some benchmark optimization problems.

Vol. 17

  xi,d (k + 1),      if ld ≤ xi,d (k + 1) ≤ ud ;      u − 0.3 · rand(0, 1) · (u − lb ), d d d xi,d (k + 1) =  if x (k + 1) > u ;  i,d d     l + 0.3 · rand(0, 1) · (u − lb ),   d d d    if xi,d (k + 1) < ld , (7)

3.1. Contracting the searching space self-adaptively For many experiment results reported suggest that too many generations do not bring the optimum better than the local one,[6] a novel technique through a technique space contraction self-adaptively to avoid exploitation excessively in redundant space was proposed as below. If a local optimum Qg is found, we define a new £ ¤ searching space D(t+1) = a(t+1) , b(t+1) centering Qg £ ¤ from the current searching space D(t) = a(t) , b(t) as below: (t+1)

ai

(t+1)

bi

³ ´ (t) (t) (t) = max xi − γ1 ci , ai , ³ ´ (t) (t) (t) = min xi + γ1 ci , bi ,

(5)

or (t+1)

ai

(t+1)

bi

³ ³ ´ ´ (t) (t) (t) (1) = max min xi − γ2 ci , ai , ai , ³ ³ ´ ´ (t) (t) (t) (1) = min max xi + γ2 ci , bi , bi , (6)

where u = (u1 , u2 , . . . , uD ) is the upper boundary and l = (l1 , l2 , . . . , lD ) is the lower boundary. In this way, the population will not cross the boundaries. Then the BR can eliminate the pneumonia of boundary-crossing and maintain the trend of crossing.

3.3. NQPSO With the techniques in Section 3, we can put a novel quantum particle swarm optimization. Algorithm 1 A Novel Quantum Particle Swarm Optimization (NQPSO) Step 0 Initialization. t = 0, D(0) = D, a(0) = a, b(0) = b, γ1 = 0.5, γ2 = 1.2. Step 1 Outer cycle termination condition judging. If the optimums wanted are found, output them and terminate, otherwise deal with the objective function with the deflection and stretching techniques. Step 2 Generate the initial set A(t) on [a(t), b(t)] randomly, evaluate the fitness and note the best Qg . Step 3 Use QPSO with BR strategy Eq.(7) and let S(k) =

M X

αi Pi

(8)

i=1

where i = 1, 2, . . . , s, pre-given contraction ratio γ1 , γ2 ∈ (0, 1). Then we use QPSO to search in the new space to get a new optimum Q0g , save the better one in Q0g and Qg .

3.2. BR strategy Though vi of PSO directs population to converge to the optimum, sometimes it may make xi out of feasible fields. When xi,d (k + 1) > ud , it suggests that the individual has a trend to cross the ud in its dimension d. To maintain the trend in some degree, we propose a novel strategy BR as below:

PM at each cycle, where each αi ∈ (0, 1) and i=1 αi = 1. Substitute Eq.(8) for Eq.(3) to get a current optimum; update the best one Qg . Step 4 Inner cycle Termination condition Judging. Given δ > 0 enough small, c(t) = (b(t) − a(t))/2, if max c(t) < δ, then x(t), M (t) are accepted, go to Step 1, otherwise go to Step 5. Step 5 Space contraction. Define a new region D(t + 1) = [a(t + 1), b(t + 1)] by Eq.(5) if rand < 0.95, otherwise define a new region D(t + 1) by Eq.(6); t = t + 1, go to Step 2. Through NQPSO, firstly it can avoid premature and find all the optimums sequentially with the help of outer cycle termination judging. Secondly NQPSO can avoid exploitation excessively in redundant space

No. 4

Parameters estimation online for Lorenz system by a novel quantum-behaved · · ·

1199

and search in the most prospective space of the feasible field and avoid premature. And lastly if QPSO in Step3 is valid enough, NQPSO will not contract the searching space, and in this sense QPSO is the special cases of NQPSO. For the function optimization, let x ∈ [−500, 500]10 , make the inner cycle 200 for Griewank and 400 for Rosenbrock function[7] and contracted the searching space 10 and 5 for each other while the single QPSO’s iterations are 2000. Ex.1 Griewank function

f (x) =

µ ¶ N N xi 1 X 2 Y xi − cos √ + 1. 4000 i=1 i i=1

(9)

4. Parameter NQPSO

Ex.2 Banana function

g(x) =

N −1 X

[100(xi+1 − x2i )2 + (xi − 1)2 ].

Fig.2. Comparisons for Rosenbrock.

(10)

i=1

We can get average optimum of 50 independent experiments as reported below (select one of 50): 0.063459642 for Rosenbrock function by NQPSO with space contracting and 0.097536425 by single QPSO as in Fig.1; 70.5869698130682 for Rosenbrock function by NQPSO with space contracting and 226.652920 by single QPSO as in Fig.2. From Fig.1 and Fig.2, it can be concluded that the strategy space contracting can enlarge the performance of NQPSO greatly by improving the feasibility of the solutions.

estimation

by

Though many approaches have been proposed on chaos control and synchronization, most of which are on the premise of knowing accurate value of the unknown parameters in chaotic system. This leads the research on estimating these unknown parameters. To explain the estimating method, we choose the famous Lorenz[1,13] chaos found in 1963, the first canonical chaotic attractor, which has just been mathematically recently confirmed to exist. Firstly a new system is constructed from the Lorenz system:     x ˜˙ = σ(y − x), x ˙ = σ(y − x),          y˜˙ = γx − xz − y,  y˙ = γx − xz − y, →  z˜˙ = xy − ˜bz,  z˙ = xy − bz,         L  L = (x, y, z), ˜ = (˜ x, y˜, z˜), (11) where σ = 10, γ = 28, b is unknown. When b = 8/3, system (11) is chaotic. And when b is unknown, we have to estimate it. Secondly the objective function is chosen as p2 = f (˜b) =

T ·h ° °2 X ° °˜ °L − L° .

(12)

t=0·h

Fig.1. Comparisons for Griewank.

Then the problems of estimation of parameters for chaotic system are transformed into that of nonlinear function optimization. And the smaller the function f is, the better parameter ˜b is. Thirdly let the Lorenz system evolve freely, choose any point as initial point after transience. For a

1200

Gao Fei et al

estimated of the ˜b, use the 4th-order Runge–Kutta method to resolve system (11) with h = 0.01 to get a discrete time³ serial ³ ´ of ³the ´Lorenz ³ system’s ´´ ˜ ˜ standard state at x ˜ b, t , y˜ b, t , z˜ ˜b, t , t = 0h, h, . . . , 300h. The objective function (12) is often multi-modal in the unknown parameters to be estimated ˜b’s definition field with complex structures. And the QPSO and NQPSO are chosen to resolve system (12). Algorithm 2 Parameter Estimation Step 1 Initialization. Randomly generate M individuals ˜b ∈ [1, 5], which are lower and upper bounds of the unknown parameter b of system (11) determined by existing experiences. Step 2 Fitness evaluation. Firstly, for each bi , find its corresponding state. Secondly compute its fitness by Eq.(12). Step 3 Use NQPSO to optimize the unknown parameter unless the appropriate parameters are found. To get a relatively veracious b, 20 experiments are done independently with randomly generated initial population in feasible field, and the final estimation for parameter b is averaged over these 20 results from experiments done. Then we obtain b = 2.66666666672199, which is very close to the real b. Figure 3 gives the result of the 20 experiments, where i is the number of experiments. To test Algorithm2’s validity we add white noises randomly in 3 [−0.1, 0.1] to the initial point input. And we obtain b = 2.6688660155659 as the final result, which is the averaged value of the 20 experiments done independently.

Vol. 17

k(xb , yb , zb ) − (x, y, z)k2 between system (xb , yb , zb ) with estimated b online in the case that initially there is noise in [−0.1, 0.1]3 . Now we can put the online estimation method. Algorithm 3 Parameter Estimation Online Step 1 Resolve system (11) and get an estimation of b through Algorithm 2 with every 30h. Step 2 Use the b to resolve system (11), and make the last point as the new initial point. Step 3 Repeat the above two steps 60 times till the appropriate parameters estimation online gained. We select some numerical methods to realize the online estimation, where TRA represents the trust region algorithms with line-search, Simplex is the simplex methods, and 0.618 is the golden section method. Figures 4 and 5 give the different results of these method online. In Fig.5, the error p between system (xb , yb , zb ) with estimated b online at each 30h and the real system (xb , yb , zb ) with noise in [−0.1, 0.1]3 initially through Algorithm 2 are given.

Fig.4. Results of b online by different methods.

Fig.3. Results of b with and without noises initially.

As Fig.3 shows, effects of the estimation with noise decrease sharply, so Algorithm 2 is limited. Thus we have to estimate online. Note the error p =

Fig.5. Corresponding errors online by different methods.

No. 4

Parameters estimation online for Lorenz system by a novel quantum-behaved · · ·

From Fig.4, we know that, though the estimation of b is good enough at the initial time, it is valid only in short time (as in first 900 h in the Fig.4 and Fig.5 for Lorenz chaos’ intrinsic characters). The error between the system with b estimated and the real Lorenz system becomes bigger along the system (11) which evolves just as Fig.5 shows, and NQPSO is much more effective than the other methods.

5. Conclusion In this paper, a NQPSO based on quantum computation is proposed to estimate the unknown parameter of the Lorenz chaos with two techniques of QPSO in three aspects: contracting the searching space self– adaptively, BR strategy and substituting the particles’ convex combination for their centre of mass. And the experiments done by five different methods for chaotic

References [1] Chen G R and L¨ u J H 2003 Dynamics of the Lorenz System Family: Analysis, Control and Synchronization (Beijing: Science Press) p66 (in Chinese) [2] Guan X P, Peng H P, Li L X and Wang Y Q 2001 Acta Phys. Sin. 50 26 (in Chinese) [3] Dai D, Ma X K, Li F C and You Y 2002 Acta Phys. Sin. 51 2459 (in Chinese) [4] Li L X, Peng H P, Yuan Y X and Wang X D 2007 Acta Phys. Sin. 56 51 (in Chinese) [5] Wang D F 2005 Acta Phys. Sin. 54 1495 (in Chinese) [6] Gao F and Tong H Q 2006 Acta Phys. Sin. 55 577 (in Chinese) [7] Yao X and Xu Y 2006 J. of Comp. Sci. and Tech. 21 1

1201

system done show that the proposed scheme are robust and have the advantages of easy programming and realization. As a matter of fact, NQPSO is a random heuristic method, so the suggested methods are valid probably not deterministically. And as the experiments shows that is valid in short time (in first 900 h) online estimation, while for the long time estimation, we have to appeal to traditionary Lyapunov means in chaos synchronization and control. Though experiments are done only for Lorenz chaos, we can apply the proposed method to other chaotic systems, for instance the Chen system and L¨ u [14] system, straightforwardly.

Acknowledgments We thank the anonymous reviewers for their constructive remarks and comments.

[8] Eberhart R and Shi Y H 2000 Proceedings of the 2000 Congress on Evolutionary Computation (Piscataway, NJ: IEEE) p84 [9] Parsopolos K E and Vrahatis M N 2002 Natural Computing 1 235 [10] Schweizer W 2001 Numerical quantum dynamics(Progress in Theoretical Chemistry and Physics) (Hingham, MA: Springer) p95 [11] Hogg Tad and Portnov Dmitriy 2000 Quantum Opt. Inform. Sci. 128 181 [12] Sun J, Xu W B and Feng B 2005 Proceedings of IEEE International Conference on Systems, Man and Cybernetics p3049 [13] Lorenz E N 1963 J. Atmos. Sci. 20 130 [14] Gao B J and Lu J A 2007 Chin. Phys. 16 666