vpp-ai-platform/skills-py/vpp_skills/bid_milp.py

"""Day-ahead bid optimization MILP (docs/05 §2.2, docs/07 D-1 08:00).

Decision: for each of 96 intervals, an offer quantity q_t (MWh) and an offer
price. Model (v1):

    maximise   Σ_t (w_t − c − DEVIATION_PENALTY) · q_t
    s.t.       q_t ≤ k · cap_t · 0.25 · u_t         (sellable share of certified capacity)
               q_t ≥ min_block · u_t                 (market block size, u_t ∈ {0,1})
               E_min ≤ Σ_t q_t ≤ E_max               (position ledger cascade, P7)

with w_t = (1−λ)·P50_t + λ·P10_t the risk-weighted price (λ = risk_aversion)
and c the marginal cost of delivering flexibility. In a uniform-price market a
seller is paid the clearing price, so the *offer price* is the cost floor c
(price-taker when c = 0) and the forecast's job is allocation: which intervals
get the bounded daily energy. Risk aversion tilts allocation away from
intervals whose price band is wide (low P10). The revenue distribution
evaluates the bid against each forecast quantile path with the clearing rule
"cleared where offer ≤ price".

Solver: HiGHS via scipy.optimize.milp — open source, deterministic, adequate
for 96 binaries. Pyomo/commercial solvers (docs/08) are a swap behind the
same contract.
"""

from __future__ import annotations

import time
from decimal import Decimal

import numpy as np
import scipy
from scipy.optimize import Bounds, LinearConstraint, milp
from scipy.sparse import csr_matrix, hstack, identity

from vpp_contracts.bid_optimization_request import BidOptimizationRequest
from vpp_contracts.bid_optimization_result import BidOptimizationResult

from . import SKILL_VERSIONS
from .numeric import INTERVALS, SCALE_MONEY, SCALE_MWH, SCALE_PRICE, curve, dsum, quantize, to_array

SKILL = "bid-optimization-milp"

# OPEN-QUESTION A4: Hubei deviation-assessment rule (band + penalty price
# mechanism) is undecided. Until it is, the objective carries a zero penalty
# per offered MWh so the term exists and can be wired to config, but it does
# not shape the bid. See docs/open-questions.md → market.deviation.
DEVIATION_PENALTY_YUAN_PER_MWH = Decimal("0")

STATUS_BY_SCIPY = {0: "OPTIMAL", 1: "TIME_LIMIT", 2: "INFEASIBLE", 3: "INFEASIBLE", 4: "ERROR"}
TOL = 1e-6


def _revenue(offer: list[str], qty: list[str], path: list[str]) -> Decimal:
    total = Decimal(0)
    for o, q, p in zip(offer, qty, path, strict=True):
        if Decimal(o) <= Decimal(p):
            total += Decimal(p) * Decimal(q)
    return total


def _fit_to_bounds(qty: np.ndarray, cap_mwh: np.ndarray, min_block: float, bounds) -> list[str]:
    """Quantize to SCALE_MWH without leaving [E_min, E_max].

    Rounding 96 values independently can move the sum by up to 96·½ulp, enough
    for the ledger to reject a bid the solver found feasible. Floor each value
    (sum can only drop), then top up the deficit — if any — on intervals that
    still have capacity headroom, all in exact decimal arithmetic.
    """
    ulp = Decimal(1).scaleb(-SCALE_MWH)
    vals = [(Decimal(repr(float(v))).quantize(ulp, rounding="ROUND_FLOOR")) for v in qty.tolist()]
    vals = [abs(v) if v == 0 else v for v in vals]
    caps = [Decimal(repr(float(c))).quantize(ulp, rounding="ROUND_FLOOR") for c in cap_mwh.tolist()]
    e_min = Decimal(bounds.daily_energy_min_mwh)
    total = sum(vals, Decimal(0))
    if total < e_min and total > 0:
        deficit = e_min - total
        for i in sorted(range(len(vals)), key=lambda j: caps[j] - vals[j], reverse=True):
            if vals[i] == 0 or deficit <= 0:
                continue  # never open a new interval below the block size
            room = caps[i] - vals[i]
            add = min(room, deficit).quantize(ulp, rounding="ROUND_CEILING")
            add = min(add, room)
            vals[i] += add
            deficit -= add
        # If the deficit could not be placed the solver bound was tight to
        # machine precision; the ledger check will judge the residual.
    return [format(v, "f") for v in vals]


