Introduction

The handling stability, running quality, and vehicle ride comfort is the criteria used by the passenger for assessing on-road vehicles. The handling attributes are determined by the way the forces at the tire contact patch influence the vehicle dynamics. The suspension is a complex multi-body system that connects and transmits the force between the vehicle body and the wheel. The forces can be due to a driver command such as a steering wheel, gas pedal or brake input, but also to road surface unevenness, aerodynamic forces, vibrations of the engine and drive-line, and non-uniformity of the tire/wheel assembly.

The desired behavior of the vehicle can be achieved by optimizing the design of the front suspension, which comprises a large number of design variables. One of the most popular suspension systems is the MacPherson strut [14,15], which is deployed in a large number of vehicles for its simple structure, good overall performance, package efficiency and relatively low cost as shown in Fig. 1. Fig. 1: MacPherson Suspension Design [14,15] In this work, we develop a Bayesian optimization strategy for determining parameters of a MacPherson suspension system that optimizes desired performance characteristics of an automobile. In contrast to the industry standard of hand designed parameters, we aim to achieve partial automation of the design process by employing a multibody dynamics simulation model and allowing the Bayesian optimization to iteratively explore the design space without human intervention.

Applications of Bayesian Optimization

Inverse problems are generally the problem of finding the properties of the model using indirect measurements such as finding an inverse shortest path and inverse minimum spanning tree problems [28]. Vito et al. [27] formulated the problem of finding the best solution as a linear inverse problem given some loss function and hypothesis space. This problem can also be solved using the regularized least squares algorithm. These types of problems are ill-posed problems [30]. They are often solved using Bayesian approaches. In the proposed approach, the design of the suspension of the vehicle is framed as the inverse problem and the inverse problem is solved using Bayesian Optimization.

Bayesian Optimization (BO) creates a surrogate model of the black-box objective function and finds the optimum of an unknown objective function from which samples can be obtained [16,4]. BO has been studied in many fields for the optimization of unknown functions. The tuning of the hyper-parameters, regularization terms, and optimization parameters in machine learning needs a lot of expert experience, but with the help of the tractable posterior distribution induced by the Gaussian Processes [23] leads to the efficient use of the information gathered from the past experiences, enabling optimal choices for the next parameters. The performance of the BO is highly dependent on the choice of acquisition function made to find the best-observed value. [18,13] have proposed convergence criteria for several acquisition functions in order to avoid unwanted evaluations.

Several methods exist in the literature for designing suspension architectures. Yu and Yu [31] formulated the problem of optimal vehicle suspension design with the quarter-car vehicle dynamic model using a genetic algorithm (GA). The two objectives of the optimization are to minimize the average suspension displacement and to minimize the maximum bouncing acceleration of the sprung mass [3] restricted to a number of constraints. An objective function is formulated as the sum of absolute magnitude of sprung mass (chassis) acceleration and absolute magnitude of average sprung mass displacement. The authors have considered some constraints in the problem formulation such as the maximum amplitude of sprung mass acceleration should not be less than 1g (9.8m/s 2 ), human beings feel motion sickness when vibration at a natural frequency is less than 1Hz, and the absolute value of the relative displacement between the sprung and unsprung mass should be less than 13 cm. Zhang et al. [32] proposed the design of suspension parameters of flexible multibody vehicle model using the ADAMS software and GA. A flexible multibody model having 44 degrees of freedom is compared with the multi-rigid bodies having 33 degrees of freedom. ADAMS is a tool used by the leading manufacturers to design parameters for complex design and optimization. ADAMS and GA carries the simulation analysis until the algorithm converges. Yao et al. [29] modeled a sedan MacPherson front suspension with ADAMS and studied the relationship between the wheel alignment parameters and tire wear when wheel hops. Afkar et al. [1] modeled a double wishbone suspension system using ADAMS and then optimized geometric parameters using a GA. Similarly, the authors in [20] optimized several hard points of the front suspension using the INSIGHT module of ADAMS to optimize toe and camber angle characteristics. Dye and Lankarani [7] used a neural network to fit tire test data under varying tire pressures and steady state conditions. Hurel et al. [9] developed a non-linear two-dimensional mathematical quarter-car model of MacPherson suspension. This model considered both the vertical motion of the chassis and the rotation and translation for the unsprung mass (wheel assem-bly). Zhang et al. [33] studied the geometric relationship between the structural parameter of the double wishbone front suspension and alignment parameters of the steering wheel. Liu et al. [10] studied the hardpoints of the double wishbone suspension of a Formula SAE car using correlation theory. Tao et al. [25] used a Gaussian process modeling method to design vehicle suspension parameters using a multi-disciplinary optimization architecture, and experiments were carried out for a front MacPherson suspension and a rear strut suspension system on a Altair Motion View vehicle model. The Gaussian process models are fitted with multidisciplinary optimization and efficiencies are computed.

