GA-BP in Thermal Fatigue Failure Prediction of Microelectronic Chips

: A thermal fatigue life prediction model of microelectronic chips based on thermal fatigue tests and solder / substrate interfacial singularity analysis from ﬁnite element method (FEM) analysis is established in this paper. To save the calculation of interfacial singular parameters of new chips for life prediction, and improve the accuracy of prediction results in actual applications, a hybrid genetic algorithm–artiﬁcial neural network (GA–ANN) strategy is utilized. The proposed algorithm combines the local searching ability of the gradient-based back propagation (BP) strategy with the global searching ability of a genetic algorithm. A series of combinations of the dimensions and thermal mechanical properties of the solder and the corresponding singularity parameters at the failure interface are used to train the proposed GA-BP network. The results of the network, together with the established life prediction model, are used to predict the thermal fatigue lives of new chips. The comparison between the network results and thermal fatigue lives recorded in experiments shows that the GA-BP strategy is a successful prediction technique.


Introduction
Microelectronic chips are becoming more popular, because their design in smaller dimensions meets the requirements of semiconductor packaging in terms of high-density and cost-effective performance. The solder layer plays an important role in the mechanical and electric connections of the chip; therefore, its performance, especially the thermomechanical properties of the solder, has become an important factor affecting the reliability of microelectronic chips [1][2][3]. During the service process, the chips are subjected to temperature cycling caused by circuit continuity and ambient temperature changes [4]. Due to the difference in the coefficients of thermal expansion (CTE) of various components, stress and strain always occur and gradually accumulate at the solder joints during temperature cycling. A severe CTE mismatch always results in strong interfacial shear and peeling stresses near the interfacial edges and corners. When the local stress level exceeds the bonding strength of the interface, an interfacial crack will initiate and propagate [5][6][7]. Therefore, it is of practical importance to model these stresses using analytical tools, so that the susceptibility to thermal mechanical failure can be predicted for chips with new geometries and material combinations, without costly trial and error.
No well recognized methodologies on how to quantify the interfacial strength and predict crack initiation and propagation have been established. The traditional method is to apply the fatigue law of the solder material itself directly to the solder joints, but this method is not adequate, since the stresses and strains are singular at the idealized interface edge or corner [8][9][10]. To overcome this concern, Hattori et al. suggested an approach for assessing the interface reliability of plastic IC packages using two stress intensity parameters that characterize the stress distributions near a bonded edge along the interface [11]. Similar approaches have been adopted by Gradin [12], Groth [13], and Reedy and Guess [14]. Generally, at the edge, stresses have the singular form σ = K/r 1−λ , from which the stress intensity factor K and singular order λ are used to configure a criterion for fatigue crack initiation or delamination. For various loadings, if this value is larger than a certain threshold value K c , crack initiation or delamination is assumed to occur. However, such a method is based on the concept of point failure, which may be difficult to observe through physical experiments, and is valid only for cases where the geometric shape of the solder joints and the loading conditions are the same, as this is used to inform the fatigue law [15,16]. That is, the empirical fatigue law may have strong limitations in application. For example, the stress intensity factor at the edge of the interface is always dependent on test conditions, assembly geometry, and the mechanical property of the materials and their interactions, as well as the impact of process parameters [17,18], and the process of obtaining accurate singular field parameters at the interface is also time consuming. To save the calculation of the stress intensity factor at the edge of the interface, soft-computing is a good alternative for handling this complex problem, as it is tolerant of imprecision and uncertainty. To date, various soft-computing methods such as ANN [19][20][21], fuzzy logic [22], genetic programming [23][24][25], and adaptive neuro-fuzzy inference systems [26] have been applied to the fatigue problems of engineering materials and structures.
In this paper, a thermal fatigue law is established to predict the service life of microelectronic chips, based on the concept that the fatigue failure initiates within a characteristic length, rather than one point from the solder/substrate interface. The parameters in the fatigue law are determined according to the thermal fatigue lives of chips from tests and interfacial singularity parameters, such as interfacial stress-strain intensity factors and singular orders, obtained from mechanical thermal FEM analysis under the same testing conditions. An attempt has been made to predict the interfacial singularity parameters and thermal fatigue life by applying GA-BP hybrid programming, in which the GA initiates the weights and improves the convergence of the BP network, and the BP network predicts the interfacial singularity parameters of new chips for thermal life prediction. It is observed that the predicted results are in good agreement with the experimental findings.

