Multi-objective optimization design for a battery pack of electric vehicle with surrogate models

2.1. Optimal Latin hypercube sampling for design of experiment

In engineering design, direct coupling of optimization algorithm with simulation model may not be efficient since iterative calculation and sensitivity analysis usually require enormously computational efforts. As an alternative, the metamodels or surrogate modeling techniques have been proven as an effective substitute of costly simulation analysis for optimization [14, 18-25]. In practice, surrogate modeling technique starts with the sampling data at some training points. Design of experiment (DoE) is an approach to address how to select training points effectively and appropriately. In this paper, the optimal Latin Hypercube sampling (OLHS) approach is employed to generate initial training points [28]. In order to find a set of training points with the highest uniformity, $L_P$ discrepancy is used as the criterion to gauge the uniformity in this study. Presume that there are $N$ factors of interest over a standard domain $C^N$. The goal is to choose a set of training points $x_S=\left\{x_1,x_2,\dots ,x_{n_S}\right\}$ such that these points are uniformly scattered on $C^N$. To seek a set of $x_S$ that maximizes the uniformity over all possible $n_S$ points on $C^N$, the centered $L_2$ discrepancy criterion is used in our study, given as [29]:

2.2. Metamodels in multi-objective sequential optimization

In the optimization process, the surrogate models, specifically, polynomial response surface (PRS) [15], Kriging (KRG) [16], radial basis function (RBF) [17] and etc, have been widely adopted as effective tools to approximate structural performances, and can largely reduce the number of costly finite element analysis (FEA) runs to improve the optimization efficiency and feasibility. Surrogate model technique starts with training data at sampling points as above-mentioned, and then constructs functional relationship between design variables and responses so as to predict the values at any points in the design space.

PRS model is one of the simplest and most popular surrogate models. As an effective alternative to FEA, PRS model has been widely adopted in design optimization [14, 23, 30]. Mathematically PRS metamodels from the first to the fourth order can be, respectively, written as:

where $x_i$ ($i=$ 1, 2,…, $t$) denotes the design variables determining response $y$; based on the least square method, estimated regression coefficients, namely, $a$, $b_i$, $c_{ii}$, $e_{ii}$, $f_{ii}$ and $d_{ij}$ can be obtained. $x_i{x}_j$ is the cross term that represents two-parameter interaction and the higher-order terms stand for the nonlinearity characteristic. In general, the order of PRS is no higher than 4.

The KRG model was originally developed for mining and geostatistical applications involving spatially and temporally correlated data. The stochastic process function $y\left(x\right)$ is used to describe the deterministic response of a system, which contains a regression model and a stochastic error [16]:

where ${\beta =[{\beta }_1,{\beta }_2,\dots ,{\beta }_p]}^T\text{,}{f\left(x\right)=\left[f_1\left(x\right),f_2\left(x\right),\dots ,f_p\left(x\right)\right]}^T\text{,}$ are the column vectors of regression parameters and basis functions, respectively; p represents the number of basis functions; $z\left(x\right)$ denotes a stochastic parameter with zero mean and variance as ${\sigma }^2$. The covariance matrix of $z\left(x\right)$ is shown as:

where $R$ is a correlation matrix defined by Gaussian correlation function $R(x_i,x_j)$ as follows:

here $n_{dv}$ is the number of design variables, ${\theta }_k$ is the unknown correlation parameter used to fit the model, and $x_i^k$ and $x_j^k$ are the $k$th component of training points $x_i$ and $x_j$, respectively.

Then, $\widehat{y}\left(x\right)$ is employed as predicted estimate of response $y\left(x\right)$, expressed as:

where ${y_s=\left[y\right(x_1),y(x_2),\dots ,y(x_{n_S}\left)\right]}^T$ is the response vector of the ${\mathrm{n}}_{\mathrm{S}}$ training points $x_s=\{x_1,x_2,\dots ,x_{n_S}\}\text{,}$ and ${F=\left[f\right(x_1),f(x_2),\dots ,f(x_{n_S}\left)\right]}^T$ is a $n_S\times p$ matrix. ${r^T\left(x\right)=[R\left(x,x_1\right),R\left(x,x_2\right),\dots ,R\left(x,x_{n_S}\right)]}^T$ is a correction vector that implies how close between training points and untrained points. $\widehat{\beta }$ is the general least square estimator formulated as follows:

The estimate to the variance of training data from the global model is:

For calculating ${\theta }_k$ in Eq. (5), the maximum likelihood estimates can be used by solving the following maximization problem over the interval ${\theta }_k>0$, as:

where both ${\widehat{\sigma }}^2$ and $\left|R\right|$ are the functions of ${\theta }_k$.

RBF model exhibits superior prediction ability and ensures good accuracy in highly nonlinear responses problem [17]. The RBF model consists of the linear combination of basis function, and the basis function in RBF is given as:

where $\xi =‖x-{\mu }^i‖$, is an independent variable representing the Euclidean distance between the design point $x$ and the observed input ${\mu }^i$. The RBF model is given as:

The accuracy of surrogate models should be assessed after they are established by adding extraneous confirmation sample points. Three metrics adopted herein are R-square, relative average absolute error (RAAE), and relative maximum absolute error (RMAE) [31], respectively, as:

where $y_i$ and ${\widehat{y}}_i$ indicate the exact function value and the corresponding surrogate value at confirmation point $i$, respectively. $\stackrel{-}{y}$ is the mean of $y_i$, $q$ is the number of the confirmation points. Generally, a larger value of R-square and a smaller value of RAAE inferring a higher accuracy for overall performance in the design space are preferred. On the other hand, a larger RMAE value denotes less accuracy in one region of the design space even though a very good global measurement can be given by R-square and RAAE. In optimization applications, the first two metrics are usually paid more attention for the global behavior concerning.

2.3. Multi-objective particle swarm optimization (MOPSO) algorithm

Particle swarm optimization (PSO), is a relatively new heuristic and derivative-free global optimum algorithm inspired by the choreography of a bird flock [32]. It starts with a set of initial random solutions to obtain an optimal solution by repeated iterations. The multi-objective particle swarm optimization (MOPSO) [27] is an extended version to PSO, which has gained considerable attention for its capacity of fast convergence and obtaining well-distributed Pareto frontier. To maintain the distribution uniformity of solutions, MOPSO uses an external archive to save the non-dominated solutions during the search procedure and clips the non-dominated solutions by a crowded degree algorithm. Finally, MOPSO selects the elite individual from the sparse regions in the external archive as the global optimal solution.

3.1. Finite element modeling

As a critical and independent part assembled to carry and protect Lithium-ion power battery set, especially considering of the complicated driving condition of electric vehicle, the battery pack requires high mechanics performance to carry out its functions. Fig. 1(a) and (b) display the global geometric model and inner structure model of an air-cooling battery pack using the commercial software Catia.

Based on the FEA code Hypermesh, this battery pack, composed of panels and beams, is modeled mainly with PSHELL shell elements with their material property of Q235 steel assigned linear elastic isotropic material model MAT1. For connectivity, REB2 and CWELD elements are used to model weld spots, as well as 1D PBAR element for the representative of battery connecting bars. For lack of detailed information about battery material and mechanics property, like other studies [4, 6-8], some simplified treatment measurements are employed in this research to assure the mechanical effect of battery on the whole pack structure. Taking account of cuboid geometry shape feature, battery sets are substituted by PSHELL elements with appropriate thickness and density. A much lower elastic modulus of the battery compared to steel is adopted here for neglecting its own mechanical response in FEA process. The entire FEA model has 117837 elements with the total mass of 0.3211 ton, which is shown in Fig. 2.

Fig. 1Global geometric model and inner structure model of the battery pack

a) Global geometric model of the battery pack

b) Inner structure model of the battery pack

Fig. 2FEA model of the battery pack and its inner battery sets

a) FEA model of the battery pack

b) FEA model of the inner battery sets

Similar to [4, 6-8], as an important load-condition in structure analysis and optimization for battery pack, restrained modal analysis is studied in this work. As the boundary condition shown in Fig. 3, 6 freedoms of each connecting position between battery pack bottom plate and the vehicle body are constrained.

Fig. 3Boundary condition of battery pack FEA model

A higher basic frequency (1st order restrained modal frequency) of battery pack can make it avoid the resonance with the vehicle [11], which is concerned most in this work. The first 20 order restrained frequencies and corresponding positions of vibration are illustrated in Table 1. Fig. 4 demonstrates the vibration mode of the first order frequency.

Besides, limited driving conditions of vehicles, such as frontal-collision and side-collision happened in a very short time, would cause super large acceleration acting on all the components [11]. For the sake of battery safety, the battery pack should not happen to crack under such transient shock, which means stress in the battery pack should not exceed the material tensile limit. Consequently, in this research, according to the regulations described in Chinese Automobile industry standard – General Requirement of Traction Battery Enclosure for Electric Vehicles (QC/T 989-2014), transient shock response of the battery pack under 40 g acceleration acted along its $x$-axis and $y$-axis in 6 ms, is concerned [12]. The loading process of acceleration excitation is shown in Fig. 5.

