Multiobjective Optimization Design for Lightweight and Crash Safety of Body-in-White Based on Entropy Weighted Grey Relational Analysis and MNSGA-II

In order to improve the lightweight level, crash safety performance and optimization design efficiency of body-in-white (BIW), this article proposes a lightweight multi-objective optimization design method for mixed-material body. The implicit parametric model of the BIW is created by using SFE-CONCEPT software, and the validity and correctness of the model are verified by tests. Material, shape and dimension parameters are introduced as design variables for design of experiments (DOEs), and 26 important design variables are screened out by combining contribution analysis with nonlinear main effects analysis. The approximate model method is used to fit the Kriging surrogate model and the RBF surrogate model, and it is found that the RBF surrogate model can better reflect the relationship between nonlinear crash performance and optimization variables. A hybrid method combined entropy weighted grey relational analysis (EGRA) with modified non-dominated sorting genetic algorithm (MNSGA-II) is proposed to carry out the lightweight multi-objective optimization of BIW in front crash and side impact, which improves the population diversity of multi-objective optimization problems and quantifies the comprehensive performance of each scheme. Comparing and analyzing the optimization platform recommending scheme, the technique for order preference by similarity to an ideal solution (TOPSIS) method preferring scheme and the EGRA method optimum scheme, it is found that EGRA method can obtain the optimal compromise scheme, and the performance improvement of the BIW are more obvious and the improvement rates are also more balanced. The results verify the feasibility of the ranking method, avoid the blindness of optimal solution selection, and establish an objective evaluation method of multi-objective optimization results. The optimization results show that the improvement rates of the BIW lightweight coefficient, the average value of the maximum acceleration of the B-pillars on both sides during the front crash, and the maximum intrusion displacement of the B-pillar chest during the side impact have reached 11.5%, 6.5%, and 6.8%, respectively. Other performance response improvement rates are also above 3.3%, the lightweight and crash safety performance are significantly improved.


I. INTRODUCTION
Lightweight technology is one of the main ways to achieve energy saving and emission reduction of cars. The The associate editor coordinating the review of this manuscript and approving it for publication was Hamed Azami .
body-in-white (BIW) accounts for 20%-30% of the car's mass, and its lightweight degree plays an important role in the lightweight of cars [1], [2]. The BIW structure is not only an important load-bearing component but also subjected to multi-directional impacts from the road surface, and the stress situation is complex. At the same time, the BIW is also an important guarantee of occupant survival space, and its lightweight design directly affects the crash safety performance of the car [3]- [5]. Therefore, the lightweight design of BIW is a multi-disciplinary and multi-objective optimization design project, and it is of great significance to establish an effective multi-objective optimization design method for lightweight and crash safety of BIW.
In the early design stage of a car, it is challenging to define complex design variables due to the limitations of mesh quality and connection relationships between components. Implicit parametric technology is used to define shape variables, which is provided by computer-aided engineering software and has great shape variation and regrinding ability [6]- [8]. Wang and Wang et al. [9], [10] established an implicit parametric model of BIW using the implicit parametric method. By comparing the prediction results with the test results, it is verified that the established implicit parametric model of BIW can be used in the design and development of the BIW conceptual design stage, and then carried out a lightweight multi-objective collaborative optimization for BIW structure. Wang and Lv et al. [11] proposed a new front-end submodule lightweight optimization method, constructed the implicit parametric coupling model of BIW, and the concept of analysis-oriented design was realized by comprehensively considering the construction efficiency and car performance of the BIW. Yao et al. [12] optimized the beam structure parameters and improved the static and dynamic performance of the body by establishing the implicit parametric model of the typical body-in-prime structure.
BIW is a typical frame structure composed of a large number of thin-walled structures, and its lightweight and crashworthiness optimization design is a typical multi-design variable nonlinear dynamic response optimization problem. Long et al. [13] designed a lightweight method to reduce the mass of the extended-range electric vehicle (E-REV) body and key components based on the rear-end crash failure analysis, and evaluated the rear-end crash safety performance of the E-REV through research. By optimizing the material layout of BIW, Park and Choi [14] reduced the mass of BIW by 45.7% by implementing the metamodel-based optimization strategy under the condition of meeting the requirements of bending and torsional stiffness. Wang et al. [15] proposed a lightweight design scheme of variable crosssection multi-material car based on crash safety, and taking a certain type of racing car frame as the research object, a racing car frame with good crash safety performance was obtained. Ou et al. [16] reduced the mass by 50.65 kg by adopting the design-driven method and considering the crashworthiness, NVH and static stiffness performance of BIW.
Gray relational analysis (GRA) is a special design of experiments (DOEs) method, which was first proposed by Deng [17] in the 1980s to determine the relationship between design sequence and ideal sequence. The relational degree can be represented by grey relational grade (GRG), with better solutions having larger GRG. Therefore, GRG can be used as an index to evaluate multi-objective problems, and GRA is used to determine the ranking order of all design sequences [18], [19]. Xiong et al. [20] combined the GRA and principal component analysis to reduce the total mass of the car body by 4.54 kg through multi-objective optimization of the side structure of the car body. Chen et al. [21] proposed an analysis and prediction method for interior noise and sound quality based on grey system theory, and determined the relationship between subjective evaluation and psychoacoustic parameters through GRA. Based on ADAMS software and car vibration model, Yarmohammadisatri et al. [22] obtained the best optimization result of car suspension through GRA. Wang et al. [23] optimized the spring stiffness and damping of the full floating cab suspension of the dump truck based on GRA, and the ride comfort of the dump truck was improved by 15.69% after optimization.
The structure design of BIW is driven by many competing criteria that must be satisfied simultaneously during the optimization process and requires the solution of complex multi-objective optimization problems. Many researchers use different evolutionary algorithms to solve multi-objective problems with a set of non-dominated optimal design solutions as the Pareto frontier [24]. One of the most popular evolutionary algorithms is the non-dominated sorting genetic algorithm (NSGA-II) proposed by Srinivas and Deb [25]. NSGA-II is a Pareto-based method that generates Pareto frontiers where non-dominated solutions perform efficiently on at least one standard [26], [27]. Since NSGA-II has some disadvantages when dealing with multi-objective optimization problems with more than two objectives, the performance of NSGA-II will be degraded [28]. To improve NSGA-II, Nariman-Zadeh et al. [29] proposed a modified non-dominated sorting genetic algorithm (MNSGA-II), which uses the ε-elimination algorithm instead of the crowding distance. Jozefowiez et al. [30] proposed a standard multi-objective evolutionary algorithm NSGA-II, which was tested on a two-objective car routing problem, and the generated Pareto set results were improved. Bian et al. [31] proposed an MNSGA-II, which improves the search efficiency of the algorithm through the interpolation and elimination strategy based on crowding distance, and uses the dynamic depth search method to balance the global and local search capabilities of the algorithm.
Although MNSGA-II can provide designers with a large number of non-dominated optimal solutions, designers still need to use engineering knowledge to independently choose the best compromise solution. In order to avoid the influence of some subjective factors, many scholars combine NSGA-II and the nearest ideal point (NIP) method or through the technique for order preference by similarity to an ideal solution (TOPSIS) for optimal design. Kilpeläinen [32] optimized the sandwich panel structure using MNSGA-II and obtained the best compromise design using NIP and TOPSIS methods. Jiang et al. [33] proposed a multi-objective optimization method combining NSGA-II and entropy weighted TOPSIS method for the lightweight design of dump truck carriage. Shojaeefard et al. [34] used MNSGA-II to determine the optimal combination of design variables for a three-cylinder engine and the TOPSIS method to determine the optimal solution between non-dominant design points. Wang et al. [35] proposed a hybrid method combining the MNSGA-II and TOPSIS methods, and obtained the best compromise solution for the lightweight design of the front subframe of the car.
In addition, it also has great application prospects by combining multiple optimization methods in the multi-objective optimization design of the BIW. Wang et al. [36] proposed an optimization method combining radial basis function (RBF) neural network model, fuzzy subtractive clustering sequential sampling method and NSGA-II, which is used to improve the calculation efficiency and accuracy of multi-objective optimization problems in engineering. Xiong et al. [37] proposed a hybrid method combining contribution analysis, RBF neural network -response surface method hybrid surrogate modeling method, multi-objective particle swarm optimization algorithm, and TOPSIS method for lightweight design of body front-end structures. Wang et al. [38] proposed an optimization method combining contribution analysis, entropy weighted TOPSIS, DOEs and GRA to optimize the design of front-end body safety parts, which significantly improved the crash safety performance of the car and reduced the mass of safety parts by 11.02%.
The above research has improved the level of lightweight and crashworthiness of BIW from the aspects of implicit parametric of BIW, lightweight and crashworthiness design, GRA or NSGA-II improvement. Based on the above research and the implicit parametric modeling technology of BIW, this article proposes a design of hybrid material body, which takes the material, shape and dimension parameters as the design variables of BIW at the same time. Combining with contribution analysis, DOEs, TOPSIS and EGRA method, MNSGA-II is used for multi-objective optimization design of BIW lightweight and crashworthiness to maximize the lightweight design of BIW and improve the comprehensive performance.