Thermal Fatigue Tests
The chipset shown in Figure 1 has nineteen pieces of chips, including two pieces of chip Q1, one piece of chip Q2, six pieces of chip Q3, four pieces of chip Q4, and six pieces of chip Q5. The structure of every chip is similar and the specific dimensions of a single chip are shown in Figure 2. The materials shown are the wafer (silicon), solder Pb-5Sn, Cu, solder SnAg3Cu0.5, Cu, insulate layer, substrate, and silica gel. The silica gel protector is cleaned and then the chipset is cut out by a wire cutter. Five types of chips are then put into a same temperature box to undergo cyclic thermal charge, as shown in Figure 3, with a temperature amplitude of 20-150 • C, 20-120 • C and 20-100 • C, respectively.
As the damage or crack in the chip cannot be observed directly, the fatigue life limit has to be defined first. To detect the fatigue failure, the electric resistance of the chips is measured by a digital resistance meter for a certain interval. Figure 4 depicts the typical electric resistance variation of chip Q2 with the thermal cycles (20-150 • C). Figure 5 gives the SEM observation of crack nucleation at the interface of the wafer and solder Pb-5Sn (interface S12), and solder SnAg3Cu0.5 and Cu (interface S45). It is found that the main crack forms first from S45. By comparing the forming of the main crack at S45 and the electronic resistance increase (RI), we define the fatigue life limit N f by RI reaching 15% [27]. Using this method, we determine the average thermal fatigue lives of five types of chips under different experimental conditions, as listed in Table 1.   As the damage or crack in the chip cannot be observed directly, the fatigue life limit has to be defined first. To detect the fatigue failure, the electric resistance of the chips is measured by a digital resistance meter for a certain interval. Figure 4 depicts the typical electric resistance variation of chip Q2 with the thermal cycles (20-150°C). Figure 5 gives the SEM observation of crack nucleation at the interface of the wafer and solder Pb-5Sn (interface S12), and solder SnAg3Cu0.5 and Cu (interface S45). It is found that the main crack forms first from S45. By comparing the forming of the main crack   As the damage or crack in the chip cannot be observed directly, the fatigue life limit has to be defined first. To detect the fatigue failure, the electric resistance of the chips is measured by a digital resistance meter for a certain interval. Figure 4 depicts the typical electric resistance variation of chip Q2 with the thermal cycles (20-150°C). Figure 5 gives the SEM observation of crack nucleation at the interface of the wafer and solder Pb-5Sn (interface S12), and solder SnAg3Cu0.5 and Cu (interface S45). It is found that the main crack forms first from S45. By comparing the forming of the main crack   As the damage or crack in the chip cannot be observed directly, the fatigue life limit has to be defined first. To detect the fatigue failure, the electric resistance of the chips is measured by a digital resistance meter for a certain interval. Figure 4 depicts the typical electric resistance variation of chip Q2 with the thermal cycles (20-150°C). Figure 5 gives the SEM observation of crack nucleation at the interface of the wafer and solder Pb-5Sn (interface S12), and solder SnAg3Cu0.5 and Cu (interface S45). It is found that the main crack forms first from S45. By comparing the forming of the main crack   at S45 and the electronic resistance increase (RI), we define the fatigue life limit Nf by RI reaching 15% [27]. Using this method, we determine the average thermal fatigue lives of five types of chips under different experimental conditions, as listed in Table 1.

FEM Model
A three dimensional numerical analysis was carried out to obtain a measure of stress and strain under cyclic thermal loading. In this paper, a new segmented constitutive model of the solder is considered. The model consists of two parts: the relationship of strain rate with stress is linear at low stress and has a hyperbolic sine at high stress, as shown in Equation (1). All constants and other material constants in the FEM analysis can be found in the previous work [27]. at S45 and the electronic resistance increase (RI), we define the fatigue life limit Nf by RI reaching 15% [27]. Using this method, we determine the average thermal fatigue lives of five types of chips under different experimental conditions, as listed in Table 1.

