Distributionally Robust Safety Under Arbitrary Uncertainties: A Safety Filtering Approach

High-speed autonomous racing under empirical disturbances and perception noise.

• Nominal policy: uses Model Predictive Path Integral (MPPI) from this paper.
• The uncertainty arises from noisy perception of distance to track boundaries.
• The nominal vehicle model used by the planner differs from unknown true vehicle dynamics obtained from Assetto Corsa Gym.
• A distributional shift arises because DRS-gatekeeper does not explicitly account for perception noise and relies on approximate vehicle dynamics.

Vehicle Model

We consider a Dallara F317 racing car with state:

$$x_t = [p^x_t, p^y_t, \psi_t, v^x_t, v^y_t, r_t, q_t]^\top$$

where \((p^x_t, p^y_t)\) and \(\psi_t\) denote global position and heading, \((v^x_t, v^y_t)\) are body-frame velocities, \(r_t\) is the yaw rate, and \(q_t\) is the steering angle. The control input is:

$$u_t = [ a_{\text{th},t}, a_{\text{br},t},q_{\text{cmd},t}]^\top$$

with throttle \(0\leq a_{\text{th},t}\leq a_{\max}=12.0\), braking \(0\leq a_{\text{br},t}\leq a_{\text{br,max}}=18.0\), and steering commands \(|q_{\text{cmd},t}|\leq q_{\max}=0.32\).

Discrete-Time Bicycle Dynamics

We model the car's dynamics using the discrete-time bicycle dynamics \(x_{t+1} = f(x_t, u_t)\) as:

$$x_{t+1} = x_t + \Delta t \underbrace{\begin{bmatrix} \bar{v}^x_t \cos \psi_t - \bar{v}^y_t \sin \psi_t \\ \bar{v}^x_t \sin \psi_t + \bar{v}^y_t \cos \psi_t \\ \bar{r}_t \\ a_{\text{long},t} + v^y_t r_t \\ \frac{1}{m}(F_{yf} \cos q_t + F_{yr}) - v^x_t r_t \\ \frac{1}{I_z}(l_f F_{yf} \cos q_t - l_r F_{yr}) \\ k_{\delta}(q_{\text{cmd}, t} - q_{t}) \end{bmatrix}}_{\dot x_t}$$

where \(\bar{v}^x_t, \bar{v}^y_t, \bar{r}_t\) are interval-averaged quantities for numerical stability, and \(\Delta t=0.05\) is the sampling time. The longitudinal acceleration is:

$$a_{\text{long},t} = a_{\text{th},t} a_{\max} - a_{\text{br},t} a_{\text{br,max}} - \frac{C_d}{m} v^x_t |v^x_t|$$

Lateral tire forces follow a linear model:

$$F_{y,\{f,r\}} = C_{\{f,r\}} \alpha_{\{f,r\}}$$

with slip angles:

$$\alpha_f = q_t - \arctan\!\left(\frac{v^y_t + l_f r_t}{v^x_t}\right) \quad \text{and} \quad \alpha_r = -\arctan\!\left(\frac{v^y_t - l_r r_t}{v^x_t}\right)$$

Model Identification and Uncertainty

To identify model parameters, we use a dataset of offline laps on the \(4.657 \, \mathrm{km}\) Barcelona-Catalunya GP circuit obtained from simulation data, where data were collected using Assetto Corsa. Model uncertainty is quantified via one-step prediction residuals against ground-truth telemetry:

$$\begin{align} \delta v^x_t &= v^x_{t+1,\text{gt}} - (v^x_t + \dot{v}^x_t \Delta t), \\ \delta v^y_t &= v^y_{t+1,\text{gt}} - (v^y_t + \dot{v}^y_t \Delta t), \\ \delta r_t &= r_{t+1,\text{gt}} - (r_t + \dot{r}_t \Delta t). \end{align}$$

We consider an F-16 aircraft with state

$$x_t = [p^N_t, p^E_t, h_t, u_t, v_t, w_t, p_t, q_t, r_t, \phi_t, \theta_t, \psi_t]^\top,$$

where \((p^N_t, p^E_t, h_t)\) and \((\phi_t,\theta_t,\psi_t)\) denote inertial position and Euler attitude, \((u_t,v_t,w_t)\) are body-frame velocities, and \((p_t,q_t,r_t)\) are body angular rates. The control input is

$$u_t = [\delta_{A,t}, \delta_{E,t}, \delta_{R,t}, \delta_{T,t}]^\top,$$

with aileron, elevator, rudder, and throttle commands. Note that \(\delta\), \(\alpha\), and \(\beta\) in the F-16 dynamics denote control inputs, angle of attack, and sideslip angle, and are distinct from DRS-gatekeeper parameters.

We model the aircraft dynamics using a discrete-time rigid-body model \(x_{t+1} = f(x_t, u_t)\) as