def optimize_bid(req: BidOptimizationRequest) -> BidOptimizationResult:
    if req.price_forecast.kind.value != "PRICE":
        raise ValueError("price_forecast.kind must be PRICE")
    if req.price_forecast.market_date != req.market_date or req.adjustable_capacity_mw.date != req.market_date:
        raise ValueError("price forecast and capacity curve must be for market_date")

    lam = float(Decimal(req.risk.risk_aversion))
    k = float(Decimal(req.risk.commitment_buffer_k))
    min_block = float(Decimal(req.risk.min_block_mwh))
    e_min = float(Decimal(req.position_bounds.daily_energy_min_mwh))
    e_max = float(Decimal(req.position_bounds.daily_energy_max_mwh))
    if not 0.0 <= lam <= 1.0:
        raise ValueError("risk_aversion must be in [0, 1]")
    if not 0.0 < k <= 1.0:
        raise ValueError("commitment_buffer_k must be in (0, 1]")
    cost = float(Decimal(req.risk.marginal_cost_yuan_per_mwh))
    if cost < 0:
        raise ValueError("marginal_cost_yuan_per_mwh must be non-negative")
    if min_block < 0 or e_min < 0 or e_max < e_min:
        raise ValueError("invalid block size or position bounds")

    q = req.price_forecast.quantiles
    p10 = to_array(v.root for v in q.p10.values)
    p50 = to_array(v.root for v in q.p50.values)
    cap_mwh = k * to_array(v.root for v in req.adjustable_capacity_mw.values) * 0.25
    if np.any(cap_mwh < 0):
        raise ValueError("adjustable capacity must be non-negative")
    w = (1.0 - lam) * p50 + lam * p10
    penalty = float(DEVIATION_PENALTY_YUAN_PER_MWH)

    n = INTERVALS
    # Variables: x = [q_0..q_95, u_0..u_95]; minimise -(w - cost - penalty)·q
    c = np.concatenate([-(w - cost - penalty), np.zeros(n)])
    integrality = np.concatenate([np.zeros(n), np.ones(n)])
    bounds = Bounds(np.zeros(2 * n), np.concatenate([cap_mwh, np.ones(n)]))

    eye = identity(n, format="csr")
    a_cap = hstack([eye, -csr_matrix(np.diag(cap_mwh))])  # q - cap·u ≤ 0
    a_blk = hstack([eye, -min_block * eye])  # q - min_block·u ≥ 0
    a_sum = csr_matrix(np.concatenate([np.ones(n), np.zeros(n)])[None, :])
    constraints = [
        LinearConstraint(a_cap, -np.inf, 0.0),
        LinearConstraint(a_blk, 0.0, np.inf),
        LinearConstraint(a_sum, e_min, e_max),
    ]

    t0 = time.perf_counter()
    res = milp(c, constraints=constraints, integrality=integrality, bounds=bounds)
    wall_ms = int(round((time.perf_counter() - t0) * 1000))
    status = STATUS_BY_SCIPY.get(int(res.status), "ERROR")

    binding: list[str] = []
    if res.x is None:
        qty = np.zeros(n)
        if e_min > cap_mwh.sum() + TOL:
            binding.append("daily_energy_min exceeds sellable energy")
        objective = None
    else:
        qty = np.clip(res.x[:n], 0.0, None)
        qty[qty < TOL] = 0.0
        total = qty.sum()
        if abs(total - e_max) < 1e-4:
            binding.append("daily_energy_max")
        if abs(total - e_min) < 1e-4:
            binding.append("daily_energy_min")
        if np.any((cap_mwh > 0) & (np.abs(qty - cap_mwh) < 1e-6)):
            binding.append("interval_capacity")
        objective = quantize(-float(res.fun), SCALE_MONEY)

    quantities = curve(_fit_to_bounds(qty, cap_mwh, min_block, req.position_bounds), req.market_date, SCALE_MWH)
    offers = curve(np.full(n, cost), req.market_date, SCALE_PRICE)
    qv, ov = quantities["values"], offers["values"]
    p10s = [v.root for v in q.p10.values]
    p50s = [v.root for v in q.p50.values]
    p90s = [v.root for v in q.p90.values]

    result = {
        "market_date": req.market_date,
        "prices_yuan_per_mwh": offers,
        "quantities_mwh": quantities,
        "daily_energy_mwh": quantize(dsum(qv), SCALE_MWH),
        "expected_revenue_yuan": quantize(_revenue(ov, qv, p50s), SCALE_MONEY),
        "revenue_distribution_yuan": {
            "p10": quantize(_revenue(ov, qv, p10s), SCALE_MONEY),
            "p50": quantize(_revenue(ov, qv, p50s), SCALE_MONEY),
            "p90": quantize(_revenue(ov, qv, p90s), SCALE_MONEY),
        },
        "position_bounds": req.position_bounds.model_dump(),
        "solver": {
            "name": "highs",
            "version": f"scipy-{scipy.__version__}",
            "status": status,
            "objective_value": objective,
            "wall_time_ms": wall_ms,
        },
        "binding_constraints": binding,
        "skill_version": SKILL_VERSIONS[SKILL],
    }
    return BidOptimizationResult.model_validate(result)