FEM Model
A three dimensional numerical analysis was carried out to obtain a measure of stress and strain under cyclic thermal loading. In this paper, a new segmented constitutive model of the solder is considered. The model consists of two parts: the relationship of strain rate with stress is linear at low stress and has a hyperbolic sine at high stress, as shown in Equation (1). All constants and other material constants in the FEM analysis can be found in the previous work [27].

FEM Model
A three dimensional numerical analysis was carried out to obtain a measure of stress and strain under cyclic thermal loading. In this paper, a new segmented constitutive model of the solder is considered. The model consists of two parts: the relationship of strain rate with stress is linear at low stress and has a hyperbolic sine at high stress, as shown in Equation (1). All constants and other material constants in the FEM analysis can be found in the previous work [27]. where σ is the equivalent stress, σ v is the linear viscous creep limit (which is the cross stress between the linear viscous and power-law creep), H is the activation energy, R is the universal gas constant, T is the absolute temperature value, and T R is 273 degrees.
An unsteady thermal conduction analysis and thermal stress analysis were carried out with the corresponding chips. The FEM model of chip Q2 is shown in Figure 6, where the total number of elements is 78,890, the number of nodes is 85,200, and the smallest element size at the interfacial singular zone is 0.00001 mm × 0.00003 mm × 0.00003 mm. The convergence analysis was conducted to ensure the accuracy of our results, as well as the computational efficiency. In the thermal conduction analysis, the surrounding temperature variation is exactly the same as that in the thermal fatigue tests. To express the stress and strain singularity at the edge of the interface between the solder and substrate, specific points at the edges S12 and S45 are selected, as shown in Figure 7. The typical stress and equivalent strain variations with time at point A and B are shown in Figures 8 and 9. The stress near the interfacial edge varies for the first several cycles, but after three cycles it is saturated. The most severe case does not occur at the initial or stable stage, but at the stage of heating. In this case, the normal stress on interface S12 is compressive, and in conclusion, the shear stress is the dominant factor of interface failure. At S45 the maximum traction σ t = (σ 2 + τ 2 ) 1/2 dominates interface failure because the normal stress is tensile.
where σ is the equivalent stress, σv is the linear viscous creep limit (which is the cross stress between the linear viscous and power-law creep), H is the activation energy, R is the universal gas constant, T is the absolute temperature value, and TR is 273 degrees. An unsteady thermal conduction analysis and thermal stress analysis were carried out with the corresponding chips. The FEM model of chip Q2 is shown in Figure 6, where the total number of elements is 78,890, the number of nodes is 85,200, and the smallest element size at the interfacial singular zone is 0.00001 mm × 0.00003 mm × 0.00003 mm. The convergence analysis was conducted to ensure the accuracy of our results, as well as the computational efficiency. In the thermal conduction analysis, the surrounding temperature variation is exactly the same as that in the thermal fatigue tests. To express the stress and strain singularity at the edge of the interface between the solder and substrate, specific points at the edges S12 and S45 are selected, as shown in Figure 7. The typical stress and equivalent strain variations with time at point A and B are shown in Figures 8 and 9. The stress near the interfacial edge varies for the first several cycles, but after three cycles it is saturated. The most severe case does not occur at the initial or stable stage, but at the stage of heating. In this case, the normal stress on interface S12 is compressive, and in conclusion, the shear stress is the dominant factor of interface failure. At S45 the maximum traction σt = (σ 2 +τ 2 ) 1/2 dominates interface failure because the normal stress is tensile.
where σ is the equivalent stress, σv is the linear viscous creep limit (which is the cross stress between the linear viscous and power-law creep), H is the activation energy, R is the universal gas constant, T is the absolute temperature value, and TR is 273 degrees. An unsteady thermal conduction analysis and thermal stress analysis were carried out with the corresponding chips. The FEM model of chip Q2 is shown in Figure 6, where the total number of elements is 78,890, the number of nodes is 85,200, and the smallest element size at the interfacial singular zone is 0.00001 mm × 0.00003 mm × 0.00003 mm. The convergence analysis was conducted to ensure the accuracy of our results, as well as the computational efficiency. In the thermal conduction analysis, the surrounding temperature variation is exactly the same as that in the thermal fatigue tests. To express the stress and strain singularity at the edge of the interface between the solder and substrate, specific points at the edges S12 and S45 are selected, as shown in Figure 7. The typical stress and equivalent strain variations with time at point A and B are shown in Figures 8 and 9. The stress near the interfacial edge varies for the first several cycles, but after three cycles it is saturated. The most severe case does not occur at the initial or stable stage, but at the stage of heating. In this case, the normal stress on interface S12 is compressive, and in conclusion, the shear stress is the dominant factor of interface failure. At S45 the maximum traction σt = (σ 2 +τ 2 ) 1/2 dominates interface failure because the normal stress is tensile.     As the initial condition is set as a zero stress-strain state, the stress and strain ranges are important for fatigue life evaluation. According to the instant singular field theory [28,29]: where Ki(t) and Kε(t) denote the stress and strain intensity factors, δi(t) and ς(t) are the stress and strain singular orders, and r denotes the distance from the singular edge. The intensity factors and singular orders of maximum stress and strain at steady state are determined numerically, as shown in Figures 10 and 11. The difference of stress singular orders between the two interfaces is not big, but the stress and strain intensity factors at S45 are larger than that at S12. The singularity of strain is stronger than that of stress, indicating the strain field may be the dominant factor of thermal fatigue failure. As shown with the above methods, the strain singular orders and intensity factors of five types of chips at the maximum state are determined and listed in Table 2.  As the initial condition is set as a zero stress-strain state, the stress and strain ranges are important for fatigue life evaluation. According to the instant singular field theory [28,29]: where Ki(t) and Kε(t) denote the stress and strain intensity factors, δi(t) and ς(t) are the stress and strain singular orders, and r denotes the distance from the singular edge. The intensity factors and singular orders of maximum stress and strain at steady state are determined numerically, as shown in Figures 10 and 11. The difference of stress singular orders between the two interfaces is not big, but the stress and strain intensity factors at S45 are larger than that at S12. The singularity of strain is stronger than that of stress, indicating the strain field may be the dominant factor of thermal fatigue failure. As shown with the above methods, the strain singular orders and intensity factors of five types of chips at the maximum state are determined and listed in Table 2. As the initial condition is set as a zero stress-strain state, the stress and strain ranges are important for fatigue life evaluation. According to the instant singular field theory [28,29]: where K i (t) and K ε (t) denote the stress and strain intensity factors, δ i (t) and ς(t) are the stress and strain singular orders, and r denotes the distance from the singular edge. The intensity factors and singular orders of maximum stress and strain at steady state are determined numerically, as shown in Figures 10 and 11. The difference of stress singular orders between the two interfaces is not big, but the stress and strain intensity factors at S45 are larger than that at S12. The singularity of strain is stronger than that of stress, indicating the strain field may be the dominant factor of thermal fatigue failure. As shown with the above methods, the strain singular orders and intensity factors of five types of chips at the maximum state are determined and listed in Table 2.

