Example: Runway Crosswind & Wind-Shift Operability

An airport decides whether a runway is usable over the coming day. Two things close it: the crosswind component of the wind (its speed times the sine of the angle between the wind and the runway heading) exceeding an aircraft limit, and abrupt wind shifts between hours. Both are driven by a wind that evolves over the day — and wind direction is a circular quantity that can wander across North and switches between two regimes. The full runnable source is examples/runway_crosswind/model.py.

This example is the process side of prophys end-to-end: a coupled stochastic wind, a circular direction with a bimodal step, and operability read off the marginal of the day’s worst crosswind — its expectation, its tail, and the exceedance probability that decides downtime.

The two runways inside a wind rose of 96 000 simulated hourly winds (4 000 days), petals by the direction the wind comes from and shaded by speed. The wind blows from around 250°, straight along runway 07/25, and almost square across runway 16/34. On the right, the probability that the worst hour of the day exceeds the 12 m/s crosswind limit: practically zero for 07/25, nearly one day in two for 16/34.

Step 1 — The wind as a coupled process

Wind speed is a mean-reverting Ornstein–Uhlenbeck process. Wind direction is a circular_autoregressive() process reverting to the prevailing heading, with a two-component Gaussian mixture innovation for the backing/veering regimes — a step distribution that is bimodal, not Gaussian. Sampling such a process is exactly what a mixture’s derived quantile makes possible (a bare mixture had no quantile before):

# Speed: mean-reverting around 8 m/s.
speed_proc = prp.ornstein_uhlenbeck(init=8.0, mean=8.0, theta=0.5, sigma=2.5, n_steps=N_HOURS)

direction_innovation = prp.Mixture(
    [prp.Gaussian(-6.0, 4.0), prp.Gaussian(6.0, 4.0)], logits=jnp.array([0.0, 0.0])
)
direction_proc = circular_autoregressive(
    init=PREVAILING_DIR_DEG,
    mean=PREVAILING_DIR_DEG,
    phi=0.85,
    innovation=direction_innovation,
    n_steps=N_HOURS,
    period=360.0,
)

The circular geometry matters: reverting toward the prevailing 250° from a direction near North follows the short arc across the 0/360 seam, and every sampled hour stays a valid heading in [0, 360).

The two are coupled into one JointProcess so gusty hours and shifting hours co-occur, then exposed as a single joint random variable whose components share the underlying draw:

wind_process = prp.JointProcess([speed_proc, direction_proc], corr_raw=jnp.array([0.35]))
wind = prp.JointRandomVariable("wind", wind_process, components=["speed", "direction"])

speed_path = prp.clip(wind["speed"], 0.0, None)  # (..., n_hours) [m/s]
direction_path = wind["direction"]  # (..., n_hours) [deg]

Twelve sampled days (thin lines), the hourly median and the band holding 90 % of 4 000 days. Top: the speed reverts to 8 m/s. Middle: the direction wanders around the prevailing 250°, some days drifting 40° away for hours. Bottom: the crosswind on runway 16/34 that follows from both. Most hours stay below the limit, but a single hour above it (dots) is enough to make the day’s worst crosswind exceed it.

Step 2 — Operability read off the wind

The instantaneous crosswind on a runway of a given heading is \(V\,\lvert\sin(\theta - \text{heading})\rvert\), and the operability question is about the worst crosswind of the day — a maximum along the time axis, observed through small anemometer noise so it is a proper UncertainAttribute:

def crosswind_component(heading_deg: float) -> prp.Expr:
    angle = prp.deg2rad(direction_path - heading_deg)
    return prp.abs_(speed_path * prp.sin(angle))


max_crosswind = {
    name: prp.UncertainAttribute(
        f"max_crosswind_{name.replace('/', '_')}",
        prp.Gaussian(mean=_time_max(crosswind_component(heading)), sigma=anemometer_sigma),
        unit="m/s",
        domain=f"worst crosswind of the day on runway {name}",
    )
    for name, heading in RUNWAYS.items()
}

The largest hour-to-hour wind shift is a circular first difference, so a swing across North is measured as the short arc rather than a spurious near-360° jump:

hourly_shift = time_diff(direction_path, geometry=Circular(360.0))  # (..., n_hours - 1) [deg]

Step 3 — Expectation, tail, and downtime probability

Compiling marginalizes the wind process by Monte Carlo. For each runway the model reports the expected worst crosswind, its Conditional Value-at-Risk (the mean of the worst 5 % of days — the number that actually governs operations), and the downtime probability that the day’s worst crosswind exceeds the limit. The exceedance is computed two ways: directly by Monte Carlo, and from a Gaussian fitted to the sampled maxima via its first-class cdf (\(P(X > \ell) = 1 - F(\ell)\)):

compiled = model.compile(mode="eval", n_samples=4000)

e = float(compiled.expectation(attr))
cvar = float(compiled.cvar(attr, alpha=0.95))
# Downtime probability P(worst crosswind > limit), by Monte Carlo ...
p_mc = 1.0 - float(compiled.prob(attr, CROSSWIND_LIMIT))
# ... and directly from a Gaussian fitted to the sampled maxima
# (first-class cdf: P(X > limit) = 1 - F(limit)).
samples = compiled.sample(attr, jax.random.PRNGKey(0), (4000,))
fitted = prp.Gaussian(float(jnp.mean(samples)), float(jnp.std(samples)) + 1e-6)
p_cdf = 1.0 - float(fitted.cdf(jnp.asarray(CROSSWIND_LIMIT)).evaluate({}))

Left: the worst crosswind of the day on each runway, with the share above the 12 m/s limit shaded and the CVaR95 marked by triangles. Right: the probability that the worst crosswind exceeds a threshold, by Monte Carlo (solid) and from the fitted Gaussian (dashed). The two agree near the middle of each distribution, which is where the limit falls for 16/34 (48.6 % against 51.1 %), but a maximum over 24 hours is right-skewed and the Gaussian understates its upper tail, by an order of magnitude and more far out.

For the modelled airport the runway aligned with the prevailing wind sees a low expected crosswind and essentially no downtime, while the near-perpendicular runway straddles the limit with about 50 % downtime — so the model picks the aligned runway.

Left: the density of the hourly direction step, the two-component mixture of backing and veering, with its 10, 50 and 90 % quantiles obtained by inverting its CDF. Right: the largest hour-to-hour shift of the day. It averages 14.6°, and on the worst 5 % of days it is 20.2°.

Results

The console output of python examples/runway_crosswind/model.py:

=== Innovation mixture (backing/veering hourly step) ===
  median step=-0.0 deg   10-90% step range=[-9.4, +9.4] deg
  (mixture.quantile is derived from its CDF by inversion)

=== Runway operability (crosswind limit 12 m/s) ===
  runway 07/25 (hdg 250): E[max xwind]= 3.9  CVaR95= 7.5  P(downtime)= 0.0% (mc) /  0.0% (cdf)  q95= 6.2 m/s
  runway 16/34 (hdg 160): E[max xwind]=12.1  CVaR95=15.4  P(downtime)=48.6% (mc) / 51.1% (cdf)  q95=14.5 m/s
  -> preferred runway (lowest expected crosswind): 07/25

=== Wind-shift risk ===
  E[largest hourly shift] = 14.59525203704834 deg
  CVaR95[largest hourly shift] = 20.167617797851562 deg

=== Export ===
exported attributes: ['max_crosswind_07_25', 'max_crosswind_16_34', 'max_wind_shift']

Primitives exercised

Category	Used for
ornstein_uhlenbeck(), circular_autoregressive()	Wind speed (mean-reverting) and direction (circular, seam-correct)
Mixture	Bimodal backing/veering hourly step — usable as an innovation via its derived quantile
JointProcess, RandomVariable	Coupling gusty and shifting hours; exposing a shared joint draw
time_diff() with Circular	Seam-safe hour-to-hour wind-shift magnitude
expectation / cvar / prob and Gaussian cdf/quantile	Expected & tail crosswind; downtime probability by Monte Carlo and by CDF

Running it

$ python examples/runway_crosswind/model.py