Optimizing damper placement in automotive mechanical suspension using Bayesian Optimization (BO)

Overview

A solution written in Java to determine the optimal damper parameters for a full-car suspension system to minimize three objectives: ride comfort (ISO 2631 standard), vibration (frequency-weighted power spectral density), and handling (tire load variation). The system must account for various road profiles (urban, highway, off-road) and damping profiles (linear, digressive, progressive). Utilizes Bayesian optimization with a surrogate model to efficiently explore the 18-dimensional parameter space.

Given Data

Vehicle Parameters

Sprung mass: \( m_s = 1500 \, \text{kg} \)
Unsprung mass per wheel: \( m_u = 50 \, \text{kg} \)
Track width: \( t = 1.5 \, \text{m} \)
Wheelbase: \( l = 2.7 \, \text{m} \)
Roll inertia: \( I_x = 500 \, \text{kg·m}^2 \)
Pitch inertia: \( I_y = 800 \, \text{kg·m}^2 \)
Tire stiffness: \( k_t = 200000 \, \text{N/m} \)
Tire damping: \( c_t = 500 \, \text{Ns/m} \)

Damper Parameters

18 parameters with bounds (e.g., front compression damping \( c_{fc} \in [2000, 8000] \, \text{Ns/m} \), blow-off threshold \( v_b \in [0.5, 1.5] \, \text{m/s} \), etc.)
Damping profiles: linear, digressive, progressive

Road Profiles

Urban: frequency \( f \in [1, 10] \, \text{Hz} \), amplitude \( A = 0.02 \, \text{m} \), noise \( \sigma = 0.01 \, \text{m} \)
Highway: \( f \in [0.1, 2] \, \text{Hz} \), \( A = 0.05 \, \text{m} \), \( \sigma = 0.005 \, \text{m} \)
Off-road: \( f \in [0.5, 15] \, \text{Hz} \), \( A = 0.1 \, \text{m} \), \( \sigma = 0.02 \, \text{m} \)

Simulation Parameters

Time step: \( \Delta t = 0.01 \, \text{s} \)
Simulation duration: \( T = 5 \, \text{s} \)
Maximum force: \( F_{\text{max}} = 10^6 \, \text{N} \)
Maximum displacement: \( z_{\text{max}} = 0.5 \, \text{m} \)

Optimization Parameters

Maximum iterations: 100
Initial random points: 30
Pareto front size: 100 solutions
Reference point for hypervolume: \( [100, 10^6, 100] \)

Assumptions

The vehicle is modeled as a full-car suspension system with 7 degrees of freedom: chassis vertical displacement (\( z \)), roll (\( \theta \)), pitch (\( \phi \)), and four unsprung mass displacements (\( z_{fl}, z_{fr}, z_{rl}, z_{rr} \)).
Road inputs are sinusoidal with Gaussian noise, varying by road profile.
The objective function is expensive to evaluate, necessitating a surrogate model for optimization.
Objectives are conflicting (e.g., improving comfort may worsen handling), requiring a Pareto front for trade-offs.
All forces and displacements are clamped to prevent non-physical values.

Solution

Step 1: Model the Full-Car Suspension System

The suspension system is modeled using Newton’s second law for the sprung and unsprung masses. The state vector is:

\[ \mathbf{s} = [z, \theta, \phi, z_{fl}, z_{fr}, z_{rl}, z_{rr}, \dot{z}, \dot{\theta}, \dot{\phi}, \dot{z}_{fl}, \dot{z}_{fr}, \dot{z}_{rl}, \dot{z}_{rr}]^T \]

1.1 Suspension and Tire Forces

For each corner \( i \in \{fl, fr, rl, rr\} \), compute the suspension displacement and velocity:

\[ z_{\text{susp},i} = z + \epsilon_{t,i} \cdot \frac{t}{2} \theta + \epsilon_{l,i} \cdot \frac{l}{2} \phi \] \[ \dot{z}_{\text{susp},i} = \dot{z} + \epsilon_{t,i} \cdot \frac{t}{2} \dot{\theta} + \epsilon_{l,i} \cdot \frac{l}{2} \dot{\phi} \]

where \( \epsilon_{t,i} = [+1, -1, +1, -1] \), \( \epsilon_{l,i} = [+1, +1, -1, -1] \) for corners \( [fl, fr, rl, rr] \).

The relative velocity across the damper is:

\[ v_{\text{rel},i} = \dot{z}_{\text{susp},i} - \dot{z}_i \]

Spring Force:

\[ F_{\text{spring},i} = k_i (z_{\text{susp},i} - z_i) \]

where \( k_i = k_f = 25000 \text{ to } 40000 \, \text{N/m} \) (front) or \( k_r \) (rear).

Damping Force:

\[ F_{\text{damp},i} = c_{\text{eff},i} \cdot (v_{\text{rel},i} \cdot r_i) / r_i \]

where \( r_i \) is the motion ratio, and \( c_{\text{eff},i} \) depends on the damping profile:

