# Comparative Performance of Surrogate-Assisted MOEAs for Geometrical Design of Pin-Fin Heat Sinks

Authors: Siwadol Kanyakam; Sujin Bureerat
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

This paper presents the comparative performance of several surrogate-assisted multiobjective evolutionary algorithms (MOEAs) for geometrical design of a pin-fin heat sink (PFHS). The surrogate-assisted MOEAs are achieved by integrating multiobjective population-based incremental learning (PBIL) with a quadratic response surface model (QRS), a radial-basis function (RBF) interpolation technique, and a Kriging (KRG) or Gaussian process model. The mixed integer/continuous multiobjective design problem of PFHS with the objective to minimise junction temperature and fan pumping power simultaneously is posed. The optimum results obtained from using the original multiobjective PBIL and the three versions of hybrid PBIL are compared. It is shown that the hybrid PBIL using KRG is the best performer. The hybrid PBILs require less number of function evaluations to surpass the original PBIL.


---

## Body

## 1. Introduction

An air-cooled heat sink is one of the most effective and popular cooling devices for electronic packages due to its high reliability, simplicity, safety, and low cost. It is operated in such a way that the heat sink is attached to a device that needs to be cooled down such as a central processing unit (CPU) chip. Then, heat dissipation is achieved by using a fan to generate air flow through the heat sink. As a result, design of such a cooling device needs to find an optimum geometry such that cooling performance is optimised. The minimisation of the temperature between the heat sink base and the electronic device, which is called junction temperature, is one of the most used design objectives for maximising thermal performance [1–5]. Nevertheless, in heat transfer design, increasing thermal performance usually results in an increase of pressure drop across the cooling device. This requires high pumping power and consequently high operating cost. Therefore, a design process of a heat sink usually has two design objectives that are junction temperature and fan pumping power for this study.The use of multiobjective optimisers for geometrical and sizing design of heat sinks has been investigated in recent years. Our previous work shows that using MOEAs for multiobjective design of pin-fin and plate-fin heat sinks is superior to a classical design approach [1–5]. Other work related to numerical simulation and optimum design of some types of pin-fin heat sinks can be found in [6–12]. Nevertheless, it is also found that the optimisation process is time-consuming due to expensive function evaluations from performing finite volume analysis (FVA); therefore, the performance enhancement of the MOEA design process is always required. The use of a radial-basis function response surface model (RSM) in combination with a strength Pareto evolutionary algorithm for plate-fin heat sink design has been studied [1]. It was found that, with the inclusion of the RBF surrogate model in the evolutionary design process, the performance of SPEA is greatly increased. In cases of pin-fin heat sink design, investigation on the use of integrated Kriging model and PBIL for solving a heat sink design problem has been made [5]. It is found that the hybrid optimiser performance is acceptable.This work is aimed at performance enhancement of a multiobjective evolutionary optimiser by incorporating evolutionary search with a number of surrogate models. Although it is well known in the past that using a surrogate model can improve search performance of MOEAs, the performance test of several surrogate models for a particular design case still needs to be studied since the nature of design problems can affect surrogate model prediction performance. The design problem is posed to find geometrical parameters of a PFHS to minimise its junction temperature and fan pumping power. The population-based incremental learning is chosen as the optimiser since it gives the best results as shown in [2, 4]. The design process has been refined from the work presented in [5]. Surrogate models employed with PBIL are a Kriging method or Gaussian process model [13], a radial-basis function [1], and a quadratic response surface model. The comparative results show that the performance of a surrogate-assisted PBIL is superior to the nonsurrogate PBIL although the former uses lower number of function evaluations. The best surrogate model for this type of design problem is KRG.


## 2. Surrogate-Assisted PBIL

