Physics-Based Path Generation for Synthetic Trajectories

Physics-based path generation constructs spatially and temporally coherent movement traces by integrating equations of motion under explicit terrain, mechanical, and operational constraints, so every synthetic point is one a real vehicle could actually have occupied. Part of Trajectory & Movement Simulation, this page covers the specific sub-problem of turning a coarse route — a set of waypoints, or the output of a stochastic router — into a dynamically feasible trajectory that honours velocity, acceleration, jerk, and curvature limits while interacting correctly with the surface beneath it. Where the parent area frames the fidelity-versus-privacy trade-off across every movement primitive, here the contract is the continuous state vector and the laws of motion that govern how it evolves; every downstream gate tests whether a generated trace ever violates them.

Problem Framing

The failure that physics-based generation exists to prevent is the geometrically smooth but physically impossible trajectory. Two engineering failures dominate.

The first is kinematic implausibility leaking into training data. A spline fitted through waypoints, or a Markov-chain routing model that samples next-hops with no dynamics, will happily produce a path that turns a 90-degree corner at 80 km/h, accelerates from rest to highway speed in one timestep, or climbs a 30-degree grade without losing velocity. The geometry looks plausible on a basemap, but an autonomy stack, a fuel-consumption model, or a sensor simulator trained on it learns a motion prior that does not exist in the physical world. The model then fails the moment it meets a real vehicle that cannot corner that fast. Enforcing $v_{max}$, $a_{max}$, $j_{max}$, and a curvature ceiling $\kappa \le 1/R_{min}$ at every integration step is what removes those impossible states before they ever reach a dataset.

The second is silent non-reproducibility of the dynamics. A physics integrator is deterministic only if its inputs are pinned: the solver, the step size, the floating-point reduction order, the terrain raster, and any random seed used for stochastic perturbation must all be fixed and recorded, or two runs of the “same” pipeline diverge after a few hundred steps and produce traces that no audit can reconstruct. Treating the seed, the integrator configuration, the coordinate reference system version, and a hash of the digital elevation model as versioned, first-class artifacts is the discipline that closes both gaps — and it rests on the same CRS contract enforcement that fixes the spatial envelope before any motion is integrated.

Prerequisites & Toolchain

Physics-based generation sits on the standard geospatial Python stack plus a numerical integration backend and a fast raster sampler for terrain. Pin to the major versions used across this site so that integration order, RNG streams, and raster interpolation behave identically across environments.

pythonCopy
# requirements.txt — pin majors; let patch releases float
numpy==1.26.*
scipy==1.11.*           # solve_ivp ODE integrators, KD-trees, sparse QP
shapely==2.0.*          # GEOS-backed geometry for path + exclusion-zone tests
pyproj==3.6.*           # CRS transforms (PROJ 9.x under the hood)
rasterio==1.3.*         # DEM / friction-map sampling with affine transforms
geopandas==0.14.*       # GeoDataFrame I/O for the route + constraint layers

Two environment variables must be set explicitly rather than inherited, or you will hit non-deterministic behaviour that masquerades as a dynamics bug:

bashCopy
# Resolve PROJ data deterministically; never rely on the ambient default
export PROJ_LIB="$(python -c 'import pyproj; print(pyproj.datadir.get_data_dir())')"
export PYTHONHASHSEED=0   # stabilize any set/dict iteration feeding node order

All integration happens in a projected, metric CRS, never in geographic degrees. Velocity, acceleration, curvature, and slope are all distance-sensitive, so integrating in EPSG:4326 distorts every quantity by latitude and makes a turning-radius limit meaningless. Reproject origin and destination from EPSG:4326 (or EPSG:3857) to a local UTM zone or a custom engineering grid at ingestion, integrate there, and reproject the assembled trace back to EPSG:4326 only at serialization. The same terrain layers — a DEM for elevation gradients, a friction raster for surface traction, and the road-network topology for lane constraints — must be version-locked and reprojected to the same working CRS so that a sampled slope at a point matches the position the solver thinks it is at.

