A full experiment: rolling forecasts, policies and extensions on M5 data¶

This notebook puts the pieces of the earlier notebooks together in one experiment on real demand. Every forecast is refitted at every review, the dated targets feed the policies through policy_schedule, and the runs include supplier rules, a callback, an unreliable supplier and a shelf life. The experiment answers two questions with numbers:

• What does each extension change? Starting from a plain run, we add one extension at a time and measure service and cost after each step.
• Which forecast and policy work best when all extensions are on? One run_comparison call replays three forecasts with two policy rules on the same demand and the same operating setting.

M5 supplies only item identifiers, dates and daily unit sales, which we use as demand. Lead time, review period, supplier behaviour, shelf life, opening stock and costs are assumptions made for this example.

In [1]:

import numpy as np
import pandas as pd

# Show every number with two decimals.
pd.options.display.float_format = "{:.2f}".format

import numpy as np import pandas as pd # Show every number with two decimals. pd.options.display.float_format = "{:.2f}".format

1. Demand: 30 M5 items¶

The prepared sample (data/m5/sample_1.csv, source and terms in data/m5/README.md) has one row per item and day, with demand in y. The simulation opens on 31 January 2016 and replays the 42 days from 1 February to 13 March. We keep the 30 items that sold on the most days in the 180 days up to the opening date, so the choice uses history only.

In [2]:

data = pd.read_csv("data/m5/sample_1.csv", parse_dates=["date"])

opening_date = pd.Timestamp("2016-01-31")
n_days = 42
run_end = opening_date + pd.Timedelta(days=n_days)

history_pool = data[
    data["date"].between(opening_date - pd.Timedelta(days=179), opening_date)
]
selling_days = (
    history_pool.assign(sold=history_pool["y"].gt(0)).groupby("unique_id")["sold"].sum()
)
items = selling_days.sort_values(ascending=False, kind="stable").index[:30].tolist()

# Everything observed for the 30 items, and the 42 days the runs replay.
observed = data.loc[data["unique_id"].isin(items), ["unique_id", "date", "y"]]
observed = observed.sort_values(["unique_id", "date"]).reset_index(drop=True)
demand = observed[observed["date"].between(opening_date, run_end, inclusive="right")]

data = pd.read_csv("data/m5/sample_1.csv", parse_dates=["date"]) opening_date = pd.Timestamp("2016-01-31") n_days = 42 run_end = opening_date + pd.Timedelta(days=n_days) history_pool = data[ data["date"].between(opening_date - pd.Timedelta(days=179), opening_date) ] selling_days = ( history_pool.assign(sold=history_pool["y"].gt(0)).groupby("unique_id")["sold"].sum() ) items = selling_days.sort_values(ascending=False, kind="stable").index[:30].tolist() # Everything observed for the 30 items, and the 42 days the runs replay. observed = data.loc[data["unique_id"].isin(items), ["unique_id", "date", "y"]] observed = observed.sort_values(["unique_id", "date"]).reset_index(drop=True) demand = observed[observed["date"].between(opening_date, run_end, inclusive="right")]

In [3]:

import matplotlib.dates as mdates
import matplotlib.pyplot as plt

DEMAND_COLOR = "#4C78A8"
RUN_SPAN_COLOR = "#F2CF5B"

# Four items across the selling-days ranking: 1st, 10th, 20th and 30th.
example_items = [items[i] for i in [0, 9, 19, 29]]
recent = observed[observed["date"].gt(opening_date - pd.Timedelta(days=42))]

fig, axes = plt.subplots(
    2,
    2,
    figsize=(11, 5.5),
    sharex=True,
    layout="constrained",
)
for ax, item in zip(axes.flat, example_items):
    rows = recent[recent["unique_id"].eq(item) & recent["date"].le(run_end)]
    ax.axvspan(
        opening_date + pd.Timedelta(days=1),
        run_end,
        color=RUN_SPAN_COLOR,
        alpha=0.14,
    )
    ax.plot(
        rows["date"],
        rows["y"],
        marker="o",
        markersize=2.5,
        linewidth=1.2,
        color=DEMAND_COLOR,
    )
    ax.set(title=item, ylabel="Units")
    ax.set_ylim(bottom=0)
    ax.grid(axis="y", alpha=0.2)
    ax.spines[["top", "right"]].set_visible(False)
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%d %b"))
fig.suptitle(
    "Daily demand of four items; the shaded span is the 42 days the runs replay"
)
plt.show()

import matplotlib.dates as mdates import matplotlib.pyplot as plt DEMAND_COLOR = "#4C78A8" RUN_SPAN_COLOR = "#F2CF5B" # Four items across the selling-days ranking: 1st, 10th, 20th and 30th. example_items = [items[i] for i in [0, 9, 19, 29]] recent = observed[observed["date"].gt(opening_date - pd.Timedelta(days=42))] fig, axes = plt.subplots( 2, 2, figsize=(11, 5.5), sharex=True, layout="constrained", ) for ax, item in zip(axes.flat, example_items): rows = recent[recent["unique_id"].eq(item) & recent["date"].le(run_end)] ax.axvspan( opening_date + pd.Timedelta(days=1), run_end, color=RUN_SPAN_COLOR, alpha=0.14, ) ax.plot( rows["date"], rows["y"], marker="o", markersize=2.5, linewidth=1.2, color=DEMAND_COLOR, ) ax.set(title=item, ylabel="Units") ax.set_ylim(bottom=0) ax.grid(axis="y", alpha=0.2) ax.spines[["top", "right"]].set_visible(False) ax.xaxis.set_major_formatter(mdates.DateFormatter("%d %b")) fig.suptitle( "Daily demand of four items; the shaded span is the 42 days the runs replay" ) plt.show()

What to notice: the 30 items differ widely. FOODS_1_085_CA_1, the top item, sells every day of the run, between 4 and 36 units. FOODS_1_051_CA_1, the 30th, sells at most 5 units a day and nothing on 20 of the 42 days. FOODS_1_064_CA_1 sells nothing at all from 17 February to 3 March. One forecast, one policy and one setting have to serve all of them.

