- Email: [email protected]
Quantum-behaved particle swarm optimization algorithm for the reconstruction of fiber Bragg grating sensor strain profiles
Optics Communications 285 (2012) 539–545
Contents lists available at SciVerse ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
Quantum-behaved particle swarm optimization algorithm for the reconstruction of fiber Bragg grating sensor strain profiles Hongbo Zou a, b,⁎, Dakai Liang a, Jie Zeng a, Lihang Feng a a b
State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, PR China College of Electric Engineering and Renewable Energy, Three Gorges University, Yichang 443002, PR China
a r t i c l e
i n f o
Article history: Received 27 June 2011 Received in revised form 17 November 2011 Accepted 21 November 2011 Available online 5 December 2011 Keywords: Fiber Bragg grating Particle swarm optimization Quantum-behaved particle swarm optimization Strain distribution
a b s t r a c t A quantum-behaved particle swarm optimization algorithm used for demodulating the strain profile along a fiber Bragg grating (FBG) from its reflection spectrum has been demonstrated. By combining the transfer matrix method for calculating the reflection spectrum of a FBG and the quantum-behaved particle swarm optimization technique, we develop a new technique for the reconstruction of FBG sensor strain profiles. This proposed technique is verified through numerical example reconstructions of FBG sensor simulated strain profile cases and numerical simulations reveal good agreement between the original and the reconstructed strain profiles. © 2011 Elsevier B.V. All rights reserved.
1. Introduction Fiber Bragg grating (FBG) sensors have all the advantages of optical fiber sensors such as small size, light weight and immunity to electromagnetic interference. Further, they possess flexible multiplexing capability because of the inherent wavelength-encoded strain or temperature response. Hence, FBG sensors have been considered excellent sensing elements for measuring strain or temperature. For a FBG sensor subjected to constant axial strain along the FBG, only the Bragg wavelength shift needs to be measured to gain the strain data because the Bragg wavelength shift is linearly proportional to the applied strain [1]. Many researches and applications have been developed based on this property [2–5]. But in practice, the axial strain applied to the FBG sensor may be non-constant such as linear, quadratic, even high-order or discontinuous. In this condition, the reflection spectrum shape of the FBG sensor is not only shifted but also distorted with multi-peaks [6]. Therefore, it is unsuitable for only measuring the Bragg wavelength shift to obtain the strain data. Fortunately, the reflection spectrum of a FBG sensor can be easily measured with a conventional optical spectrum analyzer (OSA). Thus, how to reconstruct the strain distribution along the FBG sensor from its reflection spectrum becomes a meaningful inverse problem. The simplest approach to solve this inverse problem is the intensity spectrum based (ISB) method [7], but it is only applicable to monotonic ⁎ Corresponding author at: State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, PR China. Tel.: + 86 25 84893466; fax: + 86 25 84892294. E-mail address: [email protected] (H. Zou). 0030-4018/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2011.11.069
strain distribution with increasing or decreasing profiles. Fourier transform approach [8] is also developed with the advantage of a short reconstruction time, but it is only suitable for weak FBG and requires both the intensity and phase information of reflection spectrum. Some other reconstruction methods have been reported, such as the iterative GLM method [9], the time-frequency signal analysis [10], and the layerpeeling algorithm [11] and so on. However, all these methods need both the intensity and phase information. Generally, the measurement of phase spectrum needs a more complex setup compared with that of intensity spectrum. Therefore, the intensity based methods is preferable to the methods which required the phase of reflection spectrum. Recently some evolutionary algorithms have been developed to solve this inverse problem, mainly including the genetic algorithm (GA) [12] and simulated annealing (SA) algorithm [13]. These methods only require intensity information, so they play a successful role in the distributed strain profile reconstruction. However, GA has the drawbacks of long running time and the requirements of optimizing the parameters by trials to skip the local best solutions. SA has the similar shortcomings compared with GA. Therefore, some modifications for GA and SA are reported such as chaos genetic algorithm (CGA) [14], adaptive simulated annealing (ASA) algorithm [15], simulated annealing evolutionary (SAE) algorithm [16]. The particle swarm optimization (PSO) algorithm is a member of a wider class of swarm intelligence methods used for solving global optimization problems. This method was originally proposed by Kennedy as a simulation of social behavior of bird flocking and fish schooling in 1995 [17]. Instead of using evolutionary operators to manipulate the individuals as in other evolutionary algorithms, PSO relies on the exchange of information between individuals. Each particle in
540
H. Zou et al. / Optics Communications 285 (2012) 539–545
PSO flies in search space with a velocity, which is dynamically adjusted according to its own former information. Compared with GA and other similar evolutionary algorithms, PSO has some attractive characteristics and in many cases proved to be more effective [18]. However, PSO is not a global optimization algorithm [19]. In [20–21], Sun et al. introduce quantum theory into PSO and propose a quantum-behaved PSO (QPSO) algorithm, which can be guaranteed theoretically to find optimal solution in search space. The experimental results on some widely used benchmark functions show that the QPSO works better than standard PSO [20–21] and should be a promising algorithm. In this work, QPSO algorithm was demonstrated to solve the inverse problem and compared with PSO algorithm. The basic principles are as follows. First, we use the transfer matrix method to calculate the reconstructed grating spectrums in different grating parameters. Then we use the QPSO algorithm to search a reconstructed spectrum that in agreement with the original spectrum. Thus from the parameters corresponding to this reconstructed spectrum, we can get the applied strain profile on the FBG. This paper is structured as follows. The transfer matrix method used to calculate the reflection spectrum of the FBG is presented in Section 2. The basic principle of QPSO is introduced in Section 3. The implementation of the reconstruction and the numerical simulation results are given in Section 4. Finally, the conclusions are provided in Section 5. 2. Transfer matrix method The reflection spectrum of a non-uniform FBG can be calculated by using transfer matrix method [22]. In this section, we summarized the most important mathematical equations of the transfer matrix method. This approach divides the grating uniformly into M small sections, each with uniform coupling properties. Define Ri and Si to be the field amplitudes after traversing the ith section, the propagation through this uniform section is described as:
Ri Si
¼ Fi
Ri−1 Si−1
ð1Þ
σ 6 coshðγΔzÞ−i γ sinhðγΔzÞ Fi ¼ 6 4 k i sinhðγΔzÞ γ
3 k −i sinhðγΔzÞ 7 γ 7 5 σ coshðγΔzÞ þ i sinhðγΔzÞ γ
R R ¼F 0 S S0 f 11 f 21
f 12 f 22
π νδneff λ
ð7Þ
σ¼
2π π neff þ δneff − λ Λ ðzÞ
ð8Þ
where neff is the mode effective index of refraction of the optical fiber, δn eff is the mean mode effective index of refraction variation, ν is the fringe visibility of the index variation, λ is the wavelength of FBG. In conclusion, for any given parameters L, Λ0, neff, δneff , ν, pe and ε(z), we can use the above equations to calculate the corresponding reflection spectrum. 3. Quantum-behaved particle swarm optimization algorithm PSO is very easy to be understood and implemented. In the PSO, the fitness function is a measurement of the distance between the calculated reflection intensity spectra and the target reflection intensity spectra. PSO has already been tried and tested in various standard optimization problems with excellent results [17]. However, the main disadvantage of PSO is that global convergence cannot be guaranteed [19]. To deal with this problem, concept of a global convergence guaranteed method called as quantum-behaved PSO (QPSO), was developed by Sun et al. in 2004 [20–21]. In the QPSO, the state of a particle is described by wave function Ψ(X, t), instead of position and velocity of the PSO. The dynamic behavior of the particle is widely divergent form that of the particle in the PSO systems in that the exact values of Xi and Vi cannot be determined simultaneously. We can only learn the probability of the particle's appearing in position Xi from probability density function |Ψ(X, t)| 2. The particles move according to the following Eqs: X i ðt þ 1Þ ¼ P i −β⋅jmbest i −X i ðt Þj⋅lnð1=uÞ if k≥0:5
ð9Þ
X i ðt þ 1Þ ¼ P i þ β⋅jmbest i −X i ðt Þj⋅lnð1=uÞ if kb0:5
ð10Þ
where ð2Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where γ ¼ k2 −σ 2 , k is the ‘ac’ coupling coefficient, σ is a general ‘dc’ self-coupling coefficient, Δz is the each section length. Thus, the total grating propagation mode can be expressed by:
k¼
2
F¼
where L is the length of FBG, Λ0 is the period of FBG before the ε(z) is applied, pe is the effective photo-elastic constant. Then k and σ in Eq. (2) can be expressed by:
ð3Þ ¼ F 1 F 2 ⋯F M
ð4Þ
ð5Þ
When a FBG is subject to an inhomogeneous axial strain ε(z), the equivalent period of the FBG Λ(z) [7] can be expressed as: Λ ðzÞ ¼ Λ 0 ½1 þ ð1−pe ÞεðzÞ0≤z≤L
mbest ¼
N 1X P N d¼1 g;d
ð6Þ
ð11Þ ð12Þ
In Eqs. (9)–(12), mbest is the mean best position defined as the mean of all the best positions of the population, g represents the index of the best particle among all the particles in the population, k, u and φ are random numbers distributed uniformly in [0,1], respectively. β, called contraction-expansion coefficient, is the only parameter in QPSO algorithm. It can be tuned to control the convergence speed of the algorithms. The β value is adaptively allocated as per the Eq. (13): β ¼ β max −ðβ max −β min Þ⋅t=t max
and the reflection coefficient of the FBG is given by: 2 f r ¼ 21 f 11
P i ¼ φ⋅P i;d þ ð1−φÞ⋅P g;d
ð13Þ
where βmax is the initial contraction-expansion factor value, βmin is the final contraction-expansion factor value, t is the current generation number and tmax is the maximum number of generations. The procedure for implementing the QPSO is presented by the following: Step 1. Initialization of swarm positions: Initialize a population of particles with random positions in the D-dimensional problem space using a same probability distribution function.
H. Zou et al. / Optics Communications 285 (2012) 539–545
Step 2. Evaluation of particle's fitness: Evaluate the fitness value of each particle. We define the fitness function by measuring the error between the reconstructed spectrum Rcal(λ) and the target spectrum Rtar(λ) in this work.
f ¼
H h X
2
Table 2 Strain profiles applied to the FBG sensor. Case
Strain profile (με) (z/mm)
Case1: Case2: Case3: Case4:
linear distribution quadratic distribution sine distribution discontinuous distribution
ð14Þ
Rcal ðλi Þ−Rtar ðλi Þ
541
100z + 10 5z2 + 50z + 10 500sin(0.2πz) 50z 0 ≤ z ≤ L/2 − 100z L/2 b z ≤ L
i¼1
Step 3. Step 4.
Step 5.
Step 6. Step 7.
where λi is the sampled wavelength of the FBG spectrum and the total number is H. Rcal(λ) and Rtar(λ) are closely related to the strain profile ε(z). The task of QPSO is to search the input D-dimensional space effectively to find a strain profile ε(z) along the grating axis that minimizes the fitness function. In this work, we suppose that the strain profile ε(z) is determined by a finite set of parameters αi ∈ [a, b](i = 1, . . . ,D), thus the fitness function optimization problem is converted to optimizing the D parameters. Updating of global point: Calculate the mbest using Eq. (12). Comparison to pbest (personal best): Compare each particle's fitness with the particle's pbest. If the current value is better than pbest, then set the pbest value equal to the current value and the pbest location equal to the current location in the D-dimensional space. Comparison to gbest (global best): Compare the fitness with the population's overall previous best. If the current value is better than gbest, then reset gbest to the current particle's array index and value. Updating of particles’ position: Change the position of the particles according to Eqs. (9) and (10). Repeating the evolutionary cycle: Loop to Step 2 until a stop criterion is met, a maximum number of generations adopted in this work.
4. Implementation and numerical simulation To test the proposed strain distribution reconstruction of FBG sensor, four simulated strain distribution cases are investigated: a linear distribution (case 1), a quadratic distribution (case 2), a sine distribution (case 3) and a discontinuous distribution (case 4). The parameters of the FBG sensor employed in this work are given in Table 1 and the four strain distribution cases are shown in Table 2. First the simulated spectra (as ‘measured’ spectra) for each case are calculated using the transfer matrix method, applying the specific strain profile to be reconstructed. Then we use the QPSO algorithm to reconstruct the applied strain profile. The QPSO parameters used for strain profiles reconstruction are listed in Table 3. Based on these standard parameters, the QPSO algorithm is suitable for reconstruction of different strain profiles with almost the same convergence quality. Furthermore, the performance of QPSO algorithm was demonstrated through each case and compared with that of the PSO algorithm. The PSO parameters in present work are shown in Table 4.
Table 3 Parameters of QPSO algorithm. Parameter
Value
Max contraction-expansion coefficient (βmax) Min contraction-expansion coefficient (βmin) Max generation (maxG) Population size (N) Number of grating segments (M) Number of spectral points (H) Rang of wavelength(λ)
1 0.5 100 20 20 201 2 nm
4.1. Linear strain distribution (case 1) The original and reconstructed reflection spectrums of the FBG sensor applied with a linear strain distribution (case 1) are shown in Fig. 1. Fig. 2 plots the original and reconstructed strain profiles. Fig. 3 presents the comparison of convergence processes of QPSO and PSO in the case 1. From Figs. 1 and 2, the reconstructed results of PSO and QPSO are extremely close to the original spectrum and applied strain distribution in case 1. From Fig. 3, it can be seen that though PSO shows faster convergence, better results are found with QPSO. The slow convergence speed corresponds to good global search ability, while fast speed results in good local search ability. This clearly shows that QPSO is successful in having better global search capability in comparison to PSO and hence a better optimal result is obtained using QPSO.
4.2. Quadratic strain distribution (case 2) A quadratic strain distribution (case 2) is applied to the FBG sensor; the original and reconstructed results are shown in Figs. 4 and 5. In Figs. 4 and 5, original spectrum and strain profile are plotted in dashed curve, while PSO reconstructed spectrum and strain profile are plotted in asterisk curve and QPSO reconstructed spectrum and strain profile are plotted in circle curve. Similarly, Fig. 6 presents the comparison of convergence processes of QPSO and PSO in case 2. In Fig. 6, convergence processes of PSO and QPSO are plotted in dashed and solid curve, respectively. From Figs. 4 and 5, the reconstructed results of PSO and QPSO are extremely close to the original spectrum and applied strain distribution in case 2. Fig. 6 shows that though PSO shows faster convergence, better results are found with QPSO.
Table 4 Parameters of PSO algorithm. Table 1 Parameters of FBG sensor. Parameter
Value
Length (L) Initial period (Λ0) Mode effective index of refraction (neff) Mean index variation (δneff ) Fringe visibility (ν) Effective strain-optic constant (pe)
10 mm 534.48 nm 1.45 2 × 10− 4 1 0.22
Parameter
Value
Inertia factor (w) Cognitive learning rate (c1) Social learning rate (c2) Max generation (maxG) Population size (N) Number of grating segments (M) Number of spectral points (H) Rang of wavelength(λ)
1 2 2 100 20 20 201 2 nm
542
H. Zou et al. / Optics Communications 285 (2012) 539–545
Fig. 1. Original and reconstructed reflected spectra applied to the FBG sensor in case 1.
Fig. 4. Original and reconstructed reflected spectra applied to the FBG sensor in case 2.
Fig. 2. Original and reconstructed strain profiles applied to the FBG sensor in case 1.
Fig. 5. Original and reconstructed strain profiles applied to the FBG sensor in case 2.
4.3. Sine strain distribution (case 3) A sine strain distribution (case 3) is applied to the FBG sensor. The simulation results are shown in Figs. 7, 8 and 9. In Figs. 7 and 8, the dashed curve represents the original spectrum and strain profile, the asterisk curve represents the PSO reconstructed spectrum and strain profile, and the circle curve represents the QPSO reconstructed
Fig. 3. Comparison of convergence process of QPSO and PSO in case 1.
spectrum and strain profile. In Fig. 9, the dashed and solid curves represent the convergence processes of PSO and QPSO, respectively. From Figs. 7 and 8, the reconstructed results of PSO and QPSO are very close to the original spectrum and applied strain distribution in case 3. Fig. 9 reveals that though PSO shows faster convergence, better results are found with QPSO.
Fig. 6. Comparison of convergence process of QPSO and PSO in case 2.
H. Zou et al. / Optics Communications 285 (2012) 539–545
543
Fig. 7. Original and reconstructed reflected spectra applied to the FBG sensor in case 3.
Fig. 10. Original and reconstructed reflected spectra applied to the FBG sensor in case 4.
Fig. 8. Original and reconstructed strain profiles applied to the FBG sensor in case 3.
Fig. 11. Original and reconstructed strain profiles applied to the FBG sensor in case 4.
4.4. Discontinuous strain distribution (case 4) A discontinuous strain distribution (case 4) is applied to the FBG sensor. The simulation results are shown in Figs. 10, 11 and 12. In Figs. 10 and 11, the dashed curve is the original spectrum and strain profile, the asterisk curve is the PSO reconstructed spectrum and strain profile, and the circle curve is the QPSO reconstructed spectrum
Fig. 9. Comparison of convergence process of QPSO and PSO in case 3.
and strain profile. In Fig. 12, the dashed and solid curves are the convergence processes of PSO and QPSO, respectively. From Figs. 10 and 11, the reconstructed results of PSO and QPSO are close to the original spectrum and applied strain distribution in case 4. From Fig. 12, it can be shown that though PSO shows faster convergence, better results are found with QPSO.
Fig. 12. Comparison of convergence process of QPSO and PSO in case 4.
544
H. Zou et al. / Optics Communications 285 (2012) 539–545
4.5. Result analysis The performance of PSO and QPSO are assessed from the maximum error and root mean square error (RMS error). The simulation results of four cases are summarized in Table 5. The maximum errors between the reconstructed and origin strain profiles using PSO algorithm vary from 4.526 to 23.875 με, and the RMS errors vary from 3.883 to 15.986 με. The maximum errors between the reconstructed and origin strain profiles using QPSO algorithm vary from 2.672 to 11.207 με, and the RMS errors vary from 2.135 to 8.878 με. It can be seen that the more complex the strain distribution, the larger the errors that can be induced because the complicated strain distributions lead to complicated reflected spectra from the FBG, so the difficulty of finding optimized solutions is increased. The simulation results show that there is excellent agreement between the applied strain and the reconstructed strain profiles using PSO or QPSO. Furthermore, comparing the results using PSO with that of QPSO, the maximum errors and RMS errors using QPSO are smaller than that of PSO for each case. Because QPSO has better global search capability compared with PSO. In order to investigate the effect of noise on the accuracy of strain profile reconstruction, 10% and 20% of white Gaussian noise were added to the intensity reflection spectrum of quadratic strain profile distribution (case 2). The corresponding results of strain profile reconstruction are shown in Figs. 13,14 and Table 5. The results of Fig. 13 reveal that the addition of 10% Gaussian distribution noise is acceptable. However, when the noise increases to a 20% noise variation, Fig. 14 indicates the presence of some ripples in the reconstructed strain distribution, which are due to the noisy intensity reflection spectrum. The results of Table 5 indicate that the maximum errors using PSO algorithm vary from 22.748 to 51.549 με, RMS errors vary from 11.164 to 35.890 με; the maximum errors using QPSO algorithm vary from 8.158 to 29.055 με, RMS errors vary from 6.600 to 24.565 με. It is clear that as the noise variation increases, the maximum error and RMS error both increase. The noise free case 2 indicates extremely low maximum errors and RMS errors. However, in 10% and 20% noisy cases, the maximum errors and RMS errors are higher than that in the noise free case 2. Obviously, the maximum errors and RMS errors using QPSO are lower than that of PSO. The results presented in Table 5 confirm that the PSO and QPSO algorithm are also applicable to the measurement of noisy systems, but reveal that the quality of the reconstruction results depends on the degree of noise variation in the intensity spectrum. In the same degree of noise variation in the intensity spectrum, the QPSO algorithm is superior to the PSO algorithm. 5. Conclusions
Fig. 13. Original and reconstructed strain profiles under the noisy (10% variation) intensity spectrum in case 2.
reconstructed using QPSO or PSO algorithm even in regions close to localized damage where the strain fields can be highly non-linear and discontinuous. The proposed QPSO algorithm performed consistently well on the four cases, with better results than the PSO algorithm for these cases. In terms of convergence, the simulation results show that the QPSO converges to obtain solutions closer to the good solution and presents a small error than the PSO even in the condition where random noise exists. However, this technique only provides acceptable results for noise variations of 10% or less. The number of divided sections and grating length determine the spatial resolution of the FBG sensor and here it is 0.5 mm. By increasing the number of divided sections, the spatial resolution can be improved but the complexity of calculation will also be increased. The QPSO algorithm provides a simple and effective method to measurement strain distributions which is maybe helpful to the use of FBG sensors for structural damage monitoring and internal stress mapping of material.
Acknowledgements This work was financially supported by the National Natural Science Foundation of China under Project no. 60907038, the Natural Science Foundation of Jiangsu Province under Project no. BK2009370, the Postdoctoral Science Foundation of China under Project no. 20090461116 and the Post-doctoral Scientific Research Program of Jiangsu Province under Project no. 1001010B.
This work demonstrates a strain distribution profile reconstruction strategy for FBG sensor using the QPSO algorithm and transfer matrix approach. The effectiveness and feasibility of the QPSO algorithm was tested on simulated spectrum data for four strain distribution profiles and compared with that of the PSO algorithm. Simulation results show that the applied strain profiles can be effectively Table 5 Simulation results of strain profiles reconstruction. Case
Case Case Case Case Case Case
1 2 (noise free) 2 (10% noisy) 2 (20% noisy) 3 4
Maximum strain difference (με)
Maximum error (με)
RMS error (με)
PSO
QPSO
PSO
QPSO
1000 1000 1000 1000 1000 1000
4.526 6.042 22.748 51.549 15.818 23.875
2.672 3.129 8.158 29.055 7.734 11.207
3.883 5.057 11.164 35.890 6.613 15.986
2.135 2.411 6.600 24.565 6.353 8.878
Fig. 14. Original and reconstructed strain profiles under the noisy (20% variation) intensity spectrum in case 2.
H. Zou et al. / Optics Communications 285 (2012) 539–545
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
T. Erdogan, Journal of Lightwave Technology 15 (1997) 1277. D.S. Li, Z. Zhou, J.P. Ou, Measurement 44 (2011) 722. J. Frieden, J. Cugnoni, J. Botsis, et al., Composite Structures 92 (2010) 1905. Y.B. Dai, Y.J. Liu, J.S. Leng, et al., Optics and Lasers in Engineering 47 (2009) 1028. P.M. Lam, K.T. Lau, H.Y. Ling, et al., Optics and Lasers in Engineering 47 (2009) 1049. H.Y. Ling, K.T. Lau, W. Jin, et al., Optics Communications 270 (2007) 25. S. Huang, M. Ohn, M. LeBlanc, et al., Smart Materials and Structures 7 (1998) 248. K.A. Winick, J.E. Roman, IEEE Journal of Quantum Electronics 26 (1990) 1918. E. Peral, J. Capmany, J. Marti, IEEE Journal of Quantum Electronics 32 (1996) 2078. M.A. Muriel, J. Azana, A. Carballar, Optics Letters 23 (1998) 1526. J. Skaar, L. Wang, T. Erdogan, IEEE Journal of Quantum Electronics 37 (2001) 165. A. Gill, K. Peters, M. Studer, Measurement Science and Technology 15 (2004) 1877. P. Dong, J. Azana, A.G. Kirk, Optics Communications 228 (2003) 303.
545
[14] R.X. Zhang, S.J. Zheng, Y.J. Xia, Optics Communications 281 (2008) 3476. [15] C.Z. Shi, N. Zeng, M. Zhang, et al., Optics Communications 226 (2003) 167. [16] B.F. Yun, Y.P. Wang, A.M. Li, et al., Measurement Science and Technology 16 (2005) 2425. [17] J.F. Kennedy, R.C. Eberhart, Proceedings of IEEE International Conference on Neural Network (1995) 1942. [18] T.H. Sun, Expert Systems with Applications 36 (2008) 3428. [19] F. Van den Bergh, A.P. Engelbrecht, IEEE Transactions on Evolutionary Computation 8 (2004) 225. [20] J. Sun, B. Feng, W.B. Xu, IEEE Proceedings of Congress on Evolutionary Computation (2004) 325. [21] J. Sun, W.B. Xu, B. Feng, Cybernetics and Intelligent Systems Proceedings of the 2004 IEEE Conference (2004) 111. [22] M. Yamada, K. Sakuda, Applied Optics 26 (1987) 3474.
Contents lists available at SciVerse ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
Quantum-behaved particle swarm optimization algorithm for the reconstruction of fiber Bragg grating sensor strain profiles Hongbo Zou a, b,⁎, Dakai Liang a, Jie Zeng a, Lihang Feng a a b
State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, PR China College of Electric Engineering and Renewable Energy, Three Gorges University, Yichang 443002, PR China
a r t i c l e
i n f o
Article history: Received 27 June 2011 Received in revised form 17 November 2011 Accepted 21 November 2011 Available online 5 December 2011 Keywords: Fiber Bragg grating Particle swarm optimization Quantum-behaved particle swarm optimization Strain distribution
a b s t r a c t A quantum-behaved particle swarm optimization algorithm used for demodulating the strain profile along a fiber Bragg grating (FBG) from its reflection spectrum has been demonstrated. By combining the transfer matrix method for calculating the reflection spectrum of a FBG and the quantum-behaved particle swarm optimization technique, we develop a new technique for the reconstruction of FBG sensor strain profiles. This proposed technique is verified through numerical example reconstructions of FBG sensor simulated strain profile cases and numerical simulations reveal good agreement between the original and the reconstructed strain profiles. © 2011 Elsevier B.V. All rights reserved.
1. Introduction Fiber Bragg grating (FBG) sensors have all the advantages of optical fiber sensors such as small size, light weight and immunity to electromagnetic interference. Further, they possess flexible multiplexing capability because of the inherent wavelength-encoded strain or temperature response. Hence, FBG sensors have been considered excellent sensing elements for measuring strain or temperature. For a FBG sensor subjected to constant axial strain along the FBG, only the Bragg wavelength shift needs to be measured to gain the strain data because the Bragg wavelength shift is linearly proportional to the applied strain [1]. Many researches and applications have been developed based on this property [2–5]. But in practice, the axial strain applied to the FBG sensor may be non-constant such as linear, quadratic, even high-order or discontinuous. In this condition, the reflection spectrum shape of the FBG sensor is not only shifted but also distorted with multi-peaks [6]. Therefore, it is unsuitable for only measuring the Bragg wavelength shift to obtain the strain data. Fortunately, the reflection spectrum of a FBG sensor can be easily measured with a conventional optical spectrum analyzer (OSA). Thus, how to reconstruct the strain distribution along the FBG sensor from its reflection spectrum becomes a meaningful inverse problem. The simplest approach to solve this inverse problem is the intensity spectrum based (ISB) method [7], but it is only applicable to monotonic ⁎ Corresponding author at: State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, PR China. Tel.: + 86 25 84893466; fax: + 86 25 84892294. E-mail address: [email protected] (H. Zou). 0030-4018/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2011.11.069
strain distribution with increasing or decreasing profiles. Fourier transform approach [8] is also developed with the advantage of a short reconstruction time, but it is only suitable for weak FBG and requires both the intensity and phase information of reflection spectrum. Some other reconstruction methods have been reported, such as the iterative GLM method [9], the time-frequency signal analysis [10], and the layerpeeling algorithm [11] and so on. However, all these methods need both the intensity and phase information. Generally, the measurement of phase spectrum needs a more complex setup compared with that of intensity spectrum. Therefore, the intensity based methods is preferable to the methods which required the phase of reflection spectrum. Recently some evolutionary algorithms have been developed to solve this inverse problem, mainly including the genetic algorithm (GA) [12] and simulated annealing (SA) algorithm [13]. These methods only require intensity information, so they play a successful role in the distributed strain profile reconstruction. However, GA has the drawbacks of long running time and the requirements of optimizing the parameters by trials to skip the local best solutions. SA has the similar shortcomings compared with GA. Therefore, some modifications for GA and SA are reported such as chaos genetic algorithm (CGA) [14], adaptive simulated annealing (ASA) algorithm [15], simulated annealing evolutionary (SAE) algorithm [16]. The particle swarm optimization (PSO) algorithm is a member of a wider class of swarm intelligence methods used for solving global optimization problems. This method was originally proposed by Kennedy as a simulation of social behavior of bird flocking and fish schooling in 1995 [17]. Instead of using evolutionary operators to manipulate the individuals as in other evolutionary algorithms, PSO relies on the exchange of information between individuals. Each particle in
540
H. Zou et al. / Optics Communications 285 (2012) 539–545
PSO flies in search space with a velocity, which is dynamically adjusted according to its own former information. Compared with GA and other similar evolutionary algorithms, PSO has some attractive characteristics and in many cases proved to be more effective [18]. However, PSO is not a global optimization algorithm [19]. In [20–21], Sun et al. introduce quantum theory into PSO and propose a quantum-behaved PSO (QPSO) algorithm, which can be guaranteed theoretically to find optimal solution in search space. The experimental results on some widely used benchmark functions show that the QPSO works better than standard PSO [20–21] and should be a promising algorithm. In this work, QPSO algorithm was demonstrated to solve the inverse problem and compared with PSO algorithm. The basic principles are as follows. First, we use the transfer matrix method to calculate the reconstructed grating spectrums in different grating parameters. Then we use the QPSO algorithm to search a reconstructed spectrum that in agreement with the original spectrum. Thus from the parameters corresponding to this reconstructed spectrum, we can get the applied strain profile on the FBG. This paper is structured as follows. The transfer matrix method used to calculate the reflection spectrum of the FBG is presented in Section 2. The basic principle of QPSO is introduced in Section 3. The implementation of the reconstruction and the numerical simulation results are given in Section 4. Finally, the conclusions are provided in Section 5. 2. Transfer matrix method The reflection spectrum of a non-uniform FBG can be calculated by using transfer matrix method [22]. In this section, we summarized the most important mathematical equations of the transfer matrix method. This approach divides the grating uniformly into M small sections, each with uniform coupling properties. Define Ri and Si to be the field amplitudes after traversing the ith section, the propagation through this uniform section is described as:
Ri Si
¼ Fi
Ri−1 Si−1
ð1Þ
σ 6 coshðγΔzÞ−i γ sinhðγΔzÞ Fi ¼ 6 4 k i sinhðγΔzÞ γ
3 k −i sinhðγΔzÞ 7 γ 7 5 σ coshðγΔzÞ þ i sinhðγΔzÞ γ
R R ¼F 0 S S0 f 11 f 21
f 12 f 22
π νδneff λ
ð7Þ
σ¼
2π π neff þ δneff − λ Λ ðzÞ
ð8Þ
where neff is the mode effective index of refraction of the optical fiber, δn eff is the mean mode effective index of refraction variation, ν is the fringe visibility of the index variation, λ is the wavelength of FBG. In conclusion, for any given parameters L, Λ0, neff, δneff , ν, pe and ε(z), we can use the above equations to calculate the corresponding reflection spectrum. 3. Quantum-behaved particle swarm optimization algorithm PSO is very easy to be understood and implemented. In the PSO, the fitness function is a measurement of the distance between the calculated reflection intensity spectra and the target reflection intensity spectra. PSO has already been tried and tested in various standard optimization problems with excellent results [17]. However, the main disadvantage of PSO is that global convergence cannot be guaranteed [19]. To deal with this problem, concept of a global convergence guaranteed method called as quantum-behaved PSO (QPSO), was developed by Sun et al. in 2004 [20–21]. In the QPSO, the state of a particle is described by wave function Ψ(X, t), instead of position and velocity of the PSO. The dynamic behavior of the particle is widely divergent form that of the particle in the PSO systems in that the exact values of Xi and Vi cannot be determined simultaneously. We can only learn the probability of the particle's appearing in position Xi from probability density function |Ψ(X, t)| 2. The particles move according to the following Eqs: X i ðt þ 1Þ ¼ P i −β⋅jmbest i −X i ðt Þj⋅lnð1=uÞ if k≥0:5
ð9Þ
X i ðt þ 1Þ ¼ P i þ β⋅jmbest i −X i ðt Þj⋅lnð1=uÞ if kb0:5
ð10Þ
where ð2Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where γ ¼ k2 −σ 2 , k is the ‘ac’ coupling coefficient, σ is a general ‘dc’ self-coupling coefficient, Δz is the each section length. Thus, the total grating propagation mode can be expressed by:
k¼
2
F¼
where L is the length of FBG, Λ0 is the period of FBG before the ε(z) is applied, pe is the effective photo-elastic constant. Then k and σ in Eq. (2) can be expressed by:
ð3Þ ¼ F 1 F 2 ⋯F M
ð4Þ
ð5Þ
When a FBG is subject to an inhomogeneous axial strain ε(z), the equivalent period of the FBG Λ(z) [7] can be expressed as: Λ ðzÞ ¼ Λ 0 ½1 þ ð1−pe ÞεðzÞ0≤z≤L
mbest ¼
N 1X P N d¼1 g;d
ð6Þ
ð11Þ ð12Þ
In Eqs. (9)–(12), mbest is the mean best position defined as the mean of all the best positions of the population, g represents the index of the best particle among all the particles in the population, k, u and φ are random numbers distributed uniformly in [0,1], respectively. β, called contraction-expansion coefficient, is the only parameter in QPSO algorithm. It can be tuned to control the convergence speed of the algorithms. The β value is adaptively allocated as per the Eq. (13): β ¼ β max −ðβ max −β min Þ⋅t=t max
and the reflection coefficient of the FBG is given by: 2 f r ¼ 21 f 11
P i ¼ φ⋅P i;d þ ð1−φÞ⋅P g;d
ð13Þ
where βmax is the initial contraction-expansion factor value, βmin is the final contraction-expansion factor value, t is the current generation number and tmax is the maximum number of generations. The procedure for implementing the QPSO is presented by the following: Step 1. Initialization of swarm positions: Initialize a population of particles with random positions in the D-dimensional problem space using a same probability distribution function.
H. Zou et al. / Optics Communications 285 (2012) 539–545
Step 2. Evaluation of particle's fitness: Evaluate the fitness value of each particle. We define the fitness function by measuring the error between the reconstructed spectrum Rcal(λ) and the target spectrum Rtar(λ) in this work.
f ¼
H h X
2
Table 2 Strain profiles applied to the FBG sensor. Case
Strain profile (με) (z/mm)
Case1: Case2: Case3: Case4:
linear distribution quadratic distribution sine distribution discontinuous distribution
ð14Þ
Rcal ðλi Þ−Rtar ðλi Þ
541
100z + 10 5z2 + 50z + 10 500sin(0.2πz) 50z 0 ≤ z ≤ L/2 − 100z L/2 b z ≤ L
i¼1
Step 3. Step 4.
Step 5.
Step 6. Step 7.
where λi is the sampled wavelength of the FBG spectrum and the total number is H. Rcal(λ) and Rtar(λ) are closely related to the strain profile ε(z). The task of QPSO is to search the input D-dimensional space effectively to find a strain profile ε(z) along the grating axis that minimizes the fitness function. In this work, we suppose that the strain profile ε(z) is determined by a finite set of parameters αi ∈ [a, b](i = 1, . . . ,D), thus the fitness function optimization problem is converted to optimizing the D parameters. Updating of global point: Calculate the mbest using Eq. (12). Comparison to pbest (personal best): Compare each particle's fitness with the particle's pbest. If the current value is better than pbest, then set the pbest value equal to the current value and the pbest location equal to the current location in the D-dimensional space. Comparison to gbest (global best): Compare the fitness with the population's overall previous best. If the current value is better than gbest, then reset gbest to the current particle's array index and value. Updating of particles’ position: Change the position of the particles according to Eqs. (9) and (10). Repeating the evolutionary cycle: Loop to Step 2 until a stop criterion is met, a maximum number of generations adopted in this work.
4. Implementation and numerical simulation To test the proposed strain distribution reconstruction of FBG sensor, four simulated strain distribution cases are investigated: a linear distribution (case 1), a quadratic distribution (case 2), a sine distribution (case 3) and a discontinuous distribution (case 4). The parameters of the FBG sensor employed in this work are given in Table 1 and the four strain distribution cases are shown in Table 2. First the simulated spectra (as ‘measured’ spectra) for each case are calculated using the transfer matrix method, applying the specific strain profile to be reconstructed. Then we use the QPSO algorithm to reconstruct the applied strain profile. The QPSO parameters used for strain profiles reconstruction are listed in Table 3. Based on these standard parameters, the QPSO algorithm is suitable for reconstruction of different strain profiles with almost the same convergence quality. Furthermore, the performance of QPSO algorithm was demonstrated through each case and compared with that of the PSO algorithm. The PSO parameters in present work are shown in Table 4.
Table 3 Parameters of QPSO algorithm. Parameter
Value
Max contraction-expansion coefficient (βmax) Min contraction-expansion coefficient (βmin) Max generation (maxG) Population size (N) Number of grating segments (M) Number of spectral points (H) Rang of wavelength(λ)
1 0.5 100 20 20 201 2 nm
4.1. Linear strain distribution (case 1) The original and reconstructed reflection spectrums of the FBG sensor applied with a linear strain distribution (case 1) are shown in Fig. 1. Fig. 2 plots the original and reconstructed strain profiles. Fig. 3 presents the comparison of convergence processes of QPSO and PSO in the case 1. From Figs. 1 and 2, the reconstructed results of PSO and QPSO are extremely close to the original spectrum and applied strain distribution in case 1. From Fig. 3, it can be seen that though PSO shows faster convergence, better results are found with QPSO. The slow convergence speed corresponds to good global search ability, while fast speed results in good local search ability. This clearly shows that QPSO is successful in having better global search capability in comparison to PSO and hence a better optimal result is obtained using QPSO.
4.2. Quadratic strain distribution (case 2) A quadratic strain distribution (case 2) is applied to the FBG sensor; the original and reconstructed results are shown in Figs. 4 and 5. In Figs. 4 and 5, original spectrum and strain profile are plotted in dashed curve, while PSO reconstructed spectrum and strain profile are plotted in asterisk curve and QPSO reconstructed spectrum and strain profile are plotted in circle curve. Similarly, Fig. 6 presents the comparison of convergence processes of QPSO and PSO in case 2. In Fig. 6, convergence processes of PSO and QPSO are plotted in dashed and solid curve, respectively. From Figs. 4 and 5, the reconstructed results of PSO and QPSO are extremely close to the original spectrum and applied strain distribution in case 2. Fig. 6 shows that though PSO shows faster convergence, better results are found with QPSO.
Table 4 Parameters of PSO algorithm. Table 1 Parameters of FBG sensor. Parameter
Value
Length (L) Initial period (Λ0) Mode effective index of refraction (neff) Mean index variation (δneff ) Fringe visibility (ν) Effective strain-optic constant (pe)
10 mm 534.48 nm 1.45 2 × 10− 4 1 0.22
Parameter
Value
Inertia factor (w) Cognitive learning rate (c1) Social learning rate (c2) Max generation (maxG) Population size (N) Number of grating segments (M) Number of spectral points (H) Rang of wavelength(λ)
1 2 2 100 20 20 201 2 nm
542
H. Zou et al. / Optics Communications 285 (2012) 539–545
Fig. 1. Original and reconstructed reflected spectra applied to the FBG sensor in case 1.
Fig. 4. Original and reconstructed reflected spectra applied to the FBG sensor in case 2.
Fig. 2. Original and reconstructed strain profiles applied to the FBG sensor in case 1.
Fig. 5. Original and reconstructed strain profiles applied to the FBG sensor in case 2.
4.3. Sine strain distribution (case 3) A sine strain distribution (case 3) is applied to the FBG sensor. The simulation results are shown in Figs. 7, 8 and 9. In Figs. 7 and 8, the dashed curve represents the original spectrum and strain profile, the asterisk curve represents the PSO reconstructed spectrum and strain profile, and the circle curve represents the QPSO reconstructed
Fig. 3. Comparison of convergence process of QPSO and PSO in case 1.
spectrum and strain profile. In Fig. 9, the dashed and solid curves represent the convergence processes of PSO and QPSO, respectively. From Figs. 7 and 8, the reconstructed results of PSO and QPSO are very close to the original spectrum and applied strain distribution in case 3. Fig. 9 reveals that though PSO shows faster convergence, better results are found with QPSO.
Fig. 6. Comparison of convergence process of QPSO and PSO in case 2.
H. Zou et al. / Optics Communications 285 (2012) 539–545
543
Fig. 7. Original and reconstructed reflected spectra applied to the FBG sensor in case 3.
Fig. 10. Original and reconstructed reflected spectra applied to the FBG sensor in case 4.
Fig. 8. Original and reconstructed strain profiles applied to the FBG sensor in case 3.
Fig. 11. Original and reconstructed strain profiles applied to the FBG sensor in case 4.
4.4. Discontinuous strain distribution (case 4) A discontinuous strain distribution (case 4) is applied to the FBG sensor. The simulation results are shown in Figs. 10, 11 and 12. In Figs. 10 and 11, the dashed curve is the original spectrum and strain profile, the asterisk curve is the PSO reconstructed spectrum and strain profile, and the circle curve is the QPSO reconstructed spectrum
Fig. 9. Comparison of convergence process of QPSO and PSO in case 3.
and strain profile. In Fig. 12, the dashed and solid curves are the convergence processes of PSO and QPSO, respectively. From Figs. 10 and 11, the reconstructed results of PSO and QPSO are close to the original spectrum and applied strain distribution in case 4. From Fig. 12, it can be shown that though PSO shows faster convergence, better results are found with QPSO.
Fig. 12. Comparison of convergence process of QPSO and PSO in case 4.
544
H. Zou et al. / Optics Communications 285 (2012) 539–545
4.5. Result analysis The performance of PSO and QPSO are assessed from the maximum error and root mean square error (RMS error). The simulation results of four cases are summarized in Table 5. The maximum errors between the reconstructed and origin strain profiles using PSO algorithm vary from 4.526 to 23.875 με, and the RMS errors vary from 3.883 to 15.986 με. The maximum errors between the reconstructed and origin strain profiles using QPSO algorithm vary from 2.672 to 11.207 με, and the RMS errors vary from 2.135 to 8.878 με. It can be seen that the more complex the strain distribution, the larger the errors that can be induced because the complicated strain distributions lead to complicated reflected spectra from the FBG, so the difficulty of finding optimized solutions is increased. The simulation results show that there is excellent agreement between the applied strain and the reconstructed strain profiles using PSO or QPSO. Furthermore, comparing the results using PSO with that of QPSO, the maximum errors and RMS errors using QPSO are smaller than that of PSO for each case. Because QPSO has better global search capability compared with PSO. In order to investigate the effect of noise on the accuracy of strain profile reconstruction, 10% and 20% of white Gaussian noise were added to the intensity reflection spectrum of quadratic strain profile distribution (case 2). The corresponding results of strain profile reconstruction are shown in Figs. 13,14 and Table 5. The results of Fig. 13 reveal that the addition of 10% Gaussian distribution noise is acceptable. However, when the noise increases to a 20% noise variation, Fig. 14 indicates the presence of some ripples in the reconstructed strain distribution, which are due to the noisy intensity reflection spectrum. The results of Table 5 indicate that the maximum errors using PSO algorithm vary from 22.748 to 51.549 με, RMS errors vary from 11.164 to 35.890 με; the maximum errors using QPSO algorithm vary from 8.158 to 29.055 με, RMS errors vary from 6.600 to 24.565 με. It is clear that as the noise variation increases, the maximum error and RMS error both increase. The noise free case 2 indicates extremely low maximum errors and RMS errors. However, in 10% and 20% noisy cases, the maximum errors and RMS errors are higher than that in the noise free case 2. Obviously, the maximum errors and RMS errors using QPSO are lower than that of PSO. The results presented in Table 5 confirm that the PSO and QPSO algorithm are also applicable to the measurement of noisy systems, but reveal that the quality of the reconstruction results depends on the degree of noise variation in the intensity spectrum. In the same degree of noise variation in the intensity spectrum, the QPSO algorithm is superior to the PSO algorithm. 5. Conclusions
Fig. 13. Original and reconstructed strain profiles under the noisy (10% variation) intensity spectrum in case 2.
reconstructed using QPSO or PSO algorithm even in regions close to localized damage where the strain fields can be highly non-linear and discontinuous. The proposed QPSO algorithm performed consistently well on the four cases, with better results than the PSO algorithm for these cases. In terms of convergence, the simulation results show that the QPSO converges to obtain solutions closer to the good solution and presents a small error than the PSO even in the condition where random noise exists. However, this technique only provides acceptable results for noise variations of 10% or less. The number of divided sections and grating length determine the spatial resolution of the FBG sensor and here it is 0.5 mm. By increasing the number of divided sections, the spatial resolution can be improved but the complexity of calculation will also be increased. The QPSO algorithm provides a simple and effective method to measurement strain distributions which is maybe helpful to the use of FBG sensors for structural damage monitoring and internal stress mapping of material.
Acknowledgements This work was financially supported by the National Natural Science Foundation of China under Project no. 60907038, the Natural Science Foundation of Jiangsu Province under Project no. BK2009370, the Postdoctoral Science Foundation of China under Project no. 20090461116 and the Post-doctoral Scientific Research Program of Jiangsu Province under Project no. 1001010B.
This work demonstrates a strain distribution profile reconstruction strategy for FBG sensor using the QPSO algorithm and transfer matrix approach. The effectiveness and feasibility of the QPSO algorithm was tested on simulated spectrum data for four strain distribution profiles and compared with that of the PSO algorithm. Simulation results show that the applied strain profiles can be effectively Table 5 Simulation results of strain profiles reconstruction. Case
Case Case Case Case Case Case
1 2 (noise free) 2 (10% noisy) 2 (20% noisy) 3 4
Maximum strain difference (με)
Maximum error (με)
RMS error (με)
PSO
QPSO
PSO
QPSO
1000 1000 1000 1000 1000 1000
4.526 6.042 22.748 51.549 15.818 23.875
2.672 3.129 8.158 29.055 7.734 11.207
3.883 5.057 11.164 35.890 6.613 15.986
2.135 2.411 6.600 24.565 6.353 8.878
Fig. 14. Original and reconstructed strain profiles under the noisy (20% variation) intensity spectrum in case 2.
H. Zou et al. / Optics Communications 285 (2012) 539–545
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
T. Erdogan, Journal of Lightwave Technology 15 (1997) 1277. D.S. Li, Z. Zhou, J.P. Ou, Measurement 44 (2011) 722. J. Frieden, J. Cugnoni, J. Botsis, et al., Composite Structures 92 (2010) 1905. Y.B. Dai, Y.J. Liu, J.S. Leng, et al., Optics and Lasers in Engineering 47 (2009) 1028. P.M. Lam, K.T. Lau, H.Y. Ling, et al., Optics and Lasers in Engineering 47 (2009) 1049. H.Y. Ling, K.T. Lau, W. Jin, et al., Optics Communications 270 (2007) 25. S. Huang, M. Ohn, M. LeBlanc, et al., Smart Materials and Structures 7 (1998) 248. K.A. Winick, J.E. Roman, IEEE Journal of Quantum Electronics 26 (1990) 1918. E. Peral, J. Capmany, J. Marti, IEEE Journal of Quantum Electronics 32 (1996) 2078. M.A. Muriel, J. Azana, A. Carballar, Optics Letters 23 (1998) 1526. J. Skaar, L. Wang, T. Erdogan, IEEE Journal of Quantum Electronics 37 (2001) 165. A. Gill, K. Peters, M. Studer, Measurement Science and Technology 15 (2004) 1877. P. Dong, J. Azana, A.G. Kirk, Optics Communications 228 (2003) 303.
545
[14] R.X. Zhang, S.J. Zheng, Y.J. Xia, Optics Communications 281 (2008) 3476. [15] C.Z. Shi, N. Zeng, M. Zhang, et al., Optics Communications 226 (2003) 167. [16] B.F. Yun, Y.P. Wang, A.M. Li, et al., Measurement Science and Technology 16 (2005) 2425. [17] J.F. Kennedy, R.C. Eberhart, Proceedings of IEEE International Conference on Neural Network (1995) 1942. [18] T.H. Sun, Expert Systems with Applications 36 (2008) 3428. [19] F. Van den Bergh, A.P. Engelbrecht, IEEE Transactions on Evolutionary Computation 8 (2004) 225. [20] J. Sun, B. Feng, W.B. Xu, IEEE Proceedings of Congress on Evolutionary Computation (2004) 325. [21] J. Sun, W.B. Xu, B. Feng, Cybernetics and Intelligent Systems Proceedings of the 2004 IEEE Conference (2004) 111. [22] M. Yamada, K. Sakuda, Applied Optics 26 (1987) 3474.