Thermal Fatigue Life Prediction Model
According to the FEM analysis and thermal fatigue tests, the failure from S45 is identified as strain-controlled creep failure [28,29]. Therefore, the basic strain-controlled fatigue law can be written as the Coffin-Manson formula:

Thermal Fatigue Life Prediction Model
According to the FEM analysis and thermal fatigue tests, the failure from S45 is identified as strain-controlled creep failure [28,29]. Therefore, the basic strain-controlled fatigue law can be written as the Coffin-Manson formula:

Thermal Fatigue Life Prediction Model
According to the FEM analysis and thermal fatigue tests, the failure from S45 is identified as strain-controlled creep failure [28,29]. Therefore, the basic strain-controlled fatigue law can be written as the Coffin-Manson formula: where ∆ε is the range of strain, and n and C are the constants to be determined by experimental results. Considering the strain range average within a characteristic failure length l ε along the interfacial edge and neglecting the dependence of ς(t) on the maximum and minimum strain states, we get: where ∆ε donates the average strain range within l ε .
Substituting ∆ε into Equation (3), we obtain: Here l ε is always dependent on ∆K ε and ς. If the separation in variables of l ε is taken into account, Equation (5) can be simplified as: where f (ς) = f 1/m 1 (ς) is a function of the singular order and m is a constant. Based on the experimental conditions outlined in this paper, the thermal fatigue life was determined, and the range of strain intensity factor was solved by numerical simulation. Figure 12 shows the relationship between K ε /(1−ς) and N f when considering the failure from S45. The fitting lines for different chips are almost parallel to each other. This fact means that m is really a material constant, and independent to the singular orders and intensity factors. The slope equals 1/m, so m = 1/0.292. Figure 13 shows the relationship of f (ς) with singular order ς in polynomial form. Therefore, the thermal fatigue life law for S45 is finally determined as: Electronics 2019, 7, x FOR PEER REVIEW 8 of 14 where Δε is the range of strain, and n and C are the constants to be determined by experimental results.
Considering the strain range average within a characteristic failure length lε along the interfacial edge and neglecting the dependence of ς(t) on the maximum and minimum strain states, we get: where ε Δ donates the average strain range within lε. Substituting ε Δ into Equation (3), we obtain: Here lε is always dependent on ΔKε and ς. If the separation in variables of lε is taken into account, Equation (5) can be simplified as: is a function of the singular order and m is a constant.
Based on the experimental conditions outlined in this paper, the thermal fatigue life was determined, and the range of strain intensity factor was solved by numerical simulation. Figure 12 shows the relationship between Kε/(1-ς) and Nf when considering the failure from S45. The fitting lines for different chips are almost parallel to each other. This fact means that m is really a material constant, and independent to the singular orders and intensity factors. The slope equals 1/m, so m=1/0.292. Figure 13 shows the relationship of f(ς) with singular order ς in polynomial form. Therefore, the thermal fatigue life law for S45 is finally determined as:  Theoretically, once the strain intensity factors and singular orders of new chips at the failure interface are obtained, their thermal fatigue lives can be predicted. At the same time, the calculation of thermal mechanical analysis to get the singularity parameters is not easy, and tends to be timeconsuming. Therefore, the GA-BP based method discussed below needs to be established for the necessary parameters in thermal fatigue life prediction.

