Optimization of solder joints in embedded mechatronic systems via Kriging-assisted CMA-ES algorithm

. In power electronics applications, embedded mechatronic systems (MSs) must meet the severe operating conditions and high levels of thermomechanical stress. The thermal fatigue of the solder joints remains the main mechanism leading to the rupture and a malfunction of the complete MS. It is the main failure to which the lifetime of embedded MS is often linked. Consequently, robust and inexpensive design optimization is needed to increase the number of life cycles of solder joints. This paper proposes an application of metamodel-assisted evolution strategy (MA-ES) which signi ﬁ cantly reduces the computational cost of ES induced by the expensive ﬁ nite element simulation, which is the objective function in optimization problems. The proposed method aims to couple the Kriging metamodel with the covariance matrix adaptation evolution strategy (CMA-ES). Kriging metamodel is used to replace the ﬁ nite element simulation in order to overcome the computational cost of ﬁ tness function evaluations ( ﬁ nite element model). Kriging is used together with CMA-ES and sequentially updated and its ﬁ delity (quality) is measured according to its ability in ranking of the population through approximate ranking procedure (ARP). The application of this method in the optimization of MS proves its ef ﬁ ciency and ability to avoid the problem of computational cost.


Introduction
At present, embedded mechatronic systems (MSs) are growing strongly. They play an important role in many areas. For example, the health field with innovative technology, the automotive industry, aeronautics, etc. The reliability assessment of these systems remains a major challenge. However, the failures of these systems are often caused by many extreme solicitations of thermal (temperature variation), and mechanical nature (shocks, vibration), etc. The most failures of mechatronic devices are caused by solder joints fatigue. In order to satisfy the reliability expectations, the component must be reliably optimized, the robustness of the solder joints should be assessed and their connections must be validated.
Thermal cycling test is one of the most commonly used reliability tests to assess thermo-mechanical reliability of the solder joints. The purpose of this test is to investigate the aptitude of solder joints to support mechanical stresses which is a major failure mechanism that can lead to solder joint fatigue failure [1][2][3][4][5]. The Thermo-mechanical stresses induced by alternating thermal cycles [6] are mainly induced by thermal expansion mismatch of different materials in assembly system. Therefore, significant stress can lead to plastic deformation whose accumulation can cause solder joint failure and consequently failure of the complete system. To avoid operational failures and reduce warranty costs, robust optimization and validation of critical solder joints should be carried out in the design and manufacturing process of mechatronic devices. The robustness of this validation includes a design optimization for the purpose to maximize the number of fatigue life cycle of device assembly.
Generally, in the structural design optimization or reliability-based design optimization (RBDO) [7,8], the derivative-based algorithms are the most used methods for RBDO [9], which require derivative computation. The main advantage of those methods, over others, is that they need a much smaller number of iterations to converge to an optimum. Nevertheless, only convergence towards a local minimum is guaranteed. The derivative-free algorithms, which are based solely on original fitness function (FF), are proved as powerful tools for global optimum research and therefore are widely used in real word problem such as engendering optimization. The evolution strategies (ESs) are one among powerful derivative-free algorithms, which are a popular methods for black-box optimization, where no expressions of so-called objective functions, are known and no derivatives can be computed [10]. The use of ESs in the real-life applications proves their power [11]. However, the disadvantage of ES is the large number of FF evaluations required to obtain a satisfying solution; moreover, in engineering application such as MSs, the evaluations of FF are carried out using expensive numerical finite element analysis. Consequently, computational cost produced by a large number of fitness evaluation remains the main challenge in the application of ES.
Metamodel-assisted evolution strategies (MA-ESs) were developed to overcome the expensive numerical finite element simulations. The metamodel is trained to approximate and to be used together with the expensive original FF (finite element code). Various MA-ESs techniques have been presented by several researchers [12,13]. In MES, several metamodels, such as Kriging metamodel, radial basis functions (RBF), support vector machines (SVM), K-nearest neighbor method (KNN), and artificial neural networks (ANN), can be used to approximate the original FF [14]. However, metamodel management and selection, which depend principally on the problem to be addressed, remain the two main issues in the development and application of MES.
In this work, the Kriging-assisted CMAE-ES algorithm is used to optimize the thermo-mechanical performance of solder joints by maximizing the number of life cycles. The rest of the document is structured as follows. In Section 2, CMA-ES and Kriging metamodel, which are the principal components of the algorithm, are shortly presented. Section 3 designed to the management and quality measure of metamodel. In Section 4, the Kriging-assisted CMA-ES (KA-CMA-ES) algorithm using modified approximate ranking procedure is discussed and the numerical tests are provided to show its effectiveness. Section 5 is devoted to the numerical simulation. At first, the numerical global and local model are developed to simulate the thermomechanical behavior of the solder joints [4], and the Ansys TM finite element software is a tool used to modelize and compute the inelastic train due to the thermal variation. The viscoplastic behavior of the solder joints is taken into account using Anand model. Secondly, the KA-CMA-ES algorithm is applied to the global optimization of the MSs, in order to optimize the thermomechanical performance of solder joints by maximizing the number of life cycles. In the 6th and the last section, the paper closes with a short conclusion and perspective on future work. At present, the CMA-ES is one of the best ESs in derivative free optimization. It has become a standard for continuous black box evolutionary optimization. The correlated mutation operator adopted by CMA-ES makes it a powerful optimization algorithm compared with other algorithms that use isotropic mutation. CMA-ES is based on parametrized multivariate normal distribution as a distribution model of a candidate population. At each generation of CMA-ES, l offspring candidates are generated from m parent. In the next generation, the new parent (m) is selected from the offspring candidates (l) using a (l, m)-selection. The m best candidates are selected according to their ranking based on the objective function, CMA-ES (l, m) algorithm uses two techniques namely the cumulative step size adaptation (CSA) and the covariance matrix adaptation (CMA) respectively for adapting the step-size and the covariance matrix. Specifically, let N ðm ðgÞ ; K ðgÞ Þ a multivariate Gaussian distribution with m (g) and K (g) are its mean and covariance matrix. In CMA-ES, the sequence of mean m (g) represents the favorite solution at generation g or the best estimate of the optimum. The covariance matrix K (g) at generation (g) is factorized in two parts: C g ∈ R n Â n and s g ∈ R þ , where C g represents a covariance matrix (a definite positive matrix) and the s g is a scalar factor to control the step length of the iteration. Finally, the l search points (individuals, offspring) are sampled through the following equation: with l ≥ 2 is the population size (number of offspring), and N ð0; C ðgÞ Þ is a multivariate Gaussian distribution with zero mean and covariance matrix C(g).
To define completely the iterative stochastic algorithm of CMA-ES, the three terms denoted m (g) , s (g) and C (g) need to be adapted at each generation.
After sampling process, the l offspring is evaluated and ranked based on the objective function. Finally, the m best individual is selected, and the mean m (g) is then updated by taking the weighted mean of the m individuals where ðw i Þ 1 i m are the m strictly positive and normalized weights, satisfying P m i¼1 ¼ 1. The second term to update is the step size s that uses the so called evolution path p ðgÞ s [15] for its accumulative adaptation. The s computation is given as follow: where EðjjN ð0; IÞjjÞ is the expectation of the Euclidean norm of an N (0, I) distributed random vector, which is approximated by: The last term to update is the covariance matrix C (g) that uses a combination of two mechanisms namely rank À one update and rank À mu update. The update of covariance matrix is given as: where p ðgÞ c ∈ R n is the covariance matrix evolution path [15].
The remaining parameters for CMA-ES update are given in [15].

