Beer distribution game
Sterman's beer distribution game, a supply chain whose orders swing wider at every echelon upstream.
- Name:
beer_game - Field:
economics/markets - Formalism: ordinary differential equations
- References: Sterman 1989
Factory
@register_model(
field='economics/markets',
tags=('supply chain', 'bullwhip effect', 'delay', 'stocks and flows', 'system dynamics', 'ode'),
references=('Sterman 1989',),
defaults={'until': 36.0, 'dt': 0.5},
observables=lambda system: [
*(ObservableSpec('orders', BreweryInteraction._orders, element=label) for label in ECHELONS),
*(ObservableSpec('shipments', EchelonElement._shipments, element=label) for label in ECHELONS),
ObservableSpec('customers.orders', _customer_orders),
],
)
def beer_game(
inventory: float = 12.0,
*,
demand: float = 4.0,
demand_step: float = 4.0,
step_time: float = 5.0,
theta: float = 0.36,
alpha: float = 0.26,
beta: float = 0.34,
dispatch_time: float = 1.0,
mail_delay: float = 2.0,
shipping_delay: float = 2.0,
production_delay: float = 3.0,
) -> System:
"""Return Sterman's beer distribution game, a supply chain whose orders swing wider at every echelon upstream.
A retailer, a wholesaler, a distributor, and a factory each hold an inventory of beer and a backlog of orders
waiting to ship. Each ships its backlog over the dispatch time, limited by its inventory, and the beer reaches the
next echelon downstream after the shipping delay. Orders reach the next echelon upstream after the mail delay, and
the factory brews its orders over the production delay. Customer demand steps from 4 to 8 cases a week in week 5.
Every echelon orders by Sterman's anchoring and adjustment rule, its expected demand plus `alpha` times the gap
between its desired stock `S'` and its effective inventory, the inventory less the backlog, plus `beta` times its
supply line, the beer it ordered and has not received. The expected demand adapts to incoming orders at the rate
`theta`. The defaults of `theta`, `alpha`, and `beta` are the median estimates of Sterman 1989, and `S'` is set so
the chain starts in equilibrium at 4 cases a week. The game is weekly and its delays are fixed, while this model is
a continuous approximation in which third-order delays stand in for the fixed ones. The amplitude of the orders
grows from the retailer to the factory, the bullwhip effect.
Args:
inventory:
The initial inventory of every echelon, in cases.
demand:
The customer demand, in cases a week.
demand_step:
The permanent change of the customer demand.
step_time:
The week of the change of the customer demand.
theta:
The rate at which the expected demand adapts to incoming orders.
alpha:
The fraction of the stock gap ordered per week.
beta:
The weight of the supply line in the stock gap.
dispatch_time:
The time over which an echelon ships its backlog.
mail_delay:
The average time an order takes to reach the next echelon upstream.
shipping_delay:
The average time a shipment takes to reach the next echelon downstream.
production_delay:
The average time the factory takes to brew its orders.
Returns:
system:
The model, with the stocks `<echelon>.inventory`, `<echelon>.in_transit`, `<echelon>.backlog`,
`<echelon>.expected`, and `<echelon>.supply_line` of the echelons `retailer`, `wholesaler`, `distributor`,
and `factory`, and the flows `<upstream>_to_<downstream>.shipments` and
`<upstream>_to_<downstream>.arrivals` of the links, `customers.sales`, and `brewery.starts` and
`brewery.completions`. The rates `<echelon>_backlog` and `<echelon>_expected` add the orders that reach
every echelon and adapt its expectation, and the rates `<echelon>_supply_line` keep its supply line. The
observables of the model are the orders `<echelon>.orders` and the shipments
`<echelon>.shipments` of every echelon, and the orders of the customers `customers.orders`.
Examples:
>>> from skcomplex.models.economics.markets import beer_game
>>> [element.label for element in beer_game().elements]
['retailer', 'wholesaler', 'distributor', 'factory']
"""
throughput = demand * dispatch_time
echelons: list[EchelonElement] = []
for label in ECHELONS:
pipeline = production_delay if label == 'factory' else shipping_delay
lead_time = pipeline if label == 'factory' else mail_delay + shipping_delay
supply_line = demand * lead_time
values = {
'inventory': inventory,
'in_transit': demand * pipeline,
'backlog': throughput,
'expected': demand,
'supply_line': supply_line,
'desired': inventory - throughput + beta * supply_line,
'theta': theta,
'alpha': alpha,
'beta': beta,
'dispatch_time': dispatch_time,
'mail_delay': mail_delay,
'shipping_delay': shipping_delay,
}
if label == 'retailer':
values.update(demand=demand, demand_step=demand_step, step_time=step_time)
elif label == 'factory':
values.update(production_delay=production_delay)
kind = {'retailer': RetailerElement, 'factory': FactoryElement}.get(label, EchelonElement)
echelons.append(kind(label).set_params(**values))
links = [
interaction
for downstream, upstream in pairwise(ECHELONS)
for interaction in (
*_get_link_rates(downstream, upstream),
SupplyLinkInteraction(f'{upstream}_to_{downstream}', downstream=downstream, upstream=upstream),
)
]
customers = [*_get_customer_rates(), RetailDemandInteraction('customers', retailer='retailer')]
brewery = [_get_brewery_rate(), BreweryInteraction('brewery', factory='factory')]
return System(elements=echelons, interactions=[*links, *customers, *brewery])
Components
class RetailerElement(EchelonElement):
"""Retailer of a supply chain, an echelon whose customers order their demand, which steps once.
Args:
label:
The label of the retailer.
visual:
How the element is drawn, or `None` for `default_visual_`.
inventory:
The initial inventory, in cases.
in_transit:
The initial beer on its way to the echelon.
backlog:
The initial orders waiting to ship.
expected:
The initial expected demand.
supply_line:
The initial beer ordered and not received.
desired:
The desired stock `S'`.
theta:
The rate at which the expected demand adapts to incoming orders.
alpha:
The fraction of the stock gap ordered per week.
beta:
The weight of the supply line in the stock gap.
dispatch_time:
The time over which the echelon ships its backlog.
mail_delay:
The average time an order of the echelon takes to reach its supplier.
shipping_delay:
The average time a shipment takes to reach the echelon from its supplier.
demand:
The customer demand, in cases a week.
demand_step:
The permanent change of the demand.
step_time:
The week of the change of the demand.
Examples:
>>> from skcomplex.models.economics.markets import RetailerElement
>>> RetailerElement('retailer').demand
4.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
inventory: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo), quantity='beer')] = 12.0,
in_transit: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo), quantity='beer')] = 8.0,
backlog: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
expected: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
supply_line: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 16.0,
desired: Annotated[float, ParameterSpec()] = 13.44,
theta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.36,
alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.26,
beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.34,
dispatch_time: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
mail_delay: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
shipping_delay: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
demand: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
demand_step: Annotated[float, ParameterSpec()] = 4.0,
step_time: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.0,
) -> None:
super().__init__(
label,
visual=visual,
inventory=inventory,
in_transit=in_transit,
backlog=backlog,
expected=expected,
supply_line=supply_line,
desired=desired,
theta=theta,
alpha=alpha,
beta=beta,
dispatch_time=dispatch_time,
mail_delay=mail_delay,
shipping_delay=shipping_delay,
)
self.demand = demand
self.demand_step = demand_step
self.step_time = step_time
class EchelonElement(Element):
"""Echelon of a supply chain, which ships its backlog from its inventory and orders by anchoring and adjustment.
It ships its backlog over the dispatch time, limited by its inventory. Its expected demand adapts to incoming
orders at the rate `theta`, and it orders its expected demand plus `alpha` times the gap between its desired stock
and its effective inventory, the inventory less the backlog, less `beta` times its supply line, never negative.
Its orders reach its supplier after the mail delay, and the beer its supplier ships arrives after the shipping
delay.
Args:
label:
The label of the echelon.
visual:
How the element is drawn, or `None` for `default_visual_`.
inventory:
The initial inventory, in cases.
in_transit:
The initial beer on its way to the echelon.
backlog:
The initial orders waiting to ship.
expected:
The initial expected demand.
supply_line:
The initial beer ordered and not received.
desired:
The desired stock `S'`.
theta:
The rate at which the expected demand adapts to incoming orders.
alpha:
The fraction of the stock gap ordered per week.
beta:
The weight of the supply line in the stock gap.
dispatch_time:
The time over which the echelon ships its backlog.
mail_delay:
The average time an order of the echelon takes to reach its supplier.
shipping_delay:
The average time a shipment takes to reach the echelon from its supplier.
Examples:
>>> from skcomplex.models.economics.markets import EchelonElement
>>> EchelonElement('retailer', inventory=10.0).inventory
10.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
inventory: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo), quantity='beer')] = 12.0,
in_transit: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo), quantity='beer')] = 8.0,
backlog: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
expected: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
supply_line: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 16.0,
desired: Annotated[float, ParameterSpec()] = 13.44,
theta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.36,
alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.26,
beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.34,
dispatch_time: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
mail_delay: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
shipping_delay: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
) -> None:
super().__init__(label, visual=visual)
self.inventory = inventory
self.in_transit = in_transit
self.backlog = backlog
self.expected = expected
self.supply_line = supply_line
self.desired = desired
self.theta = theta
self.alpha = alpha
self.beta = beta
self.dispatch_time = dispatch_time
self.mail_delay = mail_delay
self.shipping_delay = shipping_delay
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the shipping of the backlog and the decay of the expected demand.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
return EchelonElement._echelon_laws(s)
@staticmethod
def _shipments(echelon: SimpleNamespace) -> Expression:
"""Return the shipments of an echelon, its backlog limited by its inventory over the dispatch time."""
return minimum(echelon.backlog, echelon.inventory) / echelon.dispatch_time
@staticmethod
def _echelon_laws(s: SimpleNamespace) -> list[BaseLaw]:
"""Return the laws of an echelon, the shipping of its backlog and the decay of its expectation."""
return [
RateLaw(s.backlog, -EchelonElement._shipments(s)),
RateLaw(s.expected, -s.theta * s.expected),
]
class FactoryElement(EchelonElement):
"""Factory of a supply chain, an echelon whose orders a brewery brews over the production delay.
Args:
label:
The label of the factory.
visual:
How the element is drawn, or `None` for `default_visual_`.
inventory:
The initial inventory, in cases.
in_transit:
The initial beer on its way to the echelon.
backlog:
The initial orders waiting to ship.
expected:
The initial expected demand.
supply_line:
The initial beer ordered and not received.
desired:
The desired stock `S'`.
theta:
The rate at which the expected demand adapts to incoming orders.
alpha:
The fraction of the stock gap ordered per week.
beta:
The weight of the supply line in the stock gap.
dispatch_time:
The time over which the echelon ships its backlog.
mail_delay:
The average time an order of the echelon takes to reach its supplier.
shipping_delay:
The average time a shipment takes to reach the echelon from its supplier.
production_delay:
The average time the brewery takes to brew an order, through a third-order delay.
Examples:
>>> from skcomplex.models.economics.markets import FactoryElement
>>> FactoryElement('factory').production_delay
3.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
inventory: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo), quantity='beer')] = 12.0,
in_transit: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo), quantity='beer')] = 8.0,
backlog: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
expected: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
supply_line: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 16.0,
desired: Annotated[float, ParameterSpec()] = 13.44,
theta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.36,
alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.26,
beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.34,
dispatch_time: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
mail_delay: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
shipping_delay: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
production_delay: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 3.0,
) -> None:
super().__init__(
label,
visual=visual,
inventory=inventory,
in_transit=in_transit,
backlog=backlog,
expected=expected,
supply_line=supply_line,
desired=desired,
theta=theta,
alpha=alpha,
beta=beta,
dispatch_time=dispatch_time,
mail_delay=mail_delay,
shipping_delay=shipping_delay,
)
self.production_delay = production_delay
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
class SupplyLinkInteraction(TransferInteraction):
"""Link of a supply chain, which ships beer from the upstream echelon to the downstream echelon.
The shipments and their arrivals, after a third-order delay of the shipping delay of the downstream echelon, are
the legs of one transfer, from the inventory of the upstream echelon to the beer in transit to the downstream
echelon and on to its inventory.
Args:
label:
The label of the link.
downstream:
The path of the echelon that orders.
upstream:
The path of the echelon that ships.
Examples:
>>> from skcomplex.models.economics.markets import SupplyLinkInteraction
>>> SupplyLinkInteraction('link', downstream='retailer', upstream='wholesaler').upstream
'wholesaler'
"""
def __init__(
self,
label: str,
*,
downstream: Annotated[str, RoleSpec()],
upstream: Annotated[str, RoleSpec()],
) -> None:
super().__init__(
label,
source=upstream,
target=downstream,
state='inventory',
into='in_transit',
amount=lambda s: EchelonElement._shipments(s.upstream),
flow='shipments',
legs=[
(
'downstream.in_transit',
'downstream.inventory',
lambda s: SupplyLinkInteraction._arrivals(s.upstream, s.downstream),
'arrivals',
),
],
)
self.downstream = downstream
self.upstream = upstream
def _get_holders(self, namespace: SimpleNamespace) -> list[tuple[SimpleNamespace, SimpleNamespace]]:
"""Return the upstream echelon, which ships, and the downstream echelon, which receives."""
return [(namespace.upstream, namespace.downstream)]
@staticmethod
def _arrivals(upstream: SimpleNamespace, downstream: SimpleNamespace) -> Expression:
"""Return the arrivals at the downstream echelon, the shipments after a third-order shipping delay."""
return distributed_delay(EchelonElement._shipments(upstream), order=3, time=downstream.shipping_delay)
class RetailDemandInteraction(TransferInteraction):
"""Customers of a supply chain, who take the shipments of the retailer.
The sales are a transfer from the inventory of the retailer to the customers.
Args:
label:
The label of the customers.
retailer:
The path of the retailer.
Examples:
>>> from skcomplex.models.economics.markets import RetailDemandInteraction
>>> RetailDemandInteraction('customers', retailer='retailer').retailer
'retailer'
"""
def __init__(
self,
label: str,
*,
retailer: Annotated[str, RoleSpec()],
) -> None:
super().__init__(
label,
source=retailer,
target=retailer,
state='inventory',
into=Sink('customers'),
amount=lambda s: EchelonElement._shipments(s.retailer),
flow='sales',
)
self.retailer = retailer
def _get_holders(self, namespace: SimpleNamespace) -> list[tuple[SimpleNamespace, SimpleNamespace]]:
"""Return the retailer, whose inventory the customers take."""
return [(namespace.retailer, namespace.retailer)]
class BreweryInteraction(TransferInteraction):
"""Brewery of a supply chain, which brews the orders of the factory over its production delay.
The brewing is a transfer from the brewery into the beer in production at the factory and on to its inventory.
Args:
label:
The label of the brewery.
factory:
The path of the factory.
Examples:
>>> from skcomplex.models.economics.markets import BreweryInteraction
>>> BreweryInteraction('brewery', factory='factory').factory
'factory'
"""
def __init__(
self,
label: str,
*,
factory: Annotated[str, RoleSpec()],
) -> None:
super().__init__(
label,
source=factory,
target=factory,
state=Source('brewery'),
into='in_transit',
amount=lambda s: BreweryInteraction._orders(s.factory),
flow='starts',
legs=[('in_transit', 'inventory', lambda s: BreweryInteraction._completions(s.factory), 'completions')],
)
self.factory = factory
def _get_holders(self, namespace: SimpleNamespace) -> list[tuple[SimpleNamespace, SimpleNamespace]]:
"""Return the factory, whose beer in production and inventory the brewery fills."""
return [(namespace.factory, namespace.factory)]
@staticmethod
def _orders(echelon: SimpleNamespace) -> Expression:
"""Return the orders of an echelon by Sterman's anchoring and adjustment rule, never negative."""
effective = echelon.inventory - echelon.backlog
orders = echelon.expected + echelon.alpha * (echelon.desired - effective - echelon.beta * echelon.supply_line)
return maximum(orders, 0)
@staticmethod
def _completions(factory: SimpleNamespace) -> Expression:
"""Return the brews completed for the factory, its orders after a third-order production delay."""
return distributed_delay(BreweryInteraction._orders(factory), order=3, time=factory.production_delay)
Equations
retailer.inventory: diff(inventory_retailer(t), t) == 3*_material_delay3_0_wholesaler_to_retailer_2/shipping_delay_retailer - where(inventory_retailer(t) < backlog_retailer(t), inventory_retailer(t), backlog_retailer(t))/dispatch_time_retailer
retailer.in_transit: diff(in_transit_retailer(t), t) == -3*_material_delay3_0_wholesaler_to_retailer_2/shipping_delay_retailer + where(inventory_wholesaler(t) < backlog_wholesaler(t), inventory_wholesaler(t), backlog_wholesaler(t))/dispatch_time_wholesaler
retailer.backlog: diff(backlog_retailer(t), t) == demand_retailer - where(inventory_retailer(t) < backlog_retailer(t), inventory_retailer(t), backlog_retailer(t))/dispatch_time_retailer + where(t >= step_time_retailer, demand_step_retailer, 0)
retailer.expected: diff(expected_retailer(t), t) == (demand_retailer + where(t >= step_time_retailer, demand_step_retailer, 0))*theta_retailer - theta_retailer*expected_retailer(t)
retailer.supply_line: diff(supply_line_retailer(t), t) == -3*_material_delay3_0_retailer_supply_line_2/shipping_delay_retailer + where(0 > -(beta_retailer*supply_line_retailer(t) - desired_retailer - backlog_retailer(t) + inventory_retailer(t))*alpha_retailer + expected_retailer(t), 0, -(beta_retailer*supply_line_retailer(t) - desired_retailer - backlog_retailer(t) + inventory_retailer(t))*alpha_retailer + expected_retailer(t))
wholesaler.inventory: diff(inventory_wholesaler(t), t) == 3*_material_delay3_0_distributor_to_wholesaler_2/shipping_delay_wholesaler - where(inventory_wholesaler(t) < backlog_wholesaler(t), inventory_wholesaler(t), backlog_wholesaler(t))/dispatch_time_wholesaler
wholesaler.in_transit: diff(in_transit_wholesaler(t), t) == -3*_material_delay3_0_distributor_to_wholesaler_2/shipping_delay_wholesaler + where(inventory_distributor(t) < backlog_distributor(t), inventory_distributor(t), backlog_distributor(t))/dispatch_time_distributor
wholesaler.backlog: diff(backlog_wholesaler(t), t) == 3*_material_delay3_0_wholesaler_backlog_2/mail_delay_retailer - where(inventory_wholesaler(t) < backlog_wholesaler(t), inventory_wholesaler(t), backlog_wholesaler(t))/dispatch_time_wholesaler
wholesaler.expected: diff(expected_wholesaler(t), t) == -theta_wholesaler*expected_wholesaler(t) + 3*_material_delay3_0_wholesaler_expected_2*theta_wholesaler/mail_delay_retailer
wholesaler.supply_line: diff(supply_line_wholesaler(t), t) == -3*_material_delay3_0_wholesaler_supply_line_2/shipping_delay_wholesaler + where(0 > -(beta_wholesaler*supply_line_wholesaler(t) - desired_wholesaler - backlog_wholesaler(t) + inventory_wholesaler(t))*alpha_wholesaler + expected_wholesaler(t), 0, -(beta_wholesaler*supply_line_wholesaler(t) - desired_wholesaler - backlog_wholesaler(t) + inventory_wholesaler(t))*alpha_wholesaler + expected_wholesaler(t))
distributor.inventory: diff(inventory_distributor(t), t) == 3*_material_delay3_0_factory_to_distributor_2/shipping_delay_distributor - where(inventory_distributor(t) < backlog_distributor(t), inventory_distributor(t), backlog_distributor(t))/dispatch_time_distributor
distributor.in_transit: diff(in_transit_distributor(t), t) == -3*_material_delay3_0_factory_to_distributor_2/shipping_delay_distributor + where(inventory_factory(t) < backlog_factory(t), inventory_factory(t), backlog_factory(t))/dispatch_time_factory
distributor.backlog: diff(backlog_distributor(t), t) == 3*_material_delay3_0_distributor_backlog_2/mail_delay_wholesaler - where(inventory_distributor(t) < backlog_distributor(t), inventory_distributor(t), backlog_distributor(t))/dispatch_time_distributor
distributor.expected: diff(expected_distributor(t), t) == -theta_distributor*expected_distributor(t) + 3*_material_delay3_0_distributor_expected_2*theta_distributor/mail_delay_wholesaler
distributor.supply_line: diff(supply_line_distributor(t), t) == -3*_material_delay3_0_distributor_supply_line_2/shipping_delay_distributor + where(0 > -(beta_distributor*supply_line_distributor(t) - desired_distributor - backlog_distributor(t) + inventory_distributor(t))*alpha_distributor + expected_distributor(t), 0, -(beta_distributor*supply_line_distributor(t) - desired_distributor - backlog_distributor(t) + inventory_distributor(t))*alpha_distributor + expected_distributor(t))
factory.inventory: diff(inventory_factory(t), t) == 3*_material_delay3_0_brewery_2/production_delay_factory - where(inventory_factory(t) < backlog_factory(t), inventory_factory(t), backlog_factory(t))/dispatch_time_factory
factory.in_transit: diff(in_transit_factory(t), t) == -3*_material_delay3_0_brewery_2/production_delay_factory + where(0 > -(beta_factory*supply_line_factory(t) - desired_factory - backlog_factory(t) + inventory_factory(t))*alpha_factory + expected_factory(t), 0, -(beta_factory*supply_line_factory(t) - desired_factory - backlog_factory(t) + inventory_factory(t))*alpha_factory + expected_factory(t))
factory.backlog: diff(backlog_factory(t), t) == 3*_material_delay3_0_factory_backlog_2/mail_delay_distributor - where(inventory_factory(t) < backlog_factory(t), inventory_factory(t), backlog_factory(t))/dispatch_time_factory
factory.expected: diff(expected_factory(t), t) == -theta_factory*expected_factory(t) + 3*_material_delay3_0_factory_expected_2*theta_factory/mail_delay_distributor
factory.supply_line: diff(supply_line_factory(t), t) == -3*_material_delay3_0_factory_supply_line_2/production_delay_factory + where(0 > -(beta_factory*supply_line_factory(t) - desired_factory - backlog_factory(t) + inventory_factory(t))*alpha_factory + expected_factory(t), 0, -(beta_factory*supply_line_factory(t) - desired_factory - backlog_factory(t) + inventory_factory(t))*alpha_factory + expected_factory(t))