Core Concept: Equations of Motion Under Hard Constraints

A physics-based path generator advances a continuous state vector through time:

$$\mathbf{x}(t) = [\,\mathbf{p},\ \mathbf{v},\ \mathbf{a},\ \theta,\ \omega\,]^\top,$$

where $\mathbf{p}$ is position, $\mathbf{v}$ velocity, $\mathbf{a}$ acceleration, $\theta$ heading, and $\omega = \dot{\theta}$ the yaw rate. The state evolves under second-order equations of motion driven by a control input $\mathbf{u}$ (commanded longitudinal and lateral acceleration), and the synthesis engine enforces hard constraints at every step so the trace stays physically realisable.

Kinematic limits cap the motion profile. Maximum velocity $v_{max}$, acceleration $a_{max}$, deceleration $d_{max}$, and jerk $j_{max} = \lVert \dot{\mathbf{a}} \rVert$ thresholds prevent the discontinuous, infinite-jerk motion that betrays a synthetic origin. A finite jerk bound is what makes a generated stop look like braking rather than teleportation to zero velocity.

Friction and slope modulate available traction. The DEM yields a local grade angle $\phi$, and the effective acceleration vector is projected onto the tangent plane of the terrain, so a climb sheds velocity and a descent gains it exactly as real mechanics dictate. The friction coefficient $\mu$ sampled from the surface raster bounds the total acceleration the tyre contact patch can deliver — the friction-circle constraint $\sqrt{a_{lon}^2 + a_{lat}^2} \le \mu g$ couples longitudinal and lateral demand so that hard cornering leaves less grip for braking, just as it does on a real road.

Turning radius and curvature bound how sharply the path can bend. Path curvature $\kappa$ must satisfy

$$\kappa \le \frac{1}{R_{min}}, \qquad R_{min} = \frac{v^2}{a_{lat,\max}},$$

where $R_{min}$ follows from vehicle geometry, steering-actuation limits, and the lateral-acceleration ceiling. Because $R_{min}$ grows with $v^2$, the same steering input is feasible at low speed and infeasible at high speed — the constraint is velocity-coupled, and a generator that ignores that coupling produces the classic impossible high-speed hairpin.

Temporal discretization fixes the integration grid. A fixed step $\Delta t \in [0.01, 0.1]\,\text{s}$ keeps the integrator numerically stable and the replay deterministic. Adaptive step sizing is permitted only when paired with monotonic state interpolation back onto a fixed output grid, so timestamps stay ordered and sensor simulation does not alias. This same fixed, monotonic clock is what lets temporal synchronization for moving objects align several agents on a shared timeline downstream.

When a step would breach a constraint, the solver does not simply clip the output — it applies a projection step, finding the nearest feasible state (via a small quadratic program) or correcting the control input with a penalty term so that velocity, acceleration, and curvature stay mutually consistent and temporal continuity is never broken.

Step-by-Step Implementation

The pipeline runs in four ordered stages, each exposing deterministic hooks for validation and replay.

1. Spatial anchoring and coordinate alignment

Normalize all inputs into one metric working CRS, then attach the terrain samplers the dynamics will query. Transforming with pyproj keeps the origin and destination sub-metre accurate before any integration begins.

pythonCopy
import numpy as np
import rasterio
import pyproj

WORK_CRS = "EPSG:32633"  # UTM 33N — pick the zone covering your AOI

def to_working(lon: float, lat: float) -> tuple[float, float]:
    """Project EPSG:4326 lon/lat into the metric integration CRS."""
    tf = pyproj.Transformer.from_crs("EPSG:4326", WORK_CRS, always_xy=True)
    return tf.transform(lon, lat)