Fig. 4The vibration mode of first order frequency of the battery pack

a) Main view of the vibration mode of first order frequency

b) Inner view of the vibration mode of first order frequency

Through the FEA results displayed in Fig. 6, it is easy to find that the maximum von Mises stresses in battery pack $x$-axis and $y$-axis transient shock cases are both above 380 MPa, the tensile limit of steel. Therefore, security risk exists in this battery pack, and a redesign is necessary.

Fig. 5The loading process of acceleration excitation in 6 ms

Fig. 6FEA results of x-axis and y-axis transient shock on the battery pack

a) FEA results of $x$-axis transient shock

b) FEA results of $y$-axis transient shock

3.2. Definition of optimization problem

An excellent battery pack structure design should meet all of the mechanical performance requirements, in conjunction with lighter weight for energy-saving and mileage elongation consideration. Of all the mechanical properties, its first order restrained natural frequency ought to be as high as possible to avoid the resonance with the vehicle. Therefore, the basic frequency is maximized and at the same time the structural mass is minimized, whilst maintaining certain levels of other transient shock response indices along its $x$-axis and $y$-axis as mentioned. In this study, thicknesses of seven panels and one beam set in Fig. 7 are selected as the design variables whose ranges are all from 0.5 mm to 1.5 mm. The detailed information of these design variables are depicted in Table 2.

Fig. 7Illustration of design variables

The multi-objective size optimization problem for battery pack is thus specifically formulated as follows:

here, $F_m\left(x\right)$ and $F_f\left(x\right)$ are the total mass and first order restrained natural frequency of the battery pack, respectively; $S_{\left(x\right)}$ and $S_{\left(y\right)}$ are the global maximum von Mises stress in battery pack $x$-axis and $y$-axis transient shock cases.

Fig. 8Flowchart of the proposed MOPSO procedure

It takes a long time and a lot of hurdles to obtain FEA results used in the process of searching optimal mathematical solution especially for evolutionary algorithm often employed in multi-objective and multi-constraint optimization problem that involves large number of iterations. As an alternative, surrogate model is an effective way to formulate the relationship between the design variables and responses mathematically, which will also reduce the cost of computation. In general, for different problems, the precision of different kinds of surrogate models are not the same. In this investigation, three different surrogate models, PRS, KRG, and RBF are considered. Among these surrogates, the most accurate surrogate will be picked up to perform the multi-objective optimization design. To construct high accuracy surrogate, of many different DOE methods available, OLHS approach is selected to generate the training points in the design space here. In order to capture the global maximum dynamic stress in transient shock with nonlinear responses, 100 training points are generated first in this study [22, 24, 25, 33]. It should be noted that if the accuracy of metamodels based on the initial OLHS is unacceptable, the optimization procedure should be returned to the DOE step and more sampling points will be generated in addition to the existing sample points in which a max-min distance criterion is used for inserting the new points to the existing sample space [34]. Given the existing sample set $X_p$ ($s$ samples), a new sample set $X_c$ ($t$ samples) will be selected based on the max-min distance approach to maximize the minimum distance between any two sample points in the total sample set $X_A=X_p{\cup X}_c$, i.e.:

where $X_{Ci}\in X_C$ ($i=$ 1,…, $s$), $X_{Aj}\in X_A$ ($j=$1,…, $s+t$).

To obtain the optimal multi-objective design of the battery pack, MOPSO algorithm, which has relatively fast convergence property and well-distributed Pareto frontier, is employed in this paper. More sampling points will also be added based on the max-min distance criterion mentioned above if the Pareto Optimal Frontier needs to be enhanced. For clarification, the proposed optimization procedure is summarized in the flowchart as Fig. 8.

4.1. Selection and analysis of surrogate models

For the sake of choosing the most accurate and appropriate metamodels for multidiscipline problem in this study, the abovementioned three individual surrogates, namely, PRS (from 1st order polynomial to 4th order polynomial), RBF and KRG are constructed and used to run the optimization. Since the mass of the battery pack follows a linear relationship to the panel thicknesses, the 1st order PRS is adopted to model the weight. To assess the accuracy of these different surrogate schemes modeling the dynamic and transient shock responses of battery pack, 30 new validation points are also generated using the OLHS approach here. Three numerical estimators, namely R-square (R2), relative average absolute error (RAAE) and relative maximum absolute error (RMAE) [31] are used to measure the accuracy of these metamodels, as given in Eqs. (12). The computational results are summarized in Table. 3.