$$x_{t+1} = x_t + \Delta t \underbrace{\begin{bmatrix} \dot p^N_t \\ \dot p^E_t \\ \dot h_t \\ X_t + r_t v_t - q_t w_t - g \sin\theta_t \\ Y_t + p_t w_t - r_t u_t + g \sin\phi_t \cos\theta_t \\ Z_t + q_t u_t - p_t v_t + g \cos\phi_t \cos\theta_t \\ \dot p_t \\ \dot q_t \\ \dot r_t \\ p_t + q_t \sin\phi_t \tan\theta_t + r_t \cos\phi_t \tan\theta_t \\ q_t \cos\phi_t - r_t \sin\phi_t \\ \frac{q_t \sin\phi_t + r_t \cos\phi_t}{\cos\theta_t} \end{bmatrix}}_{\dot x_t}$$

where \(\Delta t=\frac{1}{30}\) is the sampling time and \(g=32.174\,\mathrm{ft/s^2}\). The position kinematics are

$$\dot p^N_t = u_t(c_{\theta_t}c_{\psi_t}) + v_t(s_{\phi_t}s_{\theta_t}c_{\psi_t} - c_{\phi_t}s_{\psi_t}) + w_t(c_{\phi_t}s_{\theta_t}c_{\psi_t} + s_{\phi_t}s_{\psi_t}),$$ $$\dot p^E_t = u_t(c_{\theta_t}s_{\psi_t}) + v_t(s_{\phi_t}s_{\theta_t}s_{\psi_t} + c_{\phi_t}c_{\psi_t}) + w_t(c_{\phi_t}s_{\theta_t}s_{\psi_t} - s_{\phi_t}c_{\psi_t}),$$ $$\dot h_t = u_t s_{\theta_t} - v_t s_{\phi_t}c_{\theta_t} - w_t c_{\phi_t}c_{\theta_t},$$

with \(c_{(\cdot)}=\cos(\cdot)\) and \(s_{(\cdot)}=\sin(\cdot)\).

Aerodynamic coefficients are predicted via polynomial approximation

$$[\hat C_{X,t},\hat C_{Y,t},\hat C_{Z,t},\hat C_{L,t},\hat C_{M,t},\hat C_{N,t}] = \phi(\upsilon_t)^\top W + B,$$

where

$$\upsilon_t=[\alpha_t,\beta_t,\mathrm{Mach}_t,p_t,q_t,r_t,\delta_{e,t},\delta_{a,t},\delta_{r,t}]^\top,$$

and \(\alpha_t=\arctan2(w_t,u_t), ~ \beta_t=\arcsin\!\left(\frac{v_t}{V_t}\right), ~ V_t=\sqrt{u_t^2+v_t^2+w_t^2}\). Using dynamic pressure \(q_{\infty,t}=\frac12\rho V_t^2\), wing area \(S\), span \(b\), mean chord \(\bar c\), and mass \(m\):

$$X_t=q_{\infty,t}\frac{S}{m}\hat C_{X,t}+T(\delta_{t,t}),\quad Y_t=q_{\infty,t}\frac{S}{m}\hat C_{Y,t},\quad Z_t=q_{\infty,t}\frac{S}{m}\hat C_{Z,t},$$ $$L_t=q_{\infty,t}Sb\,\hat C_{L,t},\quad M_t=q_{\infty,t}S\bar c\,\hat C_{M,t},\quad N_t=q_{\infty,t}Sb\,\hat C_{N,t},$$

and angular accelerations satisfy

$$\dot\omega_t = I^{-1}\!\left(\begin{bmatrix}L_t\\M_t\\N_t\end{bmatrix} - \omega_t \times (I\omega_t)\right),\qquad \omega_t=[p_t,q_t,r_t]^\top.$$

Control inputs are normalized and bounded by \(\delta_{a,t},\delta_{e,t},\delta_{r,t}\in[-1,1]\), \(\delta_{t,t}\in[0,1]\).

Trajectory Generation

Trajectory generation for the F-16 case study follows the optimization-based workflow described in air-racing-optimization.

The canyon elevation data is obtained from an OpenTopography USGS 10 m digital elevation model of Black Canyon of the Gunnison National Park.

Polynomial Model Identification and Uncertainty Injection

The nominal aerodynamic model is identified from generated simulator data. Each rollout provides state, action, and next-state samples, which are converted into aerodynamic coefficient targets. A regularized polynomial regression is then fit from the selected flight and control features to the six aerodynamic coefficient channels:

$$\mathbf{C}_{\mathrm{nom}} = \begin{bmatrix} C_X & C_Y & C_Z & C_L & C_M & C_N \end{bmatrix}^{\top}.$$

After fitting the nominal polynomial model, residual coefficient errors are computed on the same generated dataset. Structured residual trends are removed using context-dependent calibration features, namely flight condition. The remaining centered residuals are stored as an empirical uncertainty model. During rollout, uncertainty samples are drawn from residuals observed in nearby historical contexts, so the injected uncertainty remains dependent on the current flight regime rather than being modeled as fixed Gaussian noise.

In the gatekeeper rollout, the sampled uncertainty enters directly in coefficient space. Let

