"""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)