Kriging metamodel
Kriging metamodel is a geostatistical technique to interpolate deterministic noise free data. The theory of Kriging has been formalized by the mathematician Matheron [16], who was inspired by the work of Krige [17], Kriging has subsequently become a standard method for constructing metamodels for computer experiments [18]. It is then used to predict the value of an expensive FF (FE model). This method is also known as Gaussian process regression.
Let y(x) be defined as a function of x, with x ∈ R d and y is a vector of n observed values of y(x) on D ={ x 1 , … , x n }, a design of experiment (DOE) with dimension n * d. Kriging assumes that the function y(x) is a realization of a Gaussian process denoted by Y(x), which is given as: where h(x) is the mean of the process, and Z(x) is a Gaussian process with zero mean and covariance expressed by: with s 2 is the variance of Gaussian process and R(x (i) , x (j) ) its correlation function between any two samples x (j) and x (j) . The choice of the correlation function is an important element of Kriging. Many covariance functions are proposed in the literature [19], such as Matter, Gaussian, exponential or spherical correlation functions. The Gaussian correlation function is the most commonly used, this last allows control of both the range of influence and the smoothness of the approximation model: For any new point x, the mean and variance of prediction [18] can be respectively calculated by: whereŷðxÞ and s 2 (x) are respectively the estimated mean value and variance ofŷðxÞ, The unknown model parameters can be determined using likelihood estimates (MLEs) method [20].