In automobile manufacturing, whether ride and handling characteristics have been achieved is quantified through statistics of a kinematic curve. The design process can be time intensive as an engineer has to manually input a particular set of hard-points, pass it through the simulation tool (which can take several minutes to process), observe the kinematic curves, and then repeat the whole process to refine the design. Thus, generating one design that meets the desired kinematic constraints in current standard practice can take multiple weeks and consume significant manpower.

The main novel contributions of this paper are as follows:

-Bayesian optimization is applied to the design of the front suspension of a vehicle. -Finding the geometric locations of the hardpoints of the suspension from the desired target characteristics is formulated as a generalized-inverse problem. -We develop novel convergence criteria in the context of the generalizedinverse problem. -We verify the validity of the approach in several suspension design settings, empirically validating the predicted convergence characteristics.

2 Suspension Design using Bayesian Optimization

Bayesian Optimization

BO is a derivative free optimization approach for global optimization of an expensive black-box function f . It is a class of sequential-model based optimization algorithms that uses past evaluations of the function to find the next point to sample. For a minimization problem,

y * = arg min y∈Y f (y),(1)

if f is expensive, then there is a strong need to reduce the number of evaluations. Since the objective function is unknown, according to Bayesian theory the function is treated as random and a prior belief is placed over the function values, and as more values are observed, the posterior is updated based on the observation likelihood. These models use different types of acquisition functions to determine the next point to be evaluated based on the current posterior distribution over functions.

The surrogate model used for this optimization is a Gaussian Processes (GP). A GP is characterized by its mean µ(y), and co-variance function k(y, y ). For n data points, the function f 1;n = f (y 1 ), . . . , f (y n ) can be characterized using a multivariate Gaussian distribution as

f 1:n = N (µ (y 1:n ), K),(2)

where K is a n × n kernel matrix given by

K =    k(y 1 , y 1 ) . . . k(y 1 , y n ) . . . . . . . . . k(y n , y 1 ) . . . k(y n , y n )    ,(3)

for some positive definite kernel such as a Gaussian or Matern kernel [21].

An acquisition function proposes which points should be selected next to determine the minimizer of the function, and they trade off between exploitation and exploration. Exploitation means the acquisition function selects points where the mean of the GP is low and exploration means the acquisition function selects points where the variance of the GP is high. Various acquisition functions are used for Bayesian models, including maximum probability of improvement (MPI), expected improvement (EI), confidence bounds (CB) criteria.

Maximum Probability of Improvement evaluates f (y) at the point most likely to improve than the minimal value of f (y + ) observed so far [4] M P I(y) = P (f (y) ≥ f (y

+ ) + ξ), = Φ µ(y) − f (y + ) − ξ σ(y) ,(4)

where f (y + ) is the best observed value of the function so far, µ(y) is the posterior mean of y under the GP, σ(y) is the posterior standard deviation, and ξ is a parameter used to drive exploration, which is usually very small. Φ is the cumulative distribution function of a standard normal variable. Expected Improvement evaluates f (y) at the point where the expectation of the improvement in f (y + ) under the current estimate of the GP is the highest [22]

EI(y) = E [max{0 , f (y + ) − f (y)}], = (−µ(y) + f (y + ) + ξ) Φ −µ(y) + f (y + ) + ξ σ(y) + σ(y) φ −µ(y) + f (y + ) + ξ σ(y) , (5)

where φ is the probability density function of a standard normal variable.

The basic algorithm of Bayesian Optimization [19] is given in Algorithm 1.

Algorithm 1 Bayesian Optimization

Initial Design:

Dn init = {(yi, xi)} n init i=1 , nmax

Generalized Inverse Problems using Bayesian Optimization

