guanlan-src/app_ontology/app_ontology_guanlan/sop/fatigue_validation.py

# -*- coding: utf-8 -*-
"""fatigue_validation.py — 疲劳暴露校核: per机组 **相对 TI(湍流)疲劳暴露排序** (非绝对寿命).

裁决对象 = fleet 内 **谁暴露在更强湍流 → 谁疲劳载荷循环剂量更高** 的相对排序 —
上游 (power-curve / 选址复盘) 想知道"哪几台位置的湍流暴露显著高于同场其余台"。
区别: 这**不是绝对疲劳寿命** (无气弹标定 / 无 S-N 材料曲线) → 只给**排序**不给寿命。

来源: fushan E映射 (outputs/fushan/sop/mapping_E_del.json, gitignored 真盘) 实证:
  fore-aft 塔振 DEL (雨流, 相对) 的 **within-turbine** 回归 R2 load-only 0.286 → +TI 0.40
  (dR2_TI 0.114, TI 系数 0.344, ws 系数 0.632); inter-turbine 台间 TI暴露 vs DEL
  spearman 0.42 (p=0.064, n=18)。→ 湍流强度 TI 在载荷(风速/推力/转速/功率)之外对
  fore-aft 疲劳有 **正边际**, 且台间 TI 暴露差与 DEL 剂量差同向。

★物理 trace (第一性, no black box):
  湍流强度 TI = σ(ws)/mean(ws) = 阵风波动幅度。塔架/叶根 fore-aft 疲劳 ∝ 载荷循环
  的幅值×次数 (Palmgren-Miner 线性累积 + S-N 曲线)。高 TI → 单位时间内推力波动的
  幅值与过零次数↑ → 雨流计数出的等效循环 DEL↑。故 "TI↑ → 相对疲劳暴露↑" 是**推力
  脉动驱动的载荷循环**这一物理量, 非统计巧合。within-turbine 弹性 (消逐台加速度计
  增益污染: 同台自比) 定标该驱动强度; inter-turbine 排序定位高暴露台。

★算法 (两层, 皆相对无黑箱):
  1. **暴露排序 (per台, inter-turbine)**: per台 gen 态有效窗 TI 均值 = ti_exposure;
     fleet 稳健 z = (ti_i − median)/(1.4826·MAD) (MAD=0 退化 std); **单侧** flag 高暴露
     (z ≥ z_gate) — 疲劳风险方向 = 高 TI。rank 按 ti_exposure 降序 (1=最高)。
  2. **弹性物理校验 (fleet, within-turbine, 可选)**: 有 DEL 代理列时, 逐台 z-within 变换
     log(DEL) ~ [z-within 载荷控制] (+z-within TI) 池化 OLS → R2_load_within / R2_plusTI /
     dR2_TI / TI 系数。**消增益污染** (同台自比, 镜像 mapping_E within_turbine 口径)。
     正 dR2 + 正 TI 系数 = TI→载荷弹性成立 → 排序物理前提坐实; 非正 → 响亮 caveat。

★诚实界 (task 明列 5 条 + 3原则审, 未调参掩盖):
  (a) **相对非绝对疲劳**: 无气弹标定 / 无 S-N 材料曲线 → 给 fleet 内**排序**, 不给寿命/
      DEL 绝对值/剩余年限。DEL 是相对雨流量 (逐台加速度计增益未标定)。
  (b) **TI 是暴露代理非直接损伤**: TI 高 ≠ 该台一定先坏 (材料/结构裕度/S-N 未入模);
      TI 是"输入端载荷环境"代理, 损伤还需材料侧信息才能定寿命。
  (c) **单场单月**: fushan Oct2023 1s 单场单月验证; **跨场泛化未测** (n=1 场, 非 systemic);
      迁到别场须重验 TI→DEL 弹性符号与量级 (§0.3 [暂行], SOP n<3 未定谳)。
  (d) **需 TI 通道 (ws_std/ws_mean 或直接 ti) 或塔振 DEL**: 两者皆缺 → 排序无输入 →
      **INSUFFICIENT 弃权** (非静默全零"干净")。仅 DEL 无 TI → 排序不可出 (排序轴是 TI)。
  (e) **fleet-relative → 批次盲**: 全场一致高 TI 暴露 (选址整体差) → 稳健 z 皆≈0 → **无一
      flag** = 本臂天然盲 (承 memory relative-criterion-batch-blindness + sensor 批次盲教训)。
      批次嫌疑场须先做场型识别 / 引入绝对 TI 基准 (IEC 类别 / 邻场), 不能只靠场内相对。
  (f) **DEL 代理异质 (RV-1 逮, 第五角数据面复发)**: DEL 通道跨台常两 population — fushan del_fa
      13 台 cv≈0.4 (物理) + 5 台 cv≈4 (A04-07/A11, mean 58-148 = 非物理坏通道)。原 pooled 弹性
      (dR2 0.114/TI系数 0.344) 是干净子集 (~0.21/0.46) 与死通道噪声的**稀释均值** → per台 DEL
      cv>del_cv_max(1.5) 已剔 (elasticity.del_excluded), 但**迁场必核 DEL 通道质量** (同 sensor 派生相/
      GC9 换制/pitch 高集距: 通道质量闸=生产必备)。premise_holds in-sample dR2 恒>0 → **约等 sign(TI系数)**,
      非强前提检验 (加 0.01 dR2 下界防 near-vacuous, 仍诚实标非 gold-standard)。

★可证伪判据 (falsifiable): 若高 TI 暴露台的 within-turbine DEL~TI 弹性 **≤ 0** (TI 不驱动
  DEL, 或反向) → "TI 暴露 = 疲劳剂量代理" 的物理前提在该场**假** → 排序不可作疲劳解读
  (退化为单纯湍流环境描述)。elasticity 块 ti_coef/dR2 符号即此判据的机器检查点。

★3 First Principles 自审:
  1. 物理: TI→推力脉动→载荷循环 DEL (Palmgren-Miner + S-N), 见上 trace ✓ (物理量非拟合噱头)。
  2. 统计: n 逐台声明 (min_windows) + fleet n (min_turb); 聚合层显式 (inter=排序/intra=弹性);
     falsifiable 判据显式; confounder = 载荷 (ws/rotor/power) 入 within 控制项排掉;
     单场单月 → **不作跨场 systemic claim** (诚实界 c)。
  3. 第一性: OLS + 稳健 z (median/MAD), 无黑箱; 不诉诸 memory label (fushan 数字须 data-first 复现)。

DP 验收 = tests/sop/test_fatigue_ti_dp.py (合成 deterministic + fushan 真盘 skipif sanity)。
纯核心无 I/O 可测; per-farm 加载在调用方 (skill/report 层)。
"""
from __future__ import annotations

