Robust Design of Suspension System with Polynomial Chaos Expansion and Machine Learning

During the early development of a new vehicle project, the uncertainty of parameters should be taken into consideration because the design may be perturbed due to real components’ complexity and manufacturing tolerances. Thus, the numerical validation of critical suspension specifications, such as durability and ride comfort should be carried out with random factors. In this article a multi-objective optimization methodology is proposed which involves the specification’s robustness as one of the optimization objectives. To predict the output variation from a given set of uncertain-but-bounded parameters proposed by optimization iterations, an adaptive chaos polynomial expansion (PCE) is applied to combine a local design of experiments with global response surfaces. Furthermore, in order to reduce the additional tests required for PCE construction, a machine learning algorithm based on inter-design correlation matrix firstly classifies the current design points through data mining and clustering. Then it learns how to predict the robustness of future optimized solutions with no extra simulations. At the end of the optimization, a Pareto front between specifications and their robustness can be obtained which represents the best compromises among objectives. The optimum set on the front is classified and can serve as a reference for future design. An example of a quarter car model has been tested for which the target is to optimize the global durability based on real road excitations. The statistical distribution of the parameters such as the trajectories and speeds is also taken into account. The result shows the natural incompatibility between the durability of the chassis and the robustness of this durability. Here the term robustness does not mean “strength”, but means that the performance is less sensitive to perturbations. In addition, a stochastic sampling verifies the good robustness prediction of PCE method and machine learning, based on a greatly reduced number of tests. This example demonstrates the effectiveness of the approach, in particular its ability to save computational costs for full vehicle simulation.


Introduction
The robustness of vehicle specifications is given more and more attention in Renault because original designs during the development can be perturbed by many uncertain sources: actual road charges, manufacturing tolerances, aging of materials, etc. As a result, sometimes the validation of important specifications such as chassis' durability may take a lot of time and resources to ensure the design is robust enough and can be satisfied all through the vehicle's life cycle [1].
In this paper, the term of robustness is defined as system's ability of tolerating outside perturbations. The robustness is an opposite notion of sensitivity where the least variation of vehicle performance is searched within the random input parameters. Instead of analyzing the impact of each parameter on the final output, the robust optimization focuses on minimizing the overall variations while the statistic characters of input are predefined. Meanwhile, these two notions can easily be transformed after the result of a design of experiments is obtained.
There are several numerical methods to calculate the robustness of one given set of parameters under perturbation.
• Monte-Carlo analysis is to generate random samples in uncertain spaces according to their distributions. The method is easy to integrate and reliable but it requires a huge number of samples to eliminate the random effects of sampling. Other similar methods such as Latin Hypercubic Sampling (LHS) or orthogonal design of experiments are proposed to reduce the total number of tests but are still expensive when the simulation itself is very heavy.
• Min-Max analysis is to run simulations with the combination of upper and lower bounds of uncertain parameter intervals to estimate the worst case under perturbation [1]. It needs less tests compared to random sampling and the result is reasonable as long as the effects of parameters are linear or quasi-linear. However, the statistic characters of the systems cannot be obtained when only the bounds of interval are considered and the robustness estimation may put too much emphasis on the case of which the possibility can be neglect able.
• Analytical methods such as Taylor expansion or direct interval analysis [2] are very efficient because only partial differential equations of the systems are needed instead of a large amount of Наука техника. Т. 19, № 1 (2020) Science and Technique. V. 19, No 1 (2020) и simulations. However, the mathematical expressions in the industrial problems are difficult to obtain while simplifying a complex system into an academic model may raise the problem of representativity.
One of the targets during project development is to integrate the robustness into the optimization procedure where the optimum of vehicle specifications and their robustness are searched at the same time. It becomes even more important to reduce the additional simulation number. In this article, an approach based on polynomial chaos expansion (PCE) is applied where only a limited local design of experiments is need for each design point.
The Hermite polynomial chaos was first introduced by Wiener [3] to model stochastic response of a system under Gaussian distributed parameters. Then it has been extended into other type of orthogonal polynomial bases according to different probability distributions [4,5]. In this study Chebyshev polynomials of the second kind have been applied, which fits better real industrial uncertainties.
The target of polynomial chaos expansion is to establish a relationship between a system output and parameter inputs based on a given series of polynomials. Once the coefficient in front of each polynomial is calculated, system can be described and the robustness can be calculated. For a black box system additional tests are needed in order to decide the maximum order and important interactive terms of polynomials included in this expansion.
This paper proposes a multi-objective optimization plan with the integration of adaptive-sparse polynomial chaos expansions. The polynomial chaos expansion is calculated by a projection method which reuses response surfaces constructed in the optimization process. To further reduce the additional tests, a machine learning algorithm based on data mining is applied to justify by advance the quality of response surface before running the local design of experiments.
Section II of this paper introduces the construction of polynomial chaos expansion as well as its integration into the optimization strategy. In section III a data mining and machine learning algorithm is presented which aims to further reduce sample numbers but keep the accuracy of robust-ness prediction. Section IV demonstrates the optimization of a quarter-car example with the proposed approach.