2. Rolling forecasts: a new target at every review¶

Every run reviews stock every R = 2 days and orders with a lead time of L = 2 days, so each order protects L + R = 4 days of demand. At every review the forecast is refitted on the demand known up to the day before, and gives a target: the 95% quantile of the next four days' total demand, rounded up to whole units. With 21 reviews and 30 items, each forecast is recomputed 630 times. We compare three forecasts:

• 28-day mean: the mean daily demand of the last 28 days, and the 95% quantile of a Poisson distribution with four times that mean;
• Empirical: the 95th percentile of all 4-day demand totals in the last 56 days;
• ETS: Smooth's ES("ANN") model and its 95% upper bound for the 4-day total (cumulative=True), as in Notebook 04c.

Each target row also keeps its mean (the expected 4-day total), which the (s,S) policy uses in section 3, and end, the last day the target covers.

In [4]:

import math

from smooth import ES

lead_time, review_period = 2, 2
horizon = lead_time + review_period
probability = 0.95
origins = pd.date_range(
    opening_date,
    periods=n_days // review_period,
    freq=f"{review_period}D",
)


def poisson_quantile(probability, mean):
    """Smallest k with P(N <= k) >= probability, for N ~ Poisson(mean)."""
    k, pmf = 0, math.exp(-mean)
    cdf = pmf
    while cdf < probability:
        k += 1
        pmf *= mean / k
        cdf += pmf
    return float(k)


target_rows = []
for origin in origins:
    known = observed[observed["date"].le(origin)]
    for item, item_history in known.groupby("unique_id"):
        y = item_history["y"].astype(float).reset_index(drop=True)

        recent_mean = y.iloc[-28:].mean()
        four_day_totals = y.iloc[-56:].rolling(horizon).sum().dropna()
        model = ES("ANN")
        model.fit(y)
        ets = model.predict(
            h=horizon,
            interval="prediction",
            level=probability,
            side="upper",
            cumulative=True,
        )

        for forecast, quantile, mean in [
            (
                "28-day mean",
                poisson_quantile(probability, horizon * recent_mean),
                horizon * recent_mean,
            ),
            (
                "Empirical",
                np.quantile(four_day_totals, probability),
                four_day_totals.mean(),
            ),
            ("ETS", ets.upper.iloc[0, 0], ets.mean.iloc[0]),
        ]:
            target_rows.append(
                {
                    "unique_id": item,
                    "forecast": forecast,
                    "origin": origin,
                    "target": float(math.ceil(quantile)),
                    "mean": float(mean),
                }
            )

targets = pd.DataFrame(target_rows)
targets["end"] = targets["origin"] + pd.Timedelta(days=horizon)