A. ENTROPY WEIGHTED GREY RELATIONAL ANALYSIS
GRA is a method that uses GRG to measure the degree of approximation between an experimental sequence and an ideal sequence. Since each performance index in the test data corresponds to a different dimension, it is necessary to normalize each performance response to dimensionless data between 0 and 1 before conducting GRA to facilitate the quantitative analysis of each performance response.
Depending on the characteristics of the objective function, the normalization methods used are also different. If the objective has the characteristic of 'the larger, the better', then the normalization method can be expressed as: If the objective has the characteristic of 'the lower, the better', then the normalization method can be expressed as: If the objective has the characteristic of 'closer to a specific value, the better the performance', then the normalization method can be expressed as: where x * i (k) is the normalized value of the ith response in the kth objective function; x i (k) is the initial value of the objective function; max k x i (k) and min k x i (k) are the maximum and minimum values of the kth objective function, respectively; T is the specific value.
After grey relational generation, the corresponding grey relational coefficient (GRC) can be calculated as: is the designed experimental scheme; 0i (k) is the absolute difference between x * 0 (k) and x * i (k); max and min are the maximum and minimum values of 0i (k), respectively; ζ is the discrimination coefficient, ζ ∈ [0,1].
GRG is calculated by weighted summation of GRC of each objective function, and the calculation formula of GRG is: where n is the number of objective functions; ω k is the weight coefficient of the kth objective function, n k=1 ω k = 1. Usually, since the importance of each objective response may be different, their weights can be calculated by the information entropy representing the uncertainty degree of the random variable, and the information entropy of the kth objective function is: (6) where m is the number of responses, i = 1, 2,. . . m; n is the number of objective functions, k = 1, 2,. . . n; x ik is the normalized value of the ith response in the kth objective function.
The weight coefficient of the objective function can be calculated as: where d k is the degree of deviation of the kth objective function, d k = 1-e k . Through the above calculation, the GRG of each design solution can be obtained, and the higher GRG value represents the better comprehensive performance of the optimized solution, and vice versa, the worse comprehensive performance. VOLUME 10, 2022

B. APPROXIMATE MODEL METHOD 1) KRIGING SURROGATE MODEL
The Kriging method is one of the most commonly used estimation methods in spatial data interpolation and has been widely used in recent years to build approximate models in car engineering. The kriging surrogate model covers the global trend and local nonlinearity of the response and has high accuracy in predicting the nonlinear response. The Kriging surrogate model includes the polynomial global approximate model and the random deviation, whose expressions are: where f T (x) is a polynomial with a design vector x; β is the regression coefficient vector, β = [β 1 , β 2, . . . , β n ] T ; z(x) is the random deviation, which can be expressed as a random function with mean zero and standard deviation σ . The covariance of the random deviation z(x) is: where R is an n×n symmetric correlation matrix of order 1 on both diagonals; R(x i , x j ) is the correlation function of any two sampling points x i and x j in the n training samples. Using the Gaussian correlation function R(x i , x j ) can be expressed as: where m is the number of design variables; θ k is the correlation coefficient used to fit the approximate model; x ik and x jk are the kth elements of sample points x i and x j , respectively.