Automotive engineers use mechanical discipline models to design the suspension of the vehicle, which is computationally expensive and may involve a design cycle spanning several days or weeks. Operating characteristics of the vehicle suspension components depend upon the suspension architecture, selection of components and the locations at which the components are attached to each other and the vehicle body. The suspension geometry includes the arrangement of suspension components and the location of hard-points, which are locations at which suspension components are attached to the vehicle's body. For a given set of vehicle suspension components, the suspension geometry will determine the kinematic characteristics of the vehicle suspension which includes wheel attitude and suspension travel. Kinematic characteristics are traditionally derived from a multi-body dynamic modeling/simulation of an entire suspension system subjected to spindle input loads such as may be performed by commercially available multi-body dynamic (MBD) software. Given a set of hard-point locations and corresponding kinematic curves we learn the relationship between the hard-point locations and kinematic curves using BO. This work predicts a set of kinematic curves corresponding to hard-points for MacPherson strut suspension. Predicting the kinematic curves corresponding to hard-points can reduce the amount of time required to select hard-points for a MacPherson strut suspension design process from multiple weeks to less than one day.

Statistics of kinematic curves of the vehicle form our representation of the desired vehicle characteristics. We denote the set of possible statistics as X where X ⊆ R d . This paper addresses the design of suspension hardpoints Y where Y ⊂ R m is a bounded domain for a given target x ∈ X . It is a well defined task for the engineers to design target characteristics X from the suspension design parameters Y. Consider the equation

x = g(y).(7)

We formulated this problem as an inverse problem for g : Y → X , find a value y such that g(y) ≈ x. In general the above problem is ill-posed: there may be multiple y satisfying g(y) = x for a given x or none. The existence of the solution can be restored by solving for the minimum norm solution

g † (x) := arg min y∈Y g(y) − x 2 =:f (y) ,(8)

where g † defines a generalized inverse of g [2]. We optimize the desired suspension design parameters using Bayesian Optimization (BO) where g is set to be computed by the multi-body dynamic software MSC ADAMS [17], and Bayesian optimization is used to minimize f with respect to y.

Convergence criteria

It is well known that Bayesian optimization converges under fairly general conditions to a global optimizer over a bounded domain [18,4,5]. That the algorithm will converge towards an optimum indicates that, when the requested design characteristics x are feasible, i.e.

arg min y∈Y g(y) − x 2 = 0, (9) we may simply specify a convergence criterion that terminates when g(y

+ ) − x 2 < ε,(10)

for some prespecified tolerance ε > 0. However, it may not be the case that the requested design characteristics x lie in the image of Y under g. Consider the example where passenger comfort and acceleration must simultaneously be unrealistically high: such constraints can contradict each other. In such a case it may be that for a given tolerance,

arg min y∈Y g(y) − x 2 > ε,(11)

and the optimization will never terminate. We therefore develop an additional termination criterion in the sequel. Bayesian optimization is typically presented without termination criteria, or with the assumption that there is a fixed budget of optimization iterations [8]. However, two main optimization dependent termination strategies have been proposed: (i) thresholding the acquisition function [18,13], and (ii) terminating when two sequential parameter vectors differ by a small amount [12]. We consider here the strategy of thresholding the acquisition function, as this has been analyzed from the perspective of regret minimization, which implies a bound on the convergence rate of the optimization [18,24]. Theorem 1 ([21], Sec. 2.9). Let σ n (f (y)) denote the predictive variance of a Gaussian process regression model at y given a dataset of size n. The predictive variance of using a dataset of only the first n − 1 training points is denoted σ n−1 (f (y)):

σ n (f (y)) ≤ σ n−1 (f (y)). (12)

Furthermore, for a covariance function that is strictly positive definite (i.e. the rank of the covariance matrix is equal to its size), the variance will be strictly decreasing over a bounded domain Y.

It is well known that Bayesian optimization with expected improvement will converge to an optimum [5,Theorem 2]. This leads us to the following proposition:

Proposition 1 (Termination of Expected Improvement). The following convergence criterion applied to Algorithm 1 will terminate:

max y∈Y EI(y) ≤ ε,(13)

for a fixed ε > 0.

Proof. From [5, Theorem 2], we have that

f (y + ) − arg min y∈Y f (y),(14)

is in expectation decreasing in the number of iterations, meaning that

− µ((y)) + f (y + ) + ξ,(15)

is decreasing. Furthermore, inspecting the formula for EI(y), we see that Φ and φ are bounded. From Theorem 1, we have that σ(y) is decreasing, and the result follows.

When the maximum expected improvement is small, by definition, we believe under the posterior GP model that the solution will not improve by a large amount when observing a new point, and we are close to a minimum norm solution. Furthermore, due to the quantification of uncertainty in the Gaussian process model, we may set a threshold on the acquisition function such that this second termination criterion becomes active when it is highly probable that no solution with near-zero norm exists.