Metamodel-assisted CMA evolution strategy
In brief, a good stewardship of use of metamodel with prudence to keep it in high quality during optimization process is the main assets in the use of metamodel with evolutionary strategies algorithms. In this section, the management and quality of metamodel will be discussed, and the elements of KA-CMA-ES will be presented.

Metamodel management
In evolutionary computation, the use of metamodel to approximate FF is not as simple as one may assume. There are two main concerns behind using metamodel for fitness evaluation. Firstly, the convergence of the evolutionary algorithm to the global optimum of the original FF should be ensured. Secondly, the computational cost of evolutionary algorithm, which is the principal objectif of metamodel application, should be reduced as much as possible [12,13]. One main issue to replace the FF with a global metamodel is that is very difficult to construct an accurate metamodel that can efficiently approximate the original function, it could happen that the evolution algorithm converges to a false optimum [21]. Therefore, to avoid this problem, the meta-model and the original FF are used side by side; it means that in evolutionary computation the original FF is used to evaluate all individuals in some generation (controlled generation) or some of the individuals (controlled individual) [13]. This technique is the main issue of metamodel management, which is generally divided into three main approaches. The first is the socalled No Evolution control, where the metamodel is assumed to be of high-fidelity and no individual or generation is controlled, i.e., the original FF is not used in evolutionary computation. The second is the Fixed Evolution control, where the frequency of evolution control is fixed. To be done, there are usually two approaches, one is the individual-based control [12,22], the second is generation-based control [12,22]. The last approach is the so-called adaptive Evolution control where the evolution control stands on the quality of the metamodel.

Metamodel quality measure
In addition to reduce the computational cost, the quality or fidelity of metamodel is the main issue in the design and use of metamodels for evolutionary computation. Hence, to measure the fidelity of the metamodel with avoiding the computational cost problem, it is not advisable to use a close quantitative approximation to assess the quality of metamodel. Rather, in evolutionary computation, to be faithful, the meta-model should enable the evolutionary algorithm to select the best individuals in terms of the original FF. In [23], several quality measures have been introduced.
In CMA-ES, the evaluation of objective function is performed in each individual of the generation, and the progression of the algorithm is based on the rank of individuals (rank-based selection), which is the only information to get from the fitness evaluation. In this direction, it is normal that the quality of metamodel is measured according to its ability in ranking the population [24,25]. Fortunately, the approximate ranking procedure [24], which can be classified among adaptive evolution control mechanisms, is one of the best methods to control metamodel quality by its ability in ranking the population without using the exact ranking of the entire population.
In the original approximate ranking procedure, one individual is chosen to be evaluated by original function in each iteration, then the metamodel is updated and the parent set of individuals is chosen based on the updated metamodel. This loop is stopped in the case where all individuals in the generation have been re-evaluated or the parent set does not change. Nevertheless in the case of high-dimensional or multimodal problems, the population size l is large, the amount of information added during each iteration can result in unimportant changes, even in the case of a metamodel with poor ranking predictions. To overcome this deficiency, Huang et al. [26] propose some modifications for the original approximate ranking procedure.
The first modification is to select n init = max(1, ⌊ 0.3l ⌋) individuals, instead only one, to be reevaluated by original function to update the metamodel. Noticeably, more information will be added by this change. The second modification is that in the iteration loop of approximate ranking procedure, the batch size n b = max(1, ⌊ l/10 ⌋) which is proportional to l is used.
The main objective of approximate ranking procedure is to find out how many individuals to control in each generation. However, the individuals are selected to be evaluated by the original FF according to their approximate fitness and then added to the training set T , subsequently this operation is repeated until the metamodel selection of the parents stays unchanged in two successive iterations. In Algorithm 1, n b is a batch of individuals suggested to be evaluated in each iteration of approximate ranking procedure. The used batch size is n b =⌊ l/10 ⌋ [27].