2) RBF SURROGATE MODEL
RBF neural network is a typical pre-feedback control algorithm, which mainly includes input layer, hidden layer and output layer, and its structure is shown in Fig 1. Among them, the input layer is mainly used to introduce variable information, and the number of variables determines the dimension of the input layer; the hidden layer is the core part of the RBF neural network algorithm, which maps the input information to the output layer according to a certain mathematical relationship, and the number of neurons in the hidden layer is directly related to the RBF's learning ability, fault tolerance and approximate accuracy; the output layer outputs the corresponding structural response through the linear mapping function between the hidden layer and the output layer to complete the optimization calculation process.
The RBF method uses the radial basis function as the transfer function, and uses the radial distance r from the input sample x to the center of the hidden layer as a variable for mapping. The specific mathematical expression is as follows: where f (x) is approximate of the true response value f (x); x =[x 1 , x 2 ,. . . , x n ] is the n-dimensional input vector; n is the number of basis functions; m is the number of output responses; W mk is the connection weight between the kth hidden layer node and the mth output layer node; b m is the deviation of the mth output response; v k (x) is the radial basis function that represents the distance between sample x and the ith sample x i in the design space.

C. MODIFIED NON-DOMINATED SORTING GENETIC ALGORITHM
NSGA-II was proposed by Srinivas and Deb [25] after introducing elite retention strategy, fast non-dominated sorting method and crowding distance comparison method on the basis of the NSGA algorithm. NSGA-II algorithm avoids the defect of NSGA algorithm due to the lack of diversity of populations due to shared fitness, and improves the optimization efficiency and computational accuracy. Although the crowding distance comparison method can better maintain the diversity of the population, it is still inadequate for optimization problems with more than two objective functions. In order to further improve the population diversity of multi-objective optimization problems, a fixed threshold elimination strategy is used to replace the crowding distance comparison method in MNSGA-II, which can solve multi-objective optimization problems more reasonably. The principle of MNSGA-II is shown in Fig 2. The main steps in the MNSGA-II are presented as follows: Step 1: Initialize a random parent population P t of size N according to the range of variables and constraints of the multi-objective optimization problem, and calculate the individual fitness.
Step 2: Sort the parent individuals into several Pareto fronts according to the non-domination sorting criteria.
Step 3: Generate the offspring population Q t using the selection, crossover and mutation mechanisms embedded in the genetic algorithm.
Step 4: Merge the parent population P t and the offspring population Q t to obtain a mixed population R t with a population number of 2N , and perform non-dominated sorting on the individuals of the mixed population R t .
Step 5: According to the fixed threshold ε elimination strategy, remove the bad individuals from the mixed population R t , and then randomly generate new individuals to fill the population R t , so that the population number remains unchanged at 2N .
Step 6: Perform a non-dominated sorting operation on the population R t , selecting the best individuals in order from the lower to the higher levels until the number of populations obtained is N parent populations P t+1 .
Step 7: The parent population P t+1 is subjected to selection, crossover, mutation and fixed threshold ε elimination operations to obtain the offspring population Q t+1 .
Step 8: If the termination condition is satisfied, the iteration ends. Otherwise, the optimization is repeated from step 4 until the non-dominated solution satisfies the termination condition.

A. IMPLICIT PARAMETRIC MODELING OF BIW
According to the implicit parametric modeling technology, the implicit parametric model of the BIW is established by using SFE-CONCEPT software. In the process of implicit parametric modeling of BIW, the model is mainly created according to the ideas of influence point (IP), baseline (BL), base section (BS), beam structure, connection joint, free-form surface and structural features. Among them, IP, BL and BS are the most basic parametric elements. IP is arranged in the structure space according to the geometric characteristics of the BIW, BL and BS can be defined through IP, BS is stretched into beam structure along the direction of BL, and finally the transition connection of different beam structures is made through joints to complete the creation of parametric parts. The implicit parametric modeling process of BIW is shown in Fig 3. The assembly relationship between components in the implicit parametric model of the BIW is connected through mapping technology. Through the mapping technology, each part can realize the information transfer between geometric parameters, maintain the topological relationship between parts, and effectively avoid geometric conflicts and mesh deformation. The implicit parametric model of the BIW is shown in Fig 4.

B. IMPLICIT PARAMETRIC MODEL VALIDATION OF BIW
In order to ensure the reliability of the subsequent BIW optimization design, the implicit parameterized BIW model

1) BENDING STIFFNESS SIMULATION ANALYSIS AND TEST VERIFICATION
Apply constraints and loads to the BIW according to the standard, as shown in Fig 5. Through the interval sampling of the lower surface of the front rail, sill beam and rear rail of the BIW, the maximum vertical displacement D Lmax and D Rmax of the measuring points on the left and right sides are extracted, and the static bending stiffness of the BIW is calculated by equation (12).
where K b is the static bending stiffness; F is the resultant load, which is 6400 N; D Lmax is the maximum vertical displacement of the left measurement point; D Rmax is the maximum vertical displacement of the right measurement point. In order to verify the accuracy of the model, it is necessary to compare the simulation and test values of the vertical deformation of each measurement point of the BIW. According to the same constraint method as the simulation model, VOLUME 10, 2022  The simulated and tested values of bending stiffness of BIW are shown in Table 1. The error between the simulated and tested values is +5.96%, and the error is less than 10%, which meets the accuracy requirement. The simulation and test results of the measurement points are drawn as shown in Fig 7. It can be seen from the figure that the deformation of the left and right sides of the BIW is consistent, and the overall deformation is concave. The deformation in the middle of the sill beam is large, the closer it is to the fixed endpoint, the smaller the deformation is, and the overall deformation of the front and rear rail is small. The simulation displacement curve and the test displacement curve have good consistency in the change trend and the displacement of the measuring point, although there is a certain error, but they are within the acceptable range.

