- Email: [email protected]
Particle Swarm Optimization-based algorithms for solving inverse heat conduction problems of estimating surface heat flux
International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Particle Swarm Optimization-based algorithms for solving inverse heat conduction problems of estimating surface heat flux Fung-Bao Liu ⇑ Department of Mechanical and Automation Engineering, I-Shou University, Kaohsiung 84001, Taiwan
a r t i c l e
i n f o
Article history: Received 9 February 2011 Received in revised form 14 August 2011 Available online 21 December 2011 Keywords: Inverse analysis Particle Swarm Optimization Mutation operator
a b s t r a c t An inverse analysis of estimating a time-dependent surface heat flux for a three-dimensional heat conduction problem is presented. A global optimization method known as Particle Swarm Optimization (PSO) is employed to estimate the unknown heat flux at the inner surface of a crystal tube from the knowledge of temperature measurements obtained at the external surface. Three modifications of the PSO-based algorithm, PSO with constriction factor, PSO with time-varying acceleration of the cognitive and social coefficients, and PSO with mutation are carried out to implement the optimization process of the inverse analysis. The results show that the PSO with mutation algorithm is significantly better than other PSO-based algorithms because it can overcome the drawback of trapping in the local optimum points and obtain better inverse solutions. The effects of measurement errors, number of dimensionalities, and number of generations on the inverse solutions are also investigated. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Inverse heat conduction problems are known as back-analysis problems, i.e. the reversal of the cause-effect sequence, in the field of heat transfer analysis. An inverse problem means that some of the initial, boundary conditions or material properties are not fully specified as determined from the measured temperature profiles at some specific locations. The inverse problems in most situations are likely to be ill-posed [1]. Solutions of the inverse problem are very sensitive to the measurement errors, i.e. small errors in the measured data values can produce very large errors in solutions. In general, the uniqueness and stability of an inverse problem solution are not guaranteed. In recent years, the inverse problems have been studied extensively due to their applications in various engineering disciplines. In general, the technique approaches the inverse heat conduction problem as an optimization problem, i.e. the problem is defined as the minimization of a cost function or a fitness function measuring the distance between measurements and predictions. A variety of numerical methods and computational algorithms for solving the inverse heat conduction problem have been developed over the years. Most of the analytical or numerical methods dealing with the inverse heat conduction problems of estimating ⇑ Address: Department of Mechanical and Automation Engineering, I-Shou University, No. 1, Sec. 1, Syuecheng Rd., Dashu District, Kaohsiung 84001, Taiwan. Tel.: +886 7 6577711x3228; fax: +886 7 6578853. E-mail address: fl[email protected] 0017-9310/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2011.12.007
a single surface condition are one- or two-dimensional. Only a few works estimate the unknown sources in three-dimensional problems due to the complexity of geometry or time-consuming computation. Zueco et al. [2] studied the unsteady inverse heat transfer problem in a one-dimensional geometry subjected to a time-dependent incident heat flux by using the numerical network simulation method. Two different types of incident heat flux waveform, including triangular and rectangular, had been determined. The temperature dependency of the thermal properties, conductivity and specific heat were assumed to be polynomial functions in [2]. Huang and Wu [3] investigated the one-dimensional transient inverse heat conduction problem by using a conjugate gradient method for identifying the unknown boundary heat flux. Instead of employing the conventional parabolic differential equations for the heat transfer problem, a hyperbolic heat equation model was used and analyzed. Identification of thermal conductivity and the shape of an inclusion for a two-dimensional solid body were considered by Ardakania and Khodadad [4]. The Particle Swarm Optimization algorithm (PSO), coupled with the boundary elements method, was used to implement this inverse heat conduction problem. Qi et al. [5] applied the multi-phase PSO method to solve the inverse radiation problem. The proposed method adopted the advantages of both two-group PSO and multi-start PSO algorithms by dividing particles into multi-group to increase the diversity of the swarm and introducing different phases with different flying directions and
2063
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
Nomenclature Cp c1 c2 cl cu f fb fw Gbest k L M Mg N n Pbest,i q r r1
specific heat cognitive coefficient social coefficient lower bound of cognitive and social coefficient upper bound of cognitive and social coefficient the fitness function to be minimize the best fitness of particles the worst fitness of particles best position among all the particles in the swarm thermal conductivity length of the tube total number of particles in a swarm total number of generation dimensionality current generation previous best position of particle i heat flux radial coordinate random number in the interval [0, 1]
searching ways. The standard PSO and the stochastic PSO were also carried out to compare with the proposed method for their performance and accuracy. The effectiveness and efficiency of Particle Swarm Optimization technique in inverse heat conduction analysis were investigated by Vakili and Gadala [6]. Three variations of the PSO method, i.e. basic PSO, repulsive PSO, and complete repulsive PSO, were used to solve the boundary inverse heat conduction problem, in one, two, and three dimensions. The results showed that PSO can alleviate some of the stability problems of the classical approaches for solving the inverse heat conduction problems. For the problems studied, some variants of PSO are proven to be more efficient than other algorithms. Abboudi and Artioukhine [7] proposed a numerical approach for simultaneous estimation of two boundary conditions in a twodimensional linear heat conduction problem. The inverse numerical algorithm was based on the iterative regularization method and on the conjugate gradient method. Unknown functions were parameterized in the form of a cubic B-spline. An optimal choice of the descent parameter matrix was used in the iterative regularization method. Haghighi Golbahar et al. [8] studied two-dimensional transient inverse heat conduction problem of functionally graded materials. A combination of the finite element (FE) and differential quadrature (DQ) methods was employed for solving the direct problem. The inverse analysis using conjugate gradient method was transformed into the solution of three problems; namely, the direct, sensitivity, and adjoint problems. These three problems were also solved by the coupled DQ–FE method. The surface heat fluxes on cutting edges of cutting tools were estimated in the three-dimensional inverse heat conduction problem by Huang and Lo [9]. The steepest descent method with an adjoint equation and a general-purpose commercial code were successfully applied in determining the time-dependent local heat fluxes. Although the inverse heat conduction problems have been extensively studied for different applications in the past decades, little work has been done for the inverse numerical algorithms using heat flux measurements in the objective function [10]. Louloul and Scott [11] investigated an inverse bioengineering problem of estimating the time-dependent blood perfusion and the thermal conductance between the probe and the tissue. A methodology in which a combined parameter and function estimation was developed based on the minimization of an objective function containing computed and measured heat fluxes.
r2 ri ro Si T(r, h, z, t) Ti V nþ1 i w nþ1 Xi Y(r, h, z, t) z
random number in the interval [0, 1] inner radius of tube outer radius of tube weight factor computed temperature initial temperature of the fluid current velocity of particle i inertia weight current position of particle i measured temperatures axial coordinate along the tube
Greek symbols v constriction factor q density h angular coordinate
In this work, an inverse analysis for the recovery of time-varying surface heat flux in three-dimensional cylindrical coordinates is investigated. Transient temperature measurements at the boundary, obtained from the solution of the direct problem, served as the simulated experimental data needed as input for the inverse analysis. A fitness function which is defined by the quadratic residual between the measurements and the calculated temperatures is minimized. Three modifications of the PSO-based algorithms, namely PSO with constriction factor (PSOC), PSO with time-varying coefficient (PSOTV), and PSO with mutation operator (PSOM), are carried out to implement the optimization process of the inverse analysis. The inclusion of mutation operator in the Particle Swarm Optimization algorithm which speeds up the convergence and escape from local minimum is also examined. Numerical examples are provided to compare the effectiveness of the PSOC, PSOTV, and PSOM algorithms. 2. The direct heat transfer problem The thermal behavior of the high power helicon source used to produce plasma in the Variable Specific Impulse Magnetoplasma Rocket (VASIMR™) was investigated by Mulcahy et al. [12] numer-
r z
θ
The external surface is insulated
Heat flux q(z,θ,t) is applied at the interior surface Fig. 1. Physical model.
2064
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
ically and experimentally. The same physical model was used for the direct problem in this study as shown in Fig. 1. To illustrate the methodology for developing expressions in determining the unknown surface heat flux in a circular channel, the following cylindrical geometry with constant thermo-physical properties is considered. The initial temperature is equal to Ti. When the time t > 0, the boundary conditions at boundary z = 0, z = L, and the external surface r = ro are assumed to be insulated, while an unknown heat flux q(z, h, t) is applied at the inner surface r = ri. The dimensional mathematical formulation of this transient heat conduction problem in cylindrical coordinates is given as follows [12]:
" # @T @ 2 T 1 @T 1 @ 2 T @ 2 T ; qC p ¼ k 2 þ þ þ @t @r r @r r2 @h2 @z2
ð1aÞ
and subjected to the following boundary and initial conditions.
k
@Tðr; h; z; tÞ ¼ 0; @r
at r ¼ r o ;
0 6 z 6 L;
0 < h < 2p;
t > 0; ð1bÞ
k
@Tðr;h;z;tÞ ¼ qðh;z;tÞ; at r ¼ ri ; 0 6 z 6 L; 0 < h < 2p; t > 0; @r ð1cÞ
k
@Tðr; h; z; tÞ ¼ 0; @z
at z ¼ 0;
and z ¼ L;
ri 6 r 6 ro ;
t > 0; ð1dÞ
Tðr; z; h; tÞ ¼ T i ;
at 0 6 z 6 L;
ri 6 r 6 ro ;
0 < h < 2p;
t ¼ 0: ð1eÞ
where k, q, Cp and t are the thermal conductivity, density, heat capacity and time, respectively. The problem defined by Eqs. (1a)–(1e) is called a direct problem in which the interior temperature T(r, z, h, t) is determined when the thermal properties, initial and boundary conditions are known. The direct problem can be solved by discretizing the above equations with the finite difference method, the finite element method, or the boundary element method. In this work, the weighted Crank–Nicolson implicit finite difference method is used to discretize the Eqs. (1a)–(1e) and solve the direct problem. 3. The inverse heat conduction problem For the inverse analysis of heat conduction problem, the geometric features and/or material properties are unknown, but some of the unknown boundary data can be measured and used as additional information to estimate the unknown input parameters. In this case, the boundary heat flux at r = ri is regarded as an unknown function, but other parameters in Eqs. (1a)–(1e) are known. In addition, temperature measurements at r = ro are considered available. Let temperature measurements taken by sensors at r = ro denoted as Y(ro, h, z, t). The aim of the inverse analysis is to iteratively estimate the unknown heat flux using an optimization procedure which results a negligible difference between measurements taken at the external surface and temperatures computed from the numerical model in Eqs. (1a)–(1e). Therefore, the solution of the present inverse problem is obtained by minimizing the following fitness function f(q) formulated as a sum-of-squares of residuals: Z t f Z L Z 2p min f ðqðh;z;tÞÞ ¼ ½Tðro ;h;z;tÞ Yðr o ;h;z;tÞ2 dh dz dt: ð2Þ 0
0
0
The inverse problem is recast as an optimization problem. A variety of numerical and analytical techniques have been developed to solve the optimization problems. Among them, the gradi-
ent-based methods have the advantages of the fast convergence speed and the high accuracy of the estimated solutions. However, a good initial guess solution is needed to avoid trapping into local optimum. Another approach is the stochastic methods, i.e. genetic algorithms (GA), simulated annealing (SA), and Particle Swarm Optimization algorithms (PSO). The advantages of the stochastic methods are their capability in searching for the global optimum and no computation of the complicated gradients. However, the time-consuming calculation for obtaining optimal solutions is a drawback of the stochastic methods. 3.1. The Particle Swarm Optimization algorithm The Particle Swarm Optimization algorithm introduced by Kennedy and Eberhart [13] in 1995 is a stochastic optimization technique which draws inspiration from the social behavior of a flock of birds or the collective intelligence of a group of social insects with limited individual capabilities. The standard PSO model consists of a swarm of M particles moving in a problem search space. Each particle is a potential solution of the global optimum over a given domain D. For a N-dimensional search space, the position of the ith particle is represented as Xi = (xi1, xi2, . . . , xiN). At each generation, the new particle position is found by adding a displacement to the current position where the displacement is the particle velocity multiplied by a time step of one as shown in Eq. (3)
X nþ1 ¼ X ni þ V nþ1 : i i
ð3Þ
In Eq. (3), X nþ1 and X ni represent the current and previous positions i of the particle i, V nþ1 is the current velocity of the particle i and is i represented as V nþ1 ¼ ðmi1 ; mi2 ; . . . ; miN Þ. The velocity of each particle i is also updated at each generation and is given by
V nþ1 ¼ V ni þ c1 r 1 ðP best;i X ni Þ þ c2 r2 ðGbest X ni Þ; i V nþ1 i
ð4Þ
V ni
where and are the current and previous velocities of the particle i, respectively. Each particle maintains a memory of its previous best position, say Pbest,i = (pi1, pi2, . . . , piN), where the position giving the best fitness function value. The best one among all the particles in the swarm is represented as the global best position, say Gbest = (pg1, pg2, . . . , pgN). The new velocity in Eq. (4) can be seen as the sum of three parts. The first part of Eq. (4) represents the previous velocity and is called the momentum part. The second part of the Eq. (4) represents the tracking of best position for individual particle and is called the cognition part. The third part of the Eq. (4) represents the cooperation among particles in the swarm and is called the social component. The cognitive learning coefficient c1 and social learning coefficient c2 are known as accelerating factors, and r1 and r2 are two random numbers generated by the uniform distribution within 0 and 1. The relative sizes of these components determine their contribution to the new particle velocity. The most common setting for c1 and c2 are 2.0 for the standard PSO algorithm. 3.2. The modifications of the Particle Swarm Optimization One difficulty encountered in the standard PSO algorithm is that it could easily fall into local optimum in many optimization problems. To well control the scope of the search and improve the performance of the algorithms, many modifications on the standard PSO algorithm have been proposed. The inclusion of inertia weight in the PSO algorithm was investigated by Eberhart and Shi [14]. The inertial weight, w, plays a multiplying factor of current velocity as follows:
V nþ1 ¼ wV ni þ c1 r 1 ðPbest;i X ni Þ þ c2 r2 ðGbest X ni Þ: i
ð5Þ
2065
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
Better selection of the inertial weight provides a balance between global and local exploration and exploitation. Therefore, less iteration is needed to find the optimal solution and improve the performance of the algorithm. A constriction factor was incorporated into the PSO algorithm by Clerc’s suggestion [15] to insure the convergence of the algorithm (PSOC). The velocity term was update as follows:
V nþ1 ¼ vðV ni þ c1 r 1 ðPbest;i X ni Þ þ c2 r 2 ðGbest X ni ÞÞ: i
ð6Þ
When Clerc’s constriction method is used, the constant multiplier v is approximately equal to 0.7298 and the two coefficients c1 and c2 are 2.05. The time-varying acceleration of the cognitive and social coefficients in the algorithm (PSOTV) was investigated by Suganthan [16]. In his strategy, the upper and lower bounds of the cognitive and social coefficients were predefined and both coefficients, c1 and c2, were linearly decreased over the number of generations in the search. That is, the cognitive coefficient c1 and social coefficient c2 are dynamically changed as follows:
c1 ¼ c2 ¼ ðcu cl Þ
Mg n þ cl ; Mg
ð7Þ
where n is the current generation, Mg means the total number of generations, cu is the upper bound and cl is the lower bound of the coefficients. For the modifications of the cognitive coefficient c1 and the social coefficient c2, Venter [17] concluded that better quality solutions could be obtained with the small cognitive coefficient and large social coefficients. Furthermore, Ratnaweera et al. [18] proposed another time-varying acceleration coefficient strategy in which the cognitive coefficient was reduced linearly over the number of generations in the search while the social coefficient was linearly increased inversely. In other words, the cognitive coefficient was changed as the Eq. (7), while the social coefficient was changed as follows:
c2 ¼ ðcu cl Þ
n þ cl : Mg
ð8Þ
3.3. The Particle Swarm Optimization with mutation (PSOM)
ð9Þ
where the weight factor Si is based on the fitness function of each particle and given by
( Si ¼
fw f ðX i Þ ; fw fb
if f b –fw ;
1;
otherwise;
mij ¼
0:5 xmax r 1 ;
if r2 < 0:5;
0:5 xmax r 1 ; otherwise;
ð11Þ
where r1 and r2 are two random numbers generated with uniform distribution between 0 and 1. The value of xmax is the maximum number in the search space. In each generation, the number of particles and the number of dimensionalities applied with mutation strategy are predefined and the mutation probability can be calculated. 3.4. Computation procedure for the Particle Swarm Optimization with mutation The computational steps of the PSOM algorithm described above are given as follows: Step 1: Generate the initial particles in a swarm by randomly generating the position and velocity for each particle. Step 2: Evaluate the fitness function of each particle. Step 3: Update the Pbest,i for each particle, if its fitness is smaller than the fitness of its previous best position (Pbest,i). Step 4: Update the Gbest, if the fitness function of a particle is smaller than the fitness of the best position of all particles (Gbest). Step 5: Compute the cognitive coefficient c1 by Eq. (7). Step 6: Compute the social coefficient c2 of each particle according to Eqs. (9) and (10). Step 7: Update each particle according to Eqs. (5) and (3). Step 8: Apply mutation operator according to Eq. (11). Step 9: Repeat the loop until stopping criteria or predefined number of generations is reached. 4. Numerical examples and discussion
The method of PSO sometimes suffers from a lack of diversity among particles, which leads to a stagnation stage. A number of modifications have been proposed to overcome this drawback, with varying degrees of success [19,20]. One of the modifications is to incorporate mutation operators into the PSO algorithm. Some research works have shown that mutation operators can help the PSO algorithm to prevent the premature convergence [21,22]. In this work, a PSO algorithm with a modified social coefficient and a mutation strategy proposed by Cai et al. [23] is employed for solving the inverse heat conduction problem. To overcome the decrease of searching efficiency for the PSO algorithm, a dispersed social coefficient setting is introduced to the PSO as follows:
c2i ¼ ðcu cl Þ Si þ cl ;
is, the larger the Si is. Furthermore, a mutation strategy is introduced to enhance the ability of escaping from the local optimum to avoid premature convergence. At each generation, particle i is uniformly random selected within the whole swarm and the dimensionality j of the selected particle is also random selected to perform the mutation. Then, the velocity vij is updated as follows:
ð10Þ
where fw and fb are the worst and the best fitness of particles in the swarm. The weight factor Si is an information index to represent the difference of the particles, according to their fitness values of the current position to the best position. The better the particle
A solution method for solving the inverse heat conduction problem by the Particle Swarm Optimization algorithms are carried out to determine the unknown surface heat flux in cylindrical coordinates. In order to give some idea on the physical significance of the various variables, a waste heat flux from a helicon plasma discharge in a rocket is considered. The channel is a quartz tube with density 2200 kg/m3, conductivity 1.4 W/m K, and specific heat 670 J/kg K. The length of the tube is L = 0.3 m with inner and outer diameters being ri = 0.045 m and ro = 0.0475 m, respectively [12]. The exact function of the time-dependent heat flux on the inner surface of quartz tube is given as:
2z2 2z qðz;h;tÞ ¼ qmax 2 þ þ 0:5Þð0:7 0:3cosðhÞ ; for 0:1 6 t 6 0:35; L L ð12aÞ and
qðz; h; tÞ ¼ 0;
for 0:1 < t;
t > 0:35:
ð12bÞ
Two different functional forms of the heat flux are tested. One is a step functional form with the maximum heat flux, qmax of 100 k W/m2. The other is a ramp change in which the heat flux increases linearly from zero at t = 0.1 s to the maximum at t = 0.225 s, before decreasing linearly to zero at t = 0.35 s.
2066
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
For present study, the numbers of grid points in r, z, and h directions are taken as 10, 60, and 40, respectively. The measured temperature locations are at the grid points of the external surface. The total grid number on the external surface is the number of tangential nodes time axial nodes, say 40 60. The measured time period is 0.01 s and the total measured time is 0.45 s, i.e. there are 45 time steps. Therefore, there are a total of 108,000 temperature measurements at the external surface with correspondent unknown heat fluxes at the inner surface in this study. The measurements used as input data for the inverse analysis are obtained from the solution of the direct problem. Note that real measurements generally contain errors. The recordings of temperatures at discrete regular times are modified systematically by introducing errors from a random number generator. Thus random errors are imposed onto the calculated exact temperatures as follows:
Y ¼ Y exact þ xr;
ð13Þ
where Yexact is the solution of the direct problem with the exact boundary heat flux q(z, h, t); r is the standard deviation of the measurements; and x is a random variable within 2.576 and 2.576 for a 99% confidence bound. The inverse analysis is implemented by the Particle Swarm Optimization algorithm with different modifications. Three modified algorithms used in this work are PSO with constriction factor (PSOC), PSO with time-varying acceleration coefficient (PSOTV), and PSO with mutation (PSOM). The standard PSO is also carried out for comparison. The unknown heat flux is represented by a particle i at a position Xi with velocity Vi. The number of particles in the swarm is M = 50 for this study. The initial position and velocity are randomly generated for each particle. Therefore, the position Xi is a guessed heat flux and also a potential solution for the optimal solution. For each particle, the fitness function, i.e. Eq. (2), is evaluated in which the computed temperatures are obtained by solving the direct problem with the guessed solution Xi. The particle previous best position (Pbest,i) and the best position of all particles (Gbest) are updated according to their fitness functions. The update velocity is calculated for each particle to update the new position for the next iteration. In the PSOC algorithm, the constriction factor v is equal to 0.7298. Both the cognitive learning coefficient c1 and the social learning coefficient c2 are equal to 2.05. The inertia weight w decreases linearly from 0.9 to 0.4 and the coefficient c1 decreases from 25 to 5, while the coefficient c1 increases from 5 to 25 in the
PSOTV algorithm. The mutation rates of the particle in the swarm and the dimensionality of the selected particle are 4 and 2 22% in the PSOM algorithm respectively. Fig. 2 shows the estimated results of the step heat flux by using standard PSO, PSOC, PSOTV, and PSOM algorithms at location A(0.15, 0°). The PSO-based algorithms all have the same size of swarm and they run up to 2000 generations. The inverse solutions have the best predictions by using PSOM algorithm and the worst by using PSOC algorithm. The same conclusion can be drawn in the Fig. 3 for the estimated heat flux with ramp functional form at location C(0.05, 90°). Fig. 4 shows the comparison of fitness functions of the algorithms standard PSO, PSOC, PSOTV, and PSOM at location C(0.05, 90°) for ramp heat flux. It is observed that the values of fitness function become stable after a number of generations for the cases of standard PSO, PSOC and PSOTV algorithms. The algorithms even proceed with more generations, but improvements in the results are limited. The performance of the fitness function of PSOM is the best among the three algorithms. The prediction of the heat flux is acceptable at the 2000 generations for the PSOM algorithm. Therefore, with regard to Fig. 4 and to avoid long computational time, 2000 has been chosen as a stopping criterion. In Fig. 4, it is noted that the inverse solutions of the PSOC algorithm have the fast convergence, at about 400 generations, among the three modified algorithms. However, it is trapped into an unforeseen local optimum which provides the wrong solutions as been shown in Figs. 2 and 3. The inverse solutions of the PSOTV algorithm have been improved from the PSOC algorithm, but the improvement is not noticeable after 1200 generations and the solution is also trapped into another local optimum. The social coefficient c2 is an important parameter to control the weight of each particle associated with the best position. The modification by using the time-varying acceleration to the social coefficient c2 in the PSOTV algorithm is not good enough to prevent the premature convergence. However, the inclusion of the fitness function into the social coefficient c2 has improved the search ability for better solutions of the PSOM algorithm. Furthermore, the mutation operator in the PSOM algorithm maintains the diversity of global search capability and avoids the premature convergence as been shown in Figs. 2 and 3. The estimated results of the step heat flux at different locations, A(0.15, 0°), B(0.1, 45°), and C(0.05, 90°) using the PSOM algorithm are presented in Fig. 5. The heat fluxes have the same step functional form with different maximum values at the three different
200000 160000 120000
exact PSOC PSOTV PSOM PSO
60000
80000
heat flux (W/m2)
2
heat flux (W/m )
80000
exact PSOC PSOTV PSOM PSO
40000 0
40000
20000
0
-40000 -20000
-80000 -120000
-40000
0
0.1
0.2
0.3
0.4
0.5
0.6
time (s) Fig. 2. Comparison of estimations by standard PSO, PSOC, PSOTV, and PSOM algorithms for a heat flux with step functional form at location A(0.15, 0°).
0
0.1
0.2
0.3
0.4
0.5
time (s) Fig. 3. Comparison of estimations by standard PSO, PSOC, PSOTV, and PSOM algorithms for a heat flux with ramp functional form at location C(0.05, 90°).
2067
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
120000
1.0E+09
A-exact A-estimated B-exact B-estimated C-exact C-estimated
PSOC 100000
PSOTV PSOM
80000
PSO
heat flux (W/m2 )
fitness function f
1.0E+08
1.0E+07
1.0E+06
60000 40000 20000
1.0E+05 0
1.0E+04 0
500
1000 1500 No. of generation
2000
-20000
2500
0
0.1
0.2
0.3
0.4
0.5
time (s)
Fig. 4. Convergence of the fitness functions for a ramp heat flux by standard PSO, PSOC, PSOTV, and PSOM algorithms at location C(0.05, 90°).
Fig. 6. Inverse solutions of ramp heat fluxes by PSOM algorithm at locations A(0.15, 0°), B(0.1, 45°), and C(0.05, 90°).
150000
2
heat flux (W/m )
120000 90000
60000 30000
160000 exact 2000 1000 500
120000 heat flux (W/m2)
A-exact A-estimated B-exact B-estimated C-exact C-estimated
80000 40000 0 -40000
0
-80000
-30000
0
0
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
time (s)
time (s) Fig. 5. Inverse solutions of step heat fluxes by PSOM algorithm at locations A(0.15, 0°), B(0.1, 45°), and C(0.05, 90°).
locations. In Fig. 5, the results show that the estimations are in good agreement with the exact solutions over the time interval. The same results can be observed in Fig. 6 for the ramp heat flux. The Particle Swarm Optimization algorithm employed a stochastic search for the inverse solutions. Unlike the gradient-based method, no gradient is needed to be calculated in this case. Therefore, there is no difficulty in predicting the discontinuities at both step and ramp functional forms of heat flux. Predictions with three different number of generations, namely 500, 1000, and 2000, at location A(0.15, 0°) are illustrated in Fig. 7 for the PSOM algorithm. The results show a significant improvement in the solutions from 500 to 2000 generations. The weight factor Si of the social coefficient c2 has the ability to fine tune the coefficient in the PSOM algorithm. The mutation operator helps the algorithm to escape from the local optimum. As shown in Fig. 3, the fitness function is continuing improved over the range of this study. The noise is added to the exact temperature data with varying amplitude according to Eq. (13). Fig. 8 presents the estimated heat flux q(h, z, t) at location B(0.1, 45°) for three values of random noise r = 0, 1, and 2 in the simulated temperatures. As expected, the heat
Fig. 7. Estimated results for different number of generations by PSOM algorithm at location A(0.15, 0°).
flux fluctuation increases with the increase of noise in the measured temperatures. For the case of r = 2, the perturbations are randomly scattered around the exact profile which means that the estimation procedure is able to recover the general tendency of the time-dependent heat flux but fails to recover it precisely due to the inherent ill-posed character of inverse problems. It is noted that the characteristic of the PSOM algorithm makes it possible to obtain acceptable and stable results for different amplitude of the noises imposed on the measured temperatures. Estimated results of the heat flux by using the PSOM algorithm with different number of particles in a swarm at location A(0.15, 0°) are presented in Fig. 9. Three different number of particles in the swarm used here are M = 50, 30, and 20 with the same 2000 generations for computations. As expected, better predictions have been found for the case with the larger number of particles M = 50 in the swarm. More number of particles in the swarm can increase the convergence speed of computations. It is noticed that they all have the similar performance if the same fitness function is used as the stopping criterion. However, the case with M = 30 particles in the swarm takes about 4300 generations to reach the same fitness function as the case with 50 particles and about 6100 generations for the case of M = 20 particles in the swarm.
2068
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
160000
120000
References
2
heat flux (W/m )
PSOM algorithm to recover the inner surface heat flux without using any a priori information of the unknown transient functions.
exact σ=0 σ=1 σ=2
80000
40000
0
-40000
0
0.1
0.2
0.3
0.4
0.5
time (s) Fig. 8. Estimated results for measurements with errors by PSOM algorithm at location B(0.1, 45°).
120000 exact 100000
M=50 M=30
heat flux (W/m2)
80000
M=20
60000 40000 20000 0 -20000 0
0.1
0.2
0.3
0.4
0.5
time (s) Fig. 9. Estimated results for different number of particles in the swarm by PSOM algorithm at location A(0.15, 0°).
5. Conclusion An inverse analysis using Particle Swarm Optimization algorithms has been presented to identify different types of surface heat flux in a three dimensional heat conduction problem. The estimated results by using Particle Swarm Optimization with mutation have been compared against the Particle Swarm Optimization with the constriction factor and the Particle Swarm Optimization with the time-varying acceleration coefficient. Computational results of the PSOM algorithm have been shown to improve performance for the Particle Swarm Optimization algorithms. Simulated results with or without additive noise have shown that the PSOM algorithm allows identifying the surface heat flux with good accuracy as compared with the exact solutions. The obtained results underline the feasibility of the procedure and its capabilities for the
[1] J.V. Beck, B. Blackwell, C.R.St. Clair Jr, Inverse Heat Conduction, Wiley, New York, 1985. [2] J. Zueco, F. Alhama, C.F. Gonzalez Fernandez, Numerical nonlinear inverse problem of determining wall heat flux, Heat Mass Transfer 41 (2005) 411–418. [3] C.-H. Huang, H.-H. Wu, An inverse hyperbolic heat conduction problem in estimating surface heat flux by the conjugate gradient method, Journal of Physics D-Applied Physics 39 (2006) 4087–4096. [4] M.D. Ardakania, M. Khodadad, Identification of thermal conductivity and the shape of an inclusion using the boundary elements method and the Particle Swarm Optimization algorithm, Inverse Problems in Science and Engineering 17 (7) (2009) 855–870. [5] H. Qi, L.M. Ruan, M. Shi, W. An, H.P. Tan, Application of multi-phase particle Swarm Optimization Technique to inverse radiation problem, Journal of Quantitative Spectroscopy & Radiative Heat Transfer 109 (2008) 476–493. [6] S. Vakili, M.S. Gadala, Effectiveness and efficiency of Particle Swarm Optimization technique in inverse heat conduction analysis, Numerical Heat Transfer Part B: Fundamentals 56 (2) (2009) 119–141. [7] S. Abboudi, E. Artioukhine, Parametric study and optimal algorithm of a simultaneous estimation in two-dimensional inverse heat conduction problem, Inverse Problems in Science and Engineering 16 (4) (2008) 461–482. [8] R.M. Haghighi Golbahar, M. Eghtesad, P. Malekzadeh, D.S. Necsulescu, Twodimensional inverse heat transfer analysis of functionally graded materials in estimating time-dependent surface heat flux, Numerical Heat Transfer Part A 54 (7) (2008) 744–762. [9] C.-H. Huang, H.-C. Lo, A three-dimensional inverse problem in predicting the heat fluxes distribution in the cutting tools, Numerical Heat Transfer, Part A 48 (2005) 1009–1034. [10] J. Zhou, Y. Zhang, J.K. Chen, Z.C. Feng, Inverse estimation of surface heating condition in a three-dimensional object using conjugate gradient method, International Journal of Heat and Mass Transfer 53 (2010) 2643–2654. [11] T. Louloul, E.P. Scott, An inverse heat conduction problem with heat flux measurements, International Journal for Numerical Methods in Engineering 67 (2006) 1587–1616. [12] J.M. Mulcahy, D.J. Browne, K.T. Stanton, F.R. Chang Diaz, L.D. Cassady, D.F. Berisford, R.D. Bengtson, Heat flux estimation of a plasma rocket helicon source by solution of the inverse heat conduction problem, International Journal of Heat and Mass Transfer 52 (2009) 2343–2357. [13] J. Kennedy, R.C. Eberhart, Particle Swarm Optimization, in: Proceedings of the IEEE International Conference on Neural Networks, 1995, pp.1942–1948. [14] R.C. Eberhart, Y. Shi, Evolving artificial neural network, in: Proceedings of the 1998 International Conference on Neural networks and Brain, 1998, Beijing, PL5-PL13. [15] M. Clerc, The swarm and the queen: Towards a deterministic and adaptive Particle Swarm Optimization, in: Proceedings of the Congress on Evolutionary Computation, 1999, pp. 1951–1957. [16] P.N. Suganthan, Particle swarm optimizer with neighborhood operator, in: Proceedings of IEEE International Conference on Evolutionary Computation, 1999, pp. 1958–1962. [17] G. Venter, Particle Swarm Optimization, in: Proceedings of Thirty fourth AIAA/ ASME/ASCE/AHS/ASC Structures Dynamics and Materials Conference, 2002, pp. 22–25. [18] A. Ratnaweera, S.K. Halgamuge, H.C. Watson, Self-organizing hierarchical particle swarm optimizer with time-varying acceleration coefficients, IEEE Transactions on Evolutionary Computation 8 (3) (2004) 240–255. [19] T. Kiink, J.S. Vesterstroem, J. Riget, Particle Swam Optimization with Spatial Particle Extension, in: Proceedings of the IEEE Congress on Evolutionary Computation (CEC2002), 2002, pp. 1474–1479. [20] T. Krink, M. Lovhjerg, The life cycle model: combining Particle Swarm Optimization, genetic algorithms and hill climbers, in: Proceedings of Parallel Problem Solving from Nature VI1 (PPSN 2002), 2002, pp. 621–630. [21] H. Wang, S. Zeng, Y. Liu, W. Wang, H. Shi, G. Liu, Re-diversification Based Particle Swarm Algorithm with Cauchy Mutation, in: Advances in Computation and Intelligence-Second International Symposium, ISICA 2007, Lecture Notes in Computer Science, 2007, vol.4683 LNCS, pp.362–371. [22] M. Pant, R. Thangaraj, A. Abraham, Particle Swarm Optimization using adaptive mutation, in: Proceedings of DEXA 2008, Nineteenth International Conference on Database and Expert Systems Applications, 2008, pp 519–523. [23] X.J. Cai, Z.H. Cui, J.C. Zeng, Y. Tan, Dispersed Particle Swarm Optimization, Information Processing Letters 105 (6) (2008) 231–235.
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Particle Swarm Optimization-based algorithms for solving inverse heat conduction problems of estimating surface heat flux Fung-Bao Liu ⇑ Department of Mechanical and Automation Engineering, I-Shou University, Kaohsiung 84001, Taiwan
a r t i c l e
i n f o
Article history: Received 9 February 2011 Received in revised form 14 August 2011 Available online 21 December 2011 Keywords: Inverse analysis Particle Swarm Optimization Mutation operator
a b s t r a c t An inverse analysis of estimating a time-dependent surface heat flux for a three-dimensional heat conduction problem is presented. A global optimization method known as Particle Swarm Optimization (PSO) is employed to estimate the unknown heat flux at the inner surface of a crystal tube from the knowledge of temperature measurements obtained at the external surface. Three modifications of the PSO-based algorithm, PSO with constriction factor, PSO with time-varying acceleration of the cognitive and social coefficients, and PSO with mutation are carried out to implement the optimization process of the inverse analysis. The results show that the PSO with mutation algorithm is significantly better than other PSO-based algorithms because it can overcome the drawback of trapping in the local optimum points and obtain better inverse solutions. The effects of measurement errors, number of dimensionalities, and number of generations on the inverse solutions are also investigated. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Inverse heat conduction problems are known as back-analysis problems, i.e. the reversal of the cause-effect sequence, in the field of heat transfer analysis. An inverse problem means that some of the initial, boundary conditions or material properties are not fully specified as determined from the measured temperature profiles at some specific locations. The inverse problems in most situations are likely to be ill-posed [1]. Solutions of the inverse problem are very sensitive to the measurement errors, i.e. small errors in the measured data values can produce very large errors in solutions. In general, the uniqueness and stability of an inverse problem solution are not guaranteed. In recent years, the inverse problems have been studied extensively due to their applications in various engineering disciplines. In general, the technique approaches the inverse heat conduction problem as an optimization problem, i.e. the problem is defined as the minimization of a cost function or a fitness function measuring the distance between measurements and predictions. A variety of numerical methods and computational algorithms for solving the inverse heat conduction problem have been developed over the years. Most of the analytical or numerical methods dealing with the inverse heat conduction problems of estimating ⇑ Address: Department of Mechanical and Automation Engineering, I-Shou University, No. 1, Sec. 1, Syuecheng Rd., Dashu District, Kaohsiung 84001, Taiwan. Tel.: +886 7 6577711x3228; fax: +886 7 6578853. E-mail address: fl[email protected] 0017-9310/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2011.12.007
a single surface condition are one- or two-dimensional. Only a few works estimate the unknown sources in three-dimensional problems due to the complexity of geometry or time-consuming computation. Zueco et al. [2] studied the unsteady inverse heat transfer problem in a one-dimensional geometry subjected to a time-dependent incident heat flux by using the numerical network simulation method. Two different types of incident heat flux waveform, including triangular and rectangular, had been determined. The temperature dependency of the thermal properties, conductivity and specific heat were assumed to be polynomial functions in [2]. Huang and Wu [3] investigated the one-dimensional transient inverse heat conduction problem by using a conjugate gradient method for identifying the unknown boundary heat flux. Instead of employing the conventional parabolic differential equations for the heat transfer problem, a hyperbolic heat equation model was used and analyzed. Identification of thermal conductivity and the shape of an inclusion for a two-dimensional solid body were considered by Ardakania and Khodadad [4]. The Particle Swarm Optimization algorithm (PSO), coupled with the boundary elements method, was used to implement this inverse heat conduction problem. Qi et al. [5] applied the multi-phase PSO method to solve the inverse radiation problem. The proposed method adopted the advantages of both two-group PSO and multi-start PSO algorithms by dividing particles into multi-group to increase the diversity of the swarm and introducing different phases with different flying directions and
2063
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
Nomenclature Cp c1 c2 cl cu f fb fw Gbest k L M Mg N n Pbest,i q r r1
specific heat cognitive coefficient social coefficient lower bound of cognitive and social coefficient upper bound of cognitive and social coefficient the fitness function to be minimize the best fitness of particles the worst fitness of particles best position among all the particles in the swarm thermal conductivity length of the tube total number of particles in a swarm total number of generation dimensionality current generation previous best position of particle i heat flux radial coordinate random number in the interval [0, 1]
searching ways. The standard PSO and the stochastic PSO were also carried out to compare with the proposed method for their performance and accuracy. The effectiveness and efficiency of Particle Swarm Optimization technique in inverse heat conduction analysis were investigated by Vakili and Gadala [6]. Three variations of the PSO method, i.e. basic PSO, repulsive PSO, and complete repulsive PSO, were used to solve the boundary inverse heat conduction problem, in one, two, and three dimensions. The results showed that PSO can alleviate some of the stability problems of the classical approaches for solving the inverse heat conduction problems. For the problems studied, some variants of PSO are proven to be more efficient than other algorithms. Abboudi and Artioukhine [7] proposed a numerical approach for simultaneous estimation of two boundary conditions in a twodimensional linear heat conduction problem. The inverse numerical algorithm was based on the iterative regularization method and on the conjugate gradient method. Unknown functions were parameterized in the form of a cubic B-spline. An optimal choice of the descent parameter matrix was used in the iterative regularization method. Haghighi Golbahar et al. [8] studied two-dimensional transient inverse heat conduction problem of functionally graded materials. A combination of the finite element (FE) and differential quadrature (DQ) methods was employed for solving the direct problem. The inverse analysis using conjugate gradient method was transformed into the solution of three problems; namely, the direct, sensitivity, and adjoint problems. These three problems were also solved by the coupled DQ–FE method. The surface heat fluxes on cutting edges of cutting tools were estimated in the three-dimensional inverse heat conduction problem by Huang and Lo [9]. The steepest descent method with an adjoint equation and a general-purpose commercial code were successfully applied in determining the time-dependent local heat fluxes. Although the inverse heat conduction problems have been extensively studied for different applications in the past decades, little work has been done for the inverse numerical algorithms using heat flux measurements in the objective function [10]. Louloul and Scott [11] investigated an inverse bioengineering problem of estimating the time-dependent blood perfusion and the thermal conductance between the probe and the tissue. A methodology in which a combined parameter and function estimation was developed based on the minimization of an objective function containing computed and measured heat fluxes.
r2 ri ro Si T(r, h, z, t) Ti V nþ1 i w nþ1 Xi Y(r, h, z, t) z
random number in the interval [0, 1] inner radius of tube outer radius of tube weight factor computed temperature initial temperature of the fluid current velocity of particle i inertia weight current position of particle i measured temperatures axial coordinate along the tube
Greek symbols v constriction factor q density h angular coordinate
In this work, an inverse analysis for the recovery of time-varying surface heat flux in three-dimensional cylindrical coordinates is investigated. Transient temperature measurements at the boundary, obtained from the solution of the direct problem, served as the simulated experimental data needed as input for the inverse analysis. A fitness function which is defined by the quadratic residual between the measurements and the calculated temperatures is minimized. Three modifications of the PSO-based algorithms, namely PSO with constriction factor (PSOC), PSO with time-varying coefficient (PSOTV), and PSO with mutation operator (PSOM), are carried out to implement the optimization process of the inverse analysis. The inclusion of mutation operator in the Particle Swarm Optimization algorithm which speeds up the convergence and escape from local minimum is also examined. Numerical examples are provided to compare the effectiveness of the PSOC, PSOTV, and PSOM algorithms. 2. The direct heat transfer problem The thermal behavior of the high power helicon source used to produce plasma in the Variable Specific Impulse Magnetoplasma Rocket (VASIMR™) was investigated by Mulcahy et al. [12] numer-
r z
θ
The external surface is insulated
Heat flux q(z,θ,t) is applied at the interior surface Fig. 1. Physical model.
2064
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
ically and experimentally. The same physical model was used for the direct problem in this study as shown in Fig. 1. To illustrate the methodology for developing expressions in determining the unknown surface heat flux in a circular channel, the following cylindrical geometry with constant thermo-physical properties is considered. The initial temperature is equal to Ti. When the time t > 0, the boundary conditions at boundary z = 0, z = L, and the external surface r = ro are assumed to be insulated, while an unknown heat flux q(z, h, t) is applied at the inner surface r = ri. The dimensional mathematical formulation of this transient heat conduction problem in cylindrical coordinates is given as follows [12]:
" # @T @ 2 T 1 @T 1 @ 2 T @ 2 T ; qC p ¼ k 2 þ þ þ @t @r r @r r2 @h2 @z2
ð1aÞ
and subjected to the following boundary and initial conditions.
k
@Tðr; h; z; tÞ ¼ 0; @r
at r ¼ r o ;
0 6 z 6 L;
0 < h < 2p;
t > 0; ð1bÞ
k
@Tðr;h;z;tÞ ¼ qðh;z;tÞ; at r ¼ ri ; 0 6 z 6 L; 0 < h < 2p; t > 0; @r ð1cÞ
k
@Tðr; h; z; tÞ ¼ 0; @z
at z ¼ 0;
and z ¼ L;
ri 6 r 6 ro ;
t > 0; ð1dÞ
Tðr; z; h; tÞ ¼ T i ;
at 0 6 z 6 L;
ri 6 r 6 ro ;
0 < h < 2p;
t ¼ 0: ð1eÞ
where k, q, Cp and t are the thermal conductivity, density, heat capacity and time, respectively. The problem defined by Eqs. (1a)–(1e) is called a direct problem in which the interior temperature T(r, z, h, t) is determined when the thermal properties, initial and boundary conditions are known. The direct problem can be solved by discretizing the above equations with the finite difference method, the finite element method, or the boundary element method. In this work, the weighted Crank–Nicolson implicit finite difference method is used to discretize the Eqs. (1a)–(1e) and solve the direct problem. 3. The inverse heat conduction problem For the inverse analysis of heat conduction problem, the geometric features and/or material properties are unknown, but some of the unknown boundary data can be measured and used as additional information to estimate the unknown input parameters. In this case, the boundary heat flux at r = ri is regarded as an unknown function, but other parameters in Eqs. (1a)–(1e) are known. In addition, temperature measurements at r = ro are considered available. Let temperature measurements taken by sensors at r = ro denoted as Y(ro, h, z, t). The aim of the inverse analysis is to iteratively estimate the unknown heat flux using an optimization procedure which results a negligible difference between measurements taken at the external surface and temperatures computed from the numerical model in Eqs. (1a)–(1e). Therefore, the solution of the present inverse problem is obtained by minimizing the following fitness function f(q) formulated as a sum-of-squares of residuals: Z t f Z L Z 2p min f ðqðh;z;tÞÞ ¼ ½Tðro ;h;z;tÞ Yðr o ;h;z;tÞ2 dh dz dt: ð2Þ 0
0
0
The inverse problem is recast as an optimization problem. A variety of numerical and analytical techniques have been developed to solve the optimization problems. Among them, the gradi-
ent-based methods have the advantages of the fast convergence speed and the high accuracy of the estimated solutions. However, a good initial guess solution is needed to avoid trapping into local optimum. Another approach is the stochastic methods, i.e. genetic algorithms (GA), simulated annealing (SA), and Particle Swarm Optimization algorithms (PSO). The advantages of the stochastic methods are their capability in searching for the global optimum and no computation of the complicated gradients. However, the time-consuming calculation for obtaining optimal solutions is a drawback of the stochastic methods. 3.1. The Particle Swarm Optimization algorithm The Particle Swarm Optimization algorithm introduced by Kennedy and Eberhart [13] in 1995 is a stochastic optimization technique which draws inspiration from the social behavior of a flock of birds or the collective intelligence of a group of social insects with limited individual capabilities. The standard PSO model consists of a swarm of M particles moving in a problem search space. Each particle is a potential solution of the global optimum over a given domain D. For a N-dimensional search space, the position of the ith particle is represented as Xi = (xi1, xi2, . . . , xiN). At each generation, the new particle position is found by adding a displacement to the current position where the displacement is the particle velocity multiplied by a time step of one as shown in Eq. (3)
X nþ1 ¼ X ni þ V nþ1 : i i
ð3Þ
In Eq. (3), X nþ1 and X ni represent the current and previous positions i of the particle i, V nþ1 is the current velocity of the particle i and is i represented as V nþ1 ¼ ðmi1 ; mi2 ; . . . ; miN Þ. The velocity of each particle i is also updated at each generation and is given by
V nþ1 ¼ V ni þ c1 r 1 ðP best;i X ni Þ þ c2 r2 ðGbest X ni Þ; i V nþ1 i
ð4Þ
V ni
where and are the current and previous velocities of the particle i, respectively. Each particle maintains a memory of its previous best position, say Pbest,i = (pi1, pi2, . . . , piN), where the position giving the best fitness function value. The best one among all the particles in the swarm is represented as the global best position, say Gbest = (pg1, pg2, . . . , pgN). The new velocity in Eq. (4) can be seen as the sum of three parts. The first part of Eq. (4) represents the previous velocity and is called the momentum part. The second part of the Eq. (4) represents the tracking of best position for individual particle and is called the cognition part. The third part of the Eq. (4) represents the cooperation among particles in the swarm and is called the social component. The cognitive learning coefficient c1 and social learning coefficient c2 are known as accelerating factors, and r1 and r2 are two random numbers generated by the uniform distribution within 0 and 1. The relative sizes of these components determine their contribution to the new particle velocity. The most common setting for c1 and c2 are 2.0 for the standard PSO algorithm. 3.2. The modifications of the Particle Swarm Optimization One difficulty encountered in the standard PSO algorithm is that it could easily fall into local optimum in many optimization problems. To well control the scope of the search and improve the performance of the algorithms, many modifications on the standard PSO algorithm have been proposed. The inclusion of inertia weight in the PSO algorithm was investigated by Eberhart and Shi [14]. The inertial weight, w, plays a multiplying factor of current velocity as follows:
V nþ1 ¼ wV ni þ c1 r 1 ðPbest;i X ni Þ þ c2 r2 ðGbest X ni Þ: i
ð5Þ
2065
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
Better selection of the inertial weight provides a balance between global and local exploration and exploitation. Therefore, less iteration is needed to find the optimal solution and improve the performance of the algorithm. A constriction factor was incorporated into the PSO algorithm by Clerc’s suggestion [15] to insure the convergence of the algorithm (PSOC). The velocity term was update as follows:
V nþ1 ¼ vðV ni þ c1 r 1 ðPbest;i X ni Þ þ c2 r 2 ðGbest X ni ÞÞ: i
ð6Þ
When Clerc’s constriction method is used, the constant multiplier v is approximately equal to 0.7298 and the two coefficients c1 and c2 are 2.05. The time-varying acceleration of the cognitive and social coefficients in the algorithm (PSOTV) was investigated by Suganthan [16]. In his strategy, the upper and lower bounds of the cognitive and social coefficients were predefined and both coefficients, c1 and c2, were linearly decreased over the number of generations in the search. That is, the cognitive coefficient c1 and social coefficient c2 are dynamically changed as follows:
c1 ¼ c2 ¼ ðcu cl Þ
Mg n þ cl ; Mg
ð7Þ
where n is the current generation, Mg means the total number of generations, cu is the upper bound and cl is the lower bound of the coefficients. For the modifications of the cognitive coefficient c1 and the social coefficient c2, Venter [17] concluded that better quality solutions could be obtained with the small cognitive coefficient and large social coefficients. Furthermore, Ratnaweera et al. [18] proposed another time-varying acceleration coefficient strategy in which the cognitive coefficient was reduced linearly over the number of generations in the search while the social coefficient was linearly increased inversely. In other words, the cognitive coefficient was changed as the Eq. (7), while the social coefficient was changed as follows:
c2 ¼ ðcu cl Þ
n þ cl : Mg
ð8Þ
3.3. The Particle Swarm Optimization with mutation (PSOM)
ð9Þ
where the weight factor Si is based on the fitness function of each particle and given by
( Si ¼
fw f ðX i Þ ; fw fb
if f b –fw ;
1;
otherwise;
mij ¼
0:5 xmax r 1 ;
if r2 < 0:5;
0:5 xmax r 1 ; otherwise;
ð11Þ
where r1 and r2 are two random numbers generated with uniform distribution between 0 and 1. The value of xmax is the maximum number in the search space. In each generation, the number of particles and the number of dimensionalities applied with mutation strategy are predefined and the mutation probability can be calculated. 3.4. Computation procedure for the Particle Swarm Optimization with mutation The computational steps of the PSOM algorithm described above are given as follows: Step 1: Generate the initial particles in a swarm by randomly generating the position and velocity for each particle. Step 2: Evaluate the fitness function of each particle. Step 3: Update the Pbest,i for each particle, if its fitness is smaller than the fitness of its previous best position (Pbest,i). Step 4: Update the Gbest, if the fitness function of a particle is smaller than the fitness of the best position of all particles (Gbest). Step 5: Compute the cognitive coefficient c1 by Eq. (7). Step 6: Compute the social coefficient c2 of each particle according to Eqs. (9) and (10). Step 7: Update each particle according to Eqs. (5) and (3). Step 8: Apply mutation operator according to Eq. (11). Step 9: Repeat the loop until stopping criteria or predefined number of generations is reached. 4. Numerical examples and discussion
The method of PSO sometimes suffers from a lack of diversity among particles, which leads to a stagnation stage. A number of modifications have been proposed to overcome this drawback, with varying degrees of success [19,20]. One of the modifications is to incorporate mutation operators into the PSO algorithm. Some research works have shown that mutation operators can help the PSO algorithm to prevent the premature convergence [21,22]. In this work, a PSO algorithm with a modified social coefficient and a mutation strategy proposed by Cai et al. [23] is employed for solving the inverse heat conduction problem. To overcome the decrease of searching efficiency for the PSO algorithm, a dispersed social coefficient setting is introduced to the PSO as follows:
c2i ¼ ðcu cl Þ Si þ cl ;
is, the larger the Si is. Furthermore, a mutation strategy is introduced to enhance the ability of escaping from the local optimum to avoid premature convergence. At each generation, particle i is uniformly random selected within the whole swarm and the dimensionality j of the selected particle is also random selected to perform the mutation. Then, the velocity vij is updated as follows:
ð10Þ
where fw and fb are the worst and the best fitness of particles in the swarm. The weight factor Si is an information index to represent the difference of the particles, according to their fitness values of the current position to the best position. The better the particle
A solution method for solving the inverse heat conduction problem by the Particle Swarm Optimization algorithms are carried out to determine the unknown surface heat flux in cylindrical coordinates. In order to give some idea on the physical significance of the various variables, a waste heat flux from a helicon plasma discharge in a rocket is considered. The channel is a quartz tube with density 2200 kg/m3, conductivity 1.4 W/m K, and specific heat 670 J/kg K. The length of the tube is L = 0.3 m with inner and outer diameters being ri = 0.045 m and ro = 0.0475 m, respectively [12]. The exact function of the time-dependent heat flux on the inner surface of quartz tube is given as:
2z2 2z qðz;h;tÞ ¼ qmax 2 þ þ 0:5Þð0:7 0:3cosðhÞ ; for 0:1 6 t 6 0:35; L L ð12aÞ and
qðz; h; tÞ ¼ 0;
for 0:1 < t;
t > 0:35:
ð12bÞ
Two different functional forms of the heat flux are tested. One is a step functional form with the maximum heat flux, qmax of 100 k W/m2. The other is a ramp change in which the heat flux increases linearly from zero at t = 0.1 s to the maximum at t = 0.225 s, before decreasing linearly to zero at t = 0.35 s.
2066
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
For present study, the numbers of grid points in r, z, and h directions are taken as 10, 60, and 40, respectively. The measured temperature locations are at the grid points of the external surface. The total grid number on the external surface is the number of tangential nodes time axial nodes, say 40 60. The measured time period is 0.01 s and the total measured time is 0.45 s, i.e. there are 45 time steps. Therefore, there are a total of 108,000 temperature measurements at the external surface with correspondent unknown heat fluxes at the inner surface in this study. The measurements used as input data for the inverse analysis are obtained from the solution of the direct problem. Note that real measurements generally contain errors. The recordings of temperatures at discrete regular times are modified systematically by introducing errors from a random number generator. Thus random errors are imposed onto the calculated exact temperatures as follows:
Y ¼ Y exact þ xr;
ð13Þ
where Yexact is the solution of the direct problem with the exact boundary heat flux q(z, h, t); r is the standard deviation of the measurements; and x is a random variable within 2.576 and 2.576 for a 99% confidence bound. The inverse analysis is implemented by the Particle Swarm Optimization algorithm with different modifications. Three modified algorithms used in this work are PSO with constriction factor (PSOC), PSO with time-varying acceleration coefficient (PSOTV), and PSO with mutation (PSOM). The standard PSO is also carried out for comparison. The unknown heat flux is represented by a particle i at a position Xi with velocity Vi. The number of particles in the swarm is M = 50 for this study. The initial position and velocity are randomly generated for each particle. Therefore, the position Xi is a guessed heat flux and also a potential solution for the optimal solution. For each particle, the fitness function, i.e. Eq. (2), is evaluated in which the computed temperatures are obtained by solving the direct problem with the guessed solution Xi. The particle previous best position (Pbest,i) and the best position of all particles (Gbest) are updated according to their fitness functions. The update velocity is calculated for each particle to update the new position for the next iteration. In the PSOC algorithm, the constriction factor v is equal to 0.7298. Both the cognitive learning coefficient c1 and the social learning coefficient c2 are equal to 2.05. The inertia weight w decreases linearly from 0.9 to 0.4 and the coefficient c1 decreases from 25 to 5, while the coefficient c1 increases from 5 to 25 in the
PSOTV algorithm. The mutation rates of the particle in the swarm and the dimensionality of the selected particle are 4 and 2 22% in the PSOM algorithm respectively. Fig. 2 shows the estimated results of the step heat flux by using standard PSO, PSOC, PSOTV, and PSOM algorithms at location A(0.15, 0°). The PSO-based algorithms all have the same size of swarm and they run up to 2000 generations. The inverse solutions have the best predictions by using PSOM algorithm and the worst by using PSOC algorithm. The same conclusion can be drawn in the Fig. 3 for the estimated heat flux with ramp functional form at location C(0.05, 90°). Fig. 4 shows the comparison of fitness functions of the algorithms standard PSO, PSOC, PSOTV, and PSOM at location C(0.05, 90°) for ramp heat flux. It is observed that the values of fitness function become stable after a number of generations for the cases of standard PSO, PSOC and PSOTV algorithms. The algorithms even proceed with more generations, but improvements in the results are limited. The performance of the fitness function of PSOM is the best among the three algorithms. The prediction of the heat flux is acceptable at the 2000 generations for the PSOM algorithm. Therefore, with regard to Fig. 4 and to avoid long computational time, 2000 has been chosen as a stopping criterion. In Fig. 4, it is noted that the inverse solutions of the PSOC algorithm have the fast convergence, at about 400 generations, among the three modified algorithms. However, it is trapped into an unforeseen local optimum which provides the wrong solutions as been shown in Figs. 2 and 3. The inverse solutions of the PSOTV algorithm have been improved from the PSOC algorithm, but the improvement is not noticeable after 1200 generations and the solution is also trapped into another local optimum. The social coefficient c2 is an important parameter to control the weight of each particle associated with the best position. The modification by using the time-varying acceleration to the social coefficient c2 in the PSOTV algorithm is not good enough to prevent the premature convergence. However, the inclusion of the fitness function into the social coefficient c2 has improved the search ability for better solutions of the PSOM algorithm. Furthermore, the mutation operator in the PSOM algorithm maintains the diversity of global search capability and avoids the premature convergence as been shown in Figs. 2 and 3. The estimated results of the step heat flux at different locations, A(0.15, 0°), B(0.1, 45°), and C(0.05, 90°) using the PSOM algorithm are presented in Fig. 5. The heat fluxes have the same step functional form with different maximum values at the three different
200000 160000 120000
exact PSOC PSOTV PSOM PSO
60000
80000
heat flux (W/m2)
2
heat flux (W/m )
80000
exact PSOC PSOTV PSOM PSO
40000 0
40000
20000
0
-40000 -20000
-80000 -120000
-40000
0
0.1
0.2
0.3
0.4
0.5
0.6
time (s) Fig. 2. Comparison of estimations by standard PSO, PSOC, PSOTV, and PSOM algorithms for a heat flux with step functional form at location A(0.15, 0°).
0
0.1
0.2
0.3
0.4
0.5
time (s) Fig. 3. Comparison of estimations by standard PSO, PSOC, PSOTV, and PSOM algorithms for a heat flux with ramp functional form at location C(0.05, 90°).
2067
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
120000
1.0E+09
A-exact A-estimated B-exact B-estimated C-exact C-estimated
PSOC 100000
PSOTV PSOM
80000
PSO
heat flux (W/m2 )
fitness function f
1.0E+08
1.0E+07
1.0E+06
60000 40000 20000
1.0E+05 0
1.0E+04 0
500
1000 1500 No. of generation
2000
-20000
2500
0
0.1
0.2
0.3
0.4
0.5
time (s)
Fig. 4. Convergence of the fitness functions for a ramp heat flux by standard PSO, PSOC, PSOTV, and PSOM algorithms at location C(0.05, 90°).
Fig. 6. Inverse solutions of ramp heat fluxes by PSOM algorithm at locations A(0.15, 0°), B(0.1, 45°), and C(0.05, 90°).
150000
2
heat flux (W/m )
120000 90000
60000 30000
160000 exact 2000 1000 500
120000 heat flux (W/m2)
A-exact A-estimated B-exact B-estimated C-exact C-estimated
80000 40000 0 -40000
0
-80000
-30000
0
0
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
time (s)
time (s) Fig. 5. Inverse solutions of step heat fluxes by PSOM algorithm at locations A(0.15, 0°), B(0.1, 45°), and C(0.05, 90°).
locations. In Fig. 5, the results show that the estimations are in good agreement with the exact solutions over the time interval. The same results can be observed in Fig. 6 for the ramp heat flux. The Particle Swarm Optimization algorithm employed a stochastic search for the inverse solutions. Unlike the gradient-based method, no gradient is needed to be calculated in this case. Therefore, there is no difficulty in predicting the discontinuities at both step and ramp functional forms of heat flux. Predictions with three different number of generations, namely 500, 1000, and 2000, at location A(0.15, 0°) are illustrated in Fig. 7 for the PSOM algorithm. The results show a significant improvement in the solutions from 500 to 2000 generations. The weight factor Si of the social coefficient c2 has the ability to fine tune the coefficient in the PSOM algorithm. The mutation operator helps the algorithm to escape from the local optimum. As shown in Fig. 3, the fitness function is continuing improved over the range of this study. The noise is added to the exact temperature data with varying amplitude according to Eq. (13). Fig. 8 presents the estimated heat flux q(h, z, t) at location B(0.1, 45°) for three values of random noise r = 0, 1, and 2 in the simulated temperatures. As expected, the heat
Fig. 7. Estimated results for different number of generations by PSOM algorithm at location A(0.15, 0°).
flux fluctuation increases with the increase of noise in the measured temperatures. For the case of r = 2, the perturbations are randomly scattered around the exact profile which means that the estimation procedure is able to recover the general tendency of the time-dependent heat flux but fails to recover it precisely due to the inherent ill-posed character of inverse problems. It is noted that the characteristic of the PSOM algorithm makes it possible to obtain acceptable and stable results for different amplitude of the noises imposed on the measured temperatures. Estimated results of the heat flux by using the PSOM algorithm with different number of particles in a swarm at location A(0.15, 0°) are presented in Fig. 9. Three different number of particles in the swarm used here are M = 50, 30, and 20 with the same 2000 generations for computations. As expected, better predictions have been found for the case with the larger number of particles M = 50 in the swarm. More number of particles in the swarm can increase the convergence speed of computations. It is noticed that they all have the similar performance if the same fitness function is used as the stopping criterion. However, the case with M = 30 particles in the swarm takes about 4300 generations to reach the same fitness function as the case with 50 particles and about 6100 generations for the case of M = 20 particles in the swarm.
2068
F.-B. Liu / International Journal of Heat and Mass Transfer 55 (2012) 2062–2068
160000
120000
References
2
heat flux (W/m )
PSOM algorithm to recover the inner surface heat flux without using any a priori information of the unknown transient functions.
exact σ=0 σ=1 σ=2
80000
40000
0
-40000
0
0.1
0.2
0.3
0.4
0.5
time (s) Fig. 8. Estimated results for measurements with errors by PSOM algorithm at location B(0.1, 45°).
120000 exact 100000
M=50 M=30
heat flux (W/m2)
80000
M=20
60000 40000 20000 0 -20000 0
0.1
0.2
0.3
0.4
0.5
time (s) Fig. 9. Estimated results for different number of particles in the swarm by PSOM algorithm at location A(0.15, 0°).
5. Conclusion An inverse analysis using Particle Swarm Optimization algorithms has been presented to identify different types of surface heat flux in a three dimensional heat conduction problem. The estimated results by using Particle Swarm Optimization with mutation have been compared against the Particle Swarm Optimization with the constriction factor and the Particle Swarm Optimization with the time-varying acceleration coefficient. Computational results of the PSOM algorithm have been shown to improve performance for the Particle Swarm Optimization algorithms. Simulated results with or without additive noise have shown that the PSOM algorithm allows identifying the surface heat flux with good accuracy as compared with the exact solutions. The obtained results underline the feasibility of the procedure and its capabilities for the
[1] J.V. Beck, B. Blackwell, C.R.St. Clair Jr, Inverse Heat Conduction, Wiley, New York, 1985. [2] J. Zueco, F. Alhama, C.F. Gonzalez Fernandez, Numerical nonlinear inverse problem of determining wall heat flux, Heat Mass Transfer 41 (2005) 411–418. [3] C.-H. Huang, H.-H. Wu, An inverse hyperbolic heat conduction problem in estimating surface heat flux by the conjugate gradient method, Journal of Physics D-Applied Physics 39 (2006) 4087–4096. [4] M.D. Ardakania, M. Khodadad, Identification of thermal conductivity and the shape of an inclusion using the boundary elements method and the Particle Swarm Optimization algorithm, Inverse Problems in Science and Engineering 17 (7) (2009) 855–870. [5] H. Qi, L.M. Ruan, M. Shi, W. An, H.P. Tan, Application of multi-phase particle Swarm Optimization Technique to inverse radiation problem, Journal of Quantitative Spectroscopy & Radiative Heat Transfer 109 (2008) 476–493. [6] S. Vakili, M.S. Gadala, Effectiveness and efficiency of Particle Swarm Optimization technique in inverse heat conduction analysis, Numerical Heat Transfer Part B: Fundamentals 56 (2) (2009) 119–141. [7] S. Abboudi, E. Artioukhine, Parametric study and optimal algorithm of a simultaneous estimation in two-dimensional inverse heat conduction problem, Inverse Problems in Science and Engineering 16 (4) (2008) 461–482. [8] R.M. Haghighi Golbahar, M. Eghtesad, P. Malekzadeh, D.S. Necsulescu, Twodimensional inverse heat transfer analysis of functionally graded materials in estimating time-dependent surface heat flux, Numerical Heat Transfer Part A 54 (7) (2008) 744–762. [9] C.-H. Huang, H.-C. Lo, A three-dimensional inverse problem in predicting the heat fluxes distribution in the cutting tools, Numerical Heat Transfer, Part A 48 (2005) 1009–1034. [10] J. Zhou, Y. Zhang, J.K. Chen, Z.C. Feng, Inverse estimation of surface heating condition in a three-dimensional object using conjugate gradient method, International Journal of Heat and Mass Transfer 53 (2010) 2643–2654. [11] T. Louloul, E.P. Scott, An inverse heat conduction problem with heat flux measurements, International Journal for Numerical Methods in Engineering 67 (2006) 1587–1616. [12] J.M. Mulcahy, D.J. Browne, K.T. Stanton, F.R. Chang Diaz, L.D. Cassady, D.F. Berisford, R.D. Bengtson, Heat flux estimation of a plasma rocket helicon source by solution of the inverse heat conduction problem, International Journal of Heat and Mass Transfer 52 (2009) 2343–2357. [13] J. Kennedy, R.C. Eberhart, Particle Swarm Optimization, in: Proceedings of the IEEE International Conference on Neural Networks, 1995, pp.1942–1948. [14] R.C. Eberhart, Y. Shi, Evolving artificial neural network, in: Proceedings of the 1998 International Conference on Neural networks and Brain, 1998, Beijing, PL5-PL13. [15] M. Clerc, The swarm and the queen: Towards a deterministic and adaptive Particle Swarm Optimization, in: Proceedings of the Congress on Evolutionary Computation, 1999, pp. 1951–1957. [16] P.N. Suganthan, Particle swarm optimizer with neighborhood operator, in: Proceedings of IEEE International Conference on Evolutionary Computation, 1999, pp. 1958–1962. [17] G. Venter, Particle Swarm Optimization, in: Proceedings of Thirty fourth AIAA/ ASME/ASCE/AHS/ASC Structures Dynamics and Materials Conference, 2002, pp. 22–25. [18] A. Ratnaweera, S.K. Halgamuge, H.C. Watson, Self-organizing hierarchical particle swarm optimizer with time-varying acceleration coefficients, IEEE Transactions on Evolutionary Computation 8 (3) (2004) 240–255. [19] T. Kiink, J.S. Vesterstroem, J. Riget, Particle Swam Optimization with Spatial Particle Extension, in: Proceedings of the IEEE Congress on Evolutionary Computation (CEC2002), 2002, pp. 1474–1479. [20] T. Krink, M. Lovhjerg, The life cycle model: combining Particle Swarm Optimization, genetic algorithms and hill climbers, in: Proceedings of Parallel Problem Solving from Nature VI1 (PPSN 2002), 2002, pp. 621–630. [21] H. Wang, S. Zeng, Y. Liu, W. Wang, H. Shi, G. Liu, Re-diversification Based Particle Swarm Algorithm with Cauchy Mutation, in: Advances in Computation and Intelligence-Second International Symposium, ISICA 2007, Lecture Notes in Computer Science, 2007, vol.4683 LNCS, pp.362–371. [22] M. Pant, R. Thangaraj, A. Abraham, Particle Swarm Optimization using adaptive mutation, in: Proceedings of DEXA 2008, Nineteenth International Conference on Database and Expert Systems Applications, 2008, pp 519–523. [23] X.J. Cai, Z.H. Cui, J.C. Zeng, Y. Tan, Dispersed Particle Swarm Optimization, Information Processing Letters 105 (6) (2008) 231–235.