Robust optimization with PCE
Introduction of polynomial chaos expansion and Chebyshev polynomials of the second kind. An N-dimensional random variable vector where Γ i is a one-dimension random space, its probability density function (PDF) can be defined by ( ).
i i p x Assuming that all the random variables are independent, the overall PDF P(x) can be defined as The system at one set of parameters х an be written as a converged series of polynomial basis ( ) The inner product of two orthogonal polynomial basis can be expressed as where δ ij -Kronecker delta, δ ij = 1 when i = j and To simplify the expression, the equation (2) The form of the polynomial basis i φ is defined by the distribution types of random parameters. Chebyshev polynomial of the second kind is applied in this study. As its probability density function corresponded is Wagner semicircle distribution defined in (Fig. 1) The semicircle law is more suitable to industry problems where the parameters have more weights around the design values and will never go to infinity. The polynomials can be defined by a recurrence method: The distribution function is normalized to make sure the integration in [-1.1] equals to 1, so in fact it becomes a semi ellipse ( ) The equation (3) can be expressed as for two one-dimension polynomials In equation (4) coefficients s i are unknown and need to be calculated. For a black-box system, two non-intrusive methods exist to calculate s i .
• Regression method: the order and the terms of the expansion are assumed a priori. The coefficients are then calculated by linear regression method based on the samples around the reference point. By iterations the algorithm will decide whether to add or remove terms from the expansion until all the important terms are included [6].
• Projection method: an analytical expression has been pre-defined used to represent the blackbox system. Then based on this expression, each coefficient is calculated one by one according to orthogonal projection. Additional samples are also necessary to justify if the expansion order has been converged [7].
In the optimization process expressed in Section III, the analytical expressions for the result outputs have already been calculated. The projection method will be introduced in detail and applied. Another advantage for the projection method is that the accuracy of robustness estimation depends more on the quality of analytical expressions but is less sensitive to the number of samples compared to the regression method.
Assuming f(x) is a performance function defined in the design space, one can obtain by multiplying ( ) j φ ξ to both sides of (4): The reference [8] has proposed a decomposition method of ( ) f x and the numerator of equation (9) can be expanded into univariate and bivariate terms where k j s -coefficient for the j th order univariate polynomial term of ( ) the j th order bivariate polynomial term of ( ) It can be noticed that equations (10) and (11) require only one-or two-dimensional integrations. The higher order interaction terms are neglected because of the weak non-linearity in mechanical systems, which saves greatly the computational resources.
One determination coefficient R 2 is calculated to justify the correlation between the PCE and test results where J -order of the PCE; x I -I th sample of local experiment design; ( ) It can be seen that the quality of prediction is improved with ( ) 2 R J approaching to 1. When the PCE order is increased by iterations, the PCE is converged if the difference of 2 R between two orders is smaller than a threshold ε 1 which means including new terms has nearly zero impact on the expansion ( ) ( ) Furthermore, as 2D-integration is usually more expensive to calculate, another indicator is defined to justify if the bivariate terms of ( ) If the equation (14) is satisfied it means the contribution of j th order bivariate terms of ( ) x is negligible with respect to the expansion and for the next order this term will not be calculated. A summary of the procedure of calculation of PCE around one reference solution can be seen in Fig. 2.
Once the PCE is constructed, the mean value and standard deviation of ( ) f x can be estimated by:

Fig. 2. Summary of PCE calculation by PCE method for one reference solution
Multi-objective optimization. Chassis systems has multiple demands on durability, ride comfort and handling, etc. while each of them is often incompatible with the others. Therefore, during the design phase, the optimization of chassis system is naturally multi-objective where the best compromises are searched.
The robustness of objectives included in the optimization plan are also integrated and listed as the objectives for optimizing. This step tends to make the objectives even more incompatible because empirically the optimums are less robust compared to the less good solutions. Thus, this study aims to find the relationship between the objectives and their robustness and to propose the optimums which are less sensitive to perturbations.
The mathematical expression of multi-objective optimization can be expressed as: -minimize: -under the constraints: Instead of summing all the compositions f(x) in F(x), the optimization will treat each objective separately. As a result, the optimization plan will propose a set of compromised optimums instead of only one solution, which will form a Pareto front [1,9,10]. The definition of one Pareto optimum is one which is not dominated by any other solutions. There may exist one or more solutions which have better performances in some objectives, but they must have worse solutions in other aspects than those of Pareto optimums.
Another characteristic of industrial problems is that the systems are usually black-box models and the mathematical expression of F(x) does not exist. Therefore, a Meta-model