Proposed algorithm and numerical test 4.1 Proposed algorithm
In this study, the metamodel-assisted evolutionary strategy uses the integration of Kriging metamodel in the CMA-ES algorithm with approximate ranking procedure as a quality control method [25].
The chosen Kriging uses a Gaussian correlation function and a quadratic basis function. Thus, to build Kriging metamodel, there are n k required points, which are a number of free parameters in quadratic basis function, n k ¼ p ¼ 1 2 ðn þ 1Þðn þ 2Þ. Consequently, the approximate ranking procedure will be performed in case the number of required points is assured, otherwise, the number of points jT j in the training set is more or equal to n k .
The training data set T contain all previous evaluated points; however, it is not necessary to use all previous evaluated point as training set, because the computational time to construct a metamodel increases consistently with the number of training point jT j [28]. In the Kriging-based CMA-ES algorithm, the main objective of constructed metamodel is to predict fitness values of the population of the current generation, so it is reasonable to choose the current population as training set. For CMA-ES algorithm, where multivariate normal distribution is used, the Mahalanobis distance D 2 x i ¼ x i À m ðgÞ À Á T ðs ðgÞ À Á 2 C ðgÞ Þ À1 x i À m ðgÞ À Á is appropriate to select training data set [29]. Mahalanobis expression gives the distance between each evaluated points and current distribution mean m (g) , finally the training set can be determined as follows: Algorithm 1. Approximate Ranking Procedure [26].
given: ðz k ; x k Þ l k¼1 ; m ðgÞ ; s ðgÞ ; C ðgÞ ; t; S; fðxÞ: approximate: buildf ðxÞ based on training set T and predictf ðx k Þ; k ¼ 1; :::; l rank and determine the parent set P 1 ≡ x where S is an archive containing all previous evaluated points, and x 2 n p is the quantile function for probability p of the chi-squared distribution with n degrees of freedom. It could happen that jT j n k even if the number of training data jT j is more or equal to n k . In this case, the n k required points are selected from T by changing the quantile function x 2 n p by r n k , the n k smaller Mahalanobis distance of all points in S. Consequently, the training set S can be performed as: where r 2 ¼ maxðx 2 n ðpÞ; n k Þ. Finally, the problem was avoided and the metamodel-assisted CMA-ES can be performed. The Kriging-assisted CMA-ES is provided in Algorithm 2 .