2) TORSIONAL STIFFNESS SIMULATION ANALYSIS AND TEST VERIFICATION
For the torsional stiffness test, the constraints and load standards imposed by the BIW are shown in Fig 8, and the test conditions are shown in Fig 9. The center points of the front and rear suspension of the BIW are constrained, and by measuring the displacements of the left and right loading points, the vertical displacement of the loading point is selected as the deformation displacement of the torsional stiffness, and the static torsional stiffness of the BIW is calculated by equation (13).
where K t is the torsional stiffness; M is the loaded torque; θ is the torsional angle; D L is the vertical displacement of the left loading point; D R is the vertical displacement of the right measuring point; L is the lateral distance between the left and right front suspension loading points. The simulated and tested values of BIW torsional stiffness are shown in Table 2, and the error between the simulated and tested values is +8.06%, which meets the accuracy requirement. The simulation and test results of the measurement points are drawn as shown in Fig 10. On the whole, the curve trend of the test and simulation measuring point displacement is in good consistency, the position near the rear suspension support is close to the constraint point, so the displacement of the nearby measuring points is small, almost 0. The position near the front suspension is closer to the loading point, so the displacement of its nearby measurement points is larger. Due to the approximate symmetry of the left and right sides of the    car body, the absolute values of the displacement of the measuring points on the left and right sides are also approximately equal.

3) LOW-ORDER MODAL SIMULATION ANALYSIS AND TEST VERIFICATION
During modal simulation, the simulation frequency bandwidth is set to 25-70 Hz. To make the test frequency bandwidth completely cover the simulation frequency bandwidth, the test frequency bandwidth is set to 0-70 Hz. In order to   simulate the modal characteristics in the free state, air springs with adjustable frequency stiffness were used to support the BIW during the test, and the test conditions are shown in Fig 11. Modal testing was performed using the LMS SCADAS and LMS Test. Lab equipment systems. The low-order modal shape of the BIW is shown in Fig 12, and the corresponding modal frequencies are shown in Table 3. The first-order torsional modal error of BIW is large, and the error rate is +9.50%, but the error is still within 10%, which meets the accuracy requirements. The implicit parametric model can replace the actual BIW for modal analysis research.

C. ESTABLISHMENT OF CAR CRASH MODEL 1) BENDING STIFFNESS SIMULATION ANALYSIS AND TEST VERIFICATION
Create a car crash model by assembling an implicit parametric model of the BIW with the powertrain, chassis assembly, and closure assembly.
The creation process mainly includes: Process 1: Finite element mesh discretization of the car model through shell units and solid units, and control of the mesh quality.   Process 3: According to the assembly relationship of the car, the connection relationship of each part of the car is simulated through the corresponding connection unit.
Process 4: Use the corresponding contact algorithm and control card to control the contact parameters, time step and output results of the crash model, and reasonably simulate the crash process of the car.

2) MATERIAL PROPERTY PARAMETERS
The energy absorption and deformation patterns during a car crash are directly influenced by the material model, and a correct material model is a key to ensuring the effectiveness of the crash system. The engineering stress-strain data of the material is obtained by static and dynamic tensile tests of standard samples. Taking DP350/600 steel as an example, the stress-strain curves under four different strain rates obtained by the high-speed hydraulic servo testing machine are shown in Fig 13.

3) MODELING OF FRONT CRASH AND SIDE IMPACT
According to the regulations of E-NCAP and C-NCAP and the actual situation, the car crash time is defined as 100 ms, and the front crash and side impact models of the car are established, as shown in Fig 14. In this article, the unit size is 10 mm, and the time step is empirically defined as 7×10 −7 s. The mass and center of gravity between the car crash model and the actual car are shown in Table 4. It can be seen from table 4 that the implicit parameter model of the BIW is very consistent with the actual car condition.

D. CAR CRASH MODEL VERIFICATION 1) FRONT CRASH TEST VERIFICATION
In order to facilitate the comparison of the performance indexes before and after the lightweight, the performance indexes such as car deformation pattern, acceleration and intrusion amount are directly extracted to evaluate the safety performance of the car front crash. The rationality of the front crash model is proved by analyzing the energy curve of the car front crash and the crash deformation pattern. The correctness of the front crash model is verified by comparing and analyzing the performance indicators such as the acceleration curve of the B-pillar, the maximum deformation of the upper and lower ends of the door, and the intrusion amount of the foot pedal.

a: ENERGY CHANGE RATE VERIFICATION
During the front crash of the car, the car hits the rigid wall at a speed of 56.5 km/h. During the crash, the kinetic energy decreases, the deformation internal energy increases, and the crash ends after 65 ms, and the energy change of the car remains balanced, as shown in Fig 15. The total energy in the entire crash process increased from 158300 J to 160000 J, with an increase rate of 1.1%; the hourglass energy increased by 5000 J, accounting for 3.2% of the total energy; the increase rates of the total energy and hourglass energy were both within 5%, indicating that the front crash model rationality.

b: DEFORMATION PATTERN VERIFICATION
The deformation pattern of the car in front crash is shown in Fig 16. By comparison, it can be seen that the front-end of the body undergoes large displacement and large area deformation, while the occupant compartment and its rear body maintain good integrity, ensuring the survival space  of the occupants. The deformation patterns of the test and simulation have a good consistency.

