Approximate peak-grip ellipses
diff --git a/app/static/visual_panel.js b/app/static/visual_panel.js index 671e21a3..8ad401c4 100644 --- a/app/static/visual_panel.js +++ b/app/static/visual_panel.js @@ -247,11 +247,12 @@ ctx.stroke(); const airborne = state.scene.scalar(corner.slot_airborne, f) > 0.5; + const unavailable = corner.slot_grip_unavailable !== undefined && state.scene.scalar(corner.slot_grip_unavailable, f) > 0.5; const usage = state.scene.scalar(corner.slot_usage, f); const gx = state.scene.scalar(corner.slot_x, f); const gy = state.scene.scalar(corner.slot_y, f); - if (!airborne) { + if (!airborne && !unavailable) { // Lateral is the horizontal axis: it is the one a driver feels. const px = cx + gy * radius / extent; const py = cy - gx * radius / extent; @@ -266,7 +267,7 @@ ctx.textAlign = "center"; ctx.fillText(corner.corner.toUpperCase(), cx, cy - radius - 10 * dpr); ctx.fillText( - airborne ? "off the ground" : Math.round(usage * 100) + "% of peak", + airborne ? "off the ground" : unavailable ? "grip unavailable" : Math.round(usage * 100) + "% of peak", cx, cy + radius + 18 * dpr ); diff --git a/app/visual.py b/app/visual.py index 6fbe307c..27a5a99d 100644 --- a/app/visual.py +++ b/app/visual.py @@ -33,6 +33,7 @@ import numpy as np from visualization.sim_data import SimData, resolve_load_reference +from engines.dynpy.tire import CONTACT_LOAD_N from visualization.tire_state import axle_load_transfer, grip_usage VISUAL_RESULTS_ROOT = Path("visualization/results") @@ -290,8 +291,7 @@ def _grip_payload(data: SimData, scalar_slot: Any) -> dict[str, Any]: "corner": str(corner), "axle": str(entry.get("axle", "front")), "tir": str(coeffs.get("tir", "")), - # NaN means "off the ground", which JSON cannot carry and the - # renderer must not draw as zero, so it travels as a flag. + # Invalid grip travels as a flag because scalar slots replace nonfinite values. # ``scalar_slot`` zeroes the NaN itself; calling nan_to_num # here first would turn an infinity into 1.8e308, which is # finite enough to slip past that guard and then overflows @@ -299,7 +299,8 @@ def _grip_payload(data: SimData, scalar_slot: Any) -> dict[str, Any]: "slot_x": scalar_slot(longitudinal), "slot_y": scalar_slot(lateral), "slot_usage": scalar_slot(usage), - "slot_airborne": scalar_slot(np.isnan(usage).astype(float)), + "slot_airborne": scalar_slot((fz <= CONTACT_LOAD_N).astype(float)), + "slot_grip_unavailable": scalar_slot((~np.isfinite(usage)).astype(float)), } ) if not corners: diff --git a/engines/dynpy/tire.py b/engines/dynpy/tire.py index 1bccb085..d7ef806b 100644 --- a/engines/dynpy/tire.py +++ b/engines/dynpy/tire.py @@ -33,6 +33,24 @@ MF52_COEFFICIENTS = FORCE_COEFFICIENTS + MOMENT_COEFFICIENTS +CONTACT_LOAD_N = 1e-3 + + +def peak_friction(fz_n: ArrayLike, gamma_rad: ArrayLike, p: Mapping[str, float]) -> tuple[FloatArray, FloatArray]: + """## MF5.2 Peak Coefficients + + Return signed pure-slip peak coefficients at the fitted evaluation load. + Below FZMIN, force scales with actual load while these coefficients stay fixed. + """ + fz = np.maximum(np.asarray(fz_n, dtype=float), p.get("FZMIN", CONTACT_LOAD_N)) + gamma = np.asarray(gamma_rad, dtype=float) + nominal = p["FNOMIN"] * p.get("LFZO", 1.) + dfz = (fz - nominal) / nominal + mux = (p["PDX1"] + p["PDX2"] * dfz) * (1 - p.get("PDX3", 0.) * (gamma * p.get("LGAX", 1.)) ** 2) + muy = (p["PDY1"] + p["PDY2"] * dfz) * (1 - p.get("PDY3", 0.) * (gamma * p.get("LGAY", 1.)) ** 2) + return mux * p.get("LMUX", 1.), muy * p.get("LMUY", 1.) + + @dataclass(frozen=True) class MF52Tire: """## MF5.2 Tire @@ -130,7 +148,7 @@ def slip_for_force( curves = self.forces(fz[..., None], alpha[..., None], grid, gamma[..., None])[0] result = np.zeros(fz.shape) for index in np.ndindex(fz.shape): - if fz[index] <= 1e-3: + if fz[index] <= CONTACT_LOAD_N: continue curve = curves[index] wanted = float(target[index]) @@ -193,7 +211,7 @@ def evaluate( """ load, alpha, slip, gamma, speed = np.broadcast_arrays(fz_n, alpha_rad, kappa, camber_rad, speed_mps) fz = np.maximum(load, self.fz_min_n) - scale = np.where(load > 1e-3, load / fz, 0.0) + scale = np.where(load > CONTACT_LOAD_N, load / fz, 0.0) fx, fy = self.forces(fz, alpha, slip, gamma) mx, my, mz = self._moments_from_forces(fz, alpha, slip, gamma, speed, fx, fy) return scale * fx, scale * fy, scale * mx, -scale * my, -scale * mz @@ -219,7 +237,7 @@ def rolling_moment(self, fz_n: ArrayLike, fx_n: ArrayLike, speed_mps: ArrayLike) """ load, fx, speed = np.broadcast_arrays(fz_n, fx_n, speed_mps) fz = np.maximum(load, self.fz_min_n) - scale = np.where(load > 1e-3, load / fz, 0.0) + scale = np.where(load > CONTACT_LOAD_N, load / fz, 0.0) raw_fx = np.divide(fx, scale, out=np.zeros_like(fz, dtype=float), where=scale != 0) return -scale * self._my_raw(fz, raw_fx, speed) @@ -357,11 +375,10 @@ def forces( """ load, alpha, slip, gamma = np.broadcast_arrays(fz_n, alpha_rad, kappa, camber_rad) fz = np.maximum(load, self.fz_min_n) - scale = np.where(load > 1e-3, load / fz, 0.0) + scale = np.where(load > CONTACT_LOAD_N, load / fz, 0.0) p = self.coefficients dfz = (fz - p["FNOMIN"] * p["LFZO"]) / (p["FNOMIN"] * p["LFZO"]) - mux = (p["PDX1"] + p["PDX2"] * dfz) * (1 - p["PDX3"] * (gamma * p["LGAX"]) ** 2) * p["LMUX"] - muy = (p["PDY1"] + p["PDY2"] * dfz) * (1 - p["PDY3"] * (gamma * p["LGAY"]) ** 2) * p["LMUY"] + mux, muy = peak_friction(fz, gamma, p) cx = p["PCX1"] * p["LCX"] cy = p["PCY1"] * p["LCY"] dx, dy = mux * fz, muy * fz diff --git a/tests/test_visual_tire_state.py b/tests/test_visual_tire_state.py index a9febf80..fe226b09 100644 --- a/tests/test_visual_tire_state.py +++ b/tests/test_visual_tire_state.py @@ -116,3 +116,14 @@ def test_the_bundled_tire_file_gives_plausible_fsae_grip() -> None: mu_x, mu_y = peak_mu(np.array([400.0, 650.0, 1000.0]), 0.0, friction_coefficients(tir)) assert np.all((mu_x > 1.0) & (mu_x < 4.0)) assert np.all((mu_y > 1.0) & (mu_y < 4.0)) + + +def test_grip_has_no_friction_floor_and_keeps_lightly_loaded_contact(): + coeffs = dict(COEFFS, LMUX=.001, LMUY=.001, FZMIN=100.) + mu = peak_mu(5., 0., coeffs) + np.testing.assert_allclose(mu, peak_mu(100., 0., coeffs)) + assert 0. < mu[0] < .05 + _, _, usage = grip_usage(np.array([mu[0]*5.]), np.zeros(1), np.array([5.]), 0., coeffs) + assert usage[0] == pytest.approx(1.) + _, _, zero = grip_usage(np.ones(1), np.ones(1), np.array([500.]), 0., dict(coeffs, LMUX=0.)) + assert np.isnan(zero[0]) diff --git a/visualization/tire_state.py b/visualization/tire_state.py index 14c75f61..daf5c0ed 100644 --- a/visualization/tire_state.py +++ b/visualization/tire_state.py @@ -1,6 +1,6 @@ """Tire grip use and lateral load transfer for BobVis. -Numpy only, like :mod:`visualization.navigation`, so it tests without the +Shares peak coefficients and contact conventions with DynPy without the rendering stack. Two things are derived here from signals a run already has: Friction circles @@ -29,15 +29,15 @@ from typing import Mapping import numpy as np +from engines.dynpy.tire import CONTACT_LOAD_N, peak_friction REQUIRED_TIR_KEYS = ("FNOMIN", "PDX1", "PDX2", "PDY1", "PDY2") -OPTIONAL_TIR_KEYS = {"LFZO": 1.0, "LMUX": 1.0, "LMUY": 1.0, "PDX3": 0.0, "PDY3": 0.0} +OPTIONAL_TIR_KEYS = { + "LFZO": 1.0, "LMUX": 1.0, "LMUY": 1.0, "PDX3": 0.0, "PDY3": 0.0, + "LGAX": 1.0, "LGAY": 1.0, "FZMIN": CONTACT_LOAD_N, +} -MU_FLOOR = 0.05 -"""Smallest peak friction used, so extreme load or camber never divides by ~0.""" - -MIN_LOAD_N = 10.0 -"""Below this a tire is treated as off the ground: it has no grip to use.""" +MIN_LOAD_N = CONTACT_LOAD_N LLTD_MIN_FRACTION = 0.01 """LLTD is undefined until total transfer reaches this share of the car's weight.""" @@ -61,17 +61,8 @@ def peak_mu( fz: np.ndarray | float, gamma: np.ndarray | float, coeffs: Mapping[str, float] ) -> tuple[np.ndarray, np.ndarray]: """Longitudinal and lateral peak friction at a normal load and camber (radians).""" - fz = np.asarray(fz, dtype=float) - gamma_sq = np.asarray(gamma, dtype=float) ** 2 - fz0 = float(coeffs["FNOMIN"]) * float(coeffs.get("LFZO", 1.0)) - dfz = (fz - fz0) / fz0 - - mu_x = (coeffs["PDX1"] + coeffs["PDX2"] * dfz) * (1.0 - coeffs.get("PDX3", 0.0) * gamma_sq) - mu_y = (coeffs["PDY1"] + coeffs["PDY2"] * dfz) * (1.0 - coeffs.get("PDY3", 0.0) * gamma_sq) - # PDY1 is negative in the ISO sign convention; grip is a magnitude. - mu_x = np.maximum(np.abs(mu_x * coeffs.get("LMUX", 1.0)), MU_FLOOR) - mu_y = np.maximum(np.abs(mu_y * coeffs.get("LMUY", 1.0)), MU_FLOOR) - return mu_x, mu_y + mu_x, mu_y = peak_friction(fz, gamma, coeffs) + return np.abs(mu_x), np.abs(mu_y) def grip_usage( @@ -84,7 +75,7 @@ def grip_usage( """Forces as fractions of peak grip: ``(longitudinal, lateral, usage)``. ``usage`` is the friction-ellipse radius, 1.0 at the limit. All three are - NaN where the tire is off the ground. + NaN where the tire has no contact or no grip capacity. """ fx = np.asarray(fx, dtype=float) fy = np.asarray(fy, dtype=float) @@ -93,8 +84,10 @@ def grip_usage( loaded = fz > MIN_LOAD_N safe_fz = np.where(loaded, fz, 1.0) - longitudinal = np.where(loaded, fx / (mu_x * safe_fz), np.nan) - lateral = np.where(loaded, fy / (mu_y * safe_fz), np.nan) + longitudinal = np.full_like(fx, np.nan) + lateral = np.full_like(fy, np.nan) + np.divide(fx, mu_x * safe_fz, out=longitudinal, where=loaded & (mu_x > 0.)) + np.divide(fy, mu_y * safe_fz, out=lateral, where=loaded & (mu_y > 0.)) return longitudinal, lateral, np.hypot(longitudinal, lateral)