A particular multiobjective optimisation problem can be expressed as follows.(2.1)Minimisef={f1(x),…,fk(x)}subjecttogi(x)≤0,
where x∈[xl,xu] is a vector of n design variables and gi(x) are design inequality constraints. k≥2 objective functionsfi(x)are intended to be minimised simultaneously. Figure 1 shows an arbitrary biobjective minimisation problem. The set of feasible design solutions are plotted in the objective domain (f1 and f2 coordinates). The optimum points are located on the front as shown. This is called a nondominated or Pareto front of the design problem where the members of the Pareto front are said to be equally good in optimisation point of view.Figure 1
Pareto optimal set.PBIL is an evolutionary optimiser based upon binary search space proposed by Baluja [14] in 1994, whereas it is extended for multiobjective cases by Bureerat and Sriworamas [15]. For single objective case, the PBIL search procedure starts with an initial probability vector (Pi) in which each element in the vector determines the probability of having “1” on each column of a binary population given that each row in the population is one design solution. The samples of a probability vector are shown in Table 2. It is shown that one probability vector can result in a variety of populations, and this feature of PBIL leads to impressive population diversity reproduction in a multiobjective optimisation design process [15].Having obtained a probability vector from the binary population, the best binary design solution is then used to modify the probability vector using the following relation:(2.2)Pinew=Piold(1-LR)+biLR,
wherebiis the element of the best binary solution andLRis learning rate which, in this work, is defined as(2.3)LR=0.5+rand⋅(+0.1or-0.1),
where rand ∈[0,1] is a uniform random number.For the multiobjective version of PBIL, a set of probability vectors, in which they are called a probability matrix, are used instead of only one vector so as to produce a more diverse binary population. The search procedure starts with an initial probability matrix, and a Pareto archive for collecting nondominated solutions. A binary population can then be obtained from the union set of sub-populations created from each row of the probability matrix. The nondominated solutions are obtained by sorting the combination of the nondominated solutions from the previous generation and the members in the current population. In case that the number of nondominated solutions exceeds the predefined archive size, a normal line method [2] is employed to remove some solutions from the Pareto archive. The process is repeated until a stopping criterion is fulfilled. More details of multiobjective PBIL can be seen in [2, 15].In this work, the hybridisation of PBIL and a surrogate model is carried out in such a way that training design points and their corresponding objective values are generated using a Latin hypercube sampling (LHS) technique. Then, a surrogate model (KRG, QRS, or RBF) is built. The PBIL is used to find the nondominated front of the design problem where function evaluations are based upon a surrogate model predictor. Afterwards, a few design solutions in the nondominated front are selected by means of a clustering technique [16], whereas their actual function evaluations are performed. The selected binary solutions are set as an initial population, and they are used to generate the PBIL probability matrix. Then, the usual PBIL procedure is activated for a few iterations to improve the nondominated front. From the hybrid strategy, it is expected that the number of function evaluations is lower than that used by a nonsurrogate PBIL. The computational procedure can be separated into three phases and detailed as follows.Phase I: Initialisation(1)
Generate training points using a Latin hypercube technique.(2)
Perform actual function evaluation.(3)
Construct a surrogate model (KRG, QRS, or RBF).Phase II: Optimisation using a surrogate model(1)
Generate an initial probability matrixPij = 0.5 and a Pareto archive A={}.(2)
If a termination condition is satisfied, stop. Otherwise, go to step 3.(3)
Generate a binary populationB according to Pij.(4)
Function evaluation based on surrogate prediction.(5)
Find nondominated solutions ofA∪B, and replace A with those solutions.(6)
IfA is too large, remove some of its members using a normal line technique.(7)
Use members inA to update Pij using (2.2) where LR is computed using (2.3).(8)
Go to 2.Phase III: Optimisation by actual function evaluations(1)
Select some members fromA in the previous phase using a clustering technique, and find their actual function values.(2)
Generate an initial probability matrixPij = 0.5 and a Pareto archive A from sorting the design solution in step 1 of phase III.(3)
UpdatePij using members in A.(4)
If a termination condition is satisfied, stop. Otherwise, go to step 5.(5)
Generate a binary populationB according to Pij.(6)
Actual function evaluation based on finite volume analysis.(7)
Find nondominated solutions ofA∪B, and replace A with those solutions.(8)
IfA is too large, remove some of its members using a normal line technique.(9)
Use members inA to update Pij using (2.2) where LR is computed using (2.3).(10)
Go to 4.


## 3. Surrogate Models

Lety=f(x) be a function of a design vector x sized n×1. Given a set of design solutions (or sampling points) X=[x1,…,xN] generated by using the Latin hypercute technique and their corresponding function values Y=[y1,…,yN], a surrogate or approximation model is constructed by means of curve fitting or interpolation. Approximation models used in this study are as follows.


### 3.1. Quadratic Response Surface Model

The most commonly used polynomial surrogate model or a response surface model is of the second-order polynomial or quadratic model, which can be expressed as(3.1)y̅=β0+∑βixi+∑βixixj,
where βi for i=0,…,(n+1)(n+2)/2 are the polynomial coefficients to be determined. The coefficients can be found by using a regression or least square technique.


### 3.2. Kriging Model