c: B-PILLAR ACCELERATION VERIFICATION
In the process of front crash, the car should not only use the yield deformation of the front-end safety structure to absorb a large amount of crash energy, but also reduce the acceleration of the occupant compartment as much as possible on the premise of ensuring the survival space, so as to reduce the impact force suffering in the crash process. Since the B-pillar is in the middle of the occupant compartment, its acceleration can replace the acceleration of the entire occupant compartment, and the overlapping area between the lower end of the B-pillar and the sill beam has high stiffness and small deformation, which can eliminate external environmental interference. Therefore, this article selects the overlapping area between the lower end of the B-pillar and the sill beam as the output position of the acceleration of the occupant compartment. The arrangement of the B-pillar acceleration sensor is shown in Fig 17. The comparison results of the B-pillar acceleration simulation and test are shown in Fig 18, and the overall trend is good consistency. Comparing and analyzing the maximum acceleration peaks, the simulated and test peaks of the left B-pillar acceleration are 35.68 g and 31.80 g, respectively, with an error rate of 12.2%; the simulated and test peaks of   the right B-pillar acceleration are 35.60 g and 31.08 g, respectively, with an error rate of 14.5%. Due to the complex factors such as large deformation, nonlinearity and the simplification of some plastic parts involved in the crash process, the error between the test and simulation results is large, but the trend is basically the same, and the peak error is within 15%, which can verify the correctness of the front crash model.

d: VERIFICATION OF THE MAXIMUM DEFORMATION OF THE DOOR
In the process of the front crash, due to the A-pillar rearward invasion leading to door extrusion deformation, the occupant compartment will be damaged, and even the door will be stuck, so the door deformation is an important indicator to measure the safety performance of the front crash. The small stiffness spring units are arranged horizontally at the upper and lower positions of the left and right front door frames, and the deformation of the door is evaluated by the deformation of the spring unit, as shown in Fig 19. During the test, the final deformation of the upper and lower ends of the doors was measured using a three-coordinate instrument, and the locations of the measurement points are shown in Fig 20. The results of the simulation and test are shown in Table 5, in which the error of the lower end of the right door VOLUME 10, 2022    deformation is larger, and the error rate is 12.1%, which is mainly due to the small amount of door deformation and the large error generated by the three-coordinate instrument in the measurement, but the accuracy of the model can still meet the requirements.

e: VERIFICATION OF FOOT PEDAL INTRUSION
The intrusion of the foot pedals of the car directly reflects the intrusion of the survival space of the front occupant's legs. If the intrusion is too large, the survival space of the front legs will be reduced, causing damage to the occupants. Through the simulation results, it is found that the foot pedal area has a large amount of intrusion due to the impact of the front rail, and the maximum intrusion displacement of the left and right foot pedals is 169.8 mm and 136.7 mm, respectively, as shown in Fig 21. In the front crash test, the intrusion of the pedal area is obtained by arranging regular discrete measurement points at the pedal. The distribution of both sides has 25 measurement points, and the left and right sides of the measurement point number are arranged in inverted S shape, the arrangement position as shown in Fig 22. According to the results of the front crash test, the maximum intrusion displacement of the left pedal appears at the measuring point D5 in the pedal area, which is 145.3 mm; the maximum intrusion displacement of the right pedal appears at the measuring point D23 in the pedal area, which is 129.9 mm. Due to the installation of control devices such as the accelerator pedal, clutch pedal and brake pedal on the left pedal area, the opening design appears in the left pedal area, so the intrusion displacement is larger than that of the right pedal.
The comparison results between simulation and test are shown in Table 6, in which the maximum intrusion displacement error rates of the left and right sides are 16.9% and 5.2%, respectively. Although the maximum intrusion displacement error on the left side is large, combined with the actual situation, the main reason is that some damping plates, frontend plastics and other components of the car are simplified in the simulation process, resulting in a small attenuation of crash energy and large invasion displacement. Moreover, the simulation measuring points are mainly selected in the continuous area, while the test measuring points are mainly selected in the discrete area, which leads to the difference between the simulation maximum displacement measuring point and the test maximum displacement measuring point, so that the results will be different. But overall, the error of simulation in the test is less than 20%, so the error is still within the acceptable range.

2) SIDE IMPACT TEST VERIFICATION
The key evaluation criteria of side impact are similar to those of front crash, and the rationality of side impact model is verified by analyzing the energy curve of a car side impact. By comparing the simulation results of car deformation pattern, B-pillar intrusion displacement and B-pillar intrusion acceleration with the physical test results, the correctness of the side impact model is verified.

a: ENERGY CHANGE RATE VERIFICATION
During the side impact of the car, the mobile deformation barrier (MDB) crashed with the side of the car vertically at a speed of 50 km/h. During the crash, the kinetic energy decreased, the internal energy increased, and the energy of the crash system remains balanced as a whole, as shown   in Fig 23. The total energy in the entire crash process increased from 93300 J to 95700 J, with an increase rate of 2.57%; the hourglass energy increased by 2760 J, accounting for 2.96% of the total energy; the increase rates of the total energy and hourglass energy were both within 5%, indicating that the side impact model rationality.   beam below the B-pillar is deformed by severe creases, and the two have a high degree of agreement by comparison with the test. The maximum intrusion displacement of the simulation and test is shown in Table 7. Since the intrusion amount will rebound after the test crash, the simulated intrusion amount is larger than the test value, the maximum error rate of the intrusion amount is 11.3%. The accuracy of the side impact model is reliable.

d: B-PILLAR INTRUSION ACCELERATION VERIFICATION
The simulation and test results of the acceleration at each measuring point of the B-pillar are shown in Fig 26. On the whole, the acceleration of each measuring point in the simulation and test has a good consistency, and the local wave crest positions also have a good consistency. The peak acceleration comparison between simulation and test is shown in Table 8. The results show that the acceleration error of the head measurement point is the largest, the error rate is 19.3%, but it is still within the error range of 20%, the error of chest VOLUME 10, 2022   and abdomen is about 15%, and the error of H-point is 3.3%. Overall, although there is a certain error, the acceleration curve has good consistency, so it is considered that the side impact model is still correct.

