Multi-Objective Reliability-Based Optimization of Control Arm Using MCS and NSGA-II Coupled with Entropy Weighted GRA

Abstract

Lightweight design is one of the important ways to reduce automobile fuel consumption and exhaust emissions. At the same time, the fatigue life of automobile parts also greatly affects vehicle safety. This paper proposes a multi-objective reliability optimization method by integrating Monte Carlo simulation (MCS) with the NSGA-II algorithm coupled with entropy weighted grey relational analysis (GRA) for lightweight design of the lower control arm of automobile Macpherson suspension. The dynamic load histories of the control arm were extracted through dynamic simulations of a rigid-flexible coupling vehicle model on virtual proving ground. Then, the nominal stress method was used to predict its fatigue life. Six design variables were defined to describe the geometric dimension of the control arm, while mass and fatigue life were taken as optimization objectives. The multi-objective optimization design of the control arm was carried out based on the Kriging surrogate model and NSGA-II algorithm. Aiming at the uncertainty of design variables, the reliability constraint was added to the multi-objective optimization to improve the reliability of the fatigue life of the control arm. The optimal design of the control arm was determined from Pareto solutions by entropy weighted grey relational analysis (GRA). The optimization results show that the mass of the control arm was reduced by 4.1% and the fatigue life was increased by 215.8% while its reliability increased by 7.8%. The proposed multi-objective reliability optimization method proved to be feasible and effective for lightweight design of a suspension control arm.

1. Introduction

With the continuous growth of car ownership, automobile exhausts emit a large number of pollutants. Additionally, the fuel consumption problem is becoming increasingly serious. In order to reduce the environmental and energy problems caused by the development of the automobile industry, making an automobile lightweight has become one of the effective measures [1,2,3]. A suspension control arm has to bear and transfer various loads caused by road roughness and engine vibrations. Thus, its thickness is always set to be thicker to ensure vehicle safety, which leads to the weight of the suspension system also being greatly increased. Therefore, it is of great significance to carry out a lightweight design of the control arm on the premise of ensuring its stiffness, strength, and fatigue life [4,5,6].

Scholars have paid much attention to automobile lightweight solutions, including optimization design, the application of lightweight material, and advanced manufacturing processes. Mohd et al. [7] optimized the topology of an aluminum control arm and proposed the optimal material layout method based on load and boundary conditions, which achieved a mass reduction of 41% under the condition of satisfying the strength. Yoo et al. [8] replaced the aluminum used for the control arm with carbon fiber-reinforced polymer and carried out topology optimization. At the same time, under the conditions of multiple constraints including stiffness and durability, the weight was reduced by 30%. Zhang [9] changed the material of automobile B-pillar from metal to carbon fiber-reinforced polymer, which not only reduced the mass by 50.76%, but also greatly improved the crashworthiness. Guler et al. [10] developed automobile hinge components with PA66 GF60 glass fiber-reinforced polyamide composite instead of ordinary steel. Topology optimization was carried out under a specific load to obtain the optimal size and minimum mass. Chen et al. [11] designed a new single stamping control arm and achieved a mass reduction of 33.3% compared to the original one through topology optimization and size optimization.

For engineering optimization problems, there are usually two or more objectives that need to be achieved. Generally, these objectives are in conflict with each other. Traditional single objective optimization often focuses on one key optimization performance while ignoring the correlation with the other performance. Multi-objective optimization methods are proposed to select the Pareto solutions in many engineering fields. Wang et al. [12] designed a new bumper system composed of negative Poisson’s ratio (NPR) beam and NPR absorber, and used a multi-objective optimization algorithm to search for the optimal structural parameters. The results show that the performance of pedestrian protection and vehicle crashworthiness was improved. Fossati et al. [13] used the multi-objective optimization method to optimize the passive suspension system of a whole vehicle model running on a random road, which achieved a reduction of 21.14% for the weighted mean square root of vertical angular velocity of the driver seat. Based on the Kriging surrogate model, Tang et al. [14] used a multi-objective genetic algorithm to optimize the thickness of acoustic package material. Under the condition of ensuring that the acoustic performance of automobile acoustic package remained unchanged, the mass was reduced by 20.76%. Nikkhah et al. [15] optimized the thin-walled tube with the aim of maximizing specific energy absorption (SEA) and minimizing peak crushing force (PCF). Response surface methodology (RSM) was established to find the optimal configuration of the tube for improving the energy absorption characteristics of the thin-walled tube.