For the prediction of restrained basic frequency in this study, similar to [14], 2nd order RSM is the most accurate and appropriate surrogate model with the R2 value of 0.9902 and RAAE value of 0.0805, which indicates a good performance for overall design space. As for transient acceleration shock analysis, which has relatively high nonlinearity degrees, 3rd order RSM and RBF are the most suitable metamodels for maximum global stress prediction of battery pack under $x$-axis and $y$-axis transient shock, respectively.

4.2. Results of MOPSO

Fig. 9 depicts the final Pareto optimal frontier (POF) of the multi-objective optimization for the battery pack, where the settings of MOPSO are list in Table 4.

POF converges adequately after 100 generations. A clear convexity of the curve can be observed, which provides the alternative design space for designers. It is noted that the exploration of the Pareto sets is crucial to fully understand the solution space for the multi-objective optimization problems. The POF provides the designer with a range of optimal solutions for their further decision-making. In all these Pareto plots it is shown that the mass and restrained basic frequency strongly compete with each other, with the range of mass variation from 0.3235-0.3501 ton, alongside with the range of frequency variation from 36.09-31.88 Hz. Obviously, owing to the effect of the restraint of maximum stress response under acceleration shock on the design domain, after optimization, all the values of mass and restrained basic frequency in the Pareto optimal solutions are over and below the initial design, respectively. Specifically, if the designers pay more attention to the weight reduction, the solutions at the bottom right corner should be selected. While if the designers emphasize on restrained basic frequency, the solutions at the top left corner will be considered.

Fig. 9Pareto optimal frontier (POF) and knee point of MOPSO

Thickness variation range of these eight design variables in Pareto optimal solution set is illustrated from Fig. 10(a) to Fig. 10(h). Apparently, from the eight variables, the thicknesses of the floor are distributed relatively uniform and almost occupied the whole variable space. That is because of the floor has a large surface area compared to other parts, which makes the total battery pack mass very sensitive to this floor thickness change. With the function of battery fixing for fear of its up-and-down movement in the pack, apron set is thinner, and has a size change basically between 0.5 to 1mm, since in this study we do not consider the maximum stress response under $z$-axis acceleration transient shock. Air intake plate, air outlet plate and air outlet chamber are relatively thicker with the size no less than 1.0 mm in order to meet the stress requirement in $x$-axis acceleration shock. For the sake of improving the 1st order restrained natural frequency, cover thickness should not less than 1.0 mm. The eighth design variable, namely, the thickness of the beam set, almost keeps 1.5 mm which is the upper bound of design variables in all Pareto optimal solutions for considering the reinforcement and connection functions of these beams.

Fig. 10Thickness range of eight design variables in Pareto sets after optimization

a) Thickness range of floor in Pareto sets

b) Thickness range of apron in Pareto sets

c) Thickness range of air intake plate in Pareto sets

d) Thickness range of cover in Pareto sets

e) Thickness range of air outlet plate in Pareto sets

f) Thickness range of air outlet chamber in Pareto sets

g) Thickness range of side wall in Pareto sets

h) Thickness range of beam set in Pareto sets

To select the most apposite and the best trade-off optimization scheme, Knee Point, as a normalization to determine the accuracy improvement of Minimum Distance, mathematically given as below, allows us determining a most satisfactory solution from Pareto set [33]:

where $n$ is the number of the objective components, $f_{c\tau }$ is the $\tau$th objective value in the $c$th Pareto solution, $D$ is the distance from knee point to an “utopia point”.

Compared with the initial design, it seems that the MOPSO knee point in this study deteriorate both of the two objectives (increase the mass and reduce the restrained basic frequency), which seems that our optimization study is meaningless. Actually, as above-mentioned, considering the restraint of maximum stress response under $x$-axis and $y$-axis acceleration shock, it is inevitably to increase the total mass of the battery pack to increase its strength, but simultaneously cause the basic frequency reduction. Our optimization study is just based on this situation. As shown in Fig. 9, frequency variation range in Pareto optimal solution is not much, just between 31.88-36.09 Hz, and the lower bound 31.88 Hz is over 30 Hz which is satisfactory for automotive devices [11]. Taking account of the lightest weight optimization scheme corresponding to restrained basic frequency is 31.88 Hz, this scheme is also meaningful and useful which should be selected as an alternative alongside with the Knee Point. The validation results of the two alternative optimum schemes using FEA, and their comparison with baseline are summarized in Table 5. There is small difference between the metamodel and FEA results as shown in Table 5. Both of the two schemes satisfy all these two stress response constraints presented. Nevertheless, it is noted that the Max. stress response under y-axis acceleration shock of Knee Point scheme has been pushed to the boundary of the constraint after optimization, which could lead to an unreliable solution if any perturbation of design variables presents.