A Kriging model (also known as a Gaussian process model) used in this paper is the famous MATLAB toolbox named design and analysis of computer experiments (DACEs) [13]. The estimation of function can be thought of as the combination of global and local approximation models that is(3.2)y(x)=f̅(x)+Z(x),
wheref̅(x)is a global regression model while Z(x) is a stochastic Gaussian process with zero mean and nonzero covariance representing a localised deviation. In this work, a linear function is used for a global model, which can be expressed as(3.3)f̅=β0+∑i=1nβixi=βTf,
where β=[β0,…,βn]T, f=f(x)=[1,x1,x2,…,xn]T. The covariance of Z(x) is expressed as(3.4)Coν(Z(xp),Z(xq))=σ2R[R(xp,xq)]
for p,q=1,…,N, where R is the correlation function between any two of the N design points and R is the symmetric correlation matrix size N×N with the unity diagonal [17]. The correlation function used herein is(3.5)R(xp,xq)=exp⁡(-(xp-xq)Tθ(xp-xq)),
where θi are the unknown correlation parameters to be determined by means of the maximum likelihood method. Having found β and θ, the Kriging predictor can be achieved as(3.6)y̅(x)=f(x)Tβ+rT(x)R-1(y-Fβ),
where F=[f(x1),f(x2),…,f(xn)]T and rT(x)=[R(x,x1),R(x,x2),…,R(x,xN)]. For more details, see [13].


### 3.3. Redial-Basis Function Interpolation

The redial-basis function interpolation has been used in a wide range of applications such as integration between aerodynamic and finite element grids in aeroelastic analysis [18]. The use of such a model for surrogate-assisted evolutionary optimisation is said to be commonplace [1, 19, 20]. The approximate function can be written as(3.7)y̅(x)=∑i=1NαiK(‖x-xi‖),
where αi are the coefficients to be determined and K is a radial-basis kernel (here it is set to be linear splines). The coefficients can be found from the N sampling points as detailed in [1].


## 3.1. Quadratic Response Surface Model

The most commonly used polynomial surrogate model or a response surface model is of the second-order polynomial or quadratic model, which can be expressed as(3.1)y̅=β0+∑βixi+∑βixixj,
where βi for i=0,…,(n+1)(n+2)/2 are the polynomial coefficients to be determined. The coefficients can be found by using a regression or least square technique.


## 3.2. Kriging Model

A Kriging model (also known as a Gaussian process model) used in this paper is the famous MATLAB toolbox named design and analysis of computer experiments (DACEs) [13]. The estimation of function can be thought of as the combination of global and local approximation models that is(3.2)y(x)=f̅(x)+Z(x),
wheref̅(x)is a global regression model while Z(x) is a stochastic Gaussian process with zero mean and nonzero covariance representing a localised deviation. In this work, a linear function is used for a global model, which can be expressed as(3.3)f̅=β0+∑i=1nβixi=βTf,
where β=[β0,…,βn]T, f=f(x)=[1,x1,x2,…,xn]T. The covariance of Z(x) is expressed as(3.4)Coν(Z(xp),Z(xq))=σ2R[R(xp,xq)]
for p,q=1,…,N, where R is the correlation function between any two of the N design points and R is the symmetric correlation matrix size N×N with the unity diagonal [17]. The correlation function used herein is(3.5)R(xp,xq)=exp⁡(-(xp-xq)Tθ(xp-xq)),
where θi are the unknown correlation parameters to be determined by means of the maximum likelihood method. Having found β and θ, the Kriging predictor can be achieved as(3.6)y̅(x)=f(x)Tβ+rT(x)R-1(y-Fβ),
where F=[f(x1),f(x2),…,f(xn)]T and rT(x)=[R(x,x1),R(x,x2),…,R(x,xN)]. For more details, see [13].


## 3.3. Redial-Basis Function Interpolation

The redial-basis function interpolation has been used in a wide range of applications such as integration between aerodynamic and finite element grids in aeroelastic analysis [18]. The use of such a model for surrogate-assisted evolutionary optimisation is said to be commonplace [1, 19, 20]. The approximate function can be written as(3.7)y̅(x)=∑i=1NαiK(‖x-xi‖),
where αi are the coefficients to be determined and K is a radial-basis kernel (here it is set to be linear splines). The coefficients can be found from the N sampling points as detailed in [1].


## 4. PFHS Design

