Warehouse planning when the market turns

import numpy as np

SITES, REGIONS = 5, 8
TRANSPORT = 0.02          # EUR per pallet and km


def network(number):
    """Network `number`: distances between sites and regions, each region's
    forecast, and each site's capacity and fixed cost."""
    rng = np.random.default_rng(number)
    site_xy = rng.uniform([0, 0], [600, 800], (SITES, 2))        # km
    region_xy = rng.uniform([0, 0], [600, 800], (REGIONS, 2))
    return {
        "distance": np.linalg.norm(site_xy[:, None, :] - region_xy[None, :, :], axis=2),
        "forecast": rng.uniform(60, 140, REGIONS),      # thousand pallets a year
        "capacity": rng.uniform(200, 350, SITES),       # thousand pallets a year
        "fixed_cost": rng.uniform(800, 1200, SITES),    # thousand EUR a year
    }


networks = [network(number) for number in range(20)]

shown = networks[11]
print("site  capacity  fixed cost")
for s in range(SITES):
    print(f"{s:4}  {shown['capacity'][s]:8.0f}  {shown['fixed_cost'][s]:10.0f}")
print("\nregion  forecast")
for r in range(REGIONS):
    print(f"{r:6}  {shown['forecast'][r]:8.0f}")

calm = {"error": 0.15, "correlation": 0.0, "overflow": 40.0}
turned = {"error": 0.30, "correlation": 0.5, "overflow": 80.0}


def years(net, market, n, seed):
    """`n` simulated years of demand: one row per year, one column per region."""
    rng = np.random.default_rng(seed)
    own = rng.standard_normal((n, REGIONS))     # each region's own swing
    shared = rng.standard_normal((n, 1))        # the swing all regions share
    rho = market["correlation"]
    swing = np.sqrt(1 - rho) * own + np.sqrt(rho) * shared
    return np.maximum(net["forecast"] + market["error"] * net["forecast"] * swing, 0.0)


def yearly_cost(net, market, plan, demand):
    """The cost of `plan` in each simulated year: fixed costs, transport, overflow."""
    open_sites, serves = plan
    transport = (TRANSPORT * net["distance"][serves, np.arange(REGIONS)] * demand).sum(axis=1)
    load = np.zeros((len(demand), SITES))
    for r in range(REGIONS):
        load[:, serves[r]] += demand[:, r]
    overflow = np.maximum(load - net["capacity"] * open_sites, 0.0).sum(axis=1)
    return net["fixed_cost"] @ open_sites + transport + market["overflow"] * overflow


def real_cost(number, market, plan):
    """The plan's average yearly cost over 100,000 simulated years."""
    net = networks[number]
    return yearly_cost(net, market, plan, years(net, market, 100_000, seed=number + 2000)).mean()

import pyomo.environ as pyo


def plan_for_forecast(number, market, margin=1.0):
    """The network that is cheapest if every region's demand equals its forecast,
    proven optimal by HiGHS. With a margin below 1, each warehouse counts as
    holding only that share of its capacity. Returns the plan and the cost the
    solver promises."""
    net = networks[number]
    m = pyo.ConcreteModel()
    m.open = pyo.Var(range(SITES), domain=pyo.Binary)                    # run a warehouse at site s?
    m.serve = pyo.Var(range(SITES), range(REGIONS), domain=pyo.Binary)   # does site s serve region r?
    m.over = pyo.Var(range(SITES), domain=pyo.NonNegativeReals)          # overflow at site s
    m.once = pyo.Constraint(range(REGIONS), rule=lambda m, r: sum(m.serve[s, r] for s in range(SITES)) == 1)
    m.only_open = pyo.Constraint(range(SITES), range(REGIONS), rule=lambda m, s, r: m.serve[s, r] <= m.open[s])
    m.overflow = pyo.Constraint(range(SITES), rule=lambda m, s: m.over[s] >= sum(
        net["forecast"][r] * m.serve[s, r] for r in range(REGIONS)) - margin * net["capacity"][s] * m.open[s])
    m.cost = pyo.Objective(sense=pyo.minimize, expr=(
        sum(net["fixed_cost"][s] * m.open[s] for s in range(SITES))
        + sum(TRANSPORT * net["distance"][s, r] * net["forecast"][r] * m.serve[s, r]
              for s in range(SITES) for r in range(REGIONS))
        + market["overflow"] * sum(m.over[s] for s in range(SITES))))
    solver = pyo.SolverFactory("appsi_highs")
    solver.config.mip_gap = 0.0
    assert pyo.check_optimal_termination(solver.solve(m))
    open_sites = np.array([round(m.open[s].value) for s in range(SITES)])
    serves = np.array([max(range(SITES), key=lambda s: m.serve[s, r].value) for r in range(REGIONS)])
    return (open_sites, serves), pyo.value(m.cost)


def describe(plan):
    open_sites, serves = plan
    warehouses = ", ".join(str(s) for s in range(SITES) if open_sites[s])
    return f"warehouses at sites {warehouses}; regions 0-7 served from sites {' '.join(map(str, serves))}"


plan, promised = plan_for_forecast(11, calm)
real = real_cost(11, calm, plan)
print(describe(plan))
print(f"cost the solver promises   {promised:6,.0f}")
print(f"real cost                  {real:6,.0f}   ({real / promised - 1:+.1%})")

def median_and_worst(values):
    return f"median {np.median(values):+.1%}, worst {np.max(values):+.1%}"


above_promise = []
for number in range(20):
    plan, promised = plan_for_forecast(number, calm)
    above_promise.append(real_cost(number, calm, plan) / promised - 1)