$$\mathbf{w} = \begin{bmatrix} w_{C_X} & w_{C_Y} & w_{C_Z} & w_{C_L} & w_{C_M} & w_{C_N} \end{bmatrix}^{\top}$$

denote the sampled coefficient residual. The perturbed aerodynamic coefficients are

$$\mathbf{C} = \mathbf{C}_{\mathrm{nom}} + \mathbf{w}.$$

These perturbed coefficients are then converted into dimensional forces and moments:

$$X = (C_X + w_{C_X})\frac{\bar q S}{m} + T(\delta_t), \qquad Y = (C_Y + w_{C_Y})\frac{\bar q S}{m}, \qquad Z = (C_Z + w_{C_Z})\frac{\bar q S}{m},$$ $$L = (C_L + w_{C_L})\bar q S b, \qquad M = (C_M + w_{C_M})\bar q S \bar c, \qquad N = (C_N + w_{C_N})\bar q S b.$$

Thus, uncertainty is injected before rigid-body propagation by perturbing the aerodynamic coefficients, which then changes the body-axis force and moment channels used by the rollout dynamics.

Obstacle avoidance under uncertain actuation and model mismatch.

Dubins Model

We first introduce the Dubins model in continuous time with the state:

$$x=\begin{bmatrix} p_x & p_y & \theta & v \end{bmatrix}^\top$$

with the dynamics:

$$\begin{align} \dot{p}_x &= v \cos\theta & \quad \dot{p}_y &= v \sin\theta \\ \dot{\theta} &= (1 + \eta_\omega) u_\omega + v \xi_\omega & \quad \dot{v} &= (1 + \eta_a) u_a + v \xi_a \end{align}$$

The noise vector \(\mathbf{w}^f= [\xi_a, \xi_\omega, \eta_a, \eta_\omega]^T\) depicts model mismatch and actuator gain error with a Gaussian mixture model as the true noise distribution, from which limited empirical noise samples are available. The constrained inputs are acceleration \(u_a \in [-5.0, 5.0]\) and turn rate \(u_\omega \in [-\frac{\pi}{4}, \frac{\pi}{4}]\). A lower bound on velocity is enforced in the dynamics such that \(v \ge 10.0\), which prevents the vehicle from stopping.

Uncertainty Model

The uncertainty vector \(w=[\xi_a,\xi_w,\eta_a,\eta_w]\) is generated from a two-component Gaussian mixture:

$$w \sim 0.8\,\mathcal{N}\!\left(\mu_{\mathrm{nom}},\operatorname{diag}(\sigma_{\mathrm{nom}}^2)\right) \;+ \;0.2\,\mathcal{N}\!\left(\mu_{\mathrm{out}},\operatorname{diag}(\sigma_{\mathrm{out}}^2)\right),$$

with \(\mu_{\mathrm{nom}}=[0,0,0,0]\), \(\sigma_{\mathrm{nom}}=[0.05,0.02,0.03,0.03]\), \(\mu_{\mathrm{out}}=[-0.15,\,0.08,\,-0.10,\,0.05]\), and \(\sigma_{\mathrm{out}}=[0.02,0.01,0.05,0.05]\). A finite offline dataset (default size \(5000\)) is first sampled i.i.d. from this mixture, and the simulation then draws disturbances by uniform resampling with replacement from that dataset, i.e., from its empirical distribution \(\hat P=\frac{1}{5000}\sum_{i=1}^{5000}\delta_{w^{(i)}}\).

Safety Constraint and CBF

The safety constraint is to avoid two static circular obstacles with centers \((p^i_x, p^i_y)\) and radii \(R_i\), \(i=1,2\). This defines the safe set:

$$\mathcal{S} = \left\{x \in \mathbb{R}^4 \mid h(x) = \min_{i\in\{1,2\}}\left((p_x - p^i_x)^2 + (p_y - p^i_y)^2 - R_i^2\right) > 0\right\}$$

Since the system has relative degree two with respect to the safety constraint, a higher-order CBF (HOCBF) is required. The intermediate barrier function is defined as \(h_e(x) = \dot{h}(x) + \alpha_1 h(x)\) with \(\alpha_1 > 0\). The DR-CBF method enforces the control barrier condition over a Wasserstein ambiguity set \(\mathcal{M}_N^r\) centered on the empirical distribution of offline disturbance samples:

$$\inf_{\mathbb{P} \in \mathcal{M}_N^r} \mathbb{P} \left( \dot{h}_e(x, u, \mathbf{w}^f) + \alpha_2 h_e(x) \ge 0 \right) \ge 1 - \varepsilon$$

Implementation and Control

For implementation, we discretize the continuous-time dynamics with Euler's method with \(\Delta t = 0.05\). The nominal controller for both DRS-gatekeeper and the DR-CBF produces a turn rate command proportional to the heading error to the goal and an acceleration command proportional to the target velocity error. For the backup controller in DRS-gatekeeper, we employ an orbit controller that slows the vehicle to its minimum speed and turns away from the nearest obstacle at the maximum yaw rate, entering an orbit.