Numerical test
To comprehensively evaluate an algorithm, experimental studies are important to validate and compare the suggested algorithm with other existing algorithms over a good range of test functions. In [26], the proposed KA-CMA-ES using pre-selection (PS), KA-CMA-ES using ARP and fixed generation-based control (FGC) were studied and compared to the CMA-ES algorithm to show the efficiency of the proposed algorithm.
To evaluate proposed algorithm, the set of test functions including many functions with different characteristics was used. The test contains a set of 12 continuous functions, which present a minimization problems defined as: min fðXÞ; X ¼ ½x 1 ; x 2 ; :::; x D T s:t: X ∈ S ¼ ½X LB ; X UB & where f(x) is the FF, S is the search domain which is defined by its lower and upper bounds X LB ∈ R D and X UB ∈ R D , and D is the dimension of problem. For handling the box constraints, the re-sampling method is adopted, i.e. re-sampling any infeasible solution until it becomes feasible.
The 12 test functions are listed in Table 1, in which f 1 -f 7 are unimodal problems and f 8 -f 12 are multimodal functions. All the test functions have zero as global optimum, i.e. f(X * ) = 0, which is located at X * = [0, 0, ..., 0] excluding Rosenbrock function, whose optimum is at X * = [1, 1, ..., 1] T . In the experiments, for f 1 -f 5 , the dimensionality of search space is D = 20; for f 6 -f 11 , D = 10, and f 12 has dimension D = 5. Each test function with each dimension can be considered as an optimization problem.
The Figures 1-3 present the convergence graphs of KA-CMA-ES algorithm and that of the standard CMA-ES [26]. The improvement in performance of KA-CMA-ES is clearly illustrated. The graphs present the median performance of the 25 runs on each function. The convergence graphs also indicate that KA-CMA-ES using ARP-EI converges more quickly among the others. On the one hand, the proposed KACMA-ES algorithm using ARP-EI outperforms the KA-CMA-ES using PS and FGC. On the other hand, the efficiency of the CMA-ES is well improved.

Modelization and finite element analysis
In this study, a 256 pin PQFP microcontroller placed on printed circuit board (PCB) was analyzed [30]. The component is soldered to the PCB through lead-free solder named SAC305. The seals are soldered to the PCB by an SAC305 solder joints. For reliability assessment, the solder joints of the microcontroller are the critical elements in the system. The finite element analysis is used to quantify the lifetime of component. However, the calculation using the global model (see Fig. 4) of the card with the Algorithm 2. Kriging-assisted CMA-ES.
Initialize evolution path p ð0Þ s ¼ 0; p ð0Þ c ¼ 0, the covariance matrix C (0) = I, the step size s (0) and the selection parameters according to [15]. Initialize the mean vector m (0) to a random candidate g ← 0 while Convergence criterion is not reached do

Function
Expression Research domain Sphere  microcontroller [30] has a high computational cost. In that sense, a 3D local model has been developed with the most critical solder joint, in order to decrease the computation time (Fig. 5) [4].
The local model based on the finite element technique of the critical zone is developed, it is the solder joint most stressed by thermal loading. The local submodel (Fig. 5) is composed of a FR4 PCB, an EPOXY resin component, a  copper pin and an SAC305 solder. The boundary conditions of symmetry are introduced in the modeling and the model use the Anand model to take into account the viscoplastic behavior of the solder [31]. All finite element models of this study are developed with the finite element simulation tool Ansys mechanical [32]. Accelerated thermal cycling, recommended by JEDEC standards [6], is applied as thermal load in finite element analysis. The temperature profile applied varies between À40°C and 125°C and 30 min for low dwell and high dwell time. Devices for automotive applications are generally tested in this temperature range. Figure 6 shows the thermal cycle loading. In the FE calculation, the first step is to simulate the reflow soldering process to take into account the initial constraints. This process consists in applying a temperature profile ranging from the melting temperature of solder joint to the ambient temperature of 25°C in 150°C seconds. The second step consists of applying three thermal cycles, between À40°C and 125°C, as the loading condition for the finite element simulation.
In thermomechanical analysis, using the local model, the solder joint material is assumed to have a viscoplastic behavior. Several authors have studied the response of lead-free solder joints (SnAgCu) and proposed equations to model this response. One of the models developed is the Anand model [33] which distinguishes elastic deformations from inelastic deformations, and combines the creep and the instantaneous plastic deformations [34]. However, the Anand model expresses the material viscoplastic behavior and its equation is defined as follow: where _ e p is the inelastic strain rate, Q is the activation energy, A is the pre-exponential factor (1/s), T is the temperature, R is the universal gas constant, j stands for the materials constant, and m is the strain rate sensitivity of the stress. The internal variable s represents the resistance to plastic deformation, and its evolution equation defined as: where where h 0 represents the hardening/softening constant, a is the strain rate sensitivity of hardening or softening, n stands for the strain rate sensitivity for the saturation value of deformation resistance, * describes the saturation value of s associated with a set of given temperature and strain rate, andŝ is a coefficient. Finally, for the Anand constitutive model, there are nine parameters, A, Q/R,ŝ, h 0 , j, m, n, a and s 0 , which are needed to identify the evolution of deformation resistance in equation (15), and the strain rate in equation (14). Table 2 presents the value of these materials constants for solder joint (SnAg). Table 3 shows the material properties of local model used in this study.