print("real cost above the promise:", median_and_worst(above_promise))

MARGINS = [1.0, 0.95, 0.90, 0.85, 0.80]

for margin in MARGINS:
    plan, _ = plan_for_forecast(11, calm, margin)
    print(f"plan to {margin:4.0%} of capacity: real cost {real_cost(11, calm, plan):6,.0f}   {describe(plan)}")

cheapest = {margin: 0 for margin in MARGINS}
for number in range(20):
    costs = {margin: real_cost(number, calm, plan_for_forecast(number, calm, margin)[0]) for margin in MARGINS}
    for margin, cost in costs.items():
        if cost <= min(costs.values()) * (1 + 1e-9):      # two margins that give the same plan both count
            cheapest[margin] += 1
for margin, count in cheapest.items():
    print(f"plan to {margin:4.0%} of capacity is cheapest in {count:2} of 20 networks")

from quicopt import Client
from quicopt.stochastic import Normal, set_distribution, set_scenarios, expectation, maximum
from quicopt.structured import partition, lookup_table

client = Client()


def plan_with_quicopt(number, market):
    """The network with the lowest average cost over 1,000 years that Quicopt draws
    from the forecast and its error. Returns the plan and Quicopt's estimate of
    its cost."""
    net = networks[number]
    rho = market["correlation"]
    spread = market["error"] * net["forecast"]
    m = pyo.ConcreteModel()
    set_scenarios(m, 1000, seed=7)
    split = partition(m, "network", REGIONS, SITES)      # region r is served by site part_of(r)
    m.demand = pyo.Var(range(REGIONS))                   # each region's own swing around its forecast
    for r in range(REGIONS):
        set_distribution(m, m.demand[r], Normal(net["forecast"][r], np.sqrt(1 - rho) * spread[r]),
                         name=f"demand_{r}")
    demand = [m.demand[r] for r in range(REGIONS)]
    if rho > 0:
        m.market = pyo.Var()                             # the swing all regions share
        set_distribution(m, m.market, Normal(0.0, 1.0), name="market")
        demand = [demand[r] + float(np.sqrt(rho) * spread[r]) * m.market for r in range(REGIONS)]
    freight = lookup_table(m, "freight", (TRANSPORT * net["distance"]).tolist())   # per pallet, [site][region]
    site = [split.part_of(r) for r in range(REGIONS)]
    load = [sum(pyo.Expr_if(IF=site[r] == s, THEN=demand[r], ELSE=0.0) for r in range(REGIONS))
            for s in range(SITES)]
    cost = (sum(pyo.Expr_if(IF=split.size_of(s) >= 1, THEN=net["fixed_cost"][s], ELSE=0.0) for s in range(SITES))
            + sum(freight[site[r], r] * demand[r] for r in range(REGIONS))
            + sum(market["overflow"] * maximum(load[s] - net["capacity"][s], 0.0) for s in range(SITES)))
    m.cost = pyo.Objective(expr=expectation(cost), sense=pyo.minimize)
    result = client.solve(m)
    assert result.feasible
    found = result.structures["network"]
    return (np.array([int(n >= 1) for n in found.size_of]), np.array(found.part_of)), result.objective


plan, estimate = plan_with_quicopt(11, calm)
real = real_cost(11, calm, plan)
print(describe(plan))
print(f"Quicopt's estimate   {estimate:6,.0f}")
print(f"real cost            {real:6,.0f}   ({real / estimate - 1:+.1%})")

def compare(market):
    """What each margin costs on top of Quicopt's plan, over all 20 networks, and
    how far Quicopt's real cost lies from its estimate."""
    quicopt_real, real_off = [], []
    for number in range(20):
        plan, estimate = plan_with_quicopt(number, market)
        quicopt_real.append(real_cost(number, market, plan))
        real_off.append(quicopt_real[-1] / estimate - 1)
    extra = {}
    for margin in MARGINS:
        extra[margin] = [real_cost(number, market, plan_for_forecast(number, market, margin)[0]) / quicopt_real[number] - 1
                         for number in range(20)]
        label = "plan for the forecast" if margin == 1 else f"plan to {margin:.0%} of capacity"
        print(f"{label:24}  {median_and_worst(extra[margin])}")
    print(f"\nQuicopt's real cost against its estimate: from {min(real_off):+.1%} to {max(real_off):+.1%}")
    return extra


calm_extra = compare(calm)

plan, promised = plan_for_forecast(11, turned)
print(f"plan for the forecast:     promised {promised:6,.0f}, real cost {real_cost(11, turned, plan):6,.0f}   {describe(plan)}")
for margin in MARGINS[1:]:
    plan, _ = plan_for_forecast(11, turned, margin)
    print(f"plan to {margin:.0%} of capacity:  real cost {real_cost(11, turned, plan):6,.0f}")
plan, estimate = plan_with_quicopt(11, turned)
print(f"Quicopt:                   estimate {estimate:6,.0f}, real cost {real_cost(11, turned, plan):6,.0f}   {describe(plan)}")

turned_extra = compare(turned)

above_promise = []
for number in range(20):
    plan, promised = plan_for_forecast(number, turned)
    above_promise.append(real_cost(number, turned, plan) / promised - 1)
print("plan for the forecast, real cost above the promise:", median_and_worst(above_promise))
worst = int(np.argmax(turned_extra[1.0]))
print(f"worst network for the plan for the forecast: {worst} ({turned_extra[1.0][worst]:+.1%})")