Linear: \( c_{\text{eff},i} = c_r \) (rebound, \( v_{\text{rel},i} \geq 0 \)) or \( c_c \) (compression).
Digressive: If \( |v_{\text{rel},i}| > v_b \), \( c_{\text{eff},i} \cdot \frac{v_b}{|v_{\text{rel},i}|} \).
Progressive: If \( |v_{\text{rel},i}| > v_b \), \( c_{\text{eff},i} \cdot \sqrt{\frac{|v_{\text{rel},i}|}{v_b}} \).

Suspension Force:

\[ F_{\text{susp},i} = \text{clamp}(F_{\text{spring},i} + F_{\text{damp},i}, -F_{\text{max}}, F_{\text{max}}) \]

Tire Force:

\[ F_{\text{tire},i} = \text{clamp}(k_t (z_i - z_{\text{road},i}) + c_t \dot{z}_i, 0, F_{\text{max}}) \]

where \( z_{\text{road},i} \) is the road input.

1.2 Equations of Motion

The net forces and moments on the sprung mass are:

\[ F_z = F_{\text{susp},fl} + F_{\text{susp},fr} + F_{\text{susp},rl} + F_{\text{susp},rr} \] \[ M_{\text{roll}} = \frac{t}{2} (F_{\text{susp},fl} - F_{\text{susp},fr} + F_{\text{susp},rl} - F_{\text{susp},rr}) \] \[ M_{\text{pitch}} = \frac{l}{2} (F_{\text{susp},fl} + F_{\text{susp},fr} - F_{\text{susp},rl} - F_{\text{susp},rr}) \]

For each unsprung mass:

\[ F_{u,i} = -F_{\text{susp},i} + F_{\text{tire},i} \]

Update accelerations:

\[ \ddot{z} = \frac{F_z}{m_s}, \quad \ddot{\theta} = \frac{M_{\text{roll}}}{I_x}, \quad \ddot{\phi} = \frac{M_{\text{pitch}}}{I_y}, \quad \ddot{z}_i = \frac{F_{u,i}}{m_u} \]

1.3 State Update (Euler Integration)

\[ \dot{z}(t + \Delta t) = \dot{z}(t) + \ddot{z} \Delta t, \quad z(t + \Delta t) = z(t) + \dot{z}(t) \Delta t \]

Similarly for \( \theta, \phi, z_i \). Displacements are clamped to \( [-z_{\text{max}}, z_{\text{max}}] \).

1.4 Road Inputs

Road displacement for corner \( i \):

\[ z_{\text{road},i}(t) = A \sin(2\pi f t - \delta_i) + \mathcal{N}(0, \sigma) \]

where \( f \in [f_{\text{min}}, f_{\text{max}}] \), \( \delta_i = 0.1 \, \text{s} \) for rear wheels, and noise \( \sigma \) depends on the profile.

Step 2: Compute Objective Metrics

Simulate for \( T = 5 \, \text{s} \) (\( N = 500 \) steps) per road profile, collecting chassis acceleration \( \ddot{z} \), displacement \( z \), and tire load variations.

2.1 Comfort (ISO 2631)

Root-mean-square (RMS) acceleration:

\[ a_{\text{rms}} = \sqrt{\frac{1}{N} \sum_{i=1}^N \ddot{z}_i^2} \]

where \( \ddot{z}_i = F_z / m_s \), clamped to \( [-10^3, 10^3] \, \text{m/s}^2 \).

2.2 Vibration (Frequency-Weighted PSD)

Approximate PSD by weighting velocities:

\[ v_i = \frac{z_i - z_{i-1}}{\Delta t} \] \[ \text{PSD} = \frac{1}{N-1} \sum_{i=1}^{N-1} w(f_i) v_i^2 \]

where weight \( w(f_i) = 1.0 \) if \( f_i \in [1, 10] \, \text{Hz} \), else \( 0.1 \), and \( f_i \approx \frac{1}{i \Delta t} \).

2.3 Handling

Tire load variation:

\[ F_{\text{tire},i} = k_t z_i + c_t \dot{z}_i \] \[ \sigma_{\text{tire}} = \sqrt{\frac{1}{4} \sum_{i=1}^4 (F_{\text{tire},i} - \bar{F}_{\text{tire}})^2} \]

where \( \bar{F}_{\text{tire}} \) is the mean tire load, clamped to \( [0, 10^4] \, \text{N} \).

2.4 Objective Function

For each road profile \( p \), compute metrics \( [a_{\text{rms},p}, \text{PSD}_p, \sigma_{\text{tire},p}] \). Average across profiles:

\[ O_j = \frac{1}{3} \sum_{p=1}^3 m_{j,p} + P \]

where \( j \in \{\text{comfort}, \text{vibration}, \text{handling}\} \), scaled as:

\[ O_1 = \frac{a_{\text{rms}}}{3}, \quad O_2 = \frac{\text{PSD}}{1000}, \quad O_3 = \frac{\sigma_{\text{tire}}}{1000} \]

Penalty \( P = 0.1 \) per motion ratio \( r_i \notin [0.5, 2.0] \), where:

\[ r_i = \text{clamp}\left( \frac{1}{\sqrt{x_i^2 + y_i^2 + z_i^2} \cdot \cos(\alpha_i) + 0.1}, 0.5, 2.0 \right) \]