import numpy as np
import pandas as pd

Z_GATE = 2.0          # 单侧高暴露 flag 门 (fleet 稳健 z)
MIN_WINDOWS = 200     # per台最少有效窗 (暴露均值稳)
MIN_TURB = 5          # fleet-相对 z 最少台数
MIN_POOL_ROWS = 1000  # within-turbine 弹性池化 OLS 最少行数
TI_LO = 0.01          # TI 有效下界 (去 0/负野值)
TI_HI = 0.6           # TI 有效上界 (去传感器野值)


def _robust_z(vals: np.ndarray) -> np.ndarray:
    """稳健 z = (x − median)/(1.4826·MAD); MAD=0 退化 std(ddof=1); 全常量 → 全 0。"""
    med = float(np.median(vals))
    mad = float(np.median(np.abs(vals - med))) * 1.4826
    sd = mad if mad > 0 else float(np.std(vals, ddof=1)) if len(vals) > 1 else 0.0
    if not sd or sd <= 0:
        return np.zeros_like(vals, dtype=float)
    return (vals - med) / sd


def _zwithin(df: pd.DataFrame, col: str, turbine_col: str) -> np.ndarray:
    """逐台 z-within 变换 (x − 台均)/台std → 消逐台增益/偏置污染 (同台自比)。
    单元素/常量台 std=NaN/0 → 结果 NaN → fillna(0) (无 within 方差, 不贡献斜率)。"""
    g = df.groupby(turbine_col)[col]
    z = (df[col] - g.transform("mean")) / g.transform("std")
    return z.replace([np.inf, -np.inf], np.nan).fillna(0.0).to_numpy(float)


def _r2(y: np.ndarray, yh: np.ndarray) -> float:
    ss_res = float(np.sum((y - yh) ** 2))
    ss_tot = float(np.sum((y - y.mean()) ** 2))
    return 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan")