Case Studies

We implemented the inverse problem of finding the suspension design using BO [26] in Python 3.7.3. The LBGFS algorithm with multiple restarts was used for optimizing the acquisition function. All the experiments were conducted on a workstation with two Intel Xeon CPUs and 64G memory. MSC ADAMS 18.1 was used as the discipline model, and takes several minutes per call. The problem becomes

y * = argmin y∈Y ADAM S(y) − x 2 ,(16)

for a desired set of vehicle characteristics x. We optimize the positions of Inner Tie Rod Ball and Outer Tie Rod Ball of the MacPherson front suspension architecture. We compare the convergence of Bayesian optimization using acquisition function EI to that of the commercial HEEDS Sherpa optimization software, which is commonly used in automotive design [6,11].

A classic suspension design problem is the tuning of the bump and roll steer performance, which is strictly related to the placement of the tie rod. Bump steer is a measure of the suspension toe angle with respect to suspension vertical travel and is measured in deg/m. In roll steer the wheel travel is generated by applying a roll angle to the vehicle body. Tie rod connects the steering to the steering knuckle on each front wheel. A tie rod is made of two components: inner and outer tie rod ends. This problem is relevant as the toe attitude of the wheel is a function of its vertical motion and is related to crucial vehicle attributes such as vehicle stability, cornering performance, and steering feel. Indeed, since the position of the outer tie rod defines the lever arm to the kingpin axis, the outer tie rod location determines how the forces are transmitted from the contact patch through the suspension and ultimately to the driver operating the steering wheel.

Optimization of one hardpoint

We optimized the position of Outer Tie Rod Ball as the first case study. This hardpoint has large effect on the performance of the bump steer and roll steer. We show the result of two targets such as bump steer and roll steer from the entire set of target characteristics. We achieved the convergence criterion of ADAM S(y) − x 2 ≤ 0.001 (17) in 210 iterations. The proposed design of a single hardpoint outperforms the design of the HEEDS software (Fig. 7). The acquisition function converges as shown in Fig. 2, and the norm per iteration is shown in Figs. 4-6. We have also compared to a numerical gradient based optimization method using the MAT-LAB fmincon function (Fig. 5). Fig. 7 shows that the BO approach converges more quickly than the competing techniques. Figs. 8(a)-8(c), show the points in the design space explored using HEEDS, fmincon, and Bayesian optimization. The axes of the scatter plots are normalized in the range [0, 1] with respect to the input domain of the hardpoints. Fig. 9 show the convergence of the optimized target characteristics towards the desired bump steer using HEEDS, fmincon, and Bayesian optimization. The optimized targets are normalized such that the target performance is labeled as zero and is shown with a red line. All three methods converge to similar targets, but fmincon converges with a non-zero error.

Optimization of two hardpoints

The complexity in the design increases with the increase in the number of hardpoints to be optimized. We optimized the position of Inner Tie Rod Ball and Outer Tie Rod Ball in the second case-study. The norm as a function of the iteration number for HEEDS, fmincon, and BO is shown in Figs. 10, 11, and 12, respectively. The convergence in the BO based approach is faster compared to HEEDS, as shown in Fig. 13. The normalized scatter plots for each of the hardpoints is shown in Figs. 14 and 15. Fig. 16 shows the convergence to the target bump steer for all three methods.

Conclusion

The performance of automotive design in terms of ride comfort and vehicle dynamics has a direct impact on passenger experience and customer satisfaction. The suspension system is an important automotive system that carries the entire load of the vehicle and also responsible for a smooth ride. This paper primarily focuses on designing the geometry of the MacPherson front suspension system to satisfy desired kinematic performance metrics using the MSC ADAMS Car model. Bayesian optimization is used to optimize the coordinate values of hardpoints of the front suspension. Conventionally, optimization of suspension design has been done using HEEDS software in the automobile industry. Our results demonstrate that the proposed approach is able to optimize the parts of the suspension efficiently to enhance the comfort ride. The results of some case studies were described to highlight the performance of some commonly used algorithms with the proposed approach. It performs slightly better than the close source licensed HEEDS software commonly used in the automotive industry, and substantially better than optimization with a finite difference approximation to the gradient. Furthermore, we have two convergence criteria, based on the norm or on the acquisition function, which in contrast to HEEDS, allows us to determine when global convergence criteria are satisfied.