aps-agent/server/aps_domain/robustness.py

326 lines
13 KiB
Python
Raw Permalink Normal View History

from __future__ import annotations
import copy
import random
from typing import Any, Callable
from server.aps_domain.constraints import hard_blocking_conflicts
from server.contracts import ScheduleResult
from server.engines import get_engine
from server.engines.base import EngineParams
from server.timeutil import add_minutes, fmt_date, today0
World = dict[str, Any]
DEFAULT_MONTE_CARLO_SEED = 20260731
DEFAULT_MONTE_CARLO_TRIALS = 24
def _counter() -> Callable[[str], int]:
values: dict[str, int] = {}
def next_id(kind: str) -> int:
values[kind] = values.get(kind, 3000000) + 1
return values[kind]
return next_id
def _bounded_normal(rng: random.Random, mean: float, sigma: float, low: float, high: float) -> float:
return max(low, min(high, rng.gauss(mean, sigma)))
# ---------------- 可配置分布/相关性/置信区间(矩阵 88 行剩余项) ----------------
# distributions 结构(缺省与旧行为完全一致):
# {"quantity": {"dist": "normal"|"uniform"|"triangular", "sigma"/"low"/"high"/"mode", "mean": 1.0},
# "efficiency": {...},
# "shockProbability": 0.15}
_DEFAULT_DISTRIBUTIONS: dict[str, Any] = {
"quantity": {"dist": "normal", "mean": 1.0, "sigma": 0.10, "low": 0.80, "high": 1.20},
"efficiency": {"dist": "normal", "mean": 1.0, "sigma": 0.08, "low": 0.75, "high": 1.25},
"shockProbability": 0.15,
}
def _merged_distributions(distributions: dict[str, Any] | None) -> dict[str, Any]:
"""浅合并用户配置到默认分布表(用户只覆盖要改的维度)。"""
base = {
k: dict(v) if isinstance(v, dict) else v
for k, v in _DEFAULT_DISTRIBUTIONS.items()
}
for key, value in (distributions or {}).items():
if isinstance(value, dict) and isinstance(base.get(key), dict):
base[key].update(value)
else:
base[key] = value
return base
def _clamp(value: float, low: float, high: float) -> float:
return max(low, min(high, value))
def _sample_scale(rng: random.Random, cfg: dict[str, Any], z: float | None = None) -> float:
"""按配置分布采样一个缩放因子(z 供相关性复用同一正态流;None 时内部抽样)。"""
dist = str(cfg.get("dist") or "normal")
mean = float(cfg.get("mean") or 1.0)
low = float(cfg.get("low") or 0.0)
high = float(cfg.get("high") or 2.0)
if z is None:
z = rng.gauss(0.0, 1.0)
if dist == "uniform":
lo = float(cfg.get("low") or 0.9)
hi = float(cfg.get("high") or 1.1)
return round(rng.uniform(lo, hi), 6)
if dist == "triangular":
lo = float(cfg.get("low") or 0.85)
hi = float(cfg.get("high") or 1.15)
mode = float(cfg.get("mode") or mean)
return round(rng.triangular(lo, hi, mode), 6)
sigma = float(cfg.get("sigma") or 0.10)
return round(_clamp(mean + sigma * z, low, high), 6)
def _sample_correlated_scales(
rng: random.Random, distributions: dict[str, Any], correlation: float
) -> tuple[float, float]:
"""抽取 quantity/efficiency 缩放因子(支持相关系数 rho ∈ [-1,1])。
共用标准正态流:z1, z2 独立;e = rho*z1 + sqrt(1-rho^2)*z2。
rho=0 时退化为两个独立正态,且随机流消耗顺序与旧实现一致(种子兼容)。
"""
z1 = rng.gauss(0.0, 1.0)
z2 = rng.gauss(0.0, 1.0)
rho = max(-1.0, min(1.0, float(correlation or 0.0)))
e_z = rho * z1 + (1.0 - rho * rho) ** 0.5 * z2
q_cfg = dict(distributions.get("quantity") or _DEFAULT_DISTRIBUTIONS["quantity"])
e_cfg = dict(distributions.get("efficiency") or _DEFAULT_DISTRIBUTIONS["efficiency"])
q_dist = str(q_cfg.get("dist") or "normal")
e_dist = str(e_cfg.get("dist") or "normal")
if q_dist == "uniform":
q = rng.uniform(float(q_cfg.get("low") or 0.9), float(q_cfg.get("high") or 1.1))
elif q_dist == "triangular":
q = rng.triangular(float(q_cfg.get("low") or 0.85), float(q_cfg.get("high") or 1.15),
float(q_cfg.get("mode") or 1.0))
else:
q = _sample_scale(rng, q_cfg, z=z1)
if e_dist == "uniform":
e = rng.uniform(float(e_cfg.get("low") or 0.9), float(e_cfg.get("high") or 1.1))
elif e_dist == "triangular":
e = rng.triangular(float(e_cfg.get("low") or 0.85), float(e_cfg.get("high") or 1.15),
float(e_cfg.get("mode") or 1.0))
else:
e = _sample_scale(rng, e_cfg, z=e_z)
return q, e
def _perturb_world(
world: World,
rng: random.Random,
*,
reset_schedule_products: bool,
distributions: dict[str, Any] | None = None,
correlation: float = 0.0,
) -> tuple[World, dict[str, Any]]:
sandbox = copy.deepcopy(world)
if reset_schedule_products:
sandbox["scheduleVersions"] = []
sandbox["productionOrders"] = []
sandbox["workOrders"] = []
sandbox["conflicts"] = []
dist = _merged_distributions(distributions)
quantity_scale, efficiency_scale = _sample_correlated_scales(rng, dist, correlation)
for order in sandbox.get("salesOrders") or []:
for item in order.get("items") or []:
quantity = float(item.get("quantity") or 0.0)
if quantity > 0:
item["quantity"] = max(1, int(round(quantity * quantity_scale)))
for line in sandbox.get("lines") or []:
base = float(line.get("efficiencyFactor") or 1.0)
line["efficiencyFactor"] = round(base * efficiency_scale, 6)
capacity_shock_line = None
lines = sandbox.get("lines") or []
shock_prob = float(dist.get("shockProbability", 0.15) or 0.0)
if lines and shock_prob > 0 and rng.random() < shock_prob:
shocked = lines[rng.randrange(len(lines))]
shocked["efficiencyFactor"] = round(float(shocked.get("efficiencyFactor") or 1.0) * 0.55, 6)
capacity_shock_line = shocked.get("id")
return sandbox, {
"quantityScale": round(quantity_scale, 6),
"efficiencyScale": round(efficiency_scale, 6),
"capacityShockLineId": capacity_shock_line,
}
def _percentile(values: list[float], percentile: float) -> float:
if not values:
return 0.0
ordered = sorted(values)
position = (len(ordered) - 1) * percentile
lower = int(position)
upper = min(lower + 1, len(ordered) - 1)
fraction = position - lower
return round(ordered[lower] * (1.0 - fraction) + ordered[upper] * fraction, 4)
def run_monte_carlo(
world: World,
*,
strategy: str = "COMPREHENSIVE",
engine_type: str = "RULE",
trials: int = DEFAULT_MONTE_CARLO_TRIALS,
seed: int = DEFAULT_MONTE_CARLO_SEED,
baseline_kpi: dict[str, Any] | None = None,
planning_horizon_days: int = 14,
start_date: str | None = None,
delivery_buffer_ratio: float | None = None,
freeze_window_hours: float | None = None,
constraints: dict[str, bool] | None = None,
reset_schedule_products: bool = False,
distributions: dict[str, Any] | None = None,
correlation: float = 0.0,
confidence_level: float = 0.95,
) -> dict[str, Any]:
"""Re-solve fixed-seed demand/capacity perturbations in isolated sandboxes.
经统一 Explore 通道执行(矩阵 55 行):fn 只拿深拷贝沙盒,主干永不外泄写引用;
每 trial 的 _perturb_world 再在沙盒上做独立扰动,互不污染。
矩阵 88 行剩余项:distributions 可配置分布(normal/uniform/triangular)、
correlation 需求-效率扰动相关度(-1..1)、confidence_level 置信区间置信度。
"""
from server.aps_domain.explore_boundary import run_explore
return run_explore(world, lambda sandbox: _monte_carlo_impl(
sandbox,
strategy=strategy, engine_type=engine_type, trials=trials, seed=seed,
baseline_kpi=baseline_kpi, planning_horizon_days=planning_horizon_days,
start_date=start_date, delivery_buffer_ratio=delivery_buffer_ratio,
freeze_window_hours=freeze_window_hours, constraints=constraints,
reset_schedule_products=reset_schedule_products,
distributions=distributions, correlation=correlation,
confidence_level=confidence_level,
))
from server.aps_domain.explore_boundary import run_explore
return run_explore(world, lambda sandbox: _monte_carlo_impl(
sandbox,
strategy=strategy, engine_type=engine_type, trials=trials, seed=seed,
baseline_kpi=baseline_kpi, planning_horizon_days=planning_horizon_days,
start_date=start_date, delivery_buffer_ratio=delivery_buffer_ratio,
freeze_window_hours=freeze_window_hours, constraints=constraints,
reset_schedule_products=reset_schedule_products,
))
def _monte_carlo_impl(
world: World,
*,
strategy: str,
engine_type: str,
trials: int,
seed: int,
baseline_kpi: dict[str, Any] | None,
planning_horizon_days: int,
start_date: str | None,
delivery_buffer_ratio: float | None,
freeze_window_hours: float | None,
constraints: dict[str, bool] | None,
reset_schedule_products: bool,
distributions: dict[str, Any] | None,
correlation: float,
confidence_level: float,
) -> dict[str, Any]:
"""Monte Carlo 核心实现(在统一通道沙盒内运行)。"""
sample_count = max(1, min(500, int(trials)))
baseline = baseline_kpi or {}
baseline_tardiness = max(0.0, float(baseline.get("totalTardiness") or baseline.get("tardiness") or 0.0))
baseline_conflicts = max(0, int(baseline.get("conflictCount") or baseline.get("conflicts") or 0))
criteria = {
"maxTotalTardiness": round(baseline_tardiness + max(1.0, baseline_tardiness * 0.10), 4),
"maxConflictCount": baseline_conflicts,
"requireNoHardViolations": True,
}
rng = random.Random(int(seed))
start = start_date or fmt_date(add_minutes(today0(), 24 * 60))
outcomes: list[dict[str, Any]] = []
for trial in range(sample_count):
sandbox, perturbation = _perturb_world(
world,
rng,
reset_schedule_products=reset_schedule_products,
distributions=distributions,
correlation=correlation,
)
params = EngineParams(
orderIds=[],
engineType=engine_type,
strategyTemplate=strategy,
planningHorizonDays=planning_horizon_days,
startDate=start,
deliveryBufferRatio=delivery_buffer_ratio,
freezeWindowHours=freeze_window_hours,
name=f"mc-{seed}-{trial + 1}",
constraints=constraints or EngineParams().constraints,
)
result: ScheduleResult = get_engine(engine_type).solve(sandbox, params, _counter())
hard_rows = hard_blocking_conflicts(sandbox, result.versionId, track="fixed")
hard_ids = {row.get("id") for row in hard_rows}
critical_unmapped = [
row for row in sandbox.get("conflicts") or []
if row.get("versionId") == result.versionId
and row.get("severity") == "CRITICAL"
and row.get("id") not in hard_ids
]
hard_count = len(hard_rows) + len(critical_unmapped)
accepted = (
hard_count == 0
and float(result.totalTardiness) <= criteria["maxTotalTardiness"] + 1e-9
and int(result.conflictCount) <= criteria["maxConflictCount"]
)
outcomes.append({
"trial": trial + 1,
**perturbation,
"totalTardiness": round(float(result.totalTardiness), 4),
"conflictCount": int(result.conflictCount),
"hardViolationCount": hard_count,
"accepted": accepted,
})
accepted_count = sum(1 for outcome in outcomes if outcome["accepted"])
tardiness = [float(outcome["totalTardiness"]) for outcome in outcomes]
return {
"method": "fixed-seed-monte-carlo",
"seed": int(seed),
"trials": sample_count,
"acceptanceCriteria": criteria,
"acceptedTrials": accepted_count,
"robustness": round(accepted_count / sample_count, 6),
"distribution": {
"tardinessP50": _percentile(tardiness, 0.50),
"tardinessP90": _percentile(tardiness, 0.90),
"tardinessP95": _percentile(tardiness, 0.95),
"tardinessMax": round(max(tardiness, default=0.0), 4),
"meanConflicts": round(
sum(outcome["conflictCount"] for outcome in outcomes) / sample_count, 4
),
# 矩阵 88 行剩余项:置信区间(经验分位法)+ 均值/标准差
"tardinessMean": round(sum(tardiness) / sample_count, 4) if tardiness else 0.0,
"tardinessStd": round(
(sum((v - sum(tardiness) / sample_count) ** 2 for v in tardiness) / sample_count) ** 0.5, 4
) if tardiness else 0.0,
"confidenceLevel": float(confidence_level),
"tardinessCI": [
_percentile(tardiness, (1.0 - float(confidence_level)) / 2.0),
_percentile(tardiness, 1.0 - (1.0 - float(confidence_level)) / 2.0),
],
},
"config": {
"distributions": _merged_distributions(distributions),
"correlation": float(correlation),
"confidenceLevel": float(confidence_level),
},
"outcomes": outcomes,
}