IV. MULTI-OBJECTIVE OPTIMIZATION OF BIW A. MULTI-OBJECTIVE OPTIMIZATION PROCESS
Based on the implicit parametric model of the BIW, the safety performance of the car in front crash and side impact is comprehensively considered, and the design variables are selected by the method of contribution and main effect analysis. Combined with the comprehensive performance of the BIW, define the objective function and constraints. Combined with the approximate model technology, MNSGA-II is used for multi-objective optimization of the BIW to obtain the Pareto solution set, and the comprehensive performance ranking of the Pareto solution set is performed by EGRA to select the solution with the best comprehensive performance. The lightweight multi-objective optimization process of BIW is shown in Fig 27.

B. DESIGN VARIABLE SCREENING
In this article, a hybrid material body design is proposed for the key components in the car front crash and side impact modules, while the BIW structure is improved through optimization methods such as material, shape and dimensional optimization. The design variables contain a variety of nonlinear factors, and it is difficult to obtain the relationship between the variables and the response using the traditional sensitivity method. Therefore, the method combining contribution and nonlinear main effects analysis is used to screen the initial variables. Thirty-four variables such as material, shape and dimensional variables of important safety parts and key components of the BIW were selected as initial design variables, as shown in Table 9. Combined with the static and dynamic performance of the BIW and the performance indicators of front crash and side impact, 82 sample points were collected by the optimal Latin hypercube method. The response values of each sample point are obtained by simulation, and the contribution of each initial design variable to the performance response and the nonlinear main effect values are calculated. The contribution of each initial design variable to the overall performance response of the car body is shown in Table 10. Among them, red represents positive contribution, and blue represents negative contribution; bs, ts, bf, and tf are the body bending, torsional stiffness, and the first-order bending and first-order torsional modal frequencies, respectively; mass is the car body weight. During the front crash, al and ar are the acceleration of the left and right side B-pillars, respectively, il and ir are the intrusions of the left and right side front panels, and dl and dr are the deformation of the middle of the left and right side doors, respectively. During the side impact, ib and ab are the intrusion displacement and intrusion acceleration of the B-pillar, respectively.
When defining the initial value of the design variable, by referring to the main effect relationship between the design variable and the performance response, the initial value can be better selected to avoid the reduction in the efficiency of the optimization search caused by the initial values deviating too much or too little from the optimal values. The contribution of each initial design variable to the acceleration of the left side B-pillar during the front crash process and the main effect analysis are used as examples to illustrate that the combination of the two can achieve the screening of the initial design variables and the selection of initial values, as shown in Fig 28 and Fig 29. As can be seen from Fig 28, the design variables that have the most significant impact on the front crash safety performance are the material variables DV8 and DV15, with the steel and aluminum material optimization (DV8) for the front rail having the greatest impact with a contribution of 14.1%. The main effect relationship curves of each design variable with the acceleration of the left B-pillar are shown in Fig 29, and a positive slope of the curve indicates a positive effect, while the opposite is a negative effect. There is a strong nonlinear relationship between the acceleration of the B-pillar and the design variables during the front crash. The main effect curve of the material variable DV8 has a parabolic relationship with the acceleration of the B-pillar. The upper and lower limits of the material (aluminum alloy material) have small acceleration values. The initial value of the material (raw steel material) corresponds to a larger acceleration value. Therefore, for front crash safety performance, the material variable DV8 is the design variable with the greatest impact. By studying the contribution and main effect relationship between each initial design variable and each performance response, considering the car performance comprehensively, the optimized design variables and their initial values and variation ranges as shown in Table 11 are screened out, in which the material type parameters are discrete design variables and the structural dimension parameters are continuous design variables.

C. OBJECTIVE FUNCTION
Since the lightweight performance of the BIW and the safety performance of the car are in conflict with each other, it is necessary to achieve a balance between the two in the multiobjective optimization design. Considering that the BIW lightweight coefficient QLX combines the lightweight and car safety performance, QLX is defined as the optimization target. For the front crash safety performance, considering the high initial speed during the crash, it is necessary to reduce the crash acceleration to prevent the occupants from secondary injuries, so the average value Fac of the acceleration of the left  and right sides B-pillars is taken as the optimization target. For the side impact safety performance, since the B-pillar is directly impacted by the moving barrier during the side impact, and the middle of the B-pillar is close to the weak organ of the occupant's chest, so the B-pillar chest intrusion displacement ID is taken as the optimization target. The objective function of BIW multi-objective optimization can be defined as: where QLX is the BIW lightweight coefficient; Fac is the average value of the maximum acceleration of the B-pillars on both sides during the front crash; ID is the maximum intrusion displacement of the B-pillar chest during the side impact.

D. CONSTRAINT CONDITION
In order to better ensure the comprehensive performance of the BIW during the optimization process, in addition to the objective function, the BIW bending stiffness BS, the first-order bending frequency BF, the first-order torsion frequency TF, the intrusion displacement IDF of the front panel during the front crash, the intrusion displacement IDD L of the pedal at the driver, the intrusion displacement IDD R of the pedal at the front occupant, the average deformation of the upper part of the doors on both sides DU , the average deformation of the lower part of the doors on both sides DL, the intrusion acceleration Acc of the B-pillar chest during the side impact, the intrusion displacement ID b1 of the B-pillar head, and the intrusion displacement ID b2 of the B-pillar abdomen are used as performance constraints in the optimization process. The performance constraints of the BIW are defined as: where BS 0 , BF 0 , TF 0 , IDF 0 , IDD L0 , IDD R0 , DU 0 , DL 0 , Acc 0 , ID b10 , ID b20 are the performance metrics before optimization, respectively; BS i , BF i , TF i , IDF i , IDD Li , IDD Ri , DU i , DL i , Acc i , ID b1i , ID b2i are the performance metrics after optimization, respectively.

E. MATHEMATICAL MODEL
Considering the static and dynamic performance and crash safety performance of the BIW, the mathematical model of BIW lightweight multi-objective optimization with material, shape and dimensional parameters as design variables is as follows:

A. APPROXIMATE MODEL FIT
Using the optimal Latin hypercube method, 105 sample points were randomly taken, and the Kriging surrogate model and the RBF surrogate model were fitted, respectively. In order to verify the accuracy of the approximate model, 20 randomly selected sampling points in the design sample were used to verify the accuracy of the approximate model by cross-validation method, and the accuracy verification results of the approximate model for some response indicators are shown in Fig 30. The accuracy of the approximate model is evaluated using the coefficient of determination R 2 , the root mean square error RMSE, and the maximum absolute relative error MAPE. If the R 2 value is closer to 1 and the RMSE value is closer to 0, the overall prediction accuracy of the approximate model is higher, and when the MAPE value is smaller, the local prediction accuracy of the approximate model is higher. The specific index values of the accuracy test results of each approximate model are shown in Table 12.
By comparing the Kriging surrogate model and the RBF surrogate model, it can be seen that although the prediction accuracy of the Kriging surrogate model in terms of the BIW lightweight coefficient is higher than that of the RBF surrogate model, but the accuracy of the RBF surrogate model in front crash and side impact on nonlinear large deformation responses is higher than that of the RBF surrogate model. The coefficient of determination R 2 of the RBF surrogate model is above 0.85 and has a small root mean square error RMSE and maximum absolute value relative error MAPE, which better balances the global and local prediction capabilities of the approximate model. Therefore, after comprehensive consideration in this article, the RBF surrogate model is selected to represent the functional relationship between each design variable and the target response.

B. PARETO FRONTIER
The closed optimization environment created with Isight as the basic platform, the RBF surrogate model as the optimization object, and the MNSGA-II as the driving strategy is shown in Fig 31. The BIW surrogate model was solved by using MNSGA-II, setting the population size and the maximum number of iteration steps to 160 and 300, respectively, and the crossover probability and mutation probability to 0.86 and 0.1, respectively, to obtain the Pareto frontier consisting of 258 non-dominated optimal solutions after 48,675 iterations, as shown in Fig 32.

C. PARETO FRONTIER SOLUTION SET RANKING
It can be seen from the Pareto frontier that it is difficult to achieve the optimal performance of QLX, Fac and ID at the same time, so different compromise solutions need to be selected according to different optimization requirements. Since the three objective functions of QLX, Fac and ID have the characteristic of 'the lower, the better', the data of the objective functions corresponding to the 258 non-dominated optimal solutions in the Pareto frontier are normalized, and then the GRC of each non-dominated optimal solution are calculated. Using the entropy weight method, the weights of the three objective functions of QLX, Fac, and ID are 0.27, 0.38 and 0.35, respectively. The GRG of the Pareto frontier obtained using the EGRA method is shown in Fig 33, where the 146th non-dominated solution has the largest GRG   of 0.8275, so the 146th scheme A is considered the design scheme with the best overall performance.
In order to verify the effectiveness and feasibility of the EGRA method in the process of selecting the optimal compromise optimization scheme, this article uses the TOPSIS method to obtain the optimal compromise scheme B, and obtains the recommended scheme C with the smallest lightweight coefficient in the Isight platform. When using the TOPSIS method to rank the Pareto frontier, in order to ensure the consistency of each condition, each objective function still uses the weight obtained by the entropy weight method. Fig 34 shows the similarity of the Pareto frontier obtained by the TOPSIS method, in which the 109th scheme B has the largest similarity of 0.7915, so the scheme B is the optimal design scheme obtained by the TOPSIS method.
The performance comparison between optimization schemes A, B, and C is shown in Table 13. On the whole, the optimization schemes A, B, and C can all meet the performance requirements of the BIW. Among them, the optimization scheme A can greatly improve the performance of the BIW, and the performance improvement rate is mostly between 5% and 10%, and the comprehensive performance is better. Although the optimization scheme B is better than the scheme A in terms of the lightweight coefficient QLX and the B-pillar chest acceleration Acc in the side impact, the improvement in the crash safety performance is small and the difference between the improvement rates of various performances is more obvious. The optimization scheme C focuses on the lightweight coefficient QLX of the BIW, and achieves an improvement rate of 15.6%, but the improvement rate of other performances is small. Therefore, the optimization scheme A has better comprehensive performance, and the improvement rate among various performances is relatively balanced. Compared with the TOPSIS method, the EGRA method is simple and flexible in operation and can provide effective guidance in the process of selecting the optimal compromise scheme.

D. COMPARISON OF BIW PERFORMANCE BEFORE AND AFTER OPTIMIZATION
Although the BIW optimization scheme with better comprehensive performance was obtained based on MNSGA-II, approximate model method and EGRA method, in order to ensure the reliability of the optimized design, the design variables in the optimization scheme A were assigned to the finite element model (FEM) to obtain the optimized results.

1) BENDING AND TORSIONAL STIFFNESS COMPARISON
The comparison of bending and torsional stiffness before and after optimization is shown in Table 14. The maximum error between the FEM model and the approximate model is 2.6%, so the approximate model has a high accuracy guarantee. Through the multi-objective optimization design, the bending and torsional stiffness of the BIW are improved by 6.4% and 4.7%, respectively, which improves the basic load-bearing performance of the BIW.

2) COMPARISON OF FIRST-ORDER BENDING AND TORSIONAL MODAL
The first-order bending and torsional modal of the BIW before and after optimization are shown in Table 15. The  maximum error rate between the FEM model and the approximate model is only 3.5%, which verifies the accuracy of the approximate model. Through the optimized design, the first-order bending and first-order torsional modal of the BIW are improved by 5.1% and 6.1%, respectively, improving its dynamic performance. Fig 35 shows the intrusion displacement of the front panel before and after optimization, the maximum intrusion displacement of the front panel before and after optimization occurs at the left and right knees, where the left knee intrusion is the most obvious, the optimized front panel in the overall displacement trend with the optimization before maintaining a high degree of consistency. The comparison of the intrusion displacement of the front panel before and after optimization is shown in Table 16, where the error rate between the approximation model and the FEM model is 1.4%, indicating that the approximation model of the front panel intrusion displacement has high reliability. Through the optimized design,  the maximum intrusion displacement of the front panel is reduced from 190.5 mm to 182.7 mm, with an improvement rate of 4.1%, which greatly improves the occupant's survival space.