Solder joint fatigue life prediction model
The prediction of solder joint fatigue requires a combination of finite element methods with a thermal fatigue model, which is usually obtained from experimental data and accelerated tests. This model is used to calculate the number cycles the system can resist before failure. There are many models for predicting solder joint fatigue life that can be classified into four categories (approach): stress-based approach, plastic deformation and creep approach, energy-based approach, and damage-based approach [35]. The Coffin-Manson fatigue model, classified in the plastic deformation approach, is the best and most commonly used approach. The number of cycles N f is expressed as a function of the inelastic strain range, De in , the fatigue ductility exponent a(0.73), and the fatigue ductility coefficient u(u = 3.7) [36,37]. The form of Anand model used in the present analysis is given by the following equation: The inelastic strain range, De in , is calculated using a nonlinear finite element simulation based on the Anand model.
After the termination of the thermal cycles execution in ANSYS TM . Figure 7 indicates that the strongly deformed areas of component occurs in solder joint and the Figure 8 shows the evolution of inelastic deformation in solder joint resulting from thermal cycles execution.    (Table 4) and N f is the number of fatigue cycles. The number of fatigue cycles N f is calculated through the plastic deformation using the empirical Coffin-Manson model. In other words, maximizing N f amounts to minimize the plastic deformation De in of the solder joint due to an active thermal load cycle. Consequently, in this problem, the objective function represents the plastic deformation De in of the solder joint due to an active thermal load cycle (Tab. 4).
min De in ðfdgÞ The design optimization problem is carried out with the KA-CMA-ES algorithm, implemented with Matlab ® (Fig. 10). Table 5 shows the optimal and initial design calculated by a simple CMA-ES algorithm and KA-CMA-ES, also the results show that the objective function (inelastic strain range De in ) is minimized. Figure 11 shows the new design and the new inelastic strain distribution after the application of the design optimization procedure. The maximum inelastic strain (SMX = 0.08193) is more smaller than the inelastic strain in the initial configuration (SMX = 0.13557). The typical convergence histories of plastic deformation De in are plotted in Figure 12. This plots clearly indicate that the KA-CMA-ES converges more quickly than standard CMA-ES and give same results.

Conclusion
The robust design of complex MSs often induces an expensive optimization problem. In this work, the performance of CMA-ES has been used with the Kriging metamodel to efficiently solve expensive mechatronic optimization problem. the virtual thermo-mechanical test was performed to evaluate the reliability of solder joints of a mechatronic device. The 3D FE model takes into account the nonlinearities properties of viscoplastic behavior of the solder joints. This study illustrates the interest to use the metamodeling techniques with CMA-ES, to increase the reliability of solder joints of the mechatronic devices in order to overcome the computational cost problem. For training Kriging metamodel, the most relevant data are used as training set. The approximate ranking procedure, which is a reliable adaptive evolution control, was adopted to control the fidelity of the metamodel during the KA-CMA-ES. The application of this method in solder joints global optimization shows its efficiency in reducing the number of finite element simulations to reach an optimal design against the standard CMA-ES method. As perspective we expect to modify this algorithm to constrained optimization in order to take into account the uncertainties.