GA-BP: Theoretical Background
The genetic algorithm (GA) is a domain-independent problem-solving technique, which encodes the problem into a chromosome composed of genes and generates a chromosome representing the solution to the problem through a continuous process of selecting, crossing, and mutating the chromosome. The GA operates on the chromosome code, and it isolates the characteristics of the problem itself, so it has wide adaptability. At the same time, since the algorithm operates in parallel from multiple initial points, the search process is free from convergence to a local extremum. The GA implements the search process through groups, making it different from a singlepoint search and easy to parallelize, thus improving the efficiency of the algorithm [30,31].
The classic back propagation (BP) is a three, or more than three, layer hierarchical forward neural network. The BP algorithm consists of two parts: the forward transmission of information and the back propagation of errors. In the forward transmission, the input information from the input layer passes to the output layer through the hidden layer. If the desired output is not obtained at the output layer, the error change of the output layer is calculated and passed back along the original connection path, to modify the weight of each layer of neurons until the desired goal is reached [32]. Mathematically, the BP model works based on an error gradient descent criterion, and its calculation process is a process of obtaining a smaller error. However, it is often unable to achieve global convergence and falls into local minima when solving complex problems, resulting in an invalid learning process. Additionally, the learning algorithm converges slowly, especially near the target. Finally, the input and output nodes of the network are determined by the specific problem to be studied, but the selection of hidden nodes is often based on experience and is short of theoretical guidance.
The GA is an effective way to avoid the above problems. To maintain a certain size of population, the GA handles several points in the search space, rather than a single point as in the gradient descent method, which helps the search for the global best and is free from the local minimum [33]. This can avoid the defect of the gradient descent method when used in a BP artificial neural network, and ensure a good convergence speed [34]. Theoretically, once the strain intensity factors and singular orders of new chips at the failure interface are obtained, their thermal fatigue lives can be predicted. At the same time, the calculation of thermal mechanical analysis to get the singularity parameters is not easy, and tends to be time-consuming. Therefore, the GA-BP based method discussed below needs to be established for the necessary parameters in thermal fatigue life prediction.