A pin-fin heat sink herein is an air-cooled heat sink with square heat sink base and square pin fins as shown in Figure2. Heat dissipation is accomplished by using a fan to generate air flow at the top of the heat sink. A PFHS physical model with square base (whs) is given in Figure 2. Parameters defining a heat sink geometry include base length (whs), base thickness tb, fin width (wfin), and number of fins nf. The fins heights are varied based on the 16 input parameters (H1–H16). Figure 3 displays 7×7 control points on the fin base where the 16 fin height parameters are distributed symmetrically. The fin height can be determined by using radial-basis function interpolation. Figure 4 displays a sample of fin height distribution by using a RBF interpolation [3]. A fin-to-fin space or fin stitch (sfin) is equispaced. The heat sink is cooled by a fan of Vfair velocity (m/s) which is installed on the top of a heat sink. A uniform heat load (Q) generated at the junction of heat sink, and an electronic device is set to be 120 watt. Materials and air properties and the assumptions of heat-flow in this simulation are detailed in [1, 2].Figure 2
Physical model of pin-fin heat sink.Figure 3
Fin height distribution control points.Figure 4
Fin heights distribution.In the design process, an automated procedure is achieved by interfacing the CFD software into MATLAB. The design variables include those parameters for creating PFHS geometry and the inlet air velocity. For more details of the encoding/decoding process of the design variables, see [3].The design problem is proposed to minimise biobjective functions as fan pumping power and junction temperature of a pin-fin heat sink. The fan pumping power is used to measure the consumption of electronic power, while the junction temperature determines heat sink thermal performance. Heat sink aspect ratio and space limit of the device are set to be design constraints. The multiobjective design problem of PFHS can be written as follows.Objective functions:(4.1)min⁡x{Tj,Pfan},Tj=Ta+QRhs,Pfan=ṁaΔPρa.Design variable vector:(4.2)x=[nf,wfin,whs,Vf,H1-16,tb]T.Constraints:(4.3)Tj-355°K≤0,0.0025≤wfin≤0.005,4≤Nf≤14,0.0025≤tb≤0.005,0.0≤Hfin≤0.05,0.5≤Vf≤1.0,0.03≤whs≤0.06,0.0025<sfin<0.004,
where Tj is junction temperature, Pfan is fan pumping power, Q is heat load, Rhs is heat sink thermal resistant, ΔPis pressure drop across the heat sink, ρa is air density, andṁais air mass flow rate.The continuity, momentum, and energy equations for steady-state forced convection in a Newtonian, constant properties fluid without heat generation can be written as(4.4)∇⋅V⇀=0,V⇀⋅∇V⇀=-1ρa∇P+ν∇2V⇀,V⇀⋅∇T=α∇2T,
where V⇀is fluid velocity, T is temperature, α is thermal conductivity, and ν is kinematic viscosity.These are the governing equations for thermofluid analysis of the forced convection in the pin-fin heat sink. During the optimisation process, the three equations are handled by means of finite volume analysis, which is one of the most powerful numerical methods for computational fluid dynamics. The assumptions for the analysis are as follows.(i)
Fluid flow being laminar and steady.(ii)
Constant material thermophysical properties of both air and solid material.(iii)
Uniform inlet air velocity.(iv)
Uniform heat flux throughout the entire base plate bottom surface.The heat sink body is made of aluminium. The physical properties of the solid and fluid are given in [1, 2]. The ambient temperature is set to be 298 K.The vector of design variablesx has 21 elements. SI units are used in the design objectives and constraints unless otherwise specified. In this work, 30 sampling points are created using a LHS technique, while PBIL with 100 population size and 200 generations is performed with surrogate-based function evaluations. Then, 25 design solutions from the surrogate optimisation are taken as initial solutions for a common PBIL, which is used for further exploring the nondominated front with 25 population size and 19 generations. This means that the proposed surrogate-assisted PBIL takes 30+25×20=530 actual function evaluations. For the original PBIL without using a surrogate model, the population size is set to be 25 while the number of generations is set to be 35. This implies that the nonsurrogate PBIL uses 25×35=875 actual function evaluations.


## 5. Results and Discussion