( )
Meta F x will be constructed before the optimization iterations begin [11,12].

Meta F
x consists of several response surfaces which describe the black-box system with different combination of polynomials. The equation (17) can be replaced by as: , ; , , all the response surfaces regarding to their quality by cross validation [13]. The meta-model can then replace the black-box model and be used with the genetic optimization algorithm NSGA-II [14,15].
For the robustness objectives, the calculation of PCE is based on the Meta-model constructed in this step. As the Meta-model has several expressions to represent the system, several estimations of PCE can also be made. The final robustness is a weighted sum of these estimations according to the quality of response surfaces: when the PCE is converged to 1 ε (see (13)); i p -weight indicator which judges the quality of i th response surface by comparing it with 1.
The procedure of robust multi-objective optimization is summarized in Fig 4. By iterations, the potential optimums can be proposed by the genetic algorithm based on the objectives' response surfaces. The proposed solutions and their robustness will then be validated by real numerical simulations, of which the results are reused to improve the quality of response surfaces. If the optimization is converged, the solutions proposed will form a Pareto front.
It should be noted that although both the response surfaces and PCEs are in the form of polynomials, the domain of these polynomials are different. For response surface, the variables can cover any values in the design space while for PCE the variables are limited in the neighbourhood of one design value. PCE describes the local behavior of a specific point on the response surface. That's why the global quality of response surface has a great influence on the robustness estimation.

Application of data mining and machine learning
The calculation of PCE by projection method requires much fewer tests compared to purely random sampling methods such as Monte-Carlo, which has been shown by the example in Section IV. However, the integration of robustness into the optimization requires hundreds of PCE calculation which is still expensive even if the number of tests required to compute the robustness at one point has been greatly reduced. Two approaches have been applied at the same time to further reduce the number of samples tested for each point.
The first approach is called inherited design of experiments referenced in [13]. The simulations will stock in a data base and be reused in the local DOE for future calculation of PCE if there already exist test results in a new coming point's neighbourhood. With the enrichment of the data base, the extra number of simulations tends to be reduced.
The second approach is to exploit further the data base with a data-mining algorithm and to learn to construct the PCE without extra samples. In the procedure of PCE calculation by projection method in Section II, the local DOE is used to converge the projected terms and orders in PCE by analyzing the quality of approximation between the PCE response surfaces and the simulation results. Unlike the regression method, where the coefficients of polynomials are calculated directly from the DOE results, the projection method depend mostly on the modelling quality of meta-model. If one can predict the weighting factor in equation (20) of each response surface of meta-model on the point to be studied, there will be no need to run extra simulations for PCE calculations.
The data-mining strategy referenced in [10] which is firstly used for post-processing of Pareto front starts firstly by calculating the normalized correlation distances of different information between each pair of design points i and j: where , x ij D , According to the definition, the correlation between two designs will be good when 1.

ij D →
The inter-design matrix , x D inter-objective matrix F D and weighting factor matrix P D are symmetric matrix and can be used to analyze the correlation of each pair of points in the data base (Fig. 7, 8). Another mixed matrix M can be defined to show the correlation of both input and output between two designs: where , x c F c -coefficients of mixture defined between 0 and 1.
In order to make the matrix more readable, a bipolarization algorithm cited in [10] will used to arrange the order of design point according to their resemblance level. The algorithm starts with finding the two most different design point in М ij and grouping their neighbors based on a resemblance threshold .
s This operation is looped for the rest of non-grouped points until all the points are arranged. The algorithm schema is shown in Fig. 3.
An example of the result after bipolarization is shown in Fig 3. It is a process similar to clustering method in data mining which also regroups the existing solutions according to several features (in this study the parameters and objectives). The typologies of the database can be exploited as a post-processing to find the orientation of multiobjectives in each group. The grouping result will also serve as a base for learning in the next steps.  Fig. 10. It shows that for the resembling design points, the weighting factors tend also to be alike. It is reasonable because the points in the same group also tend to be neighborhood in the meta-model, where the modelling quality of response surfaces is similar in this area. The neighbourhood has been displaced for a different group thus the modelling quality of each response surface has also changed. Наука техника. Т. 19, № 1 (2020) Science and Technique. V. 19, No 1 (2020) и The principle of learning is that if the new design points of which the robustness is to be calculated fall into one existing group, the quality of estimation of each response surfaces can also be learned based on the P ij D information in the group. As a result, the calculation needs no extra simulation.
The algorithm of learning for a new coming design point can be summarized as following: 1) a new line/column is added into the correlation matrix , 2) the minimum correlation distance is to be found between the new point and each existing group. A pre-defined threshold will judge if the new coming point belongs to any group in the data base; 3) if no groups can be referenced. The robustness will be calculated by the local DOE method in Section II. If there exists at least one group that can include the new design point, the vector P of weighting factors can be learned from those in this group. The combined PCE estimation can be obtained without running simulations.