GA-BP: Theoretical Background
The genetic algorithm (GA) is a domain-independent problem-solving technique, which encodes the problem into a chromosome composed of genes and generates a chromosome representing the solution to the problem through a continuous process of selecting, crossing, and mutating the chromosome. The GA operates on the chromosome code, and it isolates the characteristics of the problem itself, so it has wide adaptability. At the same time, since the algorithm operates in parallel from multiple initial points, the search process is free from convergence to a local extremum. The GA implements the search process through groups, making it different from a single-point search and easy to parallelize, thus improving the efficiency of the algorithm [30,31].
The classic back propagation (BP) is a three, or more than three, layer hierarchical forward neural network. The BP algorithm consists of two parts: the forward transmission of information and the back propagation of errors. In the forward transmission, the input information from the input layer passes to the output layer through the hidden layer. If the desired output is not obtained at the output layer, the error change of the output layer is calculated and passed back along the original connection path, to modify the weight of each layer of neurons until the desired goal is reached [32]. Mathematically, the BP model works based on an error gradient descent criterion, and its calculation process is a process of obtaining a smaller error. However, it is often unable to achieve global convergence and falls into local minima when solving complex problems, resulting in an invalid learning process. Additionally, the learning algorithm converges slowly, especially near the target. Finally, the input and output nodes of the network are determined by the specific problem to be studied, but the selection of hidden nodes is often based on experience and is short of theoretical guidance.
The GA is an effective way to avoid the above problems. To maintain a certain size of population, the GA handles several points in the search space, rather than a single point as in the gradient descent method, which helps the search for the global best and is free from the local minimum [33]. This can avoid the defect of the gradient descent method when used in a BP artificial neural network, and ensure a good convergence speed [34].

GA-BP Hybrid Programming
The main idea of GA-BP hybrid programming is to optimize the weights and thresholds of the BP network based on the GA. Real coding is used to evolve the weights and thresholds of the network. Firstly, the fitness function is determined, then the evolution process (the global search operator of the algorithm, such as the selection, crossover and mutation operator) is selected. Lastly, the GA searches for the optimal weights and thresholds in the whole solution space. The algorithm flow chart of a real coding genetic algorithm combined with a neural network is shown in Figure 14. In essence, the main role of the GA is to optimize the initial weight distribution and locate some better search spaces in the solution space, which facilitates the BP algorithm's search for the optimal solution in these small solution spaces. Generally, the efficiency of the GA-BP hybrid training is superior to BP training alone.

GA-BP Hybrid Programming
The main idea of GA-BP hybrid programming is to optimize the weights and thresholds of the BP network based on the GA. Real coding is used to evolve the weights and thresholds of the network. Firstly, the fitness function is determined, then the evolution process (the global search operator of the algorithm, such as the selection, crossover and mutation operator) is selected. Lastly, the GA searches for the optimal weights and thresholds in the whole solution space. The algorithm flow chart of a real coding genetic algorithm combined with a neural network is shown in Figure 14.
In essence, the main role of the GA is to optimize the initial weight distribution and locate some better search spaces in the solution space, which facilitates the BP algorithm's search for the optimal solution in these small solution spaces. Generally, the efficiency of the GA-BP hybrid training is superior to BP training alone.

BP section
Determine the topology of BP network  Figure 14. Flow chart of a GA-BP neural network.

Variables in the Model
The main aim of the present investigation is to employ GA and BP programming to predict the singular orders and strain intensity factors under thermal fatigue cycles, and then calculate the thermal fatigue lives of microelectronic chips under different conditions according to the established thermal fatigue law. The singular orders and strain intensity factors at the failure interface are always dependent on the loading condition, geometry and thermal mechanical properties of the solder, and so forth, while the connections between them cannot be explicitly expressed. Therefore, it is userfriendly to clear these mathematical difficulties by a black box. The singular orders and strain intensity factors can be expressed in the following form:

Variables in the Model
The main aim of the present investigation is to employ GA and BP programming to predict the singular orders and strain intensity factors under thermal fatigue cycles, and then calculate the thermal fatigue lives of microelectronic chips under different conditions according to the established thermal fatigue law. The singular orders and strain intensity factors at the failure interface are always dependent on the loading condition, geometry and thermal mechanical properties of the solder, and so forth, while the connections between them cannot be explicitly expressed. Therefore, it is user-friendly to clear these mathematical difficulties by a black box. The singular orders and strain intensity factors can be expressed in the following form: K ε , ς = f (load condition, geometry, material properties, · · ·) It is possible to estimate the singular orders and strain intensity factors by a trained neural network instead of FEM analysis. However, to train an efficient neural network, it is necessary to have a huge amount of FEM analysis results as the samples for training. Once the neural network is well trained, the singular orders and strain intensity factors can be obtained easily. To reach the stated goal, four variables-elasticity of modulus E, thickness h and width d of solder SnAg3Cu0.5, and cyclic peak temperature T max -are considered as the parameters of the input layer of the GA-BP network, as shown in Figure 15. According to the different values of the input variables, 60 sets of test samples are prefabricated from the FEM simulation, from which 50 sets are selected as training data, and the remaining 10 sets are used as testing data. To achieve the efficient learning of the GA-BP programming, various parameters such as population size, number of generations, function set, mutation rate, crossover rate, and so forth need to be set in advance [35]. The optimum combination of values of the above parameters is obtained by the performance of several trials, as given in Table 3. It is possible to estimate the singular orders and strain intensity factors by a trained neural network instead of FEM analysis. However, to train an efficient neural network, it is necessary to have a huge amount of FEM analysis results as the samples for training. Once the neural network is well trained, the singular orders and strain intensity factors can be obtained easily. To reach the stated goal, four variables-elasticity of modulus E, thickness h and width d of solder SnAg3Cu0.5, and cyclic peak temperature Tmax-are considered as the parameters of the input layer of the GA-BP network, as shown in Figure 15. According to the different values of the input variables, 60 sets of test samples are prefabricated from the FEM simulation, from which 50 sets are selected as training data, and the remaining 10 sets are used as testing data. To achieve the efficient learning of the GA-BP programming, various parameters such as population size, number of generations, function set, mutation rate, crossover rate, and so forth need to be set in advance [35]. The optimum combination of values of the above parameters is obtained by the performance of several trials, as given in Table  3.

Predicted Results
The predicted singular strain orders and intensity factors of the steady state of S45 are listed in Table 4. The performance of the GA-BP models in the fatigue life prediction is statistically evaluated using two prediction score metrics calculated from the test dataset: the coefficient of correlation (R 2 ) and mean-square error (MSE). The calculated R 2 and MSE are 0.9578 and 1.8527, respectively. The trained neural network can predict the singular parameters with good accuracy. The thermal fatigue lives, based on the predicted singular orders and strain intensity factors, are calculated and compared with the tested ones, as depicted in Figure 16. The result shows that the difference between the predicted output and tested output is small, meeting the error requirement. Therefore, the method of GA-BP in predicting the thermal fatigue lives of microelectronic chips is feasible.

Predicted Results
The predicted singular strain orders and intensity factors of the steady state of S45 are listed in Table 4. The performance of the GA-BP models in the fatigue life prediction is statistically evaluated using two prediction score metrics calculated from the test dataset: the coefficient of correlation (R 2 ) and mean-square error (MSE). The calculated R 2 and MSE are 0.9578 and 1.8527, respectively. The trained neural network can predict the singular parameters with good accuracy. The thermal fatigue lives, based on the predicted singular orders and strain intensity factors, are calculated and compared with the tested ones, as depicted in Figure 16. The result shows that the difference between the predicted output and tested output is small, meeting the error requirement. Therefore, the method of GA-BP in predicting the thermal fatigue lives of microelectronic chips is feasible.

Conclusions
The present study outlines an efficient approach for the prediction of thermal fatigue of microelectronic chips under cyclic thermal loading, using GA-BP hybrid programming. From the result of the present investigation, it is observed that the prediction accuracy of the thermal fatigue life by the developed GA-BP method is in good agreement with the experimental results. The application of thermal fatigue life prediction demonstrates that the thermal fatigue of chips can be evaluated uniformly, no matter what the shapes, materials and loading types are, as long as the relevant stress-strain intensity factors and singular orders can be obtained.
However, as shown by the outcome of numerical analysis and experimental results, the thermal fatigue model involved several factors such as the thermal fatigue failure mode, the diversification of singular parameters in FEM analysis, and the determination and measurement of thermal fatigue life during experiments throughout the modeling process. This causes certain restrictions in its application. Further research based on other computational intelligence approaches, like computational intelligence aided design, should be employed to improve the prediction accuracy of the established thermal fatigue model.