Kerner-Konhäuser model
The Kerner-Konhäuser model, where viscous traffic flow forms dense moving clusters.
- Name:
kerner_konhauser - Field:
engineering/traffic - Formalism: ordinary differential equations
- References: Kerner and Konhäuser 1993; Chowdhury, Santen, and Schadschneider 2000
Factory
@register_model(
field='engineering/traffic',
tags=('traffic', 'field', 'conservation law', 'finite volume', 'jams', 'pde'),
references=('Kerner and Konhäuser 1993', 'Chowdhury, Santen, and Schadschneider 2000'),
defaults={'until': 100.0, 'dt': 0.2, 'solver': SSPRKSolver(cfl=0.15)},
)
def kerner_konhauser(
n: int = 40,
*,
length: float = 10.0,
density: float = 0.25,
bump: float = 0.1,
V0: float = 1.0,
rho_max: float = 1.0,
tau: float = 1.0,
c0: float = 0.45,
mu: float = 0.2,
) -> System:
"""Return the Kerner-Konhäuser model, where viscous traffic flow forms dense moving clusters.
The density and the flow `q = rho v` on a ring road follow `rho_t + q_x = 0` and
`q_t + (q^2 / rho + c0^2 rho)_x = (rho V(rho) - q) / τ + μ v_xx`, with the equilibrium speed
`V(rho) = V0 (1 / (1 + exp((rho / rho_max - 0.25) / 0.06)) - 3.72e-6)`. The flow starts in equilibrium with a
bump of density, which grows into a cluster at intermediate densities.
Args:
n:
The number of cells.
length:
The length of the ring road.
density:
The mean density.
bump:
The relative height of the density bump.
V0:
The scale of the equilibrium speed.
rho_max:
The jam density.
tau:
The relaxation time of the speed.
c0:
The speed of sound of the traffic pressure `c0^2 rho`.
mu:
The viscosity.
Returns:
system:
The model, with the density `road.rho` and the flow `road.q`.
Examples:
>>> from skcomplex.models.engineering.traffic import kerner_konhauser
>>> kerner_konhauser(10).elements[0].values['mu']
0.2
"""
rho = _density_bump(n, density, bump)
values = {
'rho': rho,
'q': rho * _kerner_konhauser_speed(rho, V0, rho_max),
'V0': V0,
'rho_max': rho_max,
'tau': tau,
'c0': c0,
'mu': mu,
}
space = LatticeSpace(n, spacing=length / n)
road = FieldElement(
'road',
element=KernerKonhauserElement,
space=space,
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.c0**2 * s.target.rho,
)
relaxation = RelaxationInteraction(
'relaxation',
target='road',
state='q',
goal=lambda s: s.target.rho * _get_kerner_konhauser_speed(s.target),
time=lambda s: s.target.tau,
)
viscosity = RateInteraction(
'viscosity',
target='road',
state='q',
rate=lambda s: s.target.mu * laplacian(s.target.q / s.target.rho),
)
return System(elements=[road], interactions=[density, flow, relaxation, viscosity])
Components
class KernerKonhauserElement(Element):
"""Cell of a road with the traffic density `rho` and the flow `q` of the Kerner-Konhäuser model.
The density and the flow are conserved by conservation interactions, the flow relaxes to the equilibrium flow by a
relaxation interaction, and the viscosity adds the rate `mu laplacian(q / rho)` to the flow.
"""
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,
V0: 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,
c0: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 0.3,
mu: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
) -> None:
super().__init__(label, visual=visual)
self.rho = rho
self.q = q
self.V0 = V0
self.rho_max = rho_max
self.tau = tau
self.c0 = c0
self.mu = mu
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)
class RateInteraction(Interaction):
"""Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.
On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.
A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
`sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
members, or a fractional Brownian motion of a Hurst index.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the rate changes.
state:
The name of the state of the target that the rate changes.
rate:
A function of the namespace that returns the rate, or a number.
sources:
The paths of the elements or groups the rate and the noise read but do not change.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
divisor:
The name of the state or parameter of the target that divides the rate and the noise, a function of the
namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.
process:
The Wiener processes that drive the noise, or `None` for one independent process per component and member.
index:
The process of `process` that drives the noise, or `None` for every process, one per component of a vector
state.
convention:
The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.
hurst:
The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
of the namespace that returns a parameter, or `None` for a Wiener process.
conserved:
Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
total of the state.
Raises:
ValidationError:
The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
index is outside `(0, 1)`, or the noise options cannot combine.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RateInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
>>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
noise: Callable[[SimpleNamespace], object] | float | None = None,
divisor: str | Callable[[SimpleNamespace], object] | None = None,
process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
index: int | None = None,
convention: str = 'ito',
hurst: Callable[[SimpleNamespace], object] | float | None = None,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rate = rate
self.sources = sources
self.noise = noise
self.divisor = divisor
self.process = process
self.index = index
self.convention = convention
self.hurst = hurst
self.conserved = conserved
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."""
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
'state.'
)
raise ValidationError(error_msg)
for name in ('rate', 'noise'):
value = getattr(self, name)
if not (value is None or 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)
if not (
self.divisor is None
or callable(self.divisor)
or (isinstance(self.divisor, str) and self.divisor.isidentifier())
):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
'function of the namespace of the target.'
)
raise ValidationError(error_msg)
def _check_noise(self) -> None:
"""Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
if self.convention not in CONVENTIONS:
error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
raise ValidationError(error_msg)
if not (
self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
):
error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
raise ValidationError(error_msg)
if self.index is not None and self.process is None:
error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
raise ValidationError(error_msg)
options = [
name
for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
if given
]
if len(options) > 1:
error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for a state, functions, or noise options the interaction cannot take.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is neither callable nor a number, the
convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
"""
super()._check()
self._check_functions()
self._check_noise()
def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return function(namespace) if callable(function) else function
def _broadcast(self, value: object, state: object) -> object:
"""Return a value with one component per component of a vector state, repeating a number."""
if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
return vector(SR, [value] * len(state))
if isinstance(value, list | tuple):
return vector(SR, value)
return value
def _divide(self, value: object, namespace: SimpleNamespace) -> object:
"""Return a value divided by the divisor of the target, or the value without a divisor."""
if self.divisor is None:
return value
divisor = (
self.divisor(self._get_member(namespace, 'target'))
if callable(self.divisor)
else getattr(self._get_member(namespace, 'target'), self.divisor)
)
return operator.truediv(value, divisor)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.rate, namespace)
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
return None if self.noise is None else self._evaluate(self.noise, namespace)
def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
"""Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
if self.hurst is not None:
hurst = self._evaluate(self.hurst, namespace)
return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
process = None
if self.process is not None:
handle = self._get_member(namespace, 'process').dW
process = handle if self.index is None else handle[self.index]
return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the state and its noise law, both divided by the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target and the sources.
Returns:
laws:
The rate law, and the noise law when the interaction has noise.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
laws: list[BasePart] = [RateLaw(state, rate)]
noise = self._get_noise(namespace)
if noise is not None:
laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
return laws
Equations
road.rho[i]: diff(rho_road(t), t) == 0
road.q[i]: diff(q_road(t), t) == mu_road*laplacian(q_road(t)/rho_road(t)) + (V0_road*(1/(e^(16.666666666666668*rho_road(t)/rho_max_road - 4.166666666666667) + 1) - 3.72e-06)*rho_road(t) - q_road(t))/tau_road