b: COMPARISON OF DOOR DEFORMATION
The comparison of the door deformation before and after optimization is shown in Table 17. The maximum error between the approximate and FEM model is 6.1%, the reliability of the approximate model is verified. The improvement rates of the upper and lower deformations of the left side doors were 7.5% and 12.3%, respectively, and the improvement rates of the upper and lower deformations of the right side doors were 8.9% and 9.6%, respectively, which significantly improved the front crash safety performance.

c: B-PILLAR ACCELERATION COMPARISON
The acceleration of the B-pillar on both sides before and after optimization is shown in Fig 36. The change trend of acceleration before and after optimization maintains a good consistency, and the peak acceleration appears at about 60ms, which reflects the rationality of the car crash process. Comparing the data results in Table 18, the maximum error rate between the B-pillar acceleration approximate model and the FEM model is 6.1%, and the approximate model meets the accuracy requirements. After optimization, the accelerations of the B-pillars on the left and right sides are 31.3 g and 31.9 g, respectively, and the improvement rates are 12.3% and 10.4%, respectively, compared with those before optimization, which effectively reduces the crash impact force on the occupants and improves the crash safety performance of the car.

4) SIDE IMPACT SAFETY PERFORMANCE COMPARISON a: B-PILLAR ACCELERATION COMPARISON
The acceleration of each position of the B-pillar before and after optimization is shown in Fig 37. The optimized acceleration can better match the acceleration curve before VOLUME 10, 2022 optimization, while the peak acceleration after optimization is delayed compared with that before optimization, and the change trend is relatively gentle, which greatly improves the crash safety performance of the BIW and reduces the impact force on the occupants. The changes of the peak acceleration of the B-pillar before and after optimization are shown in Table 19. The maximum error rate between the approximate and FEM model is 6.2%, and the approximate model can meet the accuracy requirements. The accelerations corresponding to the optimized head, chest, abdomen and H-point positions are 23.3 g, 39.2 g, 36.7 g and 58.2 g, respectively, and the improvement rates are 5.7%, 6.2%, 5.9% and 5.7%, respectively.

b: COMPARISON OF THE AMOUNT OF INTRUSION AT EACH POSITION OF THE B-PILLAR
The maximum intrusion amount at each position of the B-pillar before and after optimization is shown in Table 20. By comparison, it can be seen that the maximum error rate between the approximate and FEM model is 5.8%, and the approximate model can meet the accuracy requirements. The maximum intrusion amount corresponding to the optimized head, chest, abdomen and H-point positions was 109.5 mm, 179.1 mm, 220.4 mm and 199.7 mm, respectively, and the improvement rates reached 6.6%, 8.7%, 4.3% and 4.4%, respectively. Fig 38 illustrates the optimization effect of scheme A by comparing the deformation of the side impact of the car before and after optimization. It is clearly observed from the figure that the intrusion displacement of the middle of the B-pillar is smaller than that before optimization, and the deformation pattern reflects the design concept that the B-pillar effectively protects the occupant's head and chest through the middle and upper parts, while the lower part absorbs the crash energy better and reduces the secondary injury of the occupant.

VI. CONCLUSION
1) In order to enhance the variable space of the parametric model, improve the optimization efficiency of the BIW and enhance the lightweight effect, this article uses the SFE-CONCEPT software to create an implicit parametric model of the BIW. By comparing the physical test results with the simulation results, the accuracy of the BIW implicit parametric model is verified in terms of bending and torsional stiffness of the BIW, low-order modal, front crash and side impact safety to ensure the reliability of the subsequent optimization design based on the car model. It realizes the parallel and integration of CAD and CAE, as well as CAE analysis driving CAD design. 2) Taking the BIW as the optimization object, material parameters, shape parameters and dimensional parameters were introduced as design variables, and a method combining contribution analysis and nonlinear main effect analysis was proposed to screen out 26 important design variables and establish a mathematical model for BIW optimization. To address the accuracy of the approximate model, the Kriging surrogate model and the RBF surrogate model are established, and it is found that the RBF surrogate model can better reflect the relationship between nonlinear crash performance and optimization variables through comparison. In order to further improve the population diversity of the multi-objective optimization problem, a fixed threshold elimination strategy was used to replace the crowding distance comparison method in MNSGA-II, which is combined with the RBF approximation model to perform multi-objective optimization of the BIW, and 258 non-dominated schemes are obtained. 3) A method is proposed to use the EGRA method to sort the comprehensive performance of the Pareto solution set of the multi-objective optimization problem of the BIW, and the compromise scheme A with the best comprehensive performance is obtained. By comparing scheme A with the optimal compromise scheme B obtained by the TOPSIS method and the recommended scheme C with the smallest lightweight coefficient obtained from the Isight platform, it is found that the improvement of various BIW properties in scheme A is more obvious and the improvement rate is more balanced, which verifies the effectiveness and feasibility of this optimization scheme and ranking method. It avoids the blindness to the optimal solution selection, quantifies the comprehensive performance of each scheme, and establishes an objective evaluation method of the multi-objective optimal design results. 4) By comparing the comprehensive performance response of the BIW before and after optimization, the BIW lightweight coefficient QLX, the average value Fac of the maximum acceleration of the B-pillar on both sides during a front crash, and the maximum intrusion displacement ID of the B-pillar chest during a side impact, the improvement rates reached 11.5%, 6.5% and 6.8%, respectively, and the improvement rates of the rest of the performance response were above 3.3%. Therefore, the optimization design method proposed in this paper can effectively improve the car safety performance of the BIW while lightweight the BIW.