class Terrain:
    """Version-locked DEM + friction sampler, both in WORK_CRS."""
    def __init__(self, dem_path: str, friction_path: str):
        self.dem = rasterio.open(dem_path)
        self.fric = rasterio.open(friction_path)
        assert self.dem.crs.to_string() == WORK_CRS, "DEM not in working CRS"

    def slope(self, x: float, y: float) -> float:
        """Local grade angle phi (rad) from a 3x3 DEM gradient."""
        win = np.array(list(self.dem.sample([(x, y)])))[0, 0]
        # gradient computed from a cached tiled array in production
        return float(np.arctan(win))  # placeholder for precomputed gradient

    def mu(self, x: float, y: float) -> float:
        """Surface friction coefficient at a point."""
        return float(np.array(list(self.fric.sample([(x, y)])))[0, 0])

2. State initialization and boundary setup

Define the initial state $\mathbf{x}_0$, the terminal tolerances, and the spatial constraints. Exclusion zones (no-go areas and privacy buffers around sensitive infrastructure) are encoded before integration so a non-compliant state can never enter the trace — geofencing is a feasible-region constraint, not a post-hoc filter.

pythonCopy
from dataclasses import dataclass
from shapely.geometry import Point, Polygon

@dataclass(frozen=True)
class Limits:
    v_max: float = 27.0      # m/s  (~97 km/h)
    a_max: float = 3.0       # m/s^2
    d_max: float = 6.0       # m/s^2 braking
    j_max: float = 4.0       # m/s^3 jerk
    a_lat_max: float = 4.0   # m/s^2 lateral
    g: float = 9.81

@dataclass(frozen=True)
class State:
    x: float; y: float
    vx: float; vy: float
    theta: float; omega: float

def feasible(p: Point, exclusion: list[Polygon]) -> bool:
    """Reject any state inside a privacy / no-go polygon."""
    return not any(zone.contains(p) for zone in exclusion)

3. Numerical integration with constraint enforcement

The core solver advances the state with a fixed-step RK4 integrator — preferred over adaptive schemes here because bitwise-reproducible replay matters more than per-step error adaptivity. After each step, the candidate state is projected back into the feasible region defined by the kinematic limits, the friction circle, and the curvature ceiling.

pythonCopy
def derivative(s: State, u: np.ndarray, terr: Terrain, lim: Limits) -> np.ndarray:
    """State derivative under slope-projected, friction-bounded dynamics."""
    phi = terr.slope(s.x, s.y)                 # local grade (rad)
    mu = terr.mu(s.x, s.y)                      # surface friction
    a_lon, a_lat = u                            # commanded controls

    # friction-circle: clamp combined demand to mu * g
    demand = np.hypot(a_lon, a_lat)
    cap = mu * lim.g
    if demand > cap:
        a_lon *= cap / demand
        a_lat *= cap / demand

    a_lon -= lim.g * np.sin(phi)               # slope robs/adds longitudinal accel
    speed = np.hypot(s.vx, s.vy)
    omega = (a_lat / speed) if speed > 1e-3 else 0.0   # yaw rate from lateral accel
    ax = a_lon * np.cos(s.theta) - a_lat * np.sin(s.theta)
    ay = a_lon * np.sin(s.theta) + a_lat * np.cos(s.theta)
    return np.array([s.vx, s.vy, ax, ay, omega, 0.0])

def rk4_step(s: State, u, terr, lim, dt: float) -> State:
    """Deterministic fixed-step Runge-Kutta 4 advance."""
    y = np.array([s.x, s.y, s.vx, s.vy, s.theta, s.omega])
    f = lambda st: derivative(State(*st), u, terr, lim)
    k1 = f(y)
    k2 = f(y + 0.5 * dt * k1)
    k3 = f(y + 0.5 * dt * k2)
    k4 = f(y + dt * k3)
    y_next = y + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
    return project(State(*y_next), lim)        # enforce hard constraints