Step 3: Bayesian Optimization Framework

The objective function is expensive, so we use Bayesian optimization with a surrogate model to select parameters \( \mathbf{x} = [c_{fc}, c_{fr}, \ldots, \text{anti-dive}] \in \mathbb{R}^{18} \) and damping profile \( d \in \{\text{linear}, \text{digressive}, \text{progressive}\} \).

3.1 Surrogate Model

Use an ensemble of Gaussian Process (GP), Random Forest (RF), and Neural Network (NN) for each objective:

\[ \hat{O}_j(\mathbf{x}) = w_{\text{gp},j} \hat{O}_{\text{gp},j} + w_{\text{rf},j} \hat{O}_{\text{rf},j} + w_{\text{nn},j} \hat{O}_{\text{nn},j} \]

GP: RBF kernel, \( k(\mathbf{x}, \mathbf{x}') = \exp\left(-\frac{\|\mathbf{x} - \mathbf{x}'\|^2}{2}\right) \), with noise \( \sigma_n = 10^{-6} \). Predict mean and variance:

\[ \mu(\mathbf{x}) = \mathbf{k}^T (\mathbf{K} + \sigma_n^2 \mathbf{I})^{-1} \mathbf{y}, \quad \sigma^2(\mathbf{x}) = k(\mathbf{x}, \mathbf{x}) - \mathbf{k}^T (\mathbf{K} + \sigma_n^2 \mathbf{I})^{-1} \mathbf{k} \]

RF: 5 trees, max depth 8, predict mean as average and variance as tree disagreement.
NN: 2 hidden layers (32 nodes), ReLU activation, Monte Carlo dropout for variance.

Weights \( w_i \) updated every 5 points via 3-fold cross-validation error.

3.2 Acquisition Function

Select the next point by maximizing the active learning score:

\[ S(\mathbf{x}) = \alpha (V + D) + (1 - \alpha) \text{EHVI} \]

Variance: \( V = \frac{1}{3} \sum_{j=1}^3 \sigma_j^2(\mathbf{x}) / V_{\text{max}} \), \( V_{\text{max}} = 1 \).
Diversity: \( D = \min_{\mathbf{x}_i \in \text{Pareto}} \|\mathbf{x} - \mathbf{x}_i\| / D_{\text{max}} \), \( D_{\text{max}} = \sqrt{18 \cdot 100} \).
Expected Hypervolume Improvement (EHVI):

\[ \text{EHVI} = \frac{1}{100} \sum_{s=1}^{100} \left[ HV(\text{Pareto} \cup \{\mathbf{y}_s\}) - HV(\text{Pareto}) \right] \]

where \( \mathbf{y}_s = [\mu_j + \sigma_j \mathcal{N}(0,1)] \), and hypervolume is computed w.r.t. reference point \( [100, 10^6, 100] \). Include exploration term \( 0.1 \sum \sigma_j \).

Optimize \( S(\mathbf{x}) \) using Simulated Annealing with initial temperature \( T_0 = 100 \), cooling rate \( 0.95 \), and 200 iterations.

3.3 Pareto Front

Maintain a Pareto front of non-dominated solutions:

A solution \( \mathbf{y}_1 \) dominates \( \mathbf{y}_2 \) if \( y_{1,j} \leq y_{2,j} \forall j \) and \( y_{1,k} < y_{2,k} \) for some \( k \).
Limit to 100 solutions, removing least crowded points based on:

\[ d(\mathbf{s}) = \sum_{j=1}^3 \frac{y_{\text{next},j} - y_{\text{prev},j}}{r_j} \]

where \( r_j = [100, 10^6, 100] \).

3.4 Algorithm

Initialize 30 random points \( \mathbf{x}_i \), evaluate \( \mathbf{y}_i = [O_1, O_2, O_3] \), train surrogates.
For 100 iterations:
- Compute adaptive weights based on crowding entropy.
- Optimize acquisition function to select \( \mathbf{x}_{\text{next}} \).
- Evaluate \( \mathbf{y}_{\text{next}} \), update surrogates and Pareto front.
Output Pareto-optimal solutions.

Step 4: Implementation Details

Simulation: Runs Euler integration for 500 steps per road profile, averaging metrics.
Parallelization: Uses parallel objective evaluations for efficiency.
Error Handling: Invalid objectives (\( \text{NaN} \), \( \infty \), \( \geq 10^{12} \)) are retried or penalized.
Output: Prints each solution’s objectives and parameters.

Optimal SOlution

The optimization yields a Pareto front of solutions, each specifying 18 damper parameters and a damping profile. For example, a solution might be:

Comfort: \( 2.5 \, \text{m/s}^2 \), Vibration: \( 500 \, \text{(scaled units)} \), Handling: \( 300 \, \text{N} \).
Parameters: \( c_{fc} = 4000 \, \text{Ns/m} \), \( c_{fr} = 4500 \, \text{Ns/m} \), …, Profile: digressive.

The Bayesian optimization framework efficiently explores the parameter space, balancing exploration (via variance and diversity) and exploitation (via EHVI), producing a diverse set of trade-off solutions for suspension design.

Repository