Moment equations of Langevin traffic with a jam instability

The density and flux equations of Langevin traffic on a ring road, whose uniform flow breaks into jams.

• Name: traffic_density_equation
• Field: engineering/traffic
• Formalism: ordinary differential equations
• References: Bertin 2021

Factory

@register_model(
    field='engineering/traffic',
    tags=('traffic', 'field', 'jams', 'instability', 'pde'),
    references=('Bertin 2021',),
    defaults={'until': 400.0, 'dt': 4.0, 'solver': RK4Solver(dt=0.1)},
)
def traffic_density_equation(
    n: int = 100,
    *,
    length: float = 50.0,
    density: float = 0.5,
    gamma: float = 1.0,
    u0: float = 1.0,
    rho_max: float = 1.0,
    D1: float = 0.2,
    perturbation: float = 0.01,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the density and flux equations of Langevin traffic on a ring road, whose uniform flow breaks into jams.

    The density `rho` and the flux `w` of cars that relax to the desired speed `u(rho) = u0 (1 - rho / rho_max)` at
    the rate `gamma`, with the diffusion coefficient `D(rho) = D1 rho`, follow the moment equations
    `rho_t + w_x = (D rho_x)_x` and `w_t + (2 u w - rho u^2)_x = -gamma (w - rho u) + (D w_x)_x`, closed by expanding
    the velocity distribution around `u(rho)`. A long density wave of the uniform flow grows when the effective
    diffusion coefficient `D(rho) + ((rho u^2)' - (rho u)'^2) / gamma = D(rho) - rho^2 u'(rho)^2 / gamma` is negative,
    so the flow is unstable above the density `D1 gamma rho_max^2 / u0^2`, and the perturbation grows into dense jams.
    The uniform density starts with independent Gaussian perturbations of relative size
    `perturbation` in every cell, at the flux `rho u(rho)`.

    Args:
        n:
            The number of cells.

        length:
            The length of the ring road.

        density:
            The mean density.

        gamma:
            The relaxation rate of the speed.

        u0:
            The desired speed on an empty road.

        rho_max:
            The jam density, at which the desired speed vanishes.

        D1:
            The diffusion coefficient per unit density.

        perturbation:
            The relative standard deviation of the initial density perturbations.

        random_state:
            The seed or generator of the initial perturbations.

    Returns:
        system:
            The model, with the density `road.rho` and the flux `road.w`.

    Examples:
        >>> from skcomplex.models.engineering.traffic import traffic_density_equation
        >>> traffic_density_equation(10).elements[0].values['D1']
        0.2
    """
    rng = np.random.default_rng(random_state)
    rho = density * (1 + perturbation * rng.standard_normal(n))
    values = {
        'rho': rho,
        'w': rho * u0 * (1 - rho / rho_max),
        'gamma': gamma,
        'u0': u0,
        'rho_max': rho_max,
        'D1': D1,
    }
    space = LatticeSpace(n, spacing=length / n)
    road = FieldElement('road', element=TrafficMomentsElement, space=space, values=values)
    return System(elements=[road])

Components

class TrafficMomentsElement(Element):
    """Density `rho` and flow `w` of Langevin traffic on a road, the first two moments of its speeds."""

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        rho: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
        w: Annotated[object, StateSpec()] = 0.25,
        gamma: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        u0: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        rho_max: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        D1: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
    ) -> None:
        super().__init__(label, visual=visual)
        self.rho = rho
        self.w = w
        self.gamma = gamma
        self.u0 = u0
        self.rho_max = rho_max
        self.D1 = D1

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the closed moment equations of the density and the flow."""
        return [
            RateLaw(s.rho, -divergence([s.w]) + laplacian(s.D1 * s.rho**2 / 2)),
            RateLaw(
                s.w,
                -divergence(
                    [
                        2 * LangevinDrivingInteraction._target_speed(s.rho, s) * s.w
                        - s.rho * LangevinDrivingInteraction._target_speed(s.rho, s) ** 2,
                    ],
                )
                - s.gamma * (s.w - s.rho * LangevinDrivingInteraction._target_speed(s.rho, s))
                + s.D1 * (laplacian(s.rho * s.w) + s.rho * laplacian(s.w) - s.w * laplacian(s.rho)) / 2,
            ),
        ]

Equations

road.rho[i]: diff(rho_road(t), t) == laplacian(1/2*D1_road*rho_road(t)^2) - partial(w_road(t), 0)
road.w[i]: diff(w_road(t), t) == 1/2*(laplacian(w_road(t))*rho_road(t) - laplacian(rho_road(t))*w_road(t) + laplacian(rho_road(t)*w_road(t)))*D1_road - (u0_road*(rho_road(t)/rho_max_road - 1)*rho_road(t) + w_road(t))*gamma_road - partial(-u0_road^2*(rho_road(t)/rho_max_road - 1)^2*rho_road(t) - 2*u0_road*(rho_road(t)/rho_max_road - 1)*w_road(t), 0)