def project(s: State, lim: Limits) -> State:
    """Nearest-feasible projection: clamp speed and curvature."""
    speed = np.hypot(s.vx, s.vy)
    if speed > lim.v_max:                       # velocity ceiling
        scale = lim.v_max / speed
        s = State(s.x, s.y, s.vx * scale, s.vy * scale, s.theta, s.omega)
        speed = lim.v_max
    r_min = (speed ** 2) / lim.a_lat_max if speed > 1e-3 else 0.0
    kappa = abs(s.omega) / speed if speed > 1e-3 else 0.0
    if r_min > 0 and kappa > 1.0 / r_min:       # curvature ceiling
        s = State(s.x, s.y, s.vx, s.vy, s.theta,
                  np.sign(s.omega) * speed / r_min)
    return s

Every solver parameter — the step size, the limit set, any perturbation seed, and the environment hashes — is logged at this stage so the run can be reconstructed exactly.

4. Trace assembly, hashing, and handoff

Concatenate the per-step states into a trace, reproject to EPSG:4326, and stamp it with a content hash derived from the seed, CRS version, DEM hash, and solver config. The clean, deterministic trace is then handed to noise injection and stochastic drift to add realistic sensor error without contaminating the physical baseline.

pythonCopy
import hashlib, json

def trace_hash(seed: int, dem_sha: str, lim: Limits, dt: float) -> str:
    """Content-addressable id for exact replay and audit."""
    payload = json.dumps(
        {"seed": seed, "crs": WORK_CRS, "dem": dem_sha,
         "limits": lim.__dict__, "dt": dt},
        sort_keys=True,
    ).encode()
    return hashlib.sha256(payload).hexdigest()

Validation & Testing

Post-integration validation asserts that no physical law was broken anywhere in the trace. These checks belong in CI/CD integration for spatial data so a build that emits an infeasible trajectory is rejected before promotion.

pythonCopy
def validate(trace: np.ndarray, lim: Limits, dt: float) -> None:
    """trace columns: x, y, vx, vy, t. Raise on any physical violation."""
    pos = trace[:, :2]
    vel = trace[:, 2:4]
    t = trace[:, 4]

    speed = np.linalg.norm(vel, axis=1)
    assert speed.max() <= lim.v_max + 1e-6, "velocity bound exceeded"

    accel = np.diff(vel, axis=0) / dt
    a_mag = np.linalg.norm(accel, axis=1)
    assert a_mag.max() <= lim.a_max + 1e-6, "acceleration bound exceeded"

    jerk = np.linalg.norm(np.diff(accel, axis=0) / dt, axis=1)
    assert jerk.max() <= lim.j_max + 1e-6, "jerk bound exceeded"

    assert np.all(np.diff(t) > 0), "timestamps not strictly monotonic"

Three assertion families matter. Bound checks confirm $|\mathbf{v}| \le v_{max}$, $|\mathbf{a}| \le a_{max}$, and $j \le j_{max}$ across every step. Curvature continuity confirms $\Delta\theta/\Delta t$ stays within actuator limits and that no step demands a turn tighter than $R_{min}$ at its speed. Spatial fidelity cross-references the trace against the DEM and road-network layers to catch off-road drift or lane violations, and against the exclusion polygons to prove no point intersects a restricted zone. Beyond hard bounds, statistical realism — that the synthetic speed and acceleration distributions match a reference fleet — is scored with the distance metrics described in realism metrics evaluation.

Performance & Scale Considerations

Generating a fleet of long trajectories is dominated by two costs: terrain sampling and the per-step Python overhead of the integrator.

Precompute and tile the terrain. Sampling a DEM or friction raster one point at a time through rasterio is the single biggest bottleneck. Precompute the slope and friction fields into tiled, memory-mapped NumPy arrays keyed by the affine transform, and resolve a point to a tile-local index with integer arithmetic. A memory-mapped array lets many worker processes share one read-only copy of a continental DEM without duplicating it per process.

Parallelize across agents, not within a trajectory. Each agent’s integration is sequential and cheap; the embarrassingly parallel axis is the agent population. Derive each agent’s seed deterministically from one base seed so independent workers stay reproducible, then fan out:

pythonCopy
from concurrent.futures import ProcessPoolExecutor

def child_seed(base: int, agent_id: int) -> int:
    return int.from_bytes(
        hashlib.sha256(f"{base}:{agent_id}".encode()).digest()[:8], "big")

def run_fleet(base_seed: int, n_agents: int, gen) -> list:
    seeds = [child_seed(base_seed, i) for i in range(n_agents)]
    with ProcessPoolExecutor() as pool:
        return list(pool.map(gen, seeds))   # gen(seed) -> one trace

Cache intermediate states for replay. QA cycles re-integrate the same trajectories repeatedly. Persisting the per-step state to content-addressable storage keyed by the trace hash lets validation replay a trace from cache instead of recomputing it, and guarantees the replayed trace is byte-identical to the audited one. Keep vx/vy in float64; downcasting to float32 mid-integration changes the reduction order and breaks bitwise reproducibility across machines.

Failure Modes & Troubleshooting

Integrating in geographic degrees

Symptom: turning-radius and slope limits behave inconsistently across latitudes, and a path that is feasible near the equator is flagged infeasible at high latitude. Root cause: running the solver in EPSG:4326, where a degree is not a constant ground distance, so $v^2/R_{min}$ and the friction circle are computed on distorted units. Fix: reproject to a local UTM or equal-area CRS at ingestion, integrate there, and convert back to EPSG:4326 only at export.

Impossible high-speed cornering

Symptom: generated traces take sharp turns at highway speed that no vehicle could hold. Root cause: a curvature limit applied with a constant $R_{min}$ instead of the velocity-coupled $R_{min} = v^2/a_{lat,\max}$, so the bound is too loose at speed. Fix: recompute $R_{min}$ from the current speed at every step inside the projection, and add a friction-circle clamp so cornering demand and braking demand compete for the same grip.

Step-size instability over rough terrain

Symptom: velocity oscillates or blows up where the DEM has a discontinuity such as a cliff edge or a bridge deck. Root cause: a fixed RK4 step too large for a stiff slope transition. Fix: cap $\Delta t$ at 0.05 s for off-road or steep terrain, or switch to an implicit solver for the stiff segment, then resample the result back onto the fixed output grid to preserve monotonic timestamps.

Non-reproducible traces across hosts

Symptom: identical configs yield diverging trajectories on different machines after a few hundred steps. Root cause: an unseeded perturbation RNG, mixed float32/float64 arithmetic changing reduction order, or set/dict iteration feeding control order. Fix: propagate one seed into every stochastic call, keep the state in float64 throughout, set PYTHONHASHSEED=0, pin the toolchain, and record the trace hash (seed + CRS + DEM hash + solver config) in a replay registry.

Privacy leakage through over-faithful dynamics

Symptom: adversarial linkage reconstructs a real driver’s route from a synthetic trace because its acceleration signature is too distinctive. Root cause: dynamics calibrated so tightly to one real journey that the trace re-exposes an identifiable behavioural fingerprint. Fix: generate entirely outside real-world PII boundaries, apply differential privacy mechanisms to any calibration drawn from real telemetry, and track the composed $(\varepsilon, \delta)$ budget so the behavioural prior cannot be inverted back to an individual.

Frequently Asked Questions

Related

• Markov Chain Routing Models — produce the network-constrained hop sequence that physics-based generation then makes physically traversable.
• Noise Injection & Stochastic Drift — add realistic sensor error to the clean physics trace, decoupled from the dynamics layer.
• Temporal Synchronization for Moving Objects — align the fixed-step clocks of many physically generated agents onto a shared, replayable timeline.
• Adding Realistic GPS Noise to Synthetic Vehicle Trajectories — a worked drill-down on perturbing a clean trace without breaking its dynamics.
• Scoping Rules & Data Contracts — the CRS and spatial-envelope contract that fixes the working projection before any motion is integrated.