- Email: [email protected]
Optimal sizing of a solar thermal building installation using particle swarm optimization
Energy 41 (2012) 31e37
Contents lists available at ScienceDirect
Energy journal homepage: www.elsevier.com/locate/energy
Optimal sizing of a solar thermal building installation using particle swarm optimization Raffaele Bornatico a, *, Michael Pfeiffer a, b, Andreas Witzig b, Lino Guzzella a a b
Institute for Dynamic Systems and Control ETH Zurich, Sonneggstrasse 3, CH 8092 Zurich, Switzerland Vela Solaris AG, CH 8400 Winterthur, Switzerland
a r t i c l e i n f o
a b s t r a c t
Article history: Received 8 October 2010 Received in revised form 6 May 2011 Accepted 15 May 2011 Available online 30 June 2011
In recent years the domestic energy management has become a non-trivial task as the number of energy sources and system components involved have increased, and all components have to operate coordinately in order to maximize global efficiency measures. In this paper a methodology is presented for finding the optimal size of the main components for a solar thermal system where particular attention is given to the optimization framework. The use of the PSO (particle swarm optimization) algorithm is proposed and the results obtained are compared with a GA (genetic algorithm) solution. Further, the relative influence of certain system parameters on the optimal configuration is investigated by means of a sensitivity analysis where the size of the collector is shown to have the greatest influence on all main output quantities while the size of the auxiliary power unit presents a relatively small influence on the solution. Finally, it is demonstrated that the accurate sizing of the energy components is necessary to minimize the energy consumption and cost of installation, while maximizing the solar fraction. The proposed methodology is shown to successfully solve the problem 2011 Elsevier Ltd. All rights reserved.
Keywords: Particle swarm optimization Multi-objective optimization Polysun Solar combisystem Solar energy Smart building
1. Introduction Solar energy for SH (space heating) and the production of DHW (domestic hot water) has become an important factor in reducing the global CO2 emissions. Within the previous and the last United Nations climate change conference it has been analyzed what needs to be done to limit the long-term concentration of greenhouse gases in the atmosphere to 450 ppm of CO2 equivalent, in line with a 2 C increase in global temperature by 2100. Solar energy is considered as one of the most promising candidates to tackle our dependency on, and use of fossil fuels, and thus for reducing the related emissions. The IEA (International Energy Agency) with its “Task 26” on solar combisystem emphasizes that if the direct use of solar energy is to make a significant contribution to the heat supply, it is necessary that solar-heating technologies must be developed and widely applied over and beyond the sole field of DHW preparation [1,2]. Accordingly, the focus of this research has been set on a solar combisystem installation that simultaneously fulfills DHW and space heating needs. Previous publications on solar combisystems have presented analyses of
* Corresponding author. Tel.: þ41 44 632 2453; fax: þ41 44 632 1139. E-mail address: [email protected] (R. Bornatico). 0360-5442/$ e see front matter 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.energy.2011.05.026
the effects of geographic location [3] and the use of different DHW load profiles [4] on solar fraction and energy consumption of the system. This paper considers a mid-sized single-family house located in Zurich, Switzerland, and presents a methodology for finding the optimal component sizes for a solar combisystem. All simulations are carried out with Vela Solaris Polysun, which is a wellestablished software tool in the field of planning and optimization of building energy systems. The simulation kernel is based on a plug-flow simulation of the thermal system [5], and it uses statistical meteorological data as an input [6,7]. The optimization routine was implemented in MATLAB. The interface from MATLAB to the simulation kernel is called Polysun Inside. Polysun Inside is a Polysun plug-in that allows all simulation functionality to be controlled from the MATLAB environment in an automatic-iterative evaluation routine. In the specific case, a parallelization over 10 CPUs has been performed which considerably decreased the optimization time. This paper is structured as follows: In Section 2 the modeling assumptions and the simulation setup are presented. Section 3 describes the optimization framework and all parameters necessary for the optimization routine to succeed. Section 4 describes the main features of the PSO algorithm and introduces the subsequent results sections.
32
R. Bornatico et al. / Energy 41 (2012) 31e37
2. Solar combisystem model in Polysun The simulation software Polysun offers a broad range of functionalities required for the analysis and design of domestic energy systems. In the Polysun catalogs, a variety of components are available with all characteristic data and efficiency maps necessary for the hydronic and thermal analysis of such systems. Analysis and design of photovoltaic systems, solar cooling, and combinations of solar thermal and heat pump systems are also possible. Polysun calculates all relevant system parameters related to the production of heat and electricity. It also comprises the calculations for system amortization and the data required for subsidy applications. The variable step solver adapts the simulation step size to capture the transient effects down to a minimum of 4 min. It follows that sufficient accuracy of the system’s simulated characteristics and in particular of its dominant transients can be assumed. The user friendly interface of Polysun allows for the easy parameterization of the system, while the software itself is capable to output and store practically all relevant physical quantity in convenient data formats. Furthermore a mean simulation time of 1 min for a one-year simulation makes the software Polysun an ideal platform to perform the mentioned analyses [7]. The setup of the solar combisystem used in this research is depicted in Fig. 1. Leftmost is the APU (auxiliary power unit) and rightmost is the solar collector. Both components are connected to the storage tank in the center of the picture. The loads are characterized by a DHW and a space heating demand. The DHW profile used divides 200 l/d of water at 45 C into the periodic daily demand depicted in Fig. 2. The effects of using a realistic profile as opposed to a representative periodic profile are discussed in detail in Ref. [4], where it is shown that the heat demand difference between the two approaches is 0.1% and the fractional energy saving difference is just around 1.5% when using a flow regulation device. The greatest influence is shown to be caused by the reduced demand during the summer holidays in the realistic profile. Considering that with Polysun absences can be included in the simulations, it is reasonable to assume that the profile used hereafter approximates sufficiently well a real DHW load. Energy for space heating is less straightforward to determine since it depends on current ambient temperature and building
Fig. 2. Daily domestic hot water heat demand.
insulation. The heating setpoint temperature is 20 C during the day and 18 C at night. Other simulation relevant parameters of the building include the U-value of 0.45; the heated area of 150 m2; the air change rate of 0.3 h1; the air infiltration of 0.6 h1; the internal heat gains of 2 W m2; the heat gains from equipment of 240 W; the heat capacity of 500 kJ K1 m2 and the solar heat gain coefficient of 0.8. The collector is oriented toward South (0 ) with an inclination angle of 45 . The resulting simulated heat demand is shown in Fig. 3 and introduces the importance of the availability of comprehensive meteorological data.
2.1. Statistical meteorological data To a large extent, the accuracy of solar system simulations for a given location depends on the availability of realistic data of solar irradiation, humidity, etc. In Polysun, meteorological data are provided by the Meteonorm database [8] containing data based on measurements from 8055 weather stations worldwide. For any given location, the data of the closest weather stations are interpolated. The generation of yearly series of weather data utilizes stochastic models, where stored monthly mean values and Markov Transition Matrices are used to generate hourly weather data [9,10]. The resulting hourly data have the same statistical properties as the measured data (i.e. average value, variance, autocorrelation), and thus represent an accurate approximation. In Fig. 4, the outdoor temperature and global irradiance time series are shown over a one-year period.
Fig. 1. A solar combisystem designed using Polysun.
R. Bornatico et al. / Energy 41 (2012) 31e37
33
Fig. 3. Yearly heat requirements for space heating in a mid-sized single-family house.
3. Optimization framework The main components of the system as presented in the previous section are the collector, the tank, and the APU. In order for these components to be optimally dimensioned, the corresponding sizes have to populate the input parameter set that is fed to the simulation kernel. In Table 1 the input set is listed together with the corresponding minimal and maximal values used throughout this work. The extreme values for the output set shown in Table 2 do not represent a given constraint. Rather, these values have been obtained from simulations throughout this work. They are to show that the investigated range for these values is significant. The goal of the optimization routine consists in finding the parameters wopt that minimize a general function f : Rn /R that is formulated in the canonical form as follows:
min f ðwÞ w
subject to gðwÞ < 0; hðwÞ ¼ 0
(1)
The functions g and h include all constraints to the problem defining the set Q of the acceptable w˛Q. These constraints include the boundary values of w and two physical constraints on the system. First, for a given w to be acceptable, the heat requirements for DHW and space heating must be satisfied. Secondly, the temperature of the fluid flowing in the collector must not exceed the 100 C threshold for longer than 1% of the simulation time. This last condition limits the stagnation time to reasonable values, thereby ensuring the safe operation and a long life of the collector. In the presented problem a numeric optimization is necessary because an algebraic solution to the problem does not exist. Furthermore it is clear that the topographical properties of the function f play a crucial role in how the optimization algorithm converges from starting values w0 toward wopt , and also that the choice of the right cost function is determinant for the successful
convergence of the algorithm to the point of interest. Accordingly, all relevant quantities should appear in the cost function. In the case of the optimization of the combisystem, the objectives of the authors are to maximize the solar fraction, to minimize the total energy use, and to minimize the additional cost of the installation.
8 9 < X fj ðwÞ = f ðwÞ : ¼ yry ¼ wj $ : fjmax ; j
(2)
Therefore, the rather intuitive linear form (2) is chosen where the solar fraction f1 ðwÞ ¼ SFnðwÞ, the total energy use of the system f2 ðwÞ ¼ EtotðwÞ and the cost of the installation f3 ðwÞ ¼ CostðwÞ depend on the sizes of the components w. Since a scalar value of f is required by the minimization algorithm, the weighted sum, with weights wi of the normalized fi, yields a dimensionless, normalized cost function f ðwÞ. Note that in order to maximize the solar fraction by a minimization of f, w1 must be negative. An important result from a user’s point of view is the pricing of the components. This term is based on a linear pricing assumption for each component and on current oil prices. This choice is shown in Ref. [11] to produce meaningful results. The relevant values used for the installation pricing are listed in Table 3. 4. The particle swarm optimization algorithm The analogy of a population of individuals, called particles (e.g. a swarm of bees) having the common goal of finding an optimal position (e.g. the best flower in a field) is useful to picture how the algorithm converges toward a solution. This socially inspired optimization technique was first introduced by Kennedy and Eberhart [12] as a promising tool that does not requires the computation of derivatives.
Fig. 4. Meteorological data for Zurich, Switzerland. Solid: outdoor temperature, dotted: global irradiance.
34
R. Bornatico et al. / Energy 41 (2012) 31e37
Table 1 Input parameter set. Item
Symbol
Min
Max
Unit
Collector area Tank volume APU size
w1 w2 w3
1 100 5
40 3000 50
m2 l kW
Table 2 Output parameter set. Item
Symbol
Min
Max
Unit
Solar fraction Energy use Cost
f1 f2 f3
2.8% 12 988 12 988
42.7% 45 407 46 226
e kW h V
Table 3 Relevant factors for the installation pricing. Item
Value
Unit
Collector Collector Storage tank APU Heating oil Heating oil Maintenance
100 360 4.22 50 10.5 0.85 100
V V/m2 V/l V/kW kW h/l V/l V/year
Each particle i updates its position toward better fitness areas in the problem domain landscape based on a self cognitive term expressed by the best position (Pi) found so far by the particle, and the global best position encountered so far by the entire swarm (G). Starting with random initial velocities v0i and random initial positions p0i inside a closed set, the basic PSO algorithm computes the velocity and position of the i-th particle, i ¼ 1; 2; .; n for every iteration k ¼ 0; 1; 2; .kmax in the following way:
vkþ1 ¼ vki þ g1i Pi pki þ g2i G pki i
(3)
¼ pki þ vkþ1 pkþ1 i i
(4)
The random numbers g1;2 ˛½0; 1 influence the magnitude of the two vectors ðPi pki Þ and ðG pki Þ, thus influencing indirectly the . Fig. 5 depicts how the update from magnitude and direction of vkþ1 i k to kþ1 is performed for a simplified two dimensional problem with no randomness g1,2 ¼ 1. Note that every dimension represents a parameter to be optimized.
Fig. 5. Simplified PSO iteration in two dimensions.
The common PSO algorithm includes three additional parameters, namely an inertia function fk and two acceleration constants a1,2. The inertia function used hereafter is linearly decreasing with respect to iterations, thus reducing the influence of past velocities, thereby enabling the algorithm to adapt to small regions as the optimization converges. This adaptation capability is of key relevance since it makes the use of hybrid methods unnecessary. The resulting velocity equation, together with the position update, represents the core of the PSO algorithm [13,14].
h i h i þ a2 g2i G pki vkþ1 ¼ fk vki þ a1 g1i Pi pki i
(5)
pkþ1 ¼ pki þ vkþ1 i i
(6)
Further extensions of the algorithm can include natural selection considerations such as a varying population size where bad particles are killed and new particles are generated in the vicinity of the global best point. Dynamic adjustment of swarm parameters have been also proposed as a way to avoid local optima [15e17]. The work of Clerc and Kennedy [18] describes the importance of including constraints on the velocity vector. Their work shows that if no constraints are set the algorithm can become unstable. Note that all these extensions increase the computational complexity of the algorithm and should therefore be avoided, if not necessary. In the problem presented here a fixed population of 10 particles has been deemed to be sufficiently effective. Further, a linearly decreasing inertia function already accounts for stability, as the introduction of velocity limits can be shown to be a loose constraint. 5. Optimization results The results reported in this section have been based on simulation results of the software Polysun when simulating the model presented in Section 2. The converged results of the minimization of Eq. (2) by means of the PSO algorithm are presented in Section 5.1, while results obtained using a genetic algorithm follow in Section 5.2. A further investigation on the solution is presented in the form of a sensitivity analysis, which follows in the subsequent Section 6. 5.1. Particle swarm optimization: results The minimization of the fitness function by means of the PSO algorithm is shown in Fig. 6. The fitness function is plotted against particle generations where each generation consists of 10 particles. Note that points which do not satisfy the optimization constraints have been omitted. The PSO algorithm converges toward a minimal fitness value of f opt; PSO ¼ 1:278 and the interpretation of this value, with respect to its physical terms, is presented in Table 4. The initial population was randomly generated inside the boundary set
Fig. 6. Fitness function minimization with PSO.
R. Bornatico et al. / Energy 41 (2012) 31e37 Table 4 PSO optimization of fitness function terms.
35
Table 5 GA optimization of fitness function terms.
Term
Starting value
Optimized value
Change
Term
Starting value
Optimized value
Change
f SFn Etot Cost
1.435 18.55% 16 428 kW h 13 210 V
1.278 21.8% 15 806 kW h 8 963 V
10.1% þ17.5% 3.8% 32.1%
f SFn Etot Cost
1.335 20.6% 16 109 kW h 10 952 V
1.278 22.15% 15 750 kW h 9 230 V
4.2% þ7.5% 2.2% 15.7%
Q. It follows that the starting values as reported in Table 4, are averaged values of the feasible points of the initial population. It is important to judge the convergence of the algorithm by considering the spread of the particles. Fig. 7 illustrates the convergence of all particles, for all parameters considered in this study. From Figs. 6 and 7 it can be seen that after 50 generations the algorithm has converged to the vicinity of the optimal values which are found to be wopt;PSO ¼ f14:5 m2 ; 498:98 l; 8:5 kWg. 5.2. Genetic algorithm optimization: results On the same problem a standard GA has been applied where, if applicable, the same optimization parameters as in the PSO were used, specifically, a population of 10 particles and the same parameter set boundaries. In this case also, the fitness function decreases toward a minimal value of f opt; GA ¼ 1:278, a value that is consistent with the result of the PSO optimization. The comparison of all other optimal fitness function terms between PSO (Table 4) and GA (Table 5) shows that the two solutions differ by a minimal amount only. Also for the GA, the convergence of the fitness function and the system parameters is shown in Figs. 8 and 9. Fig. 8 shows that the convergence of all parameters over generation number is comparable to those obtained with the PSO optimization as the algorithm converges toward an optimal value of wopt; GA ¼ f15:1 m2 ; 514 l; 8:6 kWg. A slight difference in the optimal parameters obtained for the collector area and the tank volume can be noticed. However, this does not influence significantly the fitness function values, thus proving that both algorithms converge to a similar solution. Given
a
b
that the search space is rather big and that the starting individuals are randomly drawn inside this space, the fact that the optimal values found with GA and PSO are very similar supports the consistency of the results and the solution found in this section. 6. Parameter sensitivity analysis The results in Section 5 are specific for the considered system, i.e. a daily DHW consumption of 200 l/day at 45 C and an SH setpoint temperature of 20 and 18 C during day and night respectively. First the sensitivity of the results with respect to a variation of these parameters is explored and secondly, the same analysis is performed over the optimization parameters w1;2;3 . All results in this section are produced using the optimal PSO values opt;PSO opt;PSO and fj as nominal values when necessary. wi 6.1. Operational parameter variation Thus far, constant operational parameters have been used. Among them, the amount of DHW consumed (p1) and its temperature (p2) together with the heating setpoint temperature of the building during day (p3) and night (p4) are investigated here. Initially the DHW volume p1 is varied between {50, 100, 200, 300, 400} l/day and the simulation outcomes are linearly fitted to produce the results shown in the first column of Table 6. For instance the ratio ðvSFnÞ=ðvp1 Þ is 0.01 [%/l] indicating a solar fraction increment of 1% for every additional 100 l/day DHW consumed on top of the nominal 200 l/day. Similarly, the DHW temperature (p2) is varied between {35, 40, 45, 50, 55} C, p3 ˛ {18, 19, 20, 21, 22} C and p4 ˛ {16, 17, 18, 19, 20} C to produce the other values in Table 6. The sensitivity of the results relative to the corresponding nominal values is also a useful mean of comparison. The relative sensitivity computed using Eq. (7) is graphically shown in Fig. 10 where pnom ¼ {200 l/day, 45 C, 20 C, 18 C}. The relative variation of the daily heating setpoint temperature clearly dominates all other effects and indicates that this parameter is of crucial importance in the characterization of the system. From this figure it can be read that a 1% increase in p3 causes an additional energy consumption of around 2%.
Sij ¼
v fj =fjopt;PSO vðpi Þ=pnom i
c
Fig. 7. PSO convergence of system parameters towards 14.5 m2 (a), 499 l (b), and 8.5 kW (c).
Fig. 8. Fitness function minimization with GA.
(7)
36
R. Bornatico et al. / Energy 41 (2012) 31e37 Table 7 Absolute sensitivity of optimization parameters.
a
SFn [%] Etot [kW h] Cost [V]
w1 ½m2
w2 ½l
w3 ½kW
0.53 94.4 352.9
0.01 0.78 4.2
0.3 119.30 60.6
b
c
Fig. 9. GA convergence of system parameters towards 15.1 m2 (a), 514 l (b), and 8.6 kW (c).
Fig. 11. Relative sensitivity of the optimal solar fraction (1), total energy consumption (2) and cost of installation (3) related to a variation of collector size w1 , tank volume w2 , and rated power of the auxiliary power unit w3 .
Sij ¼ Table 6 Absolute sensitivity of operational parameters.
SFn [%] Etot [kW h] Cost [V]
p1 [l/day]
p2 [ C]
p3 [ C]
p4 [ C]
0.01 9.6 0.73
0.02 70.6 5.16
1.14 1638.5 125.66
0.24 91.3 7.69
6.2. Optimization parameter variation In contrast to the previous section where operational parameters are varied for a given installation, this section deals with varying size of collector, tank volume and APU power. Each one of these parameters is varied by {10%, 5%, þ5%, þ10%} of the respective nominal values, which correspond to the optimal PSO results wopt; PSO ¼ f14:5 m2 ; 498:98 l; 8:53 kWg. The same methodology as in the previous section is applied here and the results are presented in absolute values in Table 7 and in Fig. 11 according to Eq. (8).
Fig. 10. Relative sensitivity of the optimal solar fraction (1), total energy consumption (2) and cost of installation (3) related to a variation of DHW volume p1, DHW temperature p2, SH day setpoint temperature p3 and SH night setpoint temperature p4.
opt;PSO v fj =fj opt;PSO
vðwi Þ=wi
(8)
The relative sensitivity of the solar fraction Si1 is depicted in Fig. 11 by the bars (1) where it can be seen that an increase of the collector or tank size has a positive effect on the solar fraction. This effect is dominated in relative terms, by the solar collector with S11 ¼ 0.35. Further, bar (2) of w3 shows that an increase in APU size in the vicinity of the optimal solution causes a greater energy consumption of the system due to a shifted operating point in the APU consumption and efficiency map, whereas an increase of collector size or tank volume reduce the overall energy consumption of the system. The bars labeled (3) in Fig. 11 show that an increase in collector size has the greatest relative influence on the cost of the installation. The tank volume has a smaller influence, while the influence of the APU size is relatively low. 7. Conclusions and discussion In this paper the optimal sizing of a solar thermal system has been presented. The analysis is performed on a solar combisystem for a mid-sized single-family house in Zurich, Switzerland. While the optimization framework is in principle independent on the optimization algorithm, a detailed analysis has been carried out on the performance of the Particle Swarm Optimization algorithm when it is applied for solving this problem. The results are comparable to those obtained with the more common Genetic Algorithm. When the implementation efforts and the computational power demand are considered as well, the PSO is a slightly better choice for solving the presented problem. A collector size of 14.5 m2 together with a tank volume of 498.98 l and an APU nominal power of 8.5 kW, are the optimal sizes for the main system’s components, which lead to a solar fraction of 21.8%, a total energy use of 15 806 kW h and a cost of the installation of 8 963 V. The parameter sensitivity analysis shows that the size of the collector has the greatest influence on the solar fraction, energy use and the installation cost, while the tank volume influence is significant on the solar fraction and cost of installation. The size of the APU has a relatively small effect on the both energy use and installation cost. The variation of selected operational parameters
R. Bornatico et al. / Energy 41 (2012) 31e37
has been investigated as well, and shows that the sensitivity of the daily heating setpoint temperature of the building is by far the dominant parameter when compared to all parameters considered in this study. Future work includes an extension and further investigation of the parameter set. Furthermore, the proposed methodology is to be applied over other solar installations and/or more complex models including, for example a detailed life cycle analysis of the system. Acknowledgments This work has been partly funded by the Swiss Innovation Promotion Agency CTI, whose support is gratefully acknowledged. References [1] Sawin E, Jones A, Sterman J. Final Copenhagen accord press release 19 december 09-expanded version; 2009. [2] Suter JM, Letz T, Weiss W, Inbnit J. Solar combisystems in Austria, Denmark, Finland, France, Germany, Sweden, Switzerland, the Netherlands and the USA overview 2000; 2000. [3] Lund P. Sizing and applicability considerations of solar combisystems. Solar Energy 2005;78(1):59e71. [4] Jordan U, Vajen K. Influence of the dhw load profile on the fractional energy savings: a case study of a solar combi-system with trnsys simulations. Solar Energy 2001;69(Suppl. 6):197e208. [5] Klein S, Beckmann B, Duffie J. Trnsys, a transient system simulation program, users manual, version 16; 2006.
37
[6] Duffie J, Beckman W. Solar engineering of thermal processes. 3rd ed. Wiley Interscience; 2006. [7] Witzig A, Foradini F, Probst MCM, Roecker C. Simulation tool for architects. In: Proceedings of the 2009 CISBAT international conference; 2009. [8] Meteotest. Meteonorm database. Bern, Switzerland: Meteotest; 2000. [9] Meteotest. Meteonorm handbook, part 2: theory. Bern, Switzerland: Meteotest; 2000. [10] Aguiar RJ, Collares-Pereira M, Conde JP. Simple procedure for generating sequences of daily radiation values using a library of Markov transition matrices. Solar Energy 1988;40(3):269e79. [11] Fraisse G, Bai Y, Pierrs NL, Letz T. Comparative study of various optimization criteria for sdhws and a suggestion for a new global evaluation. Solar Energy 2009;83(2):232e45. [12] Kennedy J, Eberhart R. Particle swarm optimization. In: Proceedings of the 1995 IEEE International Conference on Neural Networks, vol. 4; 1995. p. 1942e1948. [13] Birge B. Psot e a particle swarm optimization toolbox for use with matlab. In: Proceedings of the 2003 IEEE Swarm Intelligence Symposium SIS 03; 2003. p. 182e186. [14] Bratton D, Kennedy J. Defining a standard for particle swarm optimization; 2007. 120e127. [15] Dimopoulos GG, Frangopoulos CA. Optimization of energy systems based on evolutionary and social metaphors. Energy 2008;33(2):171e179. 19th International Conference on Efficiency, Cost, Optimization, Simulation and Environmental Impact of Energy Systems e ECOS 2006. [16] Shi Y, Eberhart R. A modified particle swarm optimizer. In: Proceedings of the 1998 IEEE World Congress on Computational Intelligence; 1998. p. 69e73. [17] Fourie P, Groenwold A. The particle swarm optimization algorithm in size and shape optimization. Structural and Multidisciplinary Optimization 2002;23: 259e67. [18] Clerc M, Kennedy J. The particle swarm e explosion, stability, and convergence in a multidimensional complex space. IEEE Transactions on Evolutionary Computation 2002;6(1):58e73.
Contents lists available at ScienceDirect
Energy journal homepage: www.elsevier.com/locate/energy
Optimal sizing of a solar thermal building installation using particle swarm optimization Raffaele Bornatico a, *, Michael Pfeiffer a, b, Andreas Witzig b, Lino Guzzella a a b
Institute for Dynamic Systems and Control ETH Zurich, Sonneggstrasse 3, CH 8092 Zurich, Switzerland Vela Solaris AG, CH 8400 Winterthur, Switzerland
a r t i c l e i n f o
a b s t r a c t
Article history: Received 8 October 2010 Received in revised form 6 May 2011 Accepted 15 May 2011 Available online 30 June 2011
In recent years the domestic energy management has become a non-trivial task as the number of energy sources and system components involved have increased, and all components have to operate coordinately in order to maximize global efficiency measures. In this paper a methodology is presented for finding the optimal size of the main components for a solar thermal system where particular attention is given to the optimization framework. The use of the PSO (particle swarm optimization) algorithm is proposed and the results obtained are compared with a GA (genetic algorithm) solution. Further, the relative influence of certain system parameters on the optimal configuration is investigated by means of a sensitivity analysis where the size of the collector is shown to have the greatest influence on all main output quantities while the size of the auxiliary power unit presents a relatively small influence on the solution. Finally, it is demonstrated that the accurate sizing of the energy components is necessary to minimize the energy consumption and cost of installation, while maximizing the solar fraction. The proposed methodology is shown to successfully solve the problem 2011 Elsevier Ltd. All rights reserved.
Keywords: Particle swarm optimization Multi-objective optimization Polysun Solar combisystem Solar energy Smart building
1. Introduction Solar energy for SH (space heating) and the production of DHW (domestic hot water) has become an important factor in reducing the global CO2 emissions. Within the previous and the last United Nations climate change conference it has been analyzed what needs to be done to limit the long-term concentration of greenhouse gases in the atmosphere to 450 ppm of CO2 equivalent, in line with a 2 C increase in global temperature by 2100. Solar energy is considered as one of the most promising candidates to tackle our dependency on, and use of fossil fuels, and thus for reducing the related emissions. The IEA (International Energy Agency) with its “Task 26” on solar combisystem emphasizes that if the direct use of solar energy is to make a significant contribution to the heat supply, it is necessary that solar-heating technologies must be developed and widely applied over and beyond the sole field of DHW preparation [1,2]. Accordingly, the focus of this research has been set on a solar combisystem installation that simultaneously fulfills DHW and space heating needs. Previous publications on solar combisystems have presented analyses of
* Corresponding author. Tel.: þ41 44 632 2453; fax: þ41 44 632 1139. E-mail address: [email protected] (R. Bornatico). 0360-5442/$ e see front matter 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.energy.2011.05.026
the effects of geographic location [3] and the use of different DHW load profiles [4] on solar fraction and energy consumption of the system. This paper considers a mid-sized single-family house located in Zurich, Switzerland, and presents a methodology for finding the optimal component sizes for a solar combisystem. All simulations are carried out with Vela Solaris Polysun, which is a wellestablished software tool in the field of planning and optimization of building energy systems. The simulation kernel is based on a plug-flow simulation of the thermal system [5], and it uses statistical meteorological data as an input [6,7]. The optimization routine was implemented in MATLAB. The interface from MATLAB to the simulation kernel is called Polysun Inside. Polysun Inside is a Polysun plug-in that allows all simulation functionality to be controlled from the MATLAB environment in an automatic-iterative evaluation routine. In the specific case, a parallelization over 10 CPUs has been performed which considerably decreased the optimization time. This paper is structured as follows: In Section 2 the modeling assumptions and the simulation setup are presented. Section 3 describes the optimization framework and all parameters necessary for the optimization routine to succeed. Section 4 describes the main features of the PSO algorithm and introduces the subsequent results sections.
32
R. Bornatico et al. / Energy 41 (2012) 31e37
2. Solar combisystem model in Polysun The simulation software Polysun offers a broad range of functionalities required for the analysis and design of domestic energy systems. In the Polysun catalogs, a variety of components are available with all characteristic data and efficiency maps necessary for the hydronic and thermal analysis of such systems. Analysis and design of photovoltaic systems, solar cooling, and combinations of solar thermal and heat pump systems are also possible. Polysun calculates all relevant system parameters related to the production of heat and electricity. It also comprises the calculations for system amortization and the data required for subsidy applications. The variable step solver adapts the simulation step size to capture the transient effects down to a minimum of 4 min. It follows that sufficient accuracy of the system’s simulated characteristics and in particular of its dominant transients can be assumed. The user friendly interface of Polysun allows for the easy parameterization of the system, while the software itself is capable to output and store practically all relevant physical quantity in convenient data formats. Furthermore a mean simulation time of 1 min for a one-year simulation makes the software Polysun an ideal platform to perform the mentioned analyses [7]. The setup of the solar combisystem used in this research is depicted in Fig. 1. Leftmost is the APU (auxiliary power unit) and rightmost is the solar collector. Both components are connected to the storage tank in the center of the picture. The loads are characterized by a DHW and a space heating demand. The DHW profile used divides 200 l/d of water at 45 C into the periodic daily demand depicted in Fig. 2. The effects of using a realistic profile as opposed to a representative periodic profile are discussed in detail in Ref. [4], where it is shown that the heat demand difference between the two approaches is 0.1% and the fractional energy saving difference is just around 1.5% when using a flow regulation device. The greatest influence is shown to be caused by the reduced demand during the summer holidays in the realistic profile. Considering that with Polysun absences can be included in the simulations, it is reasonable to assume that the profile used hereafter approximates sufficiently well a real DHW load. Energy for space heating is less straightforward to determine since it depends on current ambient temperature and building
Fig. 2. Daily domestic hot water heat demand.
insulation. The heating setpoint temperature is 20 C during the day and 18 C at night. Other simulation relevant parameters of the building include the U-value of 0.45; the heated area of 150 m2; the air change rate of 0.3 h1; the air infiltration of 0.6 h1; the internal heat gains of 2 W m2; the heat gains from equipment of 240 W; the heat capacity of 500 kJ K1 m2 and the solar heat gain coefficient of 0.8. The collector is oriented toward South (0 ) with an inclination angle of 45 . The resulting simulated heat demand is shown in Fig. 3 and introduces the importance of the availability of comprehensive meteorological data.
2.1. Statistical meteorological data To a large extent, the accuracy of solar system simulations for a given location depends on the availability of realistic data of solar irradiation, humidity, etc. In Polysun, meteorological data are provided by the Meteonorm database [8] containing data based on measurements from 8055 weather stations worldwide. For any given location, the data of the closest weather stations are interpolated. The generation of yearly series of weather data utilizes stochastic models, where stored monthly mean values and Markov Transition Matrices are used to generate hourly weather data [9,10]. The resulting hourly data have the same statistical properties as the measured data (i.e. average value, variance, autocorrelation), and thus represent an accurate approximation. In Fig. 4, the outdoor temperature and global irradiance time series are shown over a one-year period.
Fig. 1. A solar combisystem designed using Polysun.
R. Bornatico et al. / Energy 41 (2012) 31e37
33
Fig. 3. Yearly heat requirements for space heating in a mid-sized single-family house.
3. Optimization framework The main components of the system as presented in the previous section are the collector, the tank, and the APU. In order for these components to be optimally dimensioned, the corresponding sizes have to populate the input parameter set that is fed to the simulation kernel. In Table 1 the input set is listed together with the corresponding minimal and maximal values used throughout this work. The extreme values for the output set shown in Table 2 do not represent a given constraint. Rather, these values have been obtained from simulations throughout this work. They are to show that the investigated range for these values is significant. The goal of the optimization routine consists in finding the parameters wopt that minimize a general function f : Rn /R that is formulated in the canonical form as follows:
min f ðwÞ w
subject to gðwÞ < 0; hðwÞ ¼ 0
(1)
The functions g and h include all constraints to the problem defining the set Q of the acceptable w˛Q. These constraints include the boundary values of w and two physical constraints on the system. First, for a given w to be acceptable, the heat requirements for DHW and space heating must be satisfied. Secondly, the temperature of the fluid flowing in the collector must not exceed the 100 C threshold for longer than 1% of the simulation time. This last condition limits the stagnation time to reasonable values, thereby ensuring the safe operation and a long life of the collector. In the presented problem a numeric optimization is necessary because an algebraic solution to the problem does not exist. Furthermore it is clear that the topographical properties of the function f play a crucial role in how the optimization algorithm converges from starting values w0 toward wopt , and also that the choice of the right cost function is determinant for the successful
convergence of the algorithm to the point of interest. Accordingly, all relevant quantities should appear in the cost function. In the case of the optimization of the combisystem, the objectives of the authors are to maximize the solar fraction, to minimize the total energy use, and to minimize the additional cost of the installation.
8 9 < X fj ðwÞ = f ðwÞ : ¼ yry ¼ wj $ : fjmax ; j
(2)
Therefore, the rather intuitive linear form (2) is chosen where the solar fraction f1 ðwÞ ¼ SFnðwÞ, the total energy use of the system f2 ðwÞ ¼ EtotðwÞ and the cost of the installation f3 ðwÞ ¼ CostðwÞ depend on the sizes of the components w. Since a scalar value of f is required by the minimization algorithm, the weighted sum, with weights wi of the normalized fi, yields a dimensionless, normalized cost function f ðwÞ. Note that in order to maximize the solar fraction by a minimization of f, w1 must be negative. An important result from a user’s point of view is the pricing of the components. This term is based on a linear pricing assumption for each component and on current oil prices. This choice is shown in Ref. [11] to produce meaningful results. The relevant values used for the installation pricing are listed in Table 3. 4. The particle swarm optimization algorithm The analogy of a population of individuals, called particles (e.g. a swarm of bees) having the common goal of finding an optimal position (e.g. the best flower in a field) is useful to picture how the algorithm converges toward a solution. This socially inspired optimization technique was first introduced by Kennedy and Eberhart [12] as a promising tool that does not requires the computation of derivatives.
Fig. 4. Meteorological data for Zurich, Switzerland. Solid: outdoor temperature, dotted: global irradiance.
34
R. Bornatico et al. / Energy 41 (2012) 31e37
Table 1 Input parameter set. Item
Symbol
Min
Max
Unit
Collector area Tank volume APU size
w1 w2 w3
1 100 5
40 3000 50
m2 l kW
Table 2 Output parameter set. Item
Symbol
Min
Max
Unit
Solar fraction Energy use Cost
f1 f2 f3
2.8% 12 988 12 988
42.7% 45 407 46 226
e kW h V
Table 3 Relevant factors for the installation pricing. Item
Value
Unit
Collector Collector Storage tank APU Heating oil Heating oil Maintenance
100 360 4.22 50 10.5 0.85 100
V V/m2 V/l V/kW kW h/l V/l V/year
Each particle i updates its position toward better fitness areas in the problem domain landscape based on a self cognitive term expressed by the best position (Pi) found so far by the particle, and the global best position encountered so far by the entire swarm (G). Starting with random initial velocities v0i and random initial positions p0i inside a closed set, the basic PSO algorithm computes the velocity and position of the i-th particle, i ¼ 1; 2; .; n for every iteration k ¼ 0; 1; 2; .kmax in the following way:
vkþ1 ¼ vki þ g1i Pi pki þ g2i G pki i
(3)
¼ pki þ vkþ1 pkþ1 i i
(4)
The random numbers g1;2 ˛½0; 1 influence the magnitude of the two vectors ðPi pki Þ and ðG pki Þ, thus influencing indirectly the . Fig. 5 depicts how the update from magnitude and direction of vkþ1 i k to kþ1 is performed for a simplified two dimensional problem with no randomness g1,2 ¼ 1. Note that every dimension represents a parameter to be optimized.
Fig. 5. Simplified PSO iteration in two dimensions.
The common PSO algorithm includes three additional parameters, namely an inertia function fk and two acceleration constants a1,2. The inertia function used hereafter is linearly decreasing with respect to iterations, thus reducing the influence of past velocities, thereby enabling the algorithm to adapt to small regions as the optimization converges. This adaptation capability is of key relevance since it makes the use of hybrid methods unnecessary. The resulting velocity equation, together with the position update, represents the core of the PSO algorithm [13,14].
h i h i þ a2 g2i G pki vkþ1 ¼ fk vki þ a1 g1i Pi pki i
(5)
pkþ1 ¼ pki þ vkþ1 i i
(6)
Further extensions of the algorithm can include natural selection considerations such as a varying population size where bad particles are killed and new particles are generated in the vicinity of the global best point. Dynamic adjustment of swarm parameters have been also proposed as a way to avoid local optima [15e17]. The work of Clerc and Kennedy [18] describes the importance of including constraints on the velocity vector. Their work shows that if no constraints are set the algorithm can become unstable. Note that all these extensions increase the computational complexity of the algorithm and should therefore be avoided, if not necessary. In the problem presented here a fixed population of 10 particles has been deemed to be sufficiently effective. Further, a linearly decreasing inertia function already accounts for stability, as the introduction of velocity limits can be shown to be a loose constraint. 5. Optimization results The results reported in this section have been based on simulation results of the software Polysun when simulating the model presented in Section 2. The converged results of the minimization of Eq. (2) by means of the PSO algorithm are presented in Section 5.1, while results obtained using a genetic algorithm follow in Section 5.2. A further investigation on the solution is presented in the form of a sensitivity analysis, which follows in the subsequent Section 6. 5.1. Particle swarm optimization: results The minimization of the fitness function by means of the PSO algorithm is shown in Fig. 6. The fitness function is plotted against particle generations where each generation consists of 10 particles. Note that points which do not satisfy the optimization constraints have been omitted. The PSO algorithm converges toward a minimal fitness value of f opt; PSO ¼ 1:278 and the interpretation of this value, with respect to its physical terms, is presented in Table 4. The initial population was randomly generated inside the boundary set
Fig. 6. Fitness function minimization with PSO.
R. Bornatico et al. / Energy 41 (2012) 31e37 Table 4 PSO optimization of fitness function terms.
35
Table 5 GA optimization of fitness function terms.
Term
Starting value
Optimized value
Change
Term
Starting value
Optimized value
Change
f SFn Etot Cost
1.435 18.55% 16 428 kW h 13 210 V
1.278 21.8% 15 806 kW h 8 963 V
10.1% þ17.5% 3.8% 32.1%
f SFn Etot Cost
1.335 20.6% 16 109 kW h 10 952 V
1.278 22.15% 15 750 kW h 9 230 V
4.2% þ7.5% 2.2% 15.7%
Q. It follows that the starting values as reported in Table 4, are averaged values of the feasible points of the initial population. It is important to judge the convergence of the algorithm by considering the spread of the particles. Fig. 7 illustrates the convergence of all particles, for all parameters considered in this study. From Figs. 6 and 7 it can be seen that after 50 generations the algorithm has converged to the vicinity of the optimal values which are found to be wopt;PSO ¼ f14:5 m2 ; 498:98 l; 8:5 kWg. 5.2. Genetic algorithm optimization: results On the same problem a standard GA has been applied where, if applicable, the same optimization parameters as in the PSO were used, specifically, a population of 10 particles and the same parameter set boundaries. In this case also, the fitness function decreases toward a minimal value of f opt; GA ¼ 1:278, a value that is consistent with the result of the PSO optimization. The comparison of all other optimal fitness function terms between PSO (Table 4) and GA (Table 5) shows that the two solutions differ by a minimal amount only. Also for the GA, the convergence of the fitness function and the system parameters is shown in Figs. 8 and 9. Fig. 8 shows that the convergence of all parameters over generation number is comparable to those obtained with the PSO optimization as the algorithm converges toward an optimal value of wopt; GA ¼ f15:1 m2 ; 514 l; 8:6 kWg. A slight difference in the optimal parameters obtained for the collector area and the tank volume can be noticed. However, this does not influence significantly the fitness function values, thus proving that both algorithms converge to a similar solution. Given
a
b
that the search space is rather big and that the starting individuals are randomly drawn inside this space, the fact that the optimal values found with GA and PSO are very similar supports the consistency of the results and the solution found in this section. 6. Parameter sensitivity analysis The results in Section 5 are specific for the considered system, i.e. a daily DHW consumption of 200 l/day at 45 C and an SH setpoint temperature of 20 and 18 C during day and night respectively. First the sensitivity of the results with respect to a variation of these parameters is explored and secondly, the same analysis is performed over the optimization parameters w1;2;3 . All results in this section are produced using the optimal PSO values opt;PSO opt;PSO and fj as nominal values when necessary. wi 6.1. Operational parameter variation Thus far, constant operational parameters have been used. Among them, the amount of DHW consumed (p1) and its temperature (p2) together with the heating setpoint temperature of the building during day (p3) and night (p4) are investigated here. Initially the DHW volume p1 is varied between {50, 100, 200, 300, 400} l/day and the simulation outcomes are linearly fitted to produce the results shown in the first column of Table 6. For instance the ratio ðvSFnÞ=ðvp1 Þ is 0.01 [%/l] indicating a solar fraction increment of 1% for every additional 100 l/day DHW consumed on top of the nominal 200 l/day. Similarly, the DHW temperature (p2) is varied between {35, 40, 45, 50, 55} C, p3 ˛ {18, 19, 20, 21, 22} C and p4 ˛ {16, 17, 18, 19, 20} C to produce the other values in Table 6. The sensitivity of the results relative to the corresponding nominal values is also a useful mean of comparison. The relative sensitivity computed using Eq. (7) is graphically shown in Fig. 10 where pnom ¼ {200 l/day, 45 C, 20 C, 18 C}. The relative variation of the daily heating setpoint temperature clearly dominates all other effects and indicates that this parameter is of crucial importance in the characterization of the system. From this figure it can be read that a 1% increase in p3 causes an additional energy consumption of around 2%.
Sij ¼
v fj =fjopt;PSO vðpi Þ=pnom i
c
Fig. 7. PSO convergence of system parameters towards 14.5 m2 (a), 499 l (b), and 8.5 kW (c).
Fig. 8. Fitness function minimization with GA.
(7)
36
R. Bornatico et al. / Energy 41 (2012) 31e37 Table 7 Absolute sensitivity of optimization parameters.
a
SFn [%] Etot [kW h] Cost [V]
w1 ½m2
w2 ½l
w3 ½kW
0.53 94.4 352.9
0.01 0.78 4.2
0.3 119.30 60.6
b
c
Fig. 9. GA convergence of system parameters towards 15.1 m2 (a), 514 l (b), and 8.6 kW (c).
Fig. 11. Relative sensitivity of the optimal solar fraction (1), total energy consumption (2) and cost of installation (3) related to a variation of collector size w1 , tank volume w2 , and rated power of the auxiliary power unit w3 .
Sij ¼ Table 6 Absolute sensitivity of operational parameters.
SFn [%] Etot [kW h] Cost [V]
p1 [l/day]
p2 [ C]
p3 [ C]
p4 [ C]
0.01 9.6 0.73
0.02 70.6 5.16
1.14 1638.5 125.66
0.24 91.3 7.69
6.2. Optimization parameter variation In contrast to the previous section where operational parameters are varied for a given installation, this section deals with varying size of collector, tank volume and APU power. Each one of these parameters is varied by {10%, 5%, þ5%, þ10%} of the respective nominal values, which correspond to the optimal PSO results wopt; PSO ¼ f14:5 m2 ; 498:98 l; 8:53 kWg. The same methodology as in the previous section is applied here and the results are presented in absolute values in Table 7 and in Fig. 11 according to Eq. (8).
Fig. 10. Relative sensitivity of the optimal solar fraction (1), total energy consumption (2) and cost of installation (3) related to a variation of DHW volume p1, DHW temperature p2, SH day setpoint temperature p3 and SH night setpoint temperature p4.
opt;PSO v fj =fj opt;PSO
vðwi Þ=wi
(8)
The relative sensitivity of the solar fraction Si1 is depicted in Fig. 11 by the bars (1) where it can be seen that an increase of the collector or tank size has a positive effect on the solar fraction. This effect is dominated in relative terms, by the solar collector with S11 ¼ 0.35. Further, bar (2) of w3 shows that an increase in APU size in the vicinity of the optimal solution causes a greater energy consumption of the system due to a shifted operating point in the APU consumption and efficiency map, whereas an increase of collector size or tank volume reduce the overall energy consumption of the system. The bars labeled (3) in Fig. 11 show that an increase in collector size has the greatest relative influence on the cost of the installation. The tank volume has a smaller influence, while the influence of the APU size is relatively low. 7. Conclusions and discussion In this paper the optimal sizing of a solar thermal system has been presented. The analysis is performed on a solar combisystem for a mid-sized single-family house in Zurich, Switzerland. While the optimization framework is in principle independent on the optimization algorithm, a detailed analysis has been carried out on the performance of the Particle Swarm Optimization algorithm when it is applied for solving this problem. The results are comparable to those obtained with the more common Genetic Algorithm. When the implementation efforts and the computational power demand are considered as well, the PSO is a slightly better choice for solving the presented problem. A collector size of 14.5 m2 together with a tank volume of 498.98 l and an APU nominal power of 8.5 kW, are the optimal sizes for the main system’s components, which lead to a solar fraction of 21.8%, a total energy use of 15 806 kW h and a cost of the installation of 8 963 V. The parameter sensitivity analysis shows that the size of the collector has the greatest influence on the solar fraction, energy use and the installation cost, while the tank volume influence is significant on the solar fraction and cost of installation. The size of the APU has a relatively small effect on the both energy use and installation cost. The variation of selected operational parameters
R. Bornatico et al. / Energy 41 (2012) 31e37
has been investigated as well, and shows that the sensitivity of the daily heating setpoint temperature of the building is by far the dominant parameter when compared to all parameters considered in this study. Future work includes an extension and further investigation of the parameter set. Furthermore, the proposed methodology is to be applied over other solar installations and/or more complex models including, for example a detailed life cycle analysis of the system. Acknowledgments This work has been partly funded by the Swiss Innovation Promotion Agency CTI, whose support is gratefully acknowledged. References [1] Sawin E, Jones A, Sterman J. Final Copenhagen accord press release 19 december 09-expanded version; 2009. [2] Suter JM, Letz T, Weiss W, Inbnit J. Solar combisystems in Austria, Denmark, Finland, France, Germany, Sweden, Switzerland, the Netherlands and the USA overview 2000; 2000. [3] Lund P. Sizing and applicability considerations of solar combisystems. Solar Energy 2005;78(1):59e71. [4] Jordan U, Vajen K. Influence of the dhw load profile on the fractional energy savings: a case study of a solar combi-system with trnsys simulations. Solar Energy 2001;69(Suppl. 6):197e208. [5] Klein S, Beckmann B, Duffie J. Trnsys, a transient system simulation program, users manual, version 16; 2006.
37
[6] Duffie J, Beckman W. Solar engineering of thermal processes. 3rd ed. Wiley Interscience; 2006. [7] Witzig A, Foradini F, Probst MCM, Roecker C. Simulation tool for architects. In: Proceedings of the 2009 CISBAT international conference; 2009. [8] Meteotest. Meteonorm database. Bern, Switzerland: Meteotest; 2000. [9] Meteotest. Meteonorm handbook, part 2: theory. Bern, Switzerland: Meteotest; 2000. [10] Aguiar RJ, Collares-Pereira M, Conde JP. Simple procedure for generating sequences of daily radiation values using a library of Markov transition matrices. Solar Energy 1988;40(3):269e79. [11] Fraisse G, Bai Y, Pierrs NL, Letz T. Comparative study of various optimization criteria for sdhws and a suggestion for a new global evaluation. Solar Energy 2009;83(2):232e45. [12] Kennedy J, Eberhart R. Particle swarm optimization. In: Proceedings of the 1995 IEEE International Conference on Neural Networks, vol. 4; 1995. p. 1942e1948. [13] Birge B. Psot e a particle swarm optimization toolbox for use with matlab. In: Proceedings of the 2003 IEEE Swarm Intelligence Symposium SIS 03; 2003. p. 182e186. [14] Bratton D, Kennedy J. Defining a standard for particle swarm optimization; 2007. 120e127. [15] Dimopoulos GG, Frangopoulos CA. Optimization of energy systems based on evolutionary and social metaphors. Energy 2008;33(2):171e179. 19th International Conference on Efficiency, Cost, Optimization, Simulation and Environmental Impact of Energy Systems e ECOS 2006. [16] Shi Y, Eberhart R. A modified particle swarm optimizer. In: Proceedings of the 1998 IEEE World Congress on Computational Intelligence; 1998. p. 69e73. [17] Fourie P, Groenwold A. The particle swarm optimization algorithm in size and shape optimization. Structural and Multidisciplinary Optimization 2002;23: 259e67. [18] Clerc M, Kennedy J. The particle swarm e explosion, stability, and convergence in a multidimensional complex space. IEEE Transactions on Evolutionary Computation 2002;6(1):58e73.