def _within_elasticity(sub: pd.DataFrame, *, del_col: str, ti_col: str,
                       load_cols: list, turbine_col: str, min_pool_rows: int,
                       del_log: bool, del_cv_max: float = 1.5) -> dict | None:
    """within-turbine 池化 OLS: z-within(logDEL) ~ z-within(载荷) (+z-within TI)。
    镜像 mapping_E_fushan_analyze within_turbine 口径 (zwin + ΔR2_TI)。None = 行数不足。

    ★DEL 通道活性闸 (RV-1 逮, 第五角数据面复发): DEL 代理跨台常两 population —
    fushan del_fa 13台 cv≈0.4(物理) + 5台 cv≈4(A04-07/A11, mean 58-148=非物理坏通道) →
    pool 会稀释 (0.114 是干净13台强信号与死5台噪声的均值)。**per台 DEL cv>del_cv_max 剔除再池化**
    (承 sensor 派生相 / GC9 换制 / pitch 高集距: 通道质量闸=生产必备, 别静默 pool 垃圾), 报 del_excluded。"""
    need = [turbine_col, del_col, ti_col] + load_cols
    d = sub[need].replace([np.inf, -np.inf], np.nan).dropna()
    d = d[d[del_col] > 0]                                    # DEL>0 (log 前提)
    # per台 DEL cv 活性闸: 剔非物理坏通道 (cv=std/|mean| ≫ 物理), pool 前
    cv = d.groupby(turbine_col)[del_col].agg(lambda x: float(x.std() / abs(x.mean())) if x.mean() else np.inf)
    kept = cv[cv <= del_cv_max].index.tolist()
    excluded = sorted(str(t) for t in cv[cv > del_cv_max].index)
    d = d[d[turbine_col].astype(str).isin([str(t) for t in kept])]
    if len(d) < min_pool_rows:
        return None
    y_raw = np.log(d[del_col].to_numpy(float)) if del_log else d[del_col].to_numpy(float)
    dd = d.assign(_y=y_raw)
    yw = _zwithin(dd, "_y", turbine_col)
    ones = np.ones(len(dd))
    load_z = [_zwithin(dd, c, turbine_col) for c in load_cols]
    ti_z = _zwithin(dd, ti_col, turbine_col)
    X_load = np.column_stack([ones] + load_z) if load_z else ones.reshape(-1, 1)
    X_ti = np.column_stack([X_load, ti_z])
    b_load, *_ = np.linalg.lstsq(X_load, yw, rcond=None)
    b_ti, *_ = np.linalg.lstsq(X_ti, yw, rcond=None)
    r2_load = _r2(yw, X_load @ b_load)
    r2_ti = _r2(yw, X_ti @ b_ti)
    ti_coef = float(b_ti[-1])
    dr2 = r2_ti - r2_load
    return {
        "R2_load_within": round(r2_load, 4),
        "R2_plusTI_within": round(r2_ti, 4),
        "dR2_TI_within": round(dr2, 4),
        "TI_coef_within": round(ti_coef, 4),
        "n_windows": int(len(dd)),
        "n_turbines": int(dd[turbine_col].nunique()),
        "load_controls": list(load_cols),
        "del_excluded": excluded,                            # cv>del_cv_max 剔的死通道台 (RV-1)
        "del_cv_max": del_cv_max,
        # premise: in-sample 加回归元 dR2 恒>0 → 近似 sign(TI_coef) (RV-1 逮); 加 dR2 幅值下界防 vacuous
        "premise_holds": bool(ti_coef > 0 and dr2 > 0.01),
        "premise_note": "premise=ti_coef>0 ∧ dR2>0.01; in-sample dR2 恒>0 故约等 sign(TI_coef), 0.01 下界防 near-vacuous (RV-1)",
    }


