Towards Validation of Autonomous Vehicles Across Scales using an Integrated Digital Twin Framework

For small and mid-scale vehicles, which usually implement an electric motor for propulsion, the front/rear/all wheels are driven by applying a torque ${{}^{i}\tau_{drive}}={{}^{i}I_{w}}*{{}^{i}\dot{\omega}_{w}}$, where ${{}^{i}I_{w}}=\frac{1}{2}*{{}^{i}m_{w}}*{{}^{i}{r_{w}}^{2}}$ represents the moment of inertia, ${}^{i}\dot{\omega}_{w}$ is the angular acceleration, ${}^{i}m_{w}$ is the mass, and ${}^{i}r_{w}$ is the radius of the $i$-th wheel. The actuation delays can also be modeled by splitting the torque profile into multiple segments based on operating conditions. For full-scale vehicles, however, the powertrain comprises an engine, transmission and differential. The engine is modeled based on its torque-speed characteristics. The engine RPM is updated smoothly based on its current value $RPM_{e}$, the idle speed $RPM_{i}$, average wheel speed $RPM_{w}$, final drive ratio $FDR$, current gear ratio $GR$, and the vehicle velocity $v$. The update can be expressed as $RPM_{e}:=\left[RPM_{i}+\left(|RPM_{w}|*FDR*GR\right)\right]_{(RPM_{e},v)}$ where, $[\mathscr{F}]_{x}$ denotes evaluation of $\mathscr{F}$ at $x$. The total torque generated by the powertrain is computed as $\tau_{\text{total}}=\left[\tau_{e}\right]_{RPM_{e}}*\left[GR\right]_{G_{\#}}*FDR*\tau*\mathscr{A}$. Here, $\tau_{e}$ is the engine torque, $\tau$ is the throttle input, and $\mathscr{A}$ is a non-linear smoothing operator which increases the vehicle acceleration based on the throttle input. The automatic transmission decides to upshift/downshift the gears based on the transmission map of a given vehicle. This keeps the engine RPM in a good operating range for a given speed: $RPM_{e}=\frac{{v_{\text{MPH}}*5280*12}}{{60*2*\pi*R_{\text{tire}}}}*FDR*GR$. It is to be noted that while shifting the gears, the total torque produced by the powertrain is set to zero to simulate the clutch disengagement. It is also noteworthy that the auto-transmission is put in neutral gear once the vehicle is in standstill condition and parking gear if handbrakes are engaged in standstill condition. Additionally, switching between drive and reverse gears requires that the vehicle first be in the neutral gear to allow this transition. The total torque $\tau_{\text{total}}$ from the drivetrain is divided to the wheels based on the drive configuration of the vehicle: $\tau_{\text{out}}=\begin{cases}\frac{\tau_{\text{total}}}{2}&\text{if FWD/RWD}\\ \frac{\tau_{\text{total}}}{4}&\text{if AWD}\end{cases}$. The torque transmitted to wheels $\tau_{w}$ is modeled by dividing the output torque $\tau_{\text{out}}$ to the left and right wheels based on the steering input. The left wheel receives a torque amounting to ${}^{L}\tau_{w}=\tau_{\text{out}}*(1-\tau_{\text{drop}}*|\delta^{-}|)$, while the right wheel receives a torque equivalent to ${}^{R}\tau_{w}=\tau_{\text{out}}*(1-\tau_{\text{drop}}*|\delta^{+}|)$. Here, $\tau_{\text{drop}}$ is the torque-drop at differential and $\delta^{\pm}$ indicates positive and negative steering angles, respectively. The value of $(\tau_{\text{drop}}*|\delta^{\pm}|)$ is clamped between $[0,0.9]$.

The driving actuators for small and mid-scale vehicles simulate braking torque by applying a holding torque in idle conditions, i.e., ${{}^{i}\tau_{\text{brake}}}={{}^{i}\tau_{\text{idle}}}$. For full-scale vehicles, the braking torque is modeled as ${{}^{i}\tau_{\text{brake}}}=\frac{{{}^{i}M}*v^{2}}{2*D_{\text{brake}}}*R_{b}$ where $R_{b}$ is the brake disk radius and $D_{\text{brake}}$ is the braking distance at 60 MPH, which can be obtained from physical vehicle tests. This braking torque is applied to the wheels based on the type of brake input: for combi-brakes, this torque is applied to all the wheels, and for handbrakes, it is applied to the rear wheels only.

The steering mechanism operates by employing a steering actuator, which applies a torque $\tau_{\text{steer}}$ to achieve the desired steering angle $\delta$ with a smooth rate $\dot{\delta}$, without exceeding the steering limits $\pm\delta_{\text{lim}}$. The rate at which the vehicle steers is governed by its speed $v$ and steering sensitivity $\kappa_{\delta}$, and is represented as $\dot{\delta}=\kappa_{\delta}+\kappa_{v}*\frac{v}{v_{\text{max}}}$. Here, $\kappa_{v}$ is the speed-dependency factor of the steering mechanism. Finally, the individual angle for left $\delta_{l}$ and right $\delta_{r}$ wheels are governed by the Ackermann steering geometry, considering the wheelbase $l$ and track width $w$ of the vehicle: $\left\{\begin{matrix}\delta_{l}=\textup{tan}^{-1}\left(\frac{2*l*\textup{tan}(\delta)}{2*l+w*\textup{tan}(\delta)}\right)\\ \delta_{r}=\textup{tan}^{-1}\left(\frac{2*l*\textup{tan}(\delta)}{2*l-w*\textup{tan}(\delta)}\right)\end{matrix}\right.$.

For small and mid-scale vehicles, the suspension force acting on each sprung mass is calculated as ${{}^{i}M}*{{}^{i}{\ddot{Z}}}+{{}^{i}B}*({{}^{i}{\dot{Z}}}-{{}^{i}{\dot{z}}})+{{}^{i}K}*({{}^{i}{Z}}-{{}^{i}{z}})$, where ${}^{i}Z$ and ${}^{i}z$ denote the displacements of the sprung and unsprung masses, respectively, and ${}^{i}B$ and ${}^{i}K$ represent the damping and spring coefficients of the $i$-th suspension. For full-scale vehicles, however, the stiffness ${{}^{i}K}={{}^{i}M}*{{}^{i}\omega_{n}^{2}}$ and damping ${}^{i}B=2*^{i}\zeta*\sqrt{{{}^{i}K}*{{}^{i}M}}$ coefficients of the suspension system are computed based on the sprung mass ${{}^{i}M}$, natural frequency ${{}^{i}\omega_{n}}$, and damping ratio ${{}^{i}\zeta}$ parameters. The point of suspension force application ${{}^{i}Z_{F}}$ is calculated based on the suspension geometry: ${{}^{i}Z_{F}}={{}^{i}Z_{\text{COM}}}-{{}^{i}Z_{w}}+{{}^{i}r_{w}}-{{}^{i}Z_{f}}$, where ${}^{i}Z_{\text{COM}}$ denotes the Z-component of vehicle’s center of mass, ${}^{i}Z_{w}$ is the Z-component of the relative transformation between each wheel and the vehicle frame (${}^{V}T_{w_{i}}$), ${}^{i}r_{w}$ is the wheel radius, and ${}^{i}Z_{f}$ is the force offset determined by the suspension geometry. Lastly, the suspension displacement ${}^{i}Z_{s}$ at any given moment can be computed as ${{}^{i}Z_{s}}=\frac{{{}^{i}M}*g}{{{}^{i}Z_{0}}*{{}^{i}K}}$, where $g$ represents the acceleration due to gravity, and ${}^{i}Z_{0}$ is the suspension’s equilibrium point. Additionally, full-scale vehicle models also have a provision to include anti-roll bars, which apply a force on the left ${{}^{L}F_{r}}=K_{r}*{{}^{R}Z}-{{}^{L}Z}$ and right ${R^{F}_{r}}=K_{r}*{{}^{L}Z}-{{}^{R}Z}$ wheels as long as they are grounded at the contact point $Z_{c}$. This force is directly proportional to the stiffness of the anti-roll bar, $K_{r}$. The left and right wheel travels are given by ${{}^{L}Z}=\frac{-{{}^{L}Z_{c}}-{{}^{L}r_{w}}}{{}^{L}Z_{s}}$ and ${{}^{R}Z}=\frac{-{{}^{R}Z_{c}}-{{}^{R}r_{w}}}{{}^{R}Z_{s}}$.

Tire forces are determined based on the friction curve for each tire $\left\{\begin{matrix}{{}^{i}F_{t_{x}}}=F(^{i}S_{x})\\ {{}^{i}F_{t_{y}}}=F(^{i}S_{y})\\ \end{matrix}\right.$, where ${}^{i}S_{x}$ and ${}^{i}S_{y}$ represent the longitudinal and lateral slips of the $i$-th tire, respectively. The friction curve is approximated using a two-piece spline, defined as $F(S)=\left\{\begin{matrix}f_{0}(S);\;\;S_{0}\leq S<S_{e}\\ f_{1}(S);\;\;S_{e}\leq S<S_{a}\\ \end{matrix}\right.$, with $f_{k}(S)=a_{k}*S^{3}+b_{k}*S^{2}+c_{k}*S+d_{k}$ as a cubic polynomial function. The first segment of the spline ranges from zero $(S_{0},F_{0})$ to an extremum point $(S_{e},F_{e})$, while the second segment ranges from the extremum point $(S_{e},F_{e})$ to an asymptote point $(S_{a},F_{a})$. Tire slip is influenced by factors including tire stiffness ${}^{i}C_{\alpha}$, steering angle $\delta$, wheel speeds ${}^{i}\omega$, suspension forces ${}^{i}F_{s}$, and rigid-body momentum ${{}^{i}P}={{}^{i}M}*{{}^{i}v}$. The longitudinal slip ${}^{i}S_{x}$ of $i$-th tire is calculated by comparing the longitudinal components of its surface velocity $v_{x}$ (i.e., the longitudinal linear velocity of the vehicle) with its angular velocity ${}^{i}\omega$: ${{}^{i}S_{x}}=\frac{{{}^{i}r}*{{}^{i}\omega}-v_{x}}{v_{x}}$. The lateral slip ${}^{i}S_{y}$ depends on the tire’s slip angle $\alpha$ and is determined by comparing the longitudinal $v_{x}$ (forward velocity) and lateral $v_{y}$ (side-slip velocity) components of the vehicle’s linear velocity: ${{}^{i}S_{y}}=\tan(\alpha)=\frac{v_{y}}{\left|v_{x}\right|}$.

Small and mid-scale vehicles are modeled with constant coefficients for linear $F_{d}$ as well as angular $T_{d}$ drags, which act directly proportional to their linear $v$ and angular $\omega$ velocities. These vehicles do not create significant downforce due to unoptimized aerodynamics, limited velocities and smaller size and mass. Full-scale vehicles, on the other hand, have been modeled to simulate variable air drag $F_{\text{aero}}$ acting on the vehicle, which is computed based on the vehicle’s operating condition: $F_{\text{aero}}=\begin{cases}F_{d_{\text{max}}}&\text{if }v\geq v_{\text{max}}\\ F_{d_{\text{idle}}}&\text{if }\tau_{\text{out}}=0\\ F_{d_{\text{rev}}}&\text{if }(v\geq v_{\text{rev}})\land(G_{\#}=-1)\land(RPM_{w}<0)\\ F_{d_{\text{idle}}}&\text{otherwise}\end{cases}$ where, $v$ is the vehicle velocity, $v_{\text{max}}$ is the vehicle’s designated top-speed, $v_{\text{rev}}$ is the vehicle’s designated maximum reverse velocity, $G_{\#}$ is the operating gear, and $RPM_{w}$ is the average wheel RPM. The downforce acting on a full-scale vehicle is modeled proportional to its velocity: $F_{\text{down}}=K_{\text{down}}*|v|$, where $K_{\text{down}}$ is the downforce coefficient.

Throttle ($\tau$) and steering ($\delta$) sensors are simulated using a simple feedback loop.

Simulated incremental encoders measure wheel rotations ${}^{i}N_{\text{ticks}}={{}^{i}PPR}*{{}^{i}CGR}*{{}^{i}N_{\text{rev}}}$, where ${}^{i}N_{\text{ticks}}$ represents the measured ticks, ${}^{i}PPR$ is the encoder resolution (pulses per revolution), ${}^{i}CGR$ is the cumulative gear ratio, and ${}^{i}N_{\text{rev}}$ represents the wheel revolutions.

Positioning systems and inertial measurement units (IMU) are simulated based on temporally coherent rigid-body transform updates of the vehicle $\{v\}$ with respect to the world $\{w\}$: ${{}^{w}\mathbf{T}_{v}}=\left[\begin{array}[]{c | c}\mathbf{R}_{3\times 3}&\mathbf{t}_{3\times 1}\\ \hline\cr\mathbf{0}_{1\times 3}&1\end{array}\right]\in SE(3)$. The positioning systems provide 3-DOF positional coordinates $\{x,y,z\}$ of the vehicle, while the IMU supplies linear accelerations $\{a_{x},a_{y},a_{z}\}$, angular velocities $\{\omega_{x},\omega_{y},\omega_{z}\}$, and 3-DOF orientation data for the vehicle, either as Euler angles $\{\phi_{x},\theta_{y},\psi_{z}\}$ or as a quaternion $\{q_{0},q_{1},q_{2},q_{3}\}$.

2D LIDAR simulation employs iterative ray-casting raycast{${}^{w}\mathbf{T}_{l}$, $\vec{\mathbf{R}}$, $r_{\text{max}}$} for each angle $\theta\in\left[\theta_{\text{min}}:\theta_{\text{res}}:\theta_{\text{max}}\right]$ at a specified update rate. Here, ${{}^{w}\mathbf{T}_{l}}={{}^{w}\mathbf{T}_{v}}*{{}^{v}\mathbf{T}_{l}}\in SE(3)$ represents the relative transformation of the LIDAR {$l$} with respect to the vehicle {$v$} and the world {$w$}, $\vec{\mathbf{R}}=\left[\cos(\theta)\;\;\sin(\theta)\;\;0\right]^{T}$ defines the direction vector of each ray-cast $R$, where $r_{\text{min}}$ and $r_{\text{max}}$ denote the minimum and maximum linear ranges, $\theta_{\text{min}}$ and $\theta_{\text{max}}$ denote the minimum and maximum angular ranges, and $\theta_{\text{res}}$ represents the angular resolution of the LIDAR, respectively. The laser scan ranges are determined by checking ray-cast hits and then applying a threshold to the minimum linear range of the LIDAR, calculated as ranges[i]$=\begin{cases}\texttt{hit.dist}&\text{ if }\texttt{ray[i].hit}\text{ and }\texttt{hit.dist}\geq r_{\text{min}}\\ \infty&\text{ otherwise}\end{cases}$, where ray.hit is a Boolean flag indicating whether a ray-cast hits any colliders in the scene, and hit.dist$=\sqrt{(x_{\text{hit}}-x_{\text{ray}})^{2}+(y_{\text{hit}}-y_{\text{ray}})^{2}+(z_{\text{hit}}-z_{\text{ray}})^{2}}$ calculates the Euclidean distance from the ray-cast source $\{x_{\text{ray}},y_{\text{ray}},z_{\text{ray}}\}$ to the hit point $\{x_{\text{hit}},y_{\text{hit}},z_{\text{hit}}\}$.

3D LIDAR simulation adopts multi-channel parallel ray-casting raycast{${}^{w}\mathbf{T}_{l}$, $\vec{\mathbf{R}}$, $r_{\text{max}}$} for each angle $\theta\in\left[\theta_{\text{min}}:\theta_{\text{res}}:\theta_{\text{max}}\right]$ and each channel $\phi\in\left[\phi_{\text{min}}:\phi_{\text{res}}:\phi_{\text{max}}\right]$ at a specified update rate, with GPU acceleration (if available). Here, ${{}^{w}\mathbf{T}_{l}}={{}^{w}\mathbf{T}_{v}}*{{}^{v}\mathbf{T}_{l}}\in SE(3)$ represents the relative transformation of the LIDAR {$l$} with respect to the vehicle {$v$} and the world {$w$}, $\vec{\mathbf{R}}=\left[\cos(\theta)*\cos(\phi)\;\;\sin(\theta)*\cos(\phi)\;\;-\sin(\phi)\right]^{T}$ defines the direction vector of each ray-cast $R$, where $r_{\text{min}}$ and $r_{\text{max}}$ denote the minimum and maximum linear ranges, $\theta_{\text{min}}$ and $\theta_{\text{max}}$ denote the minimum and maximum horizontal angular ranges, $\phi_{\text{min}}$ and $\phi_{\text{max}}$ denote the minimum and maximum vertical angular ranges, and $\theta_{\text{res}}$ and $\phi_{\text{res}}$ represent the horizontal and vertical angular resolutions of the LIDAR, respectively. The thresholded ray-cast hit coordinates $\{x_{\text{hit}},y_{\text{hit}},z_{\text{hit}}\}$, from each of the casted rays is encoded into byte arrays based on the LIDAR parameters, and given out as the point cloud data.

Simulated cameras are parameterized by their focal length $f$, sensor size $\{s_{x},s_{y}\}$, target resolution, as well as the distances to the near $N$ and far $F$ clipping planes. The viewport rendering pipeline for the simulated cameras operates in three stages. First, the camera view matrix $\mathbf{V}\in SE(3)$ is computed by obtaining the relative homogeneous transform of the camera $\{c\}$ with respect to the world $\{w\}$: $\mathbf{V}=\begin{bmatrix}r_{00}&r_{01}&r_{02}&t_{0}\\ r_{10}&r_{11}&r_{12}&t_{1}\\ r_{20}&r_{21}&r_{22}&t_{2}\\ 0&0&0&1\\ \end{bmatrix}$, where $r_{ij}$ and $t_{i}$ denote the rotational and translational components, respectively. Next, the camera projection matrix $\mathbf{P}\in\mathbb{R}^{4\times 4}$ is calculated to project world coordinates into image space coordinates: $\mathbf{P}=\begin{bmatrix}\frac{2*N}{R-L}&0&\frac{R+L}{R-L}&0\\ 0&\frac{2*N}{T-B}&\frac{T+B}{T-B}&0\\ 0&0&-\frac{F+N}{F-N}&-\frac{2*F*N}{F-N}\\ 0&0&-1&0\\ \end{bmatrix}$, where $L$, $R$, $T$, and $B$ denote the left, right, top, and bottom offsets of the sensor. The camera parameters $\{f,s_{x},s_{y}\}$ are related to the terms of the projection matrix as follows: $f=\frac{2*N}{R-L}$, $a=\frac{s_{y}}{s_{x}}$, and $\frac{f}{a}=\frac{2*N}{T-B}$. The perspective projection from the simulated camera’s viewport is given as $\mathbf{C}=\mathbf{P}*\mathbf{V}*\mathbf{W}$, where $\mathbf{C}=\left[x_{c}\;\;y_{c}\;\;z_{c}\;\;w_{c}\right]^{T}$ represents image space coordinates, and $\mathbf{W}=\left[x_{w}\;\;y_{w}\;\;z_{w}\;\;w_{w}\right]^{T}$ represents world coordinates. Finally, this camera projection is transformed into normalized device coordinates (NDC) by performing perspective division (i.e., dividing throughout by $w_{c}$), leading to a viewport projection achieved by scaling and shifting the result and then utilizing the rasterization process of the graphics API (e.g., DirectX for Windows, Metal for macOS, and Vulkan for Linux). Additionally, a post-processing step simulates non-linear lens and film effects, such as lens distortion, depth of field, exposure, ambient occlusion, contact shadows, bloom, motion blur, film grain, chromatic aberration, etc.