© H. Hamdani et al., published by EDP Sciences, 2019

1 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–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 thermo-mechanical 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.

2 Optimization technics

2.1 Covariance matrix adaptation-evolution strategy (CMA-ES)

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, λ offspring candidates are generated from μ parent. In the next generation, the new parent (μ) is selected from the offspring candidates (λ) using a (λ, μ)-selection. The μ best candidates are selected according to their ranking based on the objective function, CMA-ES (λ, μ) 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  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:  and , where C ^g represents a covariance matrix (a definite positive matrix) and the σ ^g is a scalar factor to control the step length of the iteration. Finally, the λ search points (individuals, offspring) are sampled through the following equation:(1)with λ ≥ 2 is the population size (number of offspring), and  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), σ ^(g) and C ^(g) need to be adapted at each generation.

After sampling process, the λ offspring is evaluated and ranked based on the objective function. Finally, the μ best individual is selected, and the mean m ^(g) is then updated by taking the weighted mean of the μ individuals(2) where are the μ strictly positive and normalized weights, satisfying .

The second term to update is the step size σ that uses the so called evolution path  [15] for its accumulative adaptation. The σ computation is given as follow:(3)where  is the expectation of the Euclidean norm of an N (0, I) distributed random vector, which is approximated by:(4)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:(5)where  is the covariance matrix evolution path [15].

The remaining parameters for CMA-ES update are given in [15].

2.2 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 ₁, … , 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:(6) where h(x) is the mean of the process, and Z(x) is a Gaussian process with zero mean and covariance expressed by:(7) with σ ² is the variance of Gaussian process and R(x ⁽ⁱ⁾, 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:(8)where d is the dimension of design space, θ _k (k  = 1, 2, ..., d) are unknown parameters of the correlation function, and respectively are the kth component of the sampling point x _i and x _j .

For any new point x, the mean and variance of prediction [18] can be respectively calculated by:(9) (10)where and s ²(x) are respectively the estimated mean value and variance of , h ^T  = [h _i ], i  = 1, … , k a set of basis functions (e.g. polynomial functions), β = (β ₁, ..., β _k ) the associated regression coefficients, R is the correlation matrix of D, and H the matrix corresponding to the values of h ^T (D), and r(x) = [R(x, x ₁), … , R(x, x _k )] ^T is the vector of correlation functions between the untried point x and the k samples x ⁽¹⁾, ..., x ^(k).

The unknown model parameters can be determined using likelihood estimates (MLEs) method [20].

3 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.

3.1 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 so-called 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.

3.2 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 λ 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.3λ ⌋) 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, ⌊ λ/10 ⌋) which is proportional to λ 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 , 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  =⌊ λ/10 ⌋ [27].

4 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, . Consequently, the approximate ranking procedure will be performed in case the number of required points is assured, otherwise, the number of points  in the training set is more or equal to n _k .

The training data set  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   [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 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:(11)where is an archive containing all previous evaluated points, and  is the quantile function for probability p of the chi-squared distribution with n degrees of freedom. It could happen that  even if the number of training data is more or equal to n _k . In this case, the n _k required points are selected from  by changing the quantile function by , the n _k smaller Mahalanobis distance of all points in . Consequently, the training set  can be performed as:(12)where . Finally, the problem was avoided and the metamodel-assisted CMA-ES can be performed. The Kriging-assisted CMA-ES is provided in Algorithm 2 .

4.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: 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 ₁–f ₇ are unimodal problems and f ₈–f ₁₂ 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 ₁–f ₅, the dimensionality of search space is D = 20; for f ₆–f ₁₁, D = 10, and f ₁₂ 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.

Fig. 1 Convergence graphs of f 1–f 4.

Fig. 2 Convergence graphs of f 5–f 8.

Fig. 3 Convergence graphs of f 9–f 11.

5 Numerical simulation

5.1 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 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:(13)where 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, ξ 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:(14)where(15)where h ₀ 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 ₀, ξ, m, n, a and s ₀, 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.

Fig. 4 Global model of the solder joints.

Fig. 5 Local model of the solder joints.

Fig. 6 Thermal cycle loading.

5.2 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, Δε_in, the fatigue ductility exponent α(0.73), and the fatigue ductility coefficient θ(θ  =  3.7) [36,37]. The form of Anand model used in the present analysis is given by the following equation:(16)

The inelastic strain range, , 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.

Fig. 7 Side view of plastic strain distribution in the solder joint at last loading thermal cycle.

Fig. 8 Plastic strain history.

5.3 Formulation of the global optimization algorithm for the solder joint

The design optimization problem for the solder joints aims to maximize the number of fatigue cycles of the solder joint:(17){d} = { h _1, h ₂, h ₃} are the design variables (Fig. 9) which represent three design parameters which have a major impact on fatigue life of solder joint [4], {d} ^Lb and {d} ^Ub are respectively their lower and upper bound (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  of the solder joint due to an active thermal load cycle. Consequently, in this problem, the objective function represents the plastic deformation  of the solder joint due to an active thermal load cycle (Tab. 4).(18)

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 Δε _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 Δε _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.

Fig. 9 Description of optimization parameters.

Fig. 10 Flowchart of global optimization methodology.

Figure 11 Strain inelastic distribution before (a) and after (b) optimization procedure

Fig. 12 Average plot of Δεin versus number of numerical simulation.

6 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.

Acknowledgments

The authors would like to thank PHC Toubkal/17/43: integrated action Morocco France for their financial support for the realization of this work.

References

All Tables

All Figures

	Fig. 1 Convergence graphs of f 1–f 4.
In the text

	Fig. 2 Convergence graphs of f 5–f 8.
In the text

	Fig. 3 Convergence graphs of f 9–f 11.
In the text

	Fig. 4 Global model of the solder joints.
In the text

	Fig. 5 Local model of the solder joints.
In the text

	Fig. 6 Thermal cycle loading.
In the text

	Fig. 7 Side view of plastic strain distribution in the solder joint at last loading thermal cycle.
In the text

	Fig. 8 Plastic strain history.
In the text

	Fig. 9 Description of optimization parameters.
In the text

	Fig. 10 Flowchart of global optimization methodology.
In the text

	Figure 11 Strain inelastic distribution before (a) and after (b) optimization procedure
In the text

	Fig. 12 Average plot of Δεin versus number of numerical simulation.
In the text