Reliability is the probability that a product will not fail. The failure probability of the product is closely related with the uncertainty of design parameters [16,17]. The traditional optimization design is usually performed without considering the uncertainty of input variables, which may lead to failure of the optimization results in practical application. Therefore, reliability optimization has attracted much attention recently. Lim et al. [18] proposed a reliability-based multi-objective optimization design method to optimize the engine mount. They obtained a more reliable optimization solution than the deterministic optimization based on the NSGA algorithm. Zhang et al. [19] presented a new method for hybrid reliability-based design optimization under random and interval uncertainties, and used a ten-bar truss example and a design of piezoelectric energy harvester to verify the accuracy and effectiveness of the method. Combined with the response surface method and a multi-island genetic algorithm, Li et al. [20] applied an efficient reliability-based design optimization process to optimize the reliability of lattice boom. The results showed that it not only met the reliability requirements, but also ensured safety and economy.

In this paper, a multi-objective reliability optimization method by integrating MCS with NSGA-II algorithm coupled with entropy weighted GRA is proposed, and the effectiveness of this method is verified by the lightweight optimization of the lower control arm of Macpherson suspension. Firstly, the finite element model of the control arm is established and the modal and stiffness analysis are carried out. Then, the inertial release method is used to calculate the stress distribution of the control arm under unit load. The load time histories of the control arm are extracted by dynamic simulations of the rigid-flexible coupling vehicle model. Accordingly, the nominal stress method is used to calculate the fatigue life of the control arm. On this basis, six design variables are defined for describing geometric dimension of control arm based on mesh morphing technology. The optimal Latin hypercube method is applied to generate the sample points for constructing Kriging surrogate models, which are used to describe the relationships between design variables of control arm and its structure performance including mass, fatigue life, stiffness, and mode. The multi-objective reliability optimization of the control arm is further performed by combining the NSGA-II algorithm and the Monte Carlo simulations. Finally, the entropy weighted grey relational analysis is adopted to determine the optimal solution from Pareto set for realizing lightweight design of the control arm.

2. Multi-Objective Reliability Optimization Method

2.1. Kriging Surrogate Model

Kriging is one of the most used estimation methods for spatial data interpolation, and it has been widely used to build surrogate models in automotive engineering in recent years [21,22,23]. The Kriging surrogate model covers the global trend and local nonlinearity of the response, and especially has high accuracy for predicting nonlinear response. It can analyze the trend and dynamic characteristics of known information [24,25,26].

The Kriging surrogate model usually includes a polynomial function and a random distribution. Its approximate function expression is:

$$\widehat{y}\left(x\right)=f^{\mathrm{T}}\left(x\right)\cdot \beta +z\left(x\right)$$

(1)

where $\beta ={\left[{\beta }_1,\dots ,{\beta }_p\right]}^{\mathrm{T}}$ represents the regression coefficient vector; $f^{\mathrm{T}}\left(x\right)$ represents a polynomial with design vector x; $z\left(x\right)$ represents a random distribution, which can be expressed as a random function with zero mean and standard deviation of $\sigma$.

The covariance of random distribution $z\left(x\right)$ is:

$$Cov\left[z\left(x_i\right),z\left(x_j\right)\right]={\sigma }^{2}R$$

(2)

where $R=\left[R\left(x_i,x_j\right)\right]$ is the spatial correlation equation of any two nondiagonal sample points $x_i$ and $x_j$ in the matrix, which plays an absolute role in the simulation accuracy.

Using Gaussian correlation equation to express $R\left(x_i,x_j\right)$ is:

$$R\left(x_i,x_j\right)=EXP\left(-{\displaystyle \sum _{k=1}^m{\lambda }_k{\left|x_{ik}-x_{jk}\right|}^2}\right)$$

(3)

where m is the number of design variables; ${\lambda }_k$ represents the correlation coefficient of the fitted surrogate model; and x_ik and x_jk represent the kth value of x_i and x_j, respectively.

After determining the correlation function, the estimated value $\widehat{y}\left(x\right)$ of the approximate response can be obtained. The estimated value $\widehat{\beta }$ of regression coefficient $\beta$ and the estimated value ${\widehat{\sigma }}^2$ can be expressed as:

$$\{\begin{array}{l}\widehat{y}\left(x\right)=f^{\mathrm{T}}\left(x\right)\widehat{\beta }+r^{\mathrm{T}}\left(X^{\ast }\right)R^{-1}\left(y-F\widehat{\beta }\right)\\ \widehat{\beta }={\left(F^{\mathrm{T}}{R}^{-1}F\right)}^{-1}{F}^T{R}^{-1}y\\ {\widehat{\sigma }}^2=\frac{{\left(y-F\widehat{\beta }\right)}^{\mathrm{T}}{R}^{-1}\left(y-F\widehat{\beta }\right)}{n}\end{array}$$

(4)

where $r^T\left(X^{\ast }\right)=\left[R\left(X^{\ast },x_1\right),R\left(X^{\ast },x_2\right),\dots ,R\left(X^{\ast },x_n\right)\right]$ is the correlation coefficient vector between prediction point $X^{\ast }$ and sample point x; F is the $\left(n\times p\right)$ design matrix.

The maximum likelihood estimation for parameter ${\lambda }_k$ is:

$$\underset{{\lambda }_k>0}{\max \left({\lambda }_k\right)}=-\frac{\left[n\ln \left({\widehat{\sigma }}^2\right)+\ln \left|R\right|\right]}{2}$$

(5)

2.2. Monte Carlo Simulation

Monte Carlo simulation (MCS) is a very important numerical analysis method guided by probability and statistics theory. The MCS method can be used to propagate the uncertainty in design variables to predicted responses. In MCS, random design variables are first generated according to statistical distribution. Their responses are usually predicted using surrogate models because the number of simulations is very large. Then, the probability distribution of responses can be obtained through statistical analysis. The process of MCS calculation using simple random sampling method is shown in Figure 1.

Based on MCS, the reliability of the design can be defined as:

$$R=1-P_f=1-\frac{N_F}{N}$$

(6)

where $P_f$ is the failure probability; N_F is the number of samples that does not satisfy the constraints; and N is the total samples.

2.3. Entropy Weighted Grey Relational Analysis

Grey relational analysis (GRA) is an important branch of grey system theory and has been widely applied to decision-making problems. Its basic idea is to map the discrete data to the geometric shape of the space by using the linear interpolation between adjacent points of the sequence. It provides the grey relational grade to judge the close relationship between sequences based on the distance between reference sequence and comparison sequence [27,28].

In order to eliminate the non-commensurable caused by different dimensions, it is necessary to normalize the response values before grey relational analysis. The normalization methods corresponding to different response characteristics are different. If the response value has a smaller-the-better characteristic, that is, the smaller the response value is, the better the performance is, the normalization calculation method is as follows:

$$x_i^{\ast }\left(k\right)=\frac{{\max }_k{x}_i\left(k\right)-x_i\left(k\right)}{{\max }_k{x}_i\left(k\right)-\mathrm{m}{\mathrm{in}}_k{x}_i\left(k\right)}$$

(7)

If the response value has a larger-the-better characteristic, that is, the larger the response value is, the better the performance is, the normalized calculation method can be expressed as follows:

$$x_i^{\ast }\left(k\right)=\frac{x_i\left(k\right)-\mathrm{m}{\mathrm{in}}_k{x}_i\left(k\right)}{{\max }_k{x}_i\left(k\right)-{\min }_k{x}_i\left(k\right)}$$

(8)

where $x_i\left(k\right)$ and $x_i^{\ast }\left(k\right)$ represent the original and normalized value of the ith response for the kth attribute, respectively; ${\max }_k{x}_i\left(k\right)$ and ${\min }_k{x}_i\left(k\right)$ are the maximum and minimum values of the kth attribute, respectively.

The response results after normalization are converted to values between [0, 1]. The larger the value is, the better the performance is, and 1 represents the optimal value. Define $x_0^{\ast }\left(k\right)$ = 1 as the reference sequence, and the normalized results as the comparison sequence. The deviation sequence of these two sequences is:

$${\Delta }_{0i}\left(k\right)=\left|x_0^{\ast }\left(k\right)-x_i^{\ast }\left(k\right)\right|$$

(9)

The calculation formula of grey relational coefficient between comparison sequence and reference sequence is as follows:

$$\{\begin{array}{l}\gamma \left(x_0^{\ast }\left(k\right),x_i^{\ast }\left(k\right)\right)=\frac{{\Delta }_{\min }+\zeta {\Delta }_{\max }}{{\Delta }_{0i}\left(k\right)+\zeta {\Delta }_{\max }}\\ {\Delta }_{\min }={\min }_i{\min }_k{\Delta }_{0i}\left(k\right)\\ {\Delta }_{\max }={\max }_i{\max }_k{\Delta }_{0i}\left(k\right)\end{array}$$

(10)

where ${\Delta }_{\min }$ and ${\Delta }_{\max }$ are the minimum and maximum of the deviation sequence, respectively; $\zeta \in \left[0,1\right]$ is the distinguishing coefficient, generally it equals to 0.5.

The grey relational grade is calculated by weighted summation of the grey relational coefficient of each attribute. The calculation formula is:

$${\mathsf{\Gamma }}_i={\displaystyle \sum _{k=1}^n{w}_k\gamma \left(x_0^{\ast }\left(k\right),x_i^{\ast }\left(k\right)\right)}$$

(11)

where n is the number of attributes; $w_k$ is the weight coefficient of the kth attribute, ${\displaystyle \sum _{k=1}^n{w}_k}=1$.

Generally, the relative importance of each attribute may be different, and their weights can be calculated by information entropy which denotes the degree of uncertainty of a random variable [29,30]. The projection value of the ith response for the jth attribute is formulated as:

$$p_{ij}=\frac{x_{ij}}{{\displaystyle \sum _i^n{x}_{ij}}}$$

(12)

where i = 1, ‧‧‧, n, n is the number of responses; j = 1, ‧‧‧, m, m is the number of attributes; and $x_{ij}$ represents the normalized value of the ith response for the jth attribute.

The entropy value of the jth attribute is:

$$e_j=-\frac{1}{\ln n}{\displaystyle \sum _{i=1}^n{p}_{ij}\ln p_{ij}}$$

(13)

Then, the weight coefficient of attribute can be calculated as:

$$w_j=\frac{d_j}{{\displaystyle \sum _{j=1}^m{d}_j}}$$

(14)

where $d_j=1-e_j$ represents the deviation degree of the jth attribute. Generally, a larger value of deviation degree indicates that it provides more information. Therefore, this attribute has a higher weight.

2.4. Multi-Objective Reliability Optimization

The traditional deterministic multi-objective optimization can be expressed by:

$$\{\begin{array}{l}\min f_m\left(x\right),m=1,\dots ,M\\ \mathrm{s}.\mathrm{t}{.\mathrm{g}}_k\left(x\right)\le 0,k=1,\dots ,K\\ \hspace{1em}{x}_{lb}\le x\le x_{ub}\end{array}$$

(15)

where x is the vector of design variables; $x_{lb}$ and $x_{ub}$ are the lower and upper limits of the design variable, respectively; $f_m\left(x\right)$ and ${\mathrm{g}}_k\left(x\right)$ are the optimization objectives and inequality constraints, respectively.

Deterministic optimization has been widely used in engineering design. However, it cannot consider the inevitably uncertain factors such as material properties and loading conditions, which may lead the optimal solutions near the constraint boundary to failure [31,32]. To overcome this drawback, the reliability optimization was proposed, as shown in Figure 2. The deterministic optimization point is the point optimized without reliability constraint and the reliability optimization point is the point optimized with reliability constraint.

Considering the uncertainty of design variables, reliability optimization was added to deterministic multi-objective optimization. It can be described as:

$$\{\begin{array}{l}\min f_m\left(x\right),m=1,\dots ,M\\ \mathrm{s}.\mathrm{t}.R\left[G_k\left(x,y\right)\le 0\right]\ge R_j,k=1,\dots ,K\\ \hspace{1em}{h}_j\left(x\right)\le 0,j=1,\dots ,J\end{array}$$

(16)

where $f_m\left(x\right)$ is the mth objective function; x is a vector of deterministic design variables; y is a vector of random variables; R is the reliability; R_j represents expected reliability; $G_k\left(x,y\right)$ is the probability constraint; and $h_j\left(x\right)$ is the deterministic constraint.

Accordingly, this paper proposes a multi-objective reliability optimization method based on Monte Carlo simulation and NSGA-II (elitist non-dominated sorting genetic algorithm) coupled with entropy weighted grey relational analysis. The Kriging surrogate models are employed to predict the response for improving optimization efficiency, for which the optimal Latin hypercube design is adopted to sampling. The proposed multi-objective reliability optimization procedure is shown in Figure 3.

3. Finite Element Analysis of Suspension Control Arm

3.1. Finite Element Modeling of Control Arm

For building the finite element model, the shell element was used to mesh the model of the control arm, which was discretized into 11,193 elements and 11,412 nodes. The front point and rear point of control arm were connected to the body using bushings, and the outer point was connected to the steering knuckle by a ball joint. These installation points were simulated by RBE2 elements. The finite element model of the lower control arm of the Macpherson strut was established, as shown in Figure 4.

3.2. Modal Analysis of Control Arm

The vibration performance of the suspension control arm has a great influence on the comfort, and is also related to the fatigue failure of the structure. Modal analysis was mainly used to calculate the vibration frequency and vibration mode of the structure. The modal analysis of the finite element model of the control arm was carried out under unconstrained free state, and the first six natural frequencies were extracted, as shown in

Table 1. Natural frequency of the first six modes of control arm.

Order	Natural Frequency/Hz
1	324.1
2	357.1
3	602.5
4	1003.9
5	1254.2
6	1439.4

.

3.3. Stiffness Analysis of Control Arm

The longitudinal stiffness and transverse stiffness of the control arm is related to its deformation under load, which affects the ride comfort and handling stability of the vehicle. In the stiffness analysis model of control arm, the translation in X, Y, Z direction and rotation around Y, Z axis of front point were constrained (marked as 1, 2, 3, 5, 6), and three translational degrees of freedom and rotation around X, Y axis of rear point were restrained (denoted as 1, 2, 3, 4, 5). The translation in Z direction of outer point was fixed (denoted as 3), while a load of 2000 N was applied in positive X and Y direction, respectively, so as to calculate the longitudinal stiffness and transverse stiffness of control arm. The boundary condition of stiffness analysis of control arm is shown in Figure 5, and the calculation results of stiffness are given in

Table 2. Results of stiffness analysis of control arm.

Direction	Force/N	Displacement/mm	Stiffness/N·mm−1
X+	2000	0.156	12,821
Y+	2000	0.077	25,974

.

4. Fatigue Life Analysis of Suspension Control Arm

4.1. Cumulative Fatigue Damage Theory

It is well known that materials subjected to alternating load may lead to fatigue failure even at a quite low stress level. The fatigue can be divided into high cycle fatigue and low cycle fatigue. The number of high cycle fatigue cycles is usually greater than 10⁵–10⁷ times, and the alternating load is generally less than the yield limit of the material. Low cycle fatigue refers to the number of cycles is usually less than 10⁴–10⁵ times, and the alternating load is close to or higher than the yield limit of the material [33].

The fatigue cumulative damage theory includes linear, nonlinear, and bilinear fatigue cumulative damage theory. Miner cumulative damage theory is a representative linear damage theory. It assumes that the damage of component under every stress is independent. Without considering the loading sequence, the fatigue damage can be linearly superimposed, which can be expressed as:

$$D={\displaystyle \sum _{i=1}^m\frac{n_i}{N_i}}$$

(17)

where D is the total fatigue damage, which equals 1, meaning fatigue failure; m is the number of stress levels; n_i is the number of cycles accumulated at the ith stress level; and N_i is the fatigue life at the ith stress level.

4.2. Stress Analysis of Control Arm

The stress analysis of the lower control arm was carried out using the inertia release method. The control arm mainly transfers the loads between the wheel and body through three installation points. Hence, nine load cases, three cases for each connection point, were employed to predict the fatigue life of the control arm. A force of 1 N was applied for each load case. The stress distribution of the control arm for the nine load cases were obtained through finite element analysis. The stress contour of the control arm under the unit load in X direction of the outer point is shown in Figure 6.

4.3. Load Spectrum of Control Arm

The durability road of automotive proving ground is usually used to provide a load spectrum for estimating the fatigue performance of vehicle components because of the greatly reduced test cycle. The fatigue life of the control arm was predicted according to a virtual durability road, including a pebble road of 40 m, a washboard road of 40 m, and a Belgian road of 80 m [34]. Then, the dynamic simulation of a rigid-flexible coupling dynamic vehicle model, running on the virtual durability road, was performed to extract load time histories of the control arm. The dynamic load time histories of the outer point, the front point, and the rear point of the control arm were acquired, as shown in Figure 7.

4.4. Fatigue Life Prediction of Control Arm

The number of fatigue cycles of the suspension control arm is usually greater than 10⁵, which belongs to the high cycle fatigue failure. Thus, the nominal stress method was adopted to predict its fatigue life. It was firstly necessary to determine the S-N curve, in which one axis denotes stress and the other axis is the number of cycles. The relationship between fatigue life and stress can be expressed by S-N curve, and its mathematical expression is:

$${\sigma }_a={\sigma }_f^{\prime }(2N_f{)}^b$$

(18)

where ${\sigma }_a$ represents the stress amplitude; ${\sigma }_f^{\prime }$ represents the fatigue strength coefficient; $N_f$ represents the number of load cycles; and b is fatigue strength index. The material of the control arm was QSTE420TM, its yield strength was 420 MPa, and its ultimate strength was 520 MPa. Then, the S-N curve was estimated using the material mapping, as shown in Figure 8.

The S-N curve is generally measured under the stress ratio of −1, which means the compressive stress equals the tensile stress. However, the stress ratio of control under load spectrum usually may not be −1. In a certain range, the compressive stress increases the fatigue limit while the tensile stress decreases the fatigue limit. Therefore, it is necessary to correct the effect of stress ratio on fatigue damage when the mean stress is not zero. Goodman theory can be used to correct the mean stress, that is:

$$\frac{S_a}{S_e(\mathrm{R}=-1)}+\frac{S_m}{\mathrm{UTS}}=1$$

(19)

where S_a and S_m are the stress amplitude and mean stress of the material under load spectrum, respectively; UTS is the ultimate tensile strength of the material; and S_e is the stress amplitude when the stress ratio is −1.

The load spectrum of the above nine cases were integrated with stress analysis results of the control arm under unit load to calculate the dynamic stress time histories. The Goodman method and the survival rate of 95% were adopted to conduct the S-N curve correction. The nominal stress method was then used to predict the fatigue life of the control arm, as shown in Figure 9.

It can be seen from Figure 9 that the minimum number of fatigue cycles of the control arm was 1.393 × 10⁶, which appeared near the front point connection. According to the length of established road model, the fatigue life durability mileage of the control arm was 1.393 × 10⁶ × 0.16 = 222,880 km. By considering the enhancement coefficients of the pebble road, washboard road, and Belgian road, which took the value of 5.8, 4.9, and 12.6, respectively [35], the enhancement coefficient of the durability road model was calculated as 5.8 × 0.25 + 4.9 × 0.25 + 12.6 × 0.5 = 8.975. Thus, the minimum life mileage of the control arm was 222,880 × 8.975 = 2,000,348 km. The life mileage had a large margin, and there is potential for lighter weight.

5. Multi-Objective Reliability Optimization of Suspension Control Arm

5.1. Design Variables

The parameterized finite element model of the control arm was developed based on mesh morphing technology. Five shape variables and a thickness variable (denoted as x₁, x₂, ···, x₆) were defined, as shown in Figure 10. The mesh morphing of the parameterized model of the control arm was realized by translating and scaling the control node. The value ranges of design variables are given in

Table 3. Design variables and value range of control arm.

Variables	Variable Description	Initial Value	Lower Limit	Upper Limit
x1/mm	Variation of flanging height	0	−5	5
x2	Width scaling at outer point	1	0.9	1.1
x3	Width scaling at front point	1	0.9	1.1
x4	Width scaling at rear point	1	0.9	1.1
x5	Scaling of hole diameter	1	0.9	1.1
x6/mm	Thickness	5	4	6

.

5.2. Kriging Surrogate Models of Control Arm

Surrogate model technology is a method used to approximate the relationship between a set of input variables and output variables by mathematical model. The sampling method has a great influence on the accuracy of the surrogate model. The optimal Latin hypercube design can sample uniform points throughout the design space, so that the fitting of factors and responses is more accurate [36].

Therefore, the optimal Latin hypercube sampling method was adopted to generate 41 sample points for the six design variables, which were employed to build the Kriging surrogate models of the control arm. The performance indexes were obtained through finite element analysis of the control arm. The accuracy was verified by the determination coefficient R², expressed as:

$$R^2=\frac{{\displaystyle \sum _{i=1}^n({\widehat{y}}_i-\overline{y})}}{{\displaystyle \sum _{i=1}^n(y_i-\overline{y})}}$$

(20)

where n is the number of sample points used to verify the accuracy of Kriging model; ${\widehat{y}}_i$ is the predicted value of the ith response; y_i is the simulation value of the ith response; and $\overline{y}$ is the average value.

The closer to 1 the determination coefficient R² is, the higher the prediction accuracy of the surrogate model is. Ten sample points were selected to verify the accuracy of the Kriging models. The R² of the fatigue life, mass, first-order natural frequency, longitudinal stiffness, and transverse stiffness of the control arm were 0.9808, 0.9830, 0.9821, 0.9672 and 0.9936, respectively. The Kriging surrogate model accuracy verification results of fatigue life and mass are shown in Figure 11.

5.3. Deterministic Multi-Objective Optimization of Control Arm

The control arm is expected to have sufficient life without fatigue failure, while having a small mass. Further considering the performance constraints of longitudinal stiffness, transverse stiffness, and first-order natural frequency, the deterministic multi-objective optimization of the control arm can be formulated as:

$$\{\begin{array}{l}\min (m(x),-N(x))\\ findx={\left(x_1,\dots ,x_2\right)}^{\mathrm{T}}\\ \mathrm{s}.\mathrm{t}.\left|k_x(x)-k_{x0}\right|\le 0.1k_{x0}\\ \hspace{1em} \left|k_y(x)-k_{y0}\right|\le 0.1k_{y0}\\ \hspace{1em} \left|f_t(x)-f_{t0}\right|\le 0.1f_{t0}\\ \hspace{1em} x\in (x_L,x_U)\end{array}$$

(21)

where x is the design variable vector of the control arm; $m(x)$ is the control arm mass; $N(x)$ is the minimum fatigue life; $k_x(x)$, and $k_{x0}$ are the longitudinal stiffness and its initial value, respectively, $k_{x0}$ = 12,821 N·mm⁻¹; $k_y(x)$ and $k_{y0}$ are the transverse stiffness and its initial value, respectively, $k_{y0}$ = 25,974 N·mm⁻¹; $f_t(x)$ and $f_{t0}$ are the first order natural frequency and its initial value, $f_{t0}$ = 324.1 Hz; and x_L and x_U are the lower and upper limits of design variable, respectively.

The NSGA-II algorithm with a population size of 40 and generation of 100 was applied to solve the optimization of control arm. The Kriging surrogate models were used to calculate the performance indexes in the optimization process. Then, the Pareto front was obtained, as shown in the Figure 12. The optimal design parameters were selected from the Pareto solutions by considering that the mass was minimized as far as possible while meeting the fatigue life requirements. According to engineering practice, the optimization results of design variables were corrected, as shown in

Table 4. Deterministic multi-objective optimization results.

Variables	Optimal Value	Corrected Value	Responses	Value
$x_1$/mm	−4.9711	−5.0	Mass/kg	2.88
$x_2$	0.94625	0.95	Life/cycles	1,080,000
$x_3$	1.0955	1.09	First order natural frequency/Hz	333.84
$x_4$	0.90005	0.90	Longitudinal stiffness/N·mm−1	11,765
$x_5$	0.98121	0.98	Transverse stiffness/N·mm−1	23,810
$x_6$/mm	4.9934	5.0

. The corrected design parameters were used to modify the finite element model of the control arm, and its performance indexes were obtained, also listed in Table 4. The mass of the optimized control arm decreased from 3.16 kg to 2.88 kg while its fatigue life was 1.08 × 10⁶, meeting the design requirement, as well as the other performance indexes.

5.4. Multi-Objective Reliability-Based Optimization

The fatigue life of the control arm is sensitive to the uncertainties of the six design variables. Therefore, a multi-objective reliability-based optimization of the control arm was developed based on deterministic optimization by considering the reliability constraint of fatigue life, in which the Monte Carlo simulation was applied to calculate the reliability. The deterministic optimization results of design variables were determined as the mean values of the design parameters of the control arm, which were taken as probability random variables subjected to normal distribution, as shown in

Table 5. Design variable probabilistic distributions.

Variables	Mean	Standard Deviation	Distribution
$x_1$/mm	−5.0	0.75	Normal
$x_2$	0.95	0.14	Normal
$x_3$	1.09	0.16	Normal
$x_4$	0.90	0.14	Normal
$x_5$	0.98	0.15	Normal
$x_6$/mm	5.0	0.75	Normal

. Considering the reliability of fatigue life more than 10⁶ cycles, the multi-objective reliability optimization can be formulated as:

$$\{\begin{array}{l}\min (m(x),-N(x))\\ \mathrm{find}x={\left(x_1,\dots ,x_2\right)}^{\mathrm{T}}\\ \mathrm{s}.\mathrm{t}.R\left[N\left(x\right)\ge {10}^6\right]\ge R_j\\ \hspace{1em} \left|k_x(x)-k_{x0}\right|\le 0.1k_{x0}\\ \hspace{1em} \left|k_y(x)-k_{y0}\right|\le 0.1k_{y0}\\ \hspace{1em} \left|f_t(x)-f_{t0}\right|\le 0.1f_{t0}\\ x\in (x_L,x_U)\end{array}$$

(22)

where R is the reliability; and R_j represents the expected reliability, set to 95%.

There were 111 Pareto solutions obtained by solving the reliability-based optimization based on NSGA-II algorithm and Monte Carlo simulations. Then, the entropy weighted grey relational analysis was employed to determine the optimal design of the control arm. The Pareto solutions were defined as the decision-making matrix, which were used to calculate the grey relational coefficients, and some of them are listed in

Table 6. Grey relational coefficient.

	Grey Relational Coefficient					Weight
	1	2	‧‧‧	110	111
Mass	0.6775	0.4280	‧‧‧	0.1730	0.0954	0.2837
Fatigue life	0.3534	0.5154	‧‧‧	0.0454	0.1197	0.1210
First order natural frequency	0.3678	0.7201	‧‧‧	0.0206	0.0408	0.0450
Longitudinal stiffness	0.3774	0.7262	‧‧‧	0.0227	0.0446	0.0529
Transverse stiffness	0.4079	0.5795	‧‧‧	0.0365	0.0932	0.0946
Reliability	0.3520	0.3600	‧‧‧	0.2364	0.1450	0.4027

. Then, the weights of performance indexes were obtained using the entropy method, as listed in Table 6. Finally, the grey relational grade was calculated, as shown in Figure 13.

It can be seen from Figure 13 that the grey relational grade of the third Pareto solution was the maximum, that is, this grey relational grade was the closest to the expected value. Therefore, the third Pareto solution was selected as the optimal design of the control arm after multi-objective reliability optimization. The design variables obtained by the reliability-based optimization were corrected according to practical engineering, as shown in

Table 7. Multi-objective reliability optimization results.

Variables	Optimal Values	Corrected Values
$x_1$/mm	−4.9471	−4.9
$x_2$	0.90029	0.90
$x_3$	1.09890	1.10
$x_4$	0.90091	0.90
$x_5$	0.90021	0.90
$x_6$/mm	5.3401	5.3

. The finite element model of the control arm was modified using these design parameters, and the fatigue life contour is given in Figure 14.

The comparison of the performance indexes before and after multi-objective reliability optimization of the control arm is shown in

Table 8. Comparison of optimal design results.

	Initial Value	Deterministic Optimization	Reliability Optimization
Mass/kg	3.16	2.88	3.03
Life/cycles	1.39 × 106	1.08 × 106	4.39 × 106
First order natural frequency/Hz	324.09	333.84	356.49
Longitudinal stiffness/N·mm−1	12,821	11,765	12,658
Transverse stiffness/N·mm−1	25,974	23,810	26,316
Reliability of fatigue life		0.90	0.97

. Compared with the initial steel control arm, after the deterministic optimization, the mass of optimized control arm was reduced by 8.89% from 3.16 to 2.88 kg and the life was reduced by 22.47% from 1.39 × 10⁶ to 1.08 × 10⁶, with a reliability of 90%. After the multi-objective reliability optimization of the control arm, longitudinal stiffness changed little, the first mode frequency and transverse stiffness increased, and the mass changed from 3.16 kg to 3.03 kg. The mass decreased by 4.1% while the life increased by 215.8%, with a high reliability of 97%.

6. Conclusions

This paper proposed a multi-objective reliability optimization method based on Monte Carlo simulation, NSGA-II algorithm and entropy weighted grey relational analysis, applied to perform the lightweight design of a lower control arm of McPherson suspension. The finite element model of the control arm was developed to conduct the modal and stiffness analysis. The fatigue life of the control arm was predicted by a normal stress method, using the load spectrum acquired by dynamic simulations of the rigid-flexible coupling vehicle model running on the virtual durability road. Then, the deterministic multi-objective optimization of the control arm was carried out to minimize the mass and maximize the fatigue life with the constraints of first order natural frequency and longitudinal and transverse stiffness, based on the Kriging model and NSGA-II algorithm. The deterministic optimization results showed that the mass of the control arm was reduced from 3.16 kg to 2.88 kg with a reduction of 8.89%, and the life was reduced by 22.47% from 1.39 × 10⁶ to 1.08 × 10⁶, with a reliability of 90%.

On this basis, the multi-objective reliability optimization of the control arm was conducted using the NSGA-II algorithm and the Monte Carlo simulations. The entropy weighted grey relational analysis was adopted to determine the optimal design of the control arm from the Pareto solutions. The reliability optimization results indicated that the mass of the optimized control arm was reduced by 4.1% and the life was increased by 215.8% with a high reliability of 97%. This proves that the proposed multi-objective reliability optimization procedure is effective for automotive lightweight design.