GKT non-local traffic
The gas-kinetic-based traffic model, where drivers brake on the traffic at an interaction point ahead.
- Name:
gkt_traffic - Field:
engineering/traffic - Formalism: ordinary differential equations
- References: Treiber, Hennecke, and Helbing 1999; arXiv:cond-mat/0012229
Factory
@register_model(
field='engineering/traffic',
tags=('traffic', 'field', 'conservation law', 'finite volume', 'non-local', 'stop-and-go', 'pde'),
references=('Treiber, Hennecke, and Helbing 1999', 'arXiv:cond-mat/0012229'),
defaults={'until': 1.0, 'dt': 1 / 60, 'solver': SSPRKSolver(cfl=0.4)},
)
def gkt_traffic(
n: int = 200,
*,
length: float = 10.0,
density: float = 35.0,
amplitude: float = 2.0,
V0: float = 110.0,
T: float = 1.8 / 3600,
tau: float = 35 / 3600,
rho_max: float = 160.0,
gamma: float = 1.2,
) -> System:
"""Return the gas-kinetic-based traffic model, where drivers brake on the traffic at an interaction point ahead.
The density and the flow `Q = rho V` on a ring road follow `rho_t + Q_x = 0` and
`Q_t + (Q**2 / rho + rho A(rho) V**2)_x = rho (Ve* - V) / tau`, where the dynamic equilibrium speed
Ve* = V0 (1 - (theta + theta_a) / (2 A(rho_max)) (rho_a T / (1 - rho_a / rho_max))**2 B(dV))
reads the density `rho_a` and the speed `V_a` at the interaction point `x + gamma (1 / rho_max + T V)`, with the
variances `theta = A(rho) V**2`, the variance prefactor `A(rho) = A0 + dA (tanh((rho - rho_c) / d_rho) + 1)`, the
dimensionless speed difference `dV = (V - V_a) / sqrt(theta + theta_a)`, and the Boltzmann factor
`B(d) = 2 (d phi(d) + (1 + d**2) Phi(d))` of the standard normal density `phi` and distribution `Phi`. Lengths are
in kilometres and times in hours, and the defaults are the calibrated values of Treiber, Hennecke, and Helbing. The
flow starts in equilibrium with their dipole perturbation of the density, of height `amplitude`. Homogeneous
traffic is unstable between about 24 and 51 vehicles per kilometre, where the perturbation grows into stop-and-go
waves whose downstream fronts move against the traffic.
Args:
n:
The number of cells.
length:
The length of the ring road in kilometres.
density:
The mean density in vehicles per kilometre.
amplitude:
The height of the dipole perturbation in vehicles per kilometre.
V0:
The desired speed in kilometres per hour.
T:
The safe time headway in hours.
tau:
The relaxation time of the speed in hours.
rho_max:
The maximum density in vehicles per kilometre.
gamma:
The anticipation factor.
Returns:
system:
The model, with the density `road.rho` and the flow `road.Q`.
Examples:
>>> from skcomplex.models.engineering.traffic import gkt_traffic
>>> gkt_traffic(10).elements[0].values['gamma']
1.2
"""
centers = (np.arange(n) + 0.5) * length / n
start, width, wake = 0.3 * length, 0.20125, 0.805
rho = density + amplitude * (
np.cosh((centers - start) / width) ** -2 - width / wake * np.cosh((centers - start - 1.00625) / wake) ** -2
)
values = {'V0': V0, 'T': T, 'tau': tau, 'rho_max': rho_max, 'gamma': gamma}
cell = GasKineticTrafficElement('cell', V0=V0, T=T, tau=tau, rho_max=rho_max, gamma=gamma)
road = FieldElement(
'road',
element=GasKineticTrafficElement,
space=LatticeSpace(n, spacing=length / n),
values={'rho': rho, 'Q': rho * _gkt_equilibrium_speed(rho, cell), **values},
discretization='finite-volume',
)
vehicles = ConservationInteraction('density', target='road', state='rho', flux=lambda s: s.target.Q)
momentum = ConservationInteraction(
'flow',
target='road',
state='Q',
flux=_get_gkt_momentum_flux,
source=_get_gkt_relaxation,
)
return System(elements=[road], interactions=[vehicles, momentum])
Components
class GasKineticTrafficElement(Element):
"""Cell of a road with the traffic density `rho` and the flow `Q` of the gas-kinetic-based traffic model.
The density and the flow are conserved by conservation interactions, and the flow relaxes to a dynamic equilibrium
speed that brakes on the density and the speed read ahead at the interaction point. Lengths are in kilometres and
times in hours, so `T` and `tau` are fractions of an hour.
Args:
label:
The label of the cell or the road.
visual:
How the cell is drawn, or `None` for `default_visual_`.
rho:
The initial density in vehicles per kilometre.
Q:
The initial flow in vehicles per hour.
V0:
The desired speed.
T:
The safe time headway.
tau:
The relaxation time of the speed.
rho_max:
The maximum density.
gamma:
The anticipation factor, the interaction distance in units of the safe distance `1 / rho_max + T V`.
A0:
The variance prefactor of free traffic.
dA:
The increase of the variance prefactor from free to congested traffic, halved.
rho_c:
The density of the transition of the variance prefactor.
d_rho:
The width of the transition of the variance prefactor.
Examples:
>>> from skcomplex.models.engineering.traffic import GasKineticTrafficElement
>>> GasKineticTrafficElement('cell').gamma
1.2
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
rho: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 20.0,
Q: Annotated[float, StateSpec()] = 1800.0,
V0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 110.0,
T: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.8 / 3600,
tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 35 / 3600,
rho_max: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 160.0,
gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.2,
A0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.008,
dA: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.02,
rho_c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.27 * 160,
d_rho: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.05 * 160,
) -> None:
super().__init__(label, visual=visual)
self.rho = rho
self.Q = Q
self.V0 = V0
self.T = T
self.tau = tau
self.rho_max = rho_max
self.gamma = gamma
self.A0 = A0
self.dA = dA
self.rho_c = rho_c
self.d_rho = d_rho
class ConservationInteraction(Interaction):
"""Conservation law `du/dt + divergence(flux(ns)) = source(ns)` of a state `u` of a field.
The flux is one expression for a field of one dimension, or one expression per dimension, of the states of the
field, and the interaction lowers it to a finite-volume flux law with a slope limiter and a numerical flux, or,
without a numerical flux, to the rate `-divergence(flux)` that the field discretizes by its own operators, as a
spectral or finite-difference field does. Without a source the total of the state on a periodic field stays
constant. The conservation laws of several states of one field form one system, whose numerical flux reads the
wave speeds of all of them, and a flux may read parameters that vary across the cells. A flux or a source function
takes the namespace of the interaction, with the namespace `target` of the field, the time `t`, and the symbols of
the parameters a subclass declares.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state is conserved.
state:
The name of the conserved state of the field.
flux:
A function of the namespace that returns the flux, or one flux per dimension.
source:
A function of the namespace that returns the source, or `None` for none.
limiter:
The slope limiter of the reconstruction, `'minmod'`, `'van-leer'`, or `None` for first order.
riemann:
The numerical flux of a finite-volume field, `'rusanov'` or `'godunov'`, the exact solver for one state, or
`None` for the divergence of the flux on a field of another discretization.
sources:
The paths of other elements the flux and the source read, as the namespaces `sources`. A field over the
position lattice of a target on a phase space is read broadcast along the class axes.
Raises:
ValidationError:
The state is not a name, a function is not callable, or the limiter or the numerical flux is unknown.
Examples:
>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element, FieldElement
>>> from skcomplex.core.interactions import ConservationInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, u: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.u = u
>>> space = LatticeSpace((50,), spacing=0.1, boundary='periodic')
>>> values = {'u': 1.0 + np.sin(2 * np.pi * np.arange(50) / 50)}
>>> volumes = 'finite-volume'
>>> field = FieldElement('wave', element=CellElement, space=space, values=values, discretization=volumes)
>>> burgers = ConservationInteraction('burgers', target='wave', state='u', flux=lambda s: s.target.u**2 / 2)
>>> system = System(elements=[field], interactions=[burgers])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=0.5).trajectory_
>>> round(float(trajectory.values('wave.u')[-1].sum()), 6)
50.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'u',
flux: Callable[[SimpleNamespace], object] | None = None,
source: Callable[[SimpleNamespace], object] | None = None,
limiter: str | None = 'minmod',
riemann: str | None = 'rusanov',
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.flux = flux
self.source = source
self.limiter = limiter
self.riemann = riemann
self.sources = sources
def _check(self) -> None:
"""Raise an error for a state that is not a name, a function that is not callable, or an unknown scheme.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is not callable, or the limiter or the
numerical flux is unknown.
"""
super()._check()
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = f'{type(self).__name__} {self.label!r} conserves the state {self.state!r}, which is no name.'
raise ValidationError(error_msg)
for name in ('flux', 'source'):
value = getattr(self, name)
if value is not None and not callable(value):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
if self.limiter not in LIMITERS:
error_msg = f'Unknown limiter {self.limiter!r} of {self.label!r}. Use minmod, van-leer, or None.'
raise ValidationError(error_msg)
if self.riemann not in RIEMANN_SOLVERS:
error_msg = f'Unknown numerical flux {self.riemann!r} of {self.label!r}. Use rusanov, godunov, or None.'
raise ValidationError(error_msg)
def _get_flux(self, namespace: SimpleNamespace) -> object:
"""Return the flux, or `None` for none, the hook a subclass overrides to derive it from its arguments."""
return None if self.flux is None else self.flux(namespace)
def _get_transport_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the laws of the divergence of the flux, the hook a subclass overrides to lower its own flux."""
flux = self._get_flux(namespace)
if flux is None:
return []
state = getattr(self._get_member(namespace, 'target'), self.state)
if self.riemann is None:
components = list(flux) if isinstance(flux, FreeModuleElement | list | tuple) else [flux]
return [RateLaw(state, -divergence(components))]
return [FluxLaw(state, flux, limiter=self.limiter, riemann=self.riemann)]
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the laws of the divergence of the flux and the rate of the source.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the field.
Returns:
laws:
The flux law, and the rate law of the source when there is one.
"""
laws = self._get_transport_laws(namespace)
if self.source is not None:
laws.append(RateLaw(getattr(self._get_member(namespace, 'target'), self.state), self.source(namespace)))
return laws
Equations
road.rho[i]: diff(rho_road(t), t) == 0
road.Q[i]: diff(Q_road(t), t) == -1/2*((((dA_road*(tanh(-(rho_c_road - shifted(rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road)))/d_rho_road) + 1) + A0_road)*shifted(Q_road(t)/rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road))^2 + (dA_road*(tanh(-(rho_c_road - rho_road(t))/d_rho_road) + 1) + A0_road)*Q_road(t)^2/rho_road(t)^2)*(((Q_road(t)/rho_road(t) - shifted(Q_road(t)/rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road)))^2/((dA_road*(tanh(-(rho_c_road - shifted(rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road)))/d_rho_road) + 1) + A0_road)*shifted(Q_road(t)/rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road))^2 + (dA_road*(tanh(-(rho_c_road - rho_road(t))/d_rho_road) + 1) + A0_road)*Q_road(t)^2/rho_road(t)^2) + 1)*(erf(1/2*sqrt(2)*(Q_road(t)/rho_road(t) - shifted(Q_road(t)/rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road)))/sqrt((dA_road*(tanh(-(rho_c_road - shifted(rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road)))/d_rho_road) + 1) + A0_road)*shifted(Q_road(t)/rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road))^2 + (dA_road*(tanh(-(rho_c_road - rho_road(t))/d_rho_road) + 1) + A0_road)*Q_road(t)^2/rho_road(t)^2)) + 1) + sqrt(2)*(Q_road(t)/rho_road(t) - shifted(Q_road(t)/rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road)))*e^(-1/2*(Q_road(t)/rho_road(t) - shifted(Q_road(t)/rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road)))^2/((dA_road*(tanh(-(rho_c_road - shifted(rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road)))/d_rho_road) + 1) + A0_road)*shifted(Q_road(t)/rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road))^2 + (dA_road*(tanh(-(rho_c_road - rho_road(t))/d_rho_road) + 1) + A0_road)*Q_road(t)^2/rho_road(t)^2))/(sqrt(pi)*sqrt((dA_road*(tanh(-(rho_c_road - shifted(rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road)))/d_rho_road) + 1) + A0_road)*shifted(Q_road(t)/rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road))^2 + (dA_road*(tanh(-(rho_c_road - rho_road(t))/d_rho_road) + 1) + A0_road)*Q_road(t)^2/rho_road(t)^2)))*T_road^2*shifted(rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road))^2/((dA_road*(tanh(-(rho_c_road - rho_max_road)/d_rho_road) + 1) + A0_road)*(shifted(rho_road(t), gamma_road*(T_road*Q_road(t)/rho_road(t) + 1/rho_max_road))/rho_max_road - 1)^2) - 2)*V0_road + 2*Q_road(t)/rho_road(t))*rho_road(t)/tau_road