Figure5 displays nondominated fronts obtained from the various surrogate-assisted PBILs, which are defined as PBIL-QRS, PBIL-RBF, and PBIL-KRG for hybrid PBILs using a quadratic response surface model, a radial-basis function interpolation technique, and a Kriging model respectively. It can be clearly seen that the front obtained by using PBIL-KRG is superior to the others while the second best is PBIL-RBF.Figure 5
Pareto fronts from surrogate-assisted PBILs.The best front from the hybrid approach is then compared with that obtained from using the nonsurrogate PBIL (termed as PBIL-35) as shown in Figure6. Apart from the two nondominated fronts, real heat sinks currently used in the real world are modelled where there junction temperature and fan pumping power values are calculated. The real heat sinks which have cylindrical pin fins include the following.(i)
Model 17R with 0.0635 m heat sink width, constant 0.0424 m fin height, 0.0033 m fin diameter, 0.0127 m heat sink base thickness, and14×14 fins.(ii)
Model 11M with 0.0508 m heat sink width, constant 0.0269 m fin height, 0.0033 m fin diameter, 0.0048 m heat sink base thickness, and 61 fins.(iii)
Model 14M with 0.0508 m heat sink width, constant 0.0480 m fin height, 0.0033 m fin diameter, 0.0048 m heat sink base thickness, and 61 fins.Figure 6
PBIL-KRG20 versus PBIL35.The first model is termed a dense pin fin heat sink, while the others are called a sparse pin fin heat sink. The models are simulated with an inlet air velocity 0.5 m/s and constant heat load 120 watts being applied at the bottom of the heat sink base. Based on the dominance theory, the design solutions obtained from PBIL-KRG and nonsurrogate PBIL are far superior to those real heat sinks.The zoom-in of Figure6 is given in Figure 7 so as to compare PBIL-KRG and PBIL-35. It can be seen that the results obtained from using the hybrid PBIL is overall superior to its nonsurrogate counterpart. This is confirmed by the quantitative assessment in Table 1. From the table, the hypervolume values of the Pareto fronts plotted in Figures 5 and 7 are computed where the reference point is set to be [303 K, 0.1 Watt]. Based on its definition, the larger hypervolume means the better Pareto or nondominated front. Two hybrid approaches PBIL-KRG and PBIL-RBF can surpass the original PBIL although the later uses 345 more actual function evaluations. Some selected design solutions from the PBIL-KRG front are plotted in Figure 8. It is shown that the hybrid optimiser results in a variety of PFHS geometries.Table 1
Hypervolume comparison (reference point is [303, 0.1]).
OptimiserHypervolumePBIL-KRG200.42131PBIL-RBF200.41935PBIL-QRS200.39659PBIL350.41183Table 2
Probability vectors and their corresponding populations.
Population 1Population 2Population 30010,0101,11011001,1110,00011110,0001,11010101,0000,0110Probability vectors [0.5, 0.5, 0.5, 0.5] [0.25, 0.5, 0.25, 0.5] [0.5, 0.75, 0.25, 0.75].Figure 7
PBIL-KRG20 versus PBIL35 zoom-in of Figure6.Figure 8
Selected PFHSs from the front of PBIL-KRG20.Further investigation is made to examine the accuracy of the employed surrogate models. The 30 sampling points in the previous section is used to build three surrogate models (QRS, RBF, and KRG). The second set of sampling points with 30 design solutions is generated, while its function values evaluated based on the surrogate models and finite volume analysis (actual function values) are determined. The boxplot in Figure9 shows the error as an absolute value of the difference between actual function value and one from surrogate approximation. It can be seen that the most accurate surrogate model for both junction temperature and fan pumping power is KRG while the second best is RBF. This is coincident with the optimum results obtained from the hybrid PBILs.Figure 9
Accuracy assessment of the surrogate models.


## 6. Conclusions

This paper proposes hybrid PBIL using three well-known surrogate models as a quadratic response surface model, a radial-basis function interpolation technique, and a Kriging model. The hybridisation is carried out in such a way that a Latin hypercube sampling technique is employed to generate training points for constructing a surrogate model. Then, PBIL is used to tackle the design problem based upon the surrogate model prediction. Having obtained Pareto solutions from this phase, the Pareto front is then further improved by running PBIL with actual function evaluations for a few generations.The hybrid PBILs along with the nonsurrogate PBIL are implemented on the multiobjective geometrical design of pin-fin heat sinks. The design problem is posed to find PFHS geometries while minimising its junction temperature and fan pumping power. The computational results reveal that the best performer is a surrogate-assisted PBIL using KRG. This is coincident with the accuracy investigation on predicting junction temperature and fan pumping power of the surrogate models. Design solutions obtained from performing the proposed design process are superior to those real heat sinks currently used in the real world. The proposed hybrid strategy is said to be efficient. Our future work is to study alternative hybrid strategy so that MOEA search performance for solving PFHS design is improved.



---

*Source: 534783-2012-05-13.xml*