def fatigue_ti_exposure_verdict(
    df: pd.DataFrame, *, turbine_col: str,
    ti_col: str | None = None, ws_std_col: str | None = None, ws_mean_col: str | None = None,
    del_col: str | None = None,
    ws_col: str | None = None, rotor_col: str | None = None, power_col: str | None = None,
    gen_frac_col: str | None = None, gen_frac_min: float = 0.8, gen_power_min: float = 50.0,
    z_gate: float = Z_GATE, min_windows: int = MIN_WINDOWS, min_turb: int = MIN_TURB,
    min_pool_rows: int = MIN_POOL_ROWS, ti_lo: float = TI_LO, ti_hi: float = TI_HI,
    del_log: bool = True,
) -> dict:
    """per-机组 **相对 TI 疲劳暴露排序** 裁决 (相对非绝对; 见模块 docstring 诚实界)。

    返回 {
      "verdict": "CANDIDATE_RANKING" | "INSUFFICIENT",
      "per_turbine": {tid: {"ti_exposure": TI均值, "z": fleet稳健z, "flag": 高暴露bool,
                            "rank": 1..N (1=最高TI暴露), "n": 有效窗, "note"?: 弃权/低置信}},
      "elasticity": {within-turbine TI→DEL 弹性块} | None,
      "honesty": [5 条诚实界字符串],
      "note": 总体说明,
    }
    TI 源: ti_col (直接) 优先; 否则 ws_std_col/ws_mean_col 派生 TI=σ/mean; 皆缺 → INSUFFICIENT。
    flag = 单侧高暴露 (z ≥ z_gate)。elasticity 需 del_col (+可选载荷控制 ws/rotor/power)。
    df 要求: 逐行 = 一个窗 (10min 或调用方定义); gen 掩码用 gen_frac_col≥gen_frac_min 优先,
      否则 power_col>gen_power_min, 皆无则不掩 (调用方保证运行态)。阈值/口径见模块 docstring。
    """
    HONESTY = [
        "(a) 相对非绝对疲劳: 无气弹标定/无S-N曲线 → 给fleet内排序, 不给寿命/DEL绝对值/剩余年限。",
        "(b) TI是暴露代理非直接损伤: 高TI≠该台一定先坏, 材料/结构裕度/S-N未入模。",
        "(c) 单场单月(fushan Oct2023)验证, 跨场泛化未测(n=1场); 迁场须重验TI→DEL弹性符号与量级。",
        "(d) 需TI通道(ws_std/ws_mean或ti)或塔振DEL; 两者皆缺→INSUFFICIENT弃权(非静默全零)。",
        "(e) fleet-relative→批次盲: 全场一致高TI暴露→稳健z皆≈0→无一flag; 批次嫌疑场须先场型识别/绝对TI基准。",
        "(f) DEL代理异质 (RV-1逮): DEL通道跨台常两population (fushan 13台cv≈0.4物理 + 5台cv≈4非物理坏通道); "
        "pooled弹性(0.114/0.344)是干净子集(~0.21/0.46)与死通道噪声的稀释均值 → per台cv>del_cv_max已剔(del_excluded), "
        "但迁场须核DEL通道质量; premise_holds in-sample约等sign(TI系数), 非强前提检验。",
    ]
    # NaN-safe: Arrow str backend 下 astype(str) 留 float NaN → sorted 混型崩; dropna + 显式 str()
    tids = sorted({str(x) for x in df[turbine_col].dropna().unique()})

    def _insufficient(reason: str) -> dict:
        return {
            "verdict": "INSUFFICIENT",
            "per_turbine": {t: {"ti_exposure": None, "z": None, "flag": False, "rank": None,
                                "n": 0, "note": reason} for t in tids},
            "elasticity": None, "honesty": HONESTY, "note": reason,
        }

    if turbine_col not in df.columns:
        return _insufficient("无 turbine_col; abstain (INSUFFICIENT)")

    # ---- TI 源解析 (诚实界 d): ti_col 优先, 否则 ws_std/ws_mean 派生; 皆缺 → 弃权 ----
    work = pd.DataFrame({turbine_col: df[turbine_col].astype(str).to_numpy()})
    if ti_col and ti_col in df.columns:
        work["_ti"] = pd.to_numeric(df[ti_col], errors="coerce").to_numpy(float)
        ti_src = f"直接 TI 列 '{ti_col}'"
    elif ws_std_col and ws_std_col in df.columns and ws_mean_col and ws_mean_col in df.columns:
        std = pd.to_numeric(df[ws_std_col], errors="coerce").to_numpy(float)
        mean = pd.to_numeric(df[ws_mean_col], errors="coerce").to_numpy(float)
        with np.errstate(divide="ignore", invalid="ignore"):
            work["_ti"] = np.where(mean > 0, std / mean, np.nan)
        ti_src = f"派生 TI = {ws_std_col}/{ws_mean_col}"
    else:
        return _insufficient(
            "无 TI 输入 (ti_col 或 ws_std_col+ws_mean_col 均缺) → 暴露排序无输入; abstain (INSUFFICIENT)")

    # ---- gen 掩码 (fraction 优先, 单位无关; 否则 power 阈; 皆无不掩) ----
    if gen_frac_col and gen_frac_col in df.columns:
        gen = pd.to_numeric(df[gen_frac_col], errors="coerce").to_numpy(float) >= gen_frac_min
    elif power_col and power_col in df.columns:
        gen = pd.to_numeric(df[power_col], errors="coerce").to_numpy(float) > gen_power_min
    else:
        gen = np.ones(len(df), dtype=bool)
    # 载荷控制列 (within 弹性用; 存在才带)
    for src, dst in ((ws_col, "_ws"), (rotor_col, "_rotor"), (power_col, "_power")):
        if src and src in df.columns:
            work[dst] = pd.to_numeric(df[src], errors="coerce").to_numpy(float)
    if del_col and del_col in df.columns:
        work["_del"] = pd.to_numeric(df[del_col], errors="coerce").to_numpy(float)

    valid = gen & np.isfinite(work["_ti"].to_numpy(float)) \
        & (work["_ti"].to_numpy(float) > ti_lo) & (work["_ti"].to_numpy(float) < ti_hi)
    sub = work[valid].copy()

    # ---- 层1: per台 TI 暴露 + fleet 稳健 z + flag ----
    per = {t: {"ti_exposure": None, "z": None, "flag": False, "rank": None, "n": 0} for t in tids}
    grp = sub.groupby(turbine_col)["_ti"]
    ti_mean = grp.mean()
    ti_n = grp.size()
    qualified = [t for t in tids if int(ti_n.get(t, 0)) >= min_windows]
    for t in tids:
        n_t = int(ti_n.get(t, 0))
        per[t]["n"] = n_t
        if t not in qualified:
            per[t]["note"] = f"有效窗 {n_t}<{min_windows} (疑停机/无本工况/TI野值滤尽); abstain"
        else:
            per[t]["ti_exposure"] = round(float(ti_mean[t]), 4)

    # ---- 层2: within-turbine TI→DEL 弹性 (可选, 有 DEL 才算) ----
    elasticity = None
    if "_del" in sub.columns:
        load_cols = [c for c in ("_ws", "_rotor", "_power") if c in sub.columns]
        elasticity = _within_elasticity(
            sub, del_col="_del", ti_col="_ti", load_cols=load_cols,
            turbine_col=turbine_col, min_pool_rows=min_pool_rows, del_log=del_log)
        if elasticity is None:
            elasticity = {"note": f"DEL 有效行 <{min_pool_rows} → within 弹性不可估 (仅出暴露排序)"}
    else:
        elasticity = {"note": "无 DEL 代理列 → 跳过 within-turbine 弹性物理校验 (仅出暴露排序)"}

    if len(qualified) < min_turb:
        note = (f"合格台数 {len(qualified)}<{min_turb} → fleet 稳健 z 不可估; abstain (INSUFFICIENT)。"
                f" 已报 per台 ti_exposure/n (信息), 但不可裁排序。")
        for t in qualified:
            per[t]["note"] = note
        return {"verdict": "INSUFFICIENT", "per_turbine": per, "elasticity": elasticity,
                "honesty": HONESTY, "note": note}

    # fleet 稳健 z (仅合格台参与) + 单侧 flag + rank
    q_arr = np.array([float(ti_mean[t]) for t in qualified], float)
    z_arr = _robust_z(q_arr)
    order = sorted(qualified, key=lambda t: float(ti_mean[t]), reverse=True)  # 1=最高
    rank_of = {t: i + 1 for i, t in enumerate(order)}
    for t, z in zip(qualified, z_arr):
        per[t]["z"] = round(float(z), 3)
        per[t]["flag"] = bool(z >= z_gate)
        per[t]["rank"] = rank_of[t]

    n_flag = sum(1 for t in qualified if per[t]["flag"])
    premise = elasticity.get("premise_holds") if isinstance(elasticity, dict) else None
    note = (f"相对 TI 疲劳暴露排序 (候选级, 非绝对寿命): {ti_src}; 合格 {len(qualified)}/{len(tids)} 台; "
            f"高暴露 flag {n_flag} 台 (z≥{z_gate})。")
    if premise is True:
        note += " within-turbine TI→DEL 弹性正 (dR2>0 ∧ TI系数>0) → 疲劳解读物理前提坐实。"
    elif premise is False:
        note += (" ⚠ within-turbine TI→DEL 弹性非正 → 可证伪判据触发: 本场 TI 不驱动 DEL, "
                 "排序仅为湍流环境描述, 不可作疲劳剂量解读 (诚实界 falsifiable)。")
    return {"verdict": "CANDIDATE_RANKING", "per_turbine": per, "elasticity": elasticity,
            "honesty": HONESTY, "note": note}