Example with a quarter car model
A quarter car model has been applied to demonstrate the robust multi-objective optimization strategy proposed in this paper. The model has been shown in Fig. 5: 1 m -unsprung mass which sums the mass of the wheels and a part of half suspension; 2 m -sprung mass of a quarter of the car body; 1 , k 1 c -tire stiffness and damping rate; 2 , k 2 c -suspension vertical stiffness and damping rate; 0 x -coordinate of the road profile; 1 2 , x xvertical displacements of the unsprung/sprung masses.
In this example it is assumed that the masses 1 2 , m m and the properties of tire 1 1 , k c are given and the target is to optimize the properties of suspension parameters: 2 k and 2 . , where ed F -equivalent damage force; max F -maximum force in Basquin model when total damage D = 1; N -number of cycles corresponding to D = 1; d i , n i -cumulative damage and its repetitions in the force signal in the simulation which are obtained by rainflow-counting algorithm; B -constant number of the Basquin model.
The second objective is the robustness of F ed due to the uncertainty added in the system: the variation of k 2 and 2 , 1...4, i c i = to represent manufactory tolerances and aging during usage aсcompanied with passing velocity and road amplitudes as validation process perturbations.
The non-dominated design points between two objectives are shown in Fig. 6 with the PCE method and learning algorithm. The optimum solutions in durability have relatively worse robustness due to perturbations. Tab. 1 shows the comparison between a much larger sampling size (1000 tests) with Latin Hypercubic method, the local DOE with only 8 tests and learning algorithm without tests for 3 design points. The reference point is the design point when the optimization starts. Optimized point 1 is one of the optimums orientated to robustness and optimized point 2 is one oriented to durability. It can be seen that the PCE method succeeded in produce the closed estimations as ones from a larger sampling especially for the first two points. Estimation for the design points oriented to durability is less good but it still shows the ten-dency of the relationship between the objective and its robustness. The matrix of D x , D F , M regrouped M and D P in each group are shown in Fig. 7-11 after 12 th iteration. It can be seen that the existing design points can be distributed into several small groups where the weighting vectors tends to resemble each other. This offers a good learning basis for new coming design points.
In this optimization totally 113 potential design optimums of robustness have been proposed during meta-model's construction and optimization iterations, which means millions of simulations of durability objectives need to be run if all the robustness is validated by a Monte-Carlo sampling. The size of one local DOE for PCE calculation is set to 8. If a PCE method without either inherited DOE nor learning method has been applied, about 1000 simulations would need to be run theoretically. With the integration of inherited DOE, the total number of simulations can be reduced to 802.
The learning algorithm permits to reduce further 22 % amount of simulations to 628 and the robust estimation keeps the same level. The gain is encouraging when one single simulation takes a long time.    Fig. 12 shows the comparison of the Pareto front between a robust optimization with learning and one without learning. It can be seen that there are some local differences between two strategies, but the tendency of durability objective and its robustness are close. Fig. 13 shows the number of design points suitable for learning algorithm and average additional samples for one PCE calculation according the optimization iterations. The curves have oscillated at the beginning and become stable from 7 th iteration. Most of design points are groupable from 7 th iteration which means the test data base is complete especially at the region close to the Pareto front.  Fig. 13. Evolution of number of groupable numbers (in triangles) and additional tests (in circle) per point with iteration

CONCLUSION
In this paper a robust optimization method has been proposed which aims to reduce the number of tests during the optimization and gives a global view on the relationship between the design objectives and their robustness. The integration of a learning algorithm based on data mining permits to further reduce the necessary simulations. The example shows the data mining plan is suitable for small size data bases and once it is constructed the estimation of robustness needs no more simulations.