import math from smooth import ES lead_time, review_period = 2, 2 horizon = lead_time + review_period probability = 0.95 origins = pd.date_range( opening_date, periods=n_days // review_period, freq=f"{review_period}D", ) def poisson_quantile(probability, mean): """Smallest k with P(N <= k) >= probability, for N ~ Poisson(mean).""" k, pmf = 0, math.exp(-mean) cdf = pmf while cdf < probability: k += 1 pmf *= mean / k cdf += pmf return float(k) target_rows = [] for origin in origins: known = observed[observed["date"].le(origin)] for item, item_history in known.groupby("unique_id"): y = item_history["y"].astype(float).reset_index(drop=True) recent_mean = y.iloc[-28:].mean() four_day_totals = y.iloc[-56:].rolling(horizon).sum().dropna() model = ES("ANN") model.fit(y) ets = model.predict( h=horizon, interval="prediction", level=probability, side="upper", cumulative=True, ) for forecast, quantile, mean in [ ( "28-day mean", poisson_quantile(probability, horizon * recent_mean), horizon * recent_mean, ), ( "Empirical", np.quantile(four_day_totals, probability), four_day_totals.mean(), ), ("ETS", ets.upper.iloc[0, 0], ets.mean.iloc[0]), ]: target_rows.append( { "unique_id": item, "forecast": forecast, "origin": origin, "target": float(math.ceil(quantile)), "mean": float(mean), } ) targets = pd.DataFrame(target_rows) targets["end"] = targets["origin"] + pd.Timedelta(days=horizon)

In [5]:

FORECAST_COLORS = {"28-day mean": "#E45756", "Empirical": "#54A24B", "ETS": "#4C78A8"}
target_item = example_items[0]

# The demand each target had to cover: the four days after its origin.
item_demand = demand[demand["unique_id"].eq(target_item)].set_index("date")["y"]
covered = pd.Series(
    [
        item_demand.loc[
            origin + pd.Timedelta(days=1) : origin + pd.Timedelta(days=horizon)
        ].sum()
        for origin in origins
    ],
    index=origins,
)

fig, ax = plt.subplots(figsize=(11, 4), layout="constrained")
for forecast, color in FORECAST_COLORS.items():
    rows = targets[
        targets["unique_id"].eq(target_item) & targets["forecast"].eq(forecast)
    ]
    ax.step(
        rows["origin"],
        rows["target"],
        where="post",
        color=color,
        linewidth=1.5,
        label=f"{forecast} target (mean {rows['target'].mean():.1f})",
    )
ax.plot(
    covered.index[:-1],
    covered.iloc[:-1],
    marker="o",
    markersize=4,
    linewidth=1,
    color="black",
    label="Demand in the four days after the origin",
)
ax.set(
    title=f"{target_item}: 95% targets for the next four days, refitted every two days",
    xlabel="Forecast origin",
    ylabel="Units",
)
ax.set_ylim(bottom=0)
ax.set_xticks(origins[::2])
ax.xaxis.set_major_formatter(mdates.DateFormatter("%d %b"))
ax.grid(axis="y", alpha=0.2)
ax.spines[["top", "right"]].set_visible(False)
fig.legend(loc="outside lower center", ncol=2, frameon=False)
plt.show()

FORECAST_COLORS = {"28-day mean": "#E45756", "Empirical": "#54A24B", "ETS": "#4C78A8"} target_item = example_items[0] # The demand each target had to cover: the four days after its origin. item_demand = demand[demand["unique_id"].eq(target_item)].set_index("date")["y"] covered = pd.Series( [ item_demand.loc[ origin + pd.Timedelta(days=1) : origin + pd.Timedelta(days=horizon) ].sum() for origin in origins ], index=origins, ) fig, ax = plt.subplots(figsize=(11, 4), layout="constrained") for forecast, color in FORECAST_COLORS.items(): rows = targets[ targets["unique_id"].eq(target_item) & targets["forecast"].eq(forecast) ] ax.step( rows["origin"], rows["target"], where="post", color=color, linewidth=1.5, label=f"{forecast} target (mean {rows['target'].mean():.1f})", ) ax.plot( covered.index[:-1], covered.iloc[:-1], marker="o", markersize=4, linewidth=1, color="black", label="Demand in the four days after the origin", ) ax.set( title=f"{target_item}: 95% targets for the next four days, refitted every two days", xlabel="Forecast origin", ylabel="Units", ) ax.set_ylim(bottom=0) ax.set_xticks(origins[::2]) ax.xaxis.set_major_formatter(mdates.DateFormatter("%d %b")) ax.grid(axis="y", alpha=0.2) ax.spines[["top", "right"]].set_visible(False) fig.legend(loc="outside lower center", ncol=2, frameon=False) plt.show()

What to notice: the three forecasts react differently to the same demand. For FOODS_1_085_CA_1:

• ETS follows the latest demand: its target jumps from 71 to 99 units after the busy days in early February and falls back to the 70s in late February.
• Empirical changes only when a large 4-day total enters or leaves its 56-day window: it climbs from 61 to 84 by 16 February and then stays at 84.
• 28-day mean is the lowest target (mean 62.1). The Poisson quantile assumes that demand varies as much as its mean, and this item varies more. Its target is exceeded at 6 of the 20 origins, against 4 for the other two.

Every point of these lines is a separate fit that uses only the demand known at its origin.

3. Dated targets become a policy_schedule¶

Each review's target table is fitted into its own policy. The end column dates the targets (date_column="end"), so each fitted policy knows its forecast origin and checks, when it is used, that its window matches the decision. A policy_schedule is a dict {period: fitted policy}: from that period on, the run uses the new policy. Period 0 is the first day of the run, 1 February, and its policy comes from the forecast made on 31 January; period t uses the forecast made t days after the opening date.

We fit two policy rules on the same targets:

• (R,S), OrderUpToPolicy: at every review, order up to S = target.
• (s,S), ReorderPointPolicy (Notebook 05g): at a review, order only if the inventory position is at or below s = target, and then order up to S = s + the expected demand of one review period (half the 4-day mean, rounded up). It orders less often, but more at a time.

In [6]:

from stockcast.policies import OrderUpToPolicy, ReorderPointPolicy


def fit_policy(rule, forecast, origin):
    """Fit one rule on the targets one forecast made at one origin."""
    rows = targets[targets["forecast"].eq(forecast) & targets["origin"].eq(origin)]
    if rule == "(R,S)":
        policy = OrderUpToPolicy(
            lead_time=lead_time,
            review_period=review_period,
            freq="D",
            service_level=probability,
            allow_backorders=False,
            date_column="end",
        )
        return policy.fit(rows, target_column="target")
    rows = rows.assign(
        S=rows["target"] + np.ceil(rows["mean"] * review_period / horizon)
    )
    policy = ReorderPointPolicy(
        lead_time=lead_time,
        review_period=review_period,
        freq="D",
        service_level=probability,
        allow_backorders=False,
        date_column="end",
    )
    return policy.fit(rows, reorder_point_column="target", order_up_to_column="S")


first_policies, schedules = {}, {}
for rule in ["(R,S)", "(s,S)"]:
    for forecast in FORECAST_COLORS:
        label = f"{rule} {forecast}"
        snapshots = {
            (origin - opening_date).days: fit_policy(rule, forecast, origin)
            for origin in origins
        }
        first_policies[label] = snapshots.pop(0)
        schedules[label] = snapshots

# The first policy and the first two schedule entries of one branch.
first_three = [(0, first_policies["(R,S) ETS"])] + list(schedules["(R,S) ETS"].items())[
    :2
]
pd.DataFrame(
    [
        {
            "period": period,
            "decision date": opening_date + pd.Timedelta(days=period + 1),
            "forecast origin": snapshot.get_target_metadata()["forecast_origin"][:10],
            "target end": snapshot.get_target_metadata()["target_end_date"][:10],
        }
        for period, snapshot in first_three
    ]
)

from stockcast.policies import OrderUpToPolicy, ReorderPointPolicy def fit_policy(rule, forecast, origin): """Fit one rule on the targets one forecast made at one origin.""" rows = targets[targets["forecast"].eq(forecast) & targets["origin"].eq(origin)] if rule == "(R,S)": policy = OrderUpToPolicy( lead_time=lead_time, review_period=review_period, freq="D", service_level=probability, allow_backorders=False, date_column="end", ) return policy.fit(rows, target_column="target") rows = rows.assign( S=rows["target"] + np.ceil(rows["mean"] * review_period / horizon) ) policy = ReorderPointPolicy( lead_time=lead_time, review_period=review_period, freq="D", service_level=probability, allow_backorders=False, date_column="end", ) return policy.fit(rows, reorder_point_column="target", order_up_to_column="S") first_policies, schedules = {}, {} for rule in ["(R,S)", "(s,S)"]: for forecast in FORECAST_COLORS: label = f"{rule} {forecast}" snapshots = { (origin - opening_date).days: fit_policy(rule, forecast, origin) for origin in origins } first_policies[label] = snapshots.pop(0) schedules[label] = snapshots # The first policy and the first two schedule entries of one branch. first_three = [(0, first_policies["(R,S) ETS"])] + list(schedules["(R,S) ETS"].items())[ :2 ] pd.DataFrame( [ { "period": period, "decision date": opening_date + pd.Timedelta(days=period + 1), "forecast origin": snapshot.get_target_metadata()["forecast_origin"][:10], "target end": snapshot.get_target_metadata()["target_end_date"][:10], } for period, snapshot in first_three ] )

Out[6]:

	period	decision date	forecast origin	target end
0	0	2016-02-01	2016-01-31	2016-02-04
1	2	2016-02-03	2016-02-02	2016-02-06
2	4	2016-02-05	2016-02-04	2016-02-08

What to notice: the decision on 1 February (period 0) uses the forecast made on 31 January, and its target covers 1 to 4 February. Each later policy decides one day after its forecast origin and covers the next four days. The first policy of each branch is passed as policy, and the other 20 as policy_schedule. A schedule entry under the wrong period stops the run, for example: policy forecast_origin 2016-02-02 00:00:00 (derived from the 'end' column) must equal decision information date 2016-02-04 00:00:00.

4. The baseline run¶

All runs start from the same state: each item holds four days of its mean demand over the 28 days before opening, received as one delivery on 31 January, with nothing on order. max_lead_time=4 leaves room in the pipeline for deliveries that arrive two days late (section 5).

The baseline is the plainest run: ETS targets with the (R,S) rule, and no extensions. Every order arrives in full after two days, nothing expires, and any quantity can be ordered.

In [7]:

from stockcast.core import InventoryStateDataFrame, SimulationEngine

mean_before_opening = (
    observed[
        observed["date"].gt(opening_date - pd.Timedelta(days=28))
        & observed["date"].le(opening_date)
    ]
    .groupby("unique_id")["y"]
    .mean()
)
opening_stock = pd.DataFrame(
    {
        "unique_id": items,
        "on_hand": (4 * mean_before_opening.loc[items])
        .round()
        .clip(lower=1)
        .to_numpy(),
    }
)
opening = InventoryStateDataFrame.from_observed(
    opening_stock,
    opening_date=opening_date,
    max_lead_time=4,
)

policy = first_policies["(R,S) ETS"]
schedule = schedules["(R,S) ETS"]
runs = {}
runs["Baseline"] = SimulationEngine().run(
    policy=policy,
    policy_schedule=schedule,
    demand_source=demand,
    inventory=opening,
    freq="D",
)

from stockcast.core import InventoryStateDataFrame, SimulationEngine mean_before_opening = ( observed[ observed["date"].gt(opening_date - pd.Timedelta(days=28)) & observed["date"].le(opening_date) ] .groupby("unique_id")["y"] .mean() ) opening_stock = pd.DataFrame( { "unique_id": items, "on_hand": (4 * mean_before_opening.loc[items]) .round() .clip(lower=1) .to_numpy(), } ) opening = InventoryStateDataFrame.from_observed( opening_stock, opening_date=opening_date, max_lead_time=4, ) policy = first_policies["(R,S) ETS"] schedule = schedules["(R,S) ETS"] runs = {} runs["Baseline"] = SimulationEngine().run( policy=policy, policy_schedule=schedule, demand_source=demand, inventory=opening, freq="D", )

5. Adding the extensions one at a time¶

Each step below adds one extension on top of the previous ones and repeats the run with the same policy, schedule, demand and opening state. Only the arguments of SimulationEngine().run change. Each extension has its own notebook; here we only measure what it does.

Supplier rules (Notebook 05c): each order line must be at least 4 units and a multiple of 4 (a case). mode="adjust" rounds a proposal up to the next valid quantity.

In [8]:

from stockcast.core import MinimumOrderQuantity, OrderingConstraints, OrderMultiple

supplier_rules = OrderingConstraints(
    [
        MinimumOrderQuantity(4.0, mode="adjust"),
        OrderMultiple(4.0, mode="adjust"),
    ]
)
runs["+ Supplier rules"] = SimulationEngine().run(
    policy=policy,
    policy_schedule=schedule,
    demand_source=demand,
    inventory=opening,
    freq="D",
    order_constraints=supplier_rules,
)

from stockcast.core import MinimumOrderQuantity, OrderingConstraints, OrderMultiple supplier_rules = OrderingConstraints( [ MinimumOrderQuantity(4.0, mode="adjust"), OrderMultiple(4.0, mode="adjust"), ] ) runs["+ Supplier rules"] = SimulationEngine().run( policy=policy, policy_schedule=schedule, demand_source=demand, inventory=opening, freq="D", order_constraints=supplier_rules, )

A whole-order minimum, as a callback (Notebook 08): the supplier only delivers an order of at least 150 units across all items. A constraint sees one item at a time, so this rule needs a callback: on_after_prediction sees the whole order after the policy proposes it, and sets every line to zero when the total is too small. The skipped quantities are not carried forward; the next review proposes again from the current stock.

In [9]:

from stockcast.core import OrderAdjustmentResult, SimulationCallback


class WholeOrderMinimum(SimulationCallback):
    """Hold the whole order when it adds up to less than `minimum` units."""

    def __init__(self, minimum):
        self.minimum = minimum

    def on_after_prediction(self, decision, context):
        orders = decision.get_dataframe()
        ordered = orders[orders["order_quantity"] > 0]
        if ordered.empty or ordered["order_quantity"].sum() >= self.minimum:
            return None  # nothing to hold
        held = pd.DataFrame(
            {
                "unique_id": ordered["unique_id"],
                "order_quantity": 0.0,
                "reason": "order below the whole-order minimum",
            }
        )
        return OrderAdjustmentResult(held)


whole_order_minimum = WholeOrderMinimum(minimum=150)
runs["+ Whole-order minimum"] = SimulationEngine().run(
    policy=policy,
    policy_schedule=schedule,
    demand_source=demand,
    inventory=opening,
    freq="D",
    order_constraints=supplier_rules,
    callbacks=[whole_order_minimum],
)

from stockcast.core import OrderAdjustmentResult, SimulationCallback class WholeOrderMinimum(SimulationCallback): """Hold the whole order when it adds up to less than `minimum` units.""" def __init__(self, minimum): self.minimum = minimum def on_after_prediction(self, decision, context): orders = decision.get_dataframe() ordered = orders[orders["order_quantity"] > 0] if ordered.empty or ordered["order_quantity"].sum() >= self.minimum: return None # nothing to hold held = pd.DataFrame( { "unique_id": ordered["unique_id"], "order_quantity": 0.0, "reason": "order below the whole-order minimum", } ) return OrderAdjustmentResult(held) whole_order_minimum = WholeOrderMinimum(minimum=150) runs["+ Whole-order minimum"] = SimulationEngine().run( policy=policy, policy_schedule=schedule, demand_source=demand, inventory=opening, freq="D", order_constraints=supplier_rules, callbacks=[whole_order_minimum], )

An unreliable supplier (Notebook 05f): a DeliveryOutcome decides what arrives when a delivery is due. Here 15% of deliveries arrive two days late, and 10% arrive with only half of the units; the rest never comes. The SupplyModel seeds the supplier's random draws, so they repeat with the same sequence of calls. Different ordering patterns can assign those draws to different deliveries: each branch uses the same lateness and short-fill probabilities, but the affected deliveries need not match. Stockcast warns that deliveries now arrive at times other than the policy's lead time, because the targets were sized for an on-time supplier; that is the situation we want to measure, so one line silences the warning.

In [10]:

import warnings

from stockcast.core import DeliveryOutcome, Supplier, SupplyModel

# Expected: the targets were sized for an on-time supplier (see above).
warnings.filterwarnings("ignore", message="supplier deliveries arrive at times other")


class UnreliableSupplier(DeliveryOutcome):
    """Some deliveries arrive late, some arrive short.

    With probability p_late a delivery slips late_by days; with probability
    p_short only short_fill of it arrives (whole units) and the rest never comes.
    A delivery that already slipped once arrives in full.
    """

    def __init__(self, p_late, late_by, p_short, short_fill):
        self.p_late = p_late
        self.late_by = late_by
        self.p_short = p_short
        self.short_fill = short_fill

    def resolve(self, due, context):
        draw = context.rng.random(len(due))
        first_try = (due["due_period"] == due["scheduled_due_period"]).to_numpy()
        late = first_try & (draw < self.p_late)
        short = first_try & ~late & (draw < self.p_late + self.p_short)
        quantity = due["quantity"].to_numpy()
        received = np.where(short, np.floor(self.short_fill * quantity), quantity)
        return due.assign(
            received_quantity=np.where(late, 0.0, received),
            delayed_quantity=np.where(late, quantity, 0.0),
            delay_periods=self.late_by,
        )


supply = SupplyModel(
    [
        Supplier(
            "supplier",
            lead_time=lead_time,
            delivery=UnreliableSupplier(
                p_late=0.15,
                late_by=2,
                p_short=0.10,
                short_fill=0.5,
            ),
        )
    ],
    random_seed=7,
)
runs["+ Unreliable supplier"] = SimulationEngine().run(
    policy=policy,
    policy_schedule=schedule,
    demand_source=demand,
    inventory=opening,
    freq="D",
    order_constraints=supplier_rules,
    callbacks=[whole_order_minimum],
    supply=supply,
)

import warnings from stockcast.core import DeliveryOutcome, Supplier, SupplyModel # Expected: the targets were sized for an on-time supplier (see above). warnings.filterwarnings("ignore", message="supplier deliveries arrive at times other") class UnreliableSupplier(DeliveryOutcome): """Some deliveries arrive late, some arrive short. With probability p_late a delivery slips late_by days; with probability p_short only short_fill of it arrives (whole units) and the rest never comes. A delivery that already slipped once arrives in full. """ def __init__(self, p_late, late_by, p_short, short_fill): self.p_late = p_late self.late_by = late_by self.p_short = p_short self.short_fill = short_fill def resolve(self, due, context): draw = context.rng.random(len(due)) first_try = (due["due_period"] == due["scheduled_due_period"]).to_numpy() late = first_try & (draw < self.p_late) short = first_try & ~late & (draw < self.p_late + self.p_short) quantity = due["quantity"].to_numpy() received = np.where(short, np.floor(self.short_fill * quantity), quantity) return due.assign( received_quantity=np.where(late, 0.0, received), delayed_quantity=np.where(late, quantity, 0.0), delay_periods=self.late_by, ) supply = SupplyModel( [ Supplier( "supplier", lead_time=lead_time, delivery=UnreliableSupplier( p_late=0.15, late_by=2, p_short=0.10, short_fill=0.5, ), ) ], random_seed=7, ) runs["+ Unreliable supplier"] = SimulationEngine().run( policy=policy, policy_schedule=schedule, demand_source=demand, inventory=opening, freq="D", order_constraints=supplier_rules, callbacks=[whole_order_minimum], supply=supply, )

A shelf life (Notebook 05e, Notebook 07): the ShelfLife process sells the oldest delivery first and removes units 10 days after they were received. It needs the date of the opening stock, so we pass the opening stock as one delivery received on 31 January (opening_lots).

In [11]:

from stockcast.core import ShelfLife

opening_deliveries = opening_stock.rename(columns={"on_hand": "quantity"}).assign(
    received_date=opening_date
)
shelf_life = ShelfLife(shelf_life_days=10, opening_lots=opening_deliveries)
runs["+ Shelf life"] = SimulationEngine().run(
    policy=policy,
    policy_schedule=schedule,
    demand_source=demand,
    inventory=opening,
    freq="D",
    order_constraints=supplier_rules,
    callbacks=[whole_order_minimum],
    supply=supply,
    processes=[shelf_life],
)

from stockcast.core import ShelfLife opening_deliveries = opening_stock.rename(columns={"on_hand": "quantity"}).assign( received_date=opening_date ) shelf_life = ShelfLife(shelf_life_days=10, opening_lots=opening_deliveries) runs["+ Shelf life"] = SimulationEngine().run( policy=policy, policy_schedule=schedule, demand_source=demand, inventory=opening, freq="D", order_constraints=supplier_rules, callbacks=[whole_order_minimum], supply=supply, processes=[shelf_life], )

6. What each extension changes¶

InventoryEvaluator scores every run with the same metrics and costs. The costs are assumptions: 0.02 per unit held per day, 2.00 per unit of lost sales, 0.50 per expired unit and 0.50 per order line. total_cost is the sum of the components listed in cost_components.

In [12]:

from IPython.display import display

from stockcast.evaluation import (
    InventoryEvaluator,
    avg_on_hand,
    fill_rate,
    holding_cost,
    lost_sales_units,
    ordering_cost,
    shortage_cost,
    sku_order_line_count,
    total_cost,
    waste_cost,
)

costs = {
    "cost_components": ["holding", "shortage", "waste", "ordering"],
    "holding_cost_per_unit_period": 0.02,
    "shortage_cost_per_unit": 2.00,
    "waste_cost_per_unit": 0.50,
    "order_cost_per_sku_line": 0.50,
    "order_cost_per_unit": 0.00,
}
metrics = [
    fill_rate,
    lost_sales_units,
    avg_on_hand,
    sku_order_line_count,
    holding_cost,
    shortage_cost,
    waste_cost,
    ordering_cost,
    total_cost,
]


def score(results):
    """One row of metrics per run."""
    rows = {}
    for label, result in results.items():
        evaluator = InventoryEvaluator()
        evaluator.fit(result)
        rows[label] = evaluator.evaluate(metrics, groupby=[], context=costs).iloc[0]
    return pd.DataFrame(rows).T


extension_scores = score(runs)
# Three decimals, so that fill rates close to 1 stay readable.
with pd.option_context("display.float_format", "{:.3f}".format):
    display(extension_scores)

from IPython.display import display from stockcast.evaluation import ( InventoryEvaluator, avg_on_hand, fill_rate, holding_cost, lost_sales_units, ordering_cost, shortage_cost, sku_order_line_count, total_cost, waste_cost, ) costs = { "cost_components": ["holding", "shortage", "waste", "ordering"], "holding_cost_per_unit_period": 0.02, "shortage_cost_per_unit": 2.00, "waste_cost_per_unit": 0.50, "order_cost_per_sku_line": 0.50, "order_cost_per_unit": 0.00, } metrics = [ fill_rate, lost_sales_units, avg_on_hand, sku_order_line_count, holding_cost, shortage_cost, waste_cost, ordering_cost, total_cost, ] def score(results): """One row of metrics per run.""" rows = {} for label, result in results.items(): evaluator = InventoryEvaluator() evaluator.fit(result) rows[label] = evaluator.evaluate(metrics, groupby=[], context=costs).iloc[0] return pd.DataFrame(rows).T extension_scores = score(runs) # Three decimals, so that fill rates close to 1 stay readable. with pd.option_context("display.float_format", "{:.3f}".format): display(extension_scores)

	fill_rate	lost_sales_units	avg_on_hand	sku_order_line_count	holding_cost	shortage_cost	waste_cost	ordering_cost	total_cost
Baseline	0.942	250.000	17.578	443.000	442.960	500.000	0.000	221.500	1164.460
+ Supplier rules	0.949	219.000	19.081	378.000	480.840	438.000	0.000	189.000	1107.840
+ Whole-order minimum	0.921	339.000	17.583	296.000	443.100	678.000	0.000	148.000	1269.100
+ Unreliable supplier	0.890	471.000	17.191	303.000	433.220	942.000	0.000	151.500	1526.720
+ Shelf life	0.878	522.000	14.320	324.000	360.860	1044.000	269.000	162.000	1835.860

In [13]:

COST_COLORS = {
    "holding_cost": ("Holding", "#4C78A8"),
    "shortage_cost": ("Lost sales", "#E45756"),
    "waste_cost": ("Waste", "#B279A2"),
    "ordering_cost": ("Ordering", "#F58518"),
}


def plot_costs_and_fill(scores, title):
    """Stacked cost components (left) and fill rate (right), one bar per run."""
    labels = scores.index.tolist()
    fig, (cost_ax, fill_ax) = plt.subplots(
        1,
        2,
        figsize=(12, 4.5),
        gridspec_kw={"width_ratios": [3, 2]},
        layout="constrained",
    )
    bottom = np.zeros(len(labels))
    for column, (name, color) in COST_COLORS.items():
        values = scores[column].to_numpy(dtype=float)
        cost_ax.bar(labels, values, bottom=bottom, color=color, label=name)
        bottom += values
    for x, total in enumerate(bottom):
        cost_ax.annotate(
            f"{total:,.0f}",
            (x, total),
            ha="center",
            va="bottom",
            fontsize=9,
        )
    cost_ax.set(title="Total cost by component", ylabel="Cost")
    cost_ax.legend(frameon=False, ncol=4, loc="upper left")
    cost_ax.set_ylim(0, bottom.max() * 1.18)

    fill = scores["fill_rate"].to_numpy(dtype=float)
    fill_ax.bar(labels, fill, color="#9D9D9D")
    for x, value in enumerate(fill):
        fill_ax.annotate(
            f"{value:.3f}",
            (x, value),
            ha="center",
            va="bottom",
            fontsize=9,
        )
    fill_ax.set(title="Fill rate", ylim=(0, 1.08))
    for ax in (cost_ax, fill_ax):
        ax.tick_params(axis="x", rotation=25)
        ax.spines[["top", "right"]].set_visible(False)
        ax.grid(axis="y", alpha=0.2)
    fig.suptitle(title)
    plt.show()


plot_costs_and_fill(
    extension_scores,
    "Adding one extension at a time (ETS targets, (R,S) rule)",
)

COST_COLORS = { "holding_cost": ("Holding", "#4C78A8"), "shortage_cost": ("Lost sales", "#E45756"), "waste_cost": ("Waste", "#B279A2"), "ordering_cost": ("Ordering", "#F58518"), } def plot_costs_and_fill(scores, title): """Stacked cost components (left) and fill rate (right), one bar per run.""" labels = scores.index.tolist() fig, (cost_ax, fill_ax) = plt.subplots( 1, 2, figsize=(12, 4.5), gridspec_kw={"width_ratios": [3, 2]}, layout="constrained", ) bottom = np.zeros(len(labels)) for column, (name, color) in COST_COLORS.items(): values = scores[column].to_numpy(dtype=float) cost_ax.bar(labels, values, bottom=bottom, color=color, label=name) bottom += values for x, total in enumerate(bottom): cost_ax.annotate( f"{total:,.0f}", (x, total), ha="center", va="bottom", fontsize=9, ) cost_ax.set(title="Total cost by component", ylabel="Cost") cost_ax.legend(frameon=False, ncol=4, loc="upper left") cost_ax.set_ylim(0, bottom.max() * 1.18) fill = scores["fill_rate"].to_numpy(dtype=float) fill_ax.bar(labels, fill, color="#9D9D9D") for x, value in enumerate(fill): fill_ax.annotate( f"{value:.3f}", (x, value), ha="center", va="bottom", fontsize=9, ) fill_ax.set(title="Fill rate", ylim=(0, 1.08)) for ax in (cost_ax, fill_ax): ax.tick_params(axis="x", rotation=25) ax.spines[["top", "right"]].set_visible(False) ax.grid(axis="y", alpha=0.2) fig.suptitle(title) plt.show() plot_costs_and_fill( extension_scores, "Adding one extension at a time (ETS targets, (R,S) rule)", )

What to notice:

• Supplier rules improve this run: rounding every line up to a case adds stock (average on hand 17.58 → 19.08), so lost sales fall from 250 to 219 units, and fewer, larger lines cut the order lines from 443 to 378. Total cost: 1,164 → 1,108.
• The whole-order minimum holds some reviews entirely: order lines fall to 296, lost sales rise to 339 units, and the fill rate drops to 0.921.
• The unreliable supplier raises lost sales to 471 units (fill rate 0.890). The targets were sized for an on-time supplier, so a late or short delivery leaves the shelf short until the next one arrives.
• The shelf life adds waste (269) and lowers the average stock to 14.32, because expired units leave the shelf before they can be sold. The fill rate falls to 0.878.

Together, the extensions raise the total cost from 1,164 to 1,836 (+58%) and lower the fill rate from 0.942 to 0.878. From the whole-order minimum on, lost sales are the largest cost. Each step is measured on top of the earlier ones, so the size of a step depends on the order in which the extensions are added.

7. One item through the full setting¶

The last run has every extension on. The event table records, for each item and day:

• what the policy proposed (requested_order_quantity) and what was placed after the callback and the supplier rules (order_quantity);
• what arrived (received_units) and what a short delivery left out (supplier_shortfall_units);
• what expired (expired_units) and what demand was lost (shortage_units).

The figure follows one item through the 42 days.

In [14]:

TRAIL_ITEM = "FOODS_1_013_CA_1"

result = runs["+ Shelf life"]
events = result.to_event_frame()
trail = events[events["unique_id"].eq(TRAIL_ITEM)].set_index("date")
orders = result.to_order_frame()
late = orders[orders["unique_id"].eq(TRAIL_ITEM) & orders["delayed_quantity"].gt(0)]


def positive(column):
    """The column on the days it is above zero."""
    return trail[column].where(trail[column] > 0).dropna()


fig, (stock_ax, order_ax) = plt.subplots(
    2,
    1,
    figsize=(11, 6.5),
    sharex=True,
    layout="constrained",
)
stock_ax.step(
    trail.index,
    trail["ending_on_hand"],
    where="mid",
    color="#4C78A8",
    linewidth=1.5,
    label="Stock at the end of the day",
)
stock_ax.plot(
    trail.index,
    trail["demand"],
    marker="o",
    markersize=2.5,
    linewidth=1,
    color="#9D9D9D",
    label="Demand",
)
order_markers = [
    (stock_ax, positive("expired_units"), "X", "#B279A2", "Expired"),
    (stock_ax, positive("shortage_units"), "o", "#E45756", "Lost sales"),
    (
        order_ax,
        positive("requested_order_quantity"),
        "_",
        "black",
        "Proposed by the policy",
    ),
    (order_ax, positive("order_quantity"), "^", "#F58518", "Placed"),
    (order_ax, positive("received_units"), "v", "#54A24B", "Received"),
    (
        order_ax,
        late.set_index("due_date")["delayed_quantity"],
        "D",
        "#9D755D",
        "Due but late",
    ),
    (order_ax, positive("supplier_shortfall_units"), "x", "#E45756", "Never delivered"),
]
for ax, values, marker, color, label in order_markers:
    ax.plot(
        values.index,
        values,
        linestyle="none",
        marker=marker,
        markersize=12 if marker == "_" else 7,
        markeredgewidth=2,
        color=color,
        label=label,
        clip_on=False,
    )
stock_ax.set(title=f"{TRAIL_ITEM}: stock, expiry and lost sales", ylabel="Units")
order_ax.set(title="Orders and deliveries", ylabel="Units")
for ax in (stock_ax, order_ax):
    ax.set_ylim(bottom=0)
    ax.grid(axis="y", alpha=0.2)
    ax.spines[["top", "right"]].set_visible(False)
order_ax.xaxis.set_major_locator(mdates.WeekdayLocator(byweekday=mdates.MO))
order_ax.xaxis.set_major_formatter(mdates.DateFormatter("%d %b"))
fig.legend(loc="outside lower center", ncol=5, frameon=False)
plt.show()

TRAIL_ITEM = "FOODS_1_013_CA_1" result = runs["+ Shelf life"] events = result.to_event_frame() trail = events[events["unique_id"].eq(TRAIL_ITEM)].set_index("date") orders = result.to_order_frame() late = orders[orders["unique_id"].eq(TRAIL_ITEM) & orders["delayed_quantity"].gt(0)] def positive(column): """The column on the days it is above zero.""" return trail[column].where(trail[column] > 0).dropna() fig, (stock_ax, order_ax) = plt.subplots( 2, 1, figsize=(11, 6.5), sharex=True, layout="constrained", ) stock_ax.step( trail.index, trail["ending_on_hand"], where="mid", color="#4C78A8", linewidth=1.5, label="Stock at the end of the day", ) stock_ax.plot( trail.index, trail["demand"], marker="o", markersize=2.5, linewidth=1, color="#9D9D9D", label="Demand", ) order_markers = [ (stock_ax, positive("expired_units"), "X", "#B279A2", "Expired"), (stock_ax, positive("shortage_units"), "o", "#E45756", "Lost sales"), ( order_ax, positive("requested_order_quantity"), "_", "black", "Proposed by the policy", ), (order_ax, positive("order_quantity"), "^", "#F58518", "Placed"), (order_ax, positive("received_units"), "v", "#54A24B", "Received"), ( order_ax, late.set_index("due_date")["delayed_quantity"], "D", "#9D755D", "Due but late", ), (order_ax, positive("supplier_shortfall_units"), "x", "#E45756", "Never delivered"), ] for ax, values, marker, color, label in order_markers: ax.plot( values.index, values, linestyle="none", marker=marker, markersize=12 if marker == "_" else 7, markeredgewidth=2, color=color, label=label, clip_on=False, ) stock_ax.set(title=f"{TRAIL_ITEM}: stock, expiry and lost sales", ylabel="Units") order_ax.set(title="Orders and deliveries", ylabel="Units") for ax in (stock_ax, order_ax): ax.set_ylim(bottom=0) ax.grid(axis="y", alpha=0.2) ax.spines[["top", "right"]].set_visible(False) order_ax.xaxis.set_major_locator(mdates.WeekdayLocator(byweekday=mdates.MO)) order_ax.xaxis.set_major_formatter(mdates.DateFormatter("%d %b")) fig.legend(loc="outside lower center", ncol=5, frameon=False) plt.show()

What to notice: every lost unit of this item (28 in total) follows one of the extensions:

• Late delivery: the 8 units ordered on 1 February were due on 3 February and arrived on 5 February. The case rule had already rounded the policy's 7 up to 8.
• Expiry: the item sold nothing from 2 to 13 February, so 11 of the 12 units received on 5 February expired on 15 February, ten days later. The shelf was empty, and 9 units were lost on 15 and 16 February before the order placed on 15 February arrived.
• Held orders: the policy proposed 1, 5 and 9 units on 17 February, 23 February and 10 March, and the whole-order minimum held all three (black dashes with no orange triangle). After the hold on 10 March, 12 units were lost on 11 and 12 March.
• Short deliveries: on 23 February only 2 of 4 units arrived, and on 10 March only 6 of 12.

to_event_frame() gives one row per item and day. to_order_frame() gives one row per order and delivery, with its status, so the late delivery can be found there.

8. Forecasts and policies with every extension on¶

Now the operating setting is fixed and the forecast and policy change. run_comparison takes the six first policies as {label: policy} and their schedules as {label: schedule}, together with the extensions. Every branch replays the same demand from its own copy of the opening state; the callback and the shelf life are reset between branches, and the supplier draws come from the same seed. The branches therefore differ in their targets and their rule, and in which deliveries happen to be late or short: a branch that orders differently uses the supplier's draws for different deliveries.

In [15]:

comparison = SimulationEngine().run_comparison(
    first_policies,
    demand,
    opening,
    freq="D",
    policy_schedules=schedules,
    order_constraints=supplier_rules,
    callbacks=[whole_order_minimum],
    supply=supply,
    processes=[shelf_life],
)
comparison_scores = score({label: comparison[label] for label in first_policies})
with pd.option_context("display.float_format", "{:.3f}".format):
    display(comparison_scores)

comparison = SimulationEngine().run_comparison( first_policies, demand, opening, freq="D", policy_schedules=schedules, order_constraints=supplier_rules, callbacks=[whole_order_minimum], supply=supply, processes=[shelf_life], ) comparison_scores = score({label: comparison[label] for label in first_policies}) with pd.option_context("display.float_format", "{:.3f}".format): display(comparison_scores)

	fill_rate	lost_sales_units	avg_on_hand	sku_order_line_count	holding_cost	shortage_cost	waste_cost	ordering_cost	total_cost
(R,S) 28-day mean	0.798	864.000	8.915	332.000	224.660	1728.000	132.000	166.000	2250.660
(R,S) Empirical	0.904	411.000	15.182	379.000	382.580	822.000	319.500	189.500	1713.580
(R,S) ETS	0.878	522.000	14.320	324.000	360.860	1044.000	269.000	162.000	1835.860
(s,S) 28-day mean	0.892	461.000	13.290	256.000	334.920	922.000	212.000	128.000	1596.920
(s,S) Empirical	0.931	296.000	20.186	272.000	508.680	592.000	386.500	136.000	1623.180
(s,S) ETS	0.944	241.000	19.361	277.000	487.900	482.000	361.500	138.500	1469.900

In [16]:

plot_costs_and_fill(
    comparison_scores,
    "Three forecasts and two policy rules, every extension on",
)

plot_costs_and_fill( comparison_scores, "Three forecasts and two policy rules, every extension on", )

What to notice:

• The (R,S) ETS branch is the last run of section 5 (total cost 1,836, fill rate 0.878), repeated inside the comparison with the same result.
• The forecast matters most for (R,S). With the 28-day mean, the fill rate is only 0.798 and lost sales cost 1,728: the Poisson target is too low for these items. It also holds the least stock (8.92) and wastes the least (132). The empirical targets keep the most (R,S) stock (15.18) and give the best (R,S) fill rate, 0.904.
• (s,S) lowers the cost for every forecast: 2,251 → 1,597 with the 28-day mean, 1,714 → 1,623 with Empirical, and 1,836 → 1,470 with ETS. It orders up to a higher level, s plus half of the 4-day mean, so it holds more stock and loses fewer sales. It also places fewer order lines (for example 324 → 277 with ETS).
• (s,S) with ETS is the cheapest branch and has the highest fill rate: 1,470 and 0.944. Compared with (R,S) ETS, it holds more stock (19.36 against 14.32) and wastes more (361.50 against 269), but loses less than half as many sales (241 units against 522).

The ranking depends on the costs. At 1.00 per lost unit instead of 2.00, each lost-sales cost halves, and (s,S) with the 28-day mean becomes the cheapest (1,597 − 461 = 1,136, against 1,470 − 241 = 1,229 for (s,S) ETS). These results come from one six-week window of 30 items. Notebook 10 runs the (s,S) ETS setting as a daily production job.