Payne model
The Payne model, where the speed of traffic relaxes to the Greenshields speed and anticipates density.
- Name:
payne - Field:
engineering/traffic - Formalism: ordinary differential equations
- References: Payne 1971; Chowdhury, Santen, and Schadschneider 2000
Factory
@register_model(
field='engineering/traffic',
tags=('traffic', 'field', 'conservation law', 'finite volume', 'pde'),
references=('Payne 1971', 'Chowdhury, Santen, and Schadschneider 2000'),
defaults={'until': 1.0, 'dt': 0.05, 'solver': SSPRKSolver()},
)
def payne(
n: int = 200,
*,
density: float = 0.3,
bump: float = 0.2,
vmax: float = 1.0,
rho_max: float = 1.0,
tau: float = 0.5,
) -> System:
"""Return the Payne model, where the speed of traffic relaxes to the Greenshields speed and anticipates density.
The density and the flow `q = rho v` on a ring road of length 1 follow `rho_t + q_x = 0` and
`q_t + (q^2 / rho + P(rho))_x = (rho V(rho) - q) / τ`, with `V(rho) = vmax (1 - rho / rho_max)` and the
anticipation pressure `P(rho) = vmax rho / (2 τ rho_max)`. This is the speed equation
`v_t + v v_x = (V - v) / τ + V'(rho) rho_x / (2 τ rho)` in conservation form. The flow starts in equilibrium
with a bump of density.
Args:
n:
The number of cells.
density:
The mean density.
bump:
The relative height of the density bump.
vmax:
The speed on an empty road.
rho_max:
The jam density.
tau:
The relaxation time of the speed.
Returns:
system:
The model, with the density `road.rho` and the flow `road.q`.
Examples:
>>> from skcomplex.models.engineering.traffic import payne
>>> payne(10).elements[0].values['tau']
0.5
"""
rho = _density_bump(n, density, bump)
values = {'rho': rho, 'q': rho * vmax * (1 - rho / rho_max), 'vmax': vmax, 'rho_max': rho_max, 'tau': tau}
road = FieldElement(
'road',
element=PayneElement,
space=LatticeSpace(n, spacing=1 / n),
values=values,
discretization='finite-volume',
)
density = ConservationInteraction('density', target='road', state='rho', flux=lambda s: s.target.q)
flow = ConservationInteraction(
'flow',
target='road',
state='q',
flux=lambda s: s.target.q**2 / s.target.rho
+ s.target.vmax * s.target.rho / (2 * s.target.tau * s.target.rho_max),
)
relaxation = RelaxationInteraction(
'relaxation',
target='road',
state='q',
goal=lambda s: s.target.rho * s.target.vmax * (1 - s.target.rho / s.target.rho_max),
time=lambda s: s.target.tau,
)
return System(elements=[road], interactions=[density, flow, relaxation])
Components
class PayneElement(Element):
"""Cell of a road with the traffic density `rho` and the flow `q` of the Payne model.
The density and the flow are conserved by conservation interactions, and the flow relaxes to the Greenshields flow
by a relaxation interaction.
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
rho: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
q: Annotated[object, StateSpec()] = 0.21,
vmax: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
rho_max: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
tau: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
) -> None:
super().__init__(label, visual=visual)
self.rho = rho
self.q = q
self.vmax = vmax
self.rho_max = rho_max
self.tau = tau
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
class RelaxationInteraction(RateInteraction):
"""Relaxation `(goal(ns) - x) / time(ns)` of a state `x` of a target towards a goal, with optional noise.
With a constant goal and time and an additive noise of amplitude `sigma`, the state follows the Ornstein-Uhlenbeck
process `dx = (goal - x) / time dt + sigma dW`, whose stationary variance is `sigma**2 time / 2`. Without noise it
is first-order smoothing, or adaptive expectation, towards the goal. A goal or a time function takes the namespace
of the interaction, as for `RateInteraction`, and a number stands for a constant. A goal that is a number applies
to every component of a vector state.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state relaxes.
state:
The name of the state of the target that relaxes.
goal:
A function of the namespace that returns the goal, or a number.
time:
A function of the namespace that returns the relaxation time, or a positive number.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
Raises:
ValidationError:
The state is not a name, or a function is neither callable nor a number.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RelaxationInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> smoothing = RelaxationInteraction('smoothing', target='belief', state='n', goal=10.0, time=0.5)
>>> system = System(elements=[CompartmentElement('belief', n=0.0)], interactions=[smoothing])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=5.0).trajectory_
>>> round(float(trajectory.values('belief.n')[-1]), 3)
10.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
goal: Callable[[SimpleNamespace], object] | float = 0.0,
time: Callable[[SimpleNamespace], object] | float = 1.0,
noise: Callable[[SimpleNamespace], object] | float | None = None,
) -> None:
super().__init__(label, target=target, state=state, noise=noise)
self.goal = goal
self.time = time
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
super()._check_functions()
for name in ('goal', 'time'):
value = getattr(self, name)
if not (callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
def _get_goal(self, namespace: SimpleNamespace) -> object:
"""Return the goal, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.goal, namespace)
def _get_time(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation time, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.time, namespace)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation `(goal - x) / time` of the state.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The rate of the state.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
goal = self._broadcast(self._get_goal(namespace), state)
return (goal - state) / self._get_time(namespace)
Equations
road.rho[i]: diff(rho_road(t), t) == 0
road.q[i]: diff(q_road(t), t) == -(vmax_road*(rho_road(t)/rho_max_road - 1)*rho_road(t) + q_road(t))/tau_road