guanlan-src/app_algorithmModel/app_algorithmModel_guanlan/perf/trend.py

# -*- coding: utf-8 -*-
"""A翼 月度滚动趋势 + 欠发归因 router (windscada M10; SOP 1.4 + 1.3c).
纪律:
  - 2026-01 = 限电制度断点 → 禁跨断点单斜率 (trend-across-gap: 断点期 level shift 更可能是工况非退化);
    分段: 2025段斜率 / 2026段斜率 / 断点Δ(2026均值−2025H2均值) 三件分开报;
  - dev = 同月同风档配对 (只正常发电态行, 已剥限电);
  - 渐进退化判据: 段内斜率持续同号 + 端段对比 (斜率只作形态描述, SOP slope-dilution 纪律);
  - 归因 router (1.3c): 候选台逐轴排除 → 桨角(M7已排)/转速跟踪(λ代理)/刹车拖曳(温度轴)/形态(阶跃vs渐进vs恒低), 余判 INSUFFICIENT."""
import numpy as np, pandas as pd, pathlib
from src.windscada.config import farm

BREAK = '2026-01'
# 断点依据 (数据内证, 2026-08-26 外审后补): 逐月限电时长占比 2025=间歇(3.2%~66.6%波动, 夏季高)
# vs 2026=常态化(逐月 55%~73% 持续) → 2026 段内部环境同质, 分段有效; 表述用"限电常态化"非"制度切换"(制度归因需外部文件)

def monthly_dev(cfg=None):
    cfg = cfg or farm()
    df = pd.read_parquet(pathlib.Path(cfg['store']) / 'pc_monthly_bins.parquet')
    base = df.groupby(['month', 'wb']).apply(lambda g: np.average(g.p, weights=g.n), include_groups=False).rename('fleet_p')
    df = df.join(base, on=['month', 'wb'])
    def dev(g):
        if g.n.sum() < 300: return np.nan     # 月样本门: <50h 正常发电行 → 不出数 (重限电月饥饿防噪)
        w = g.fleet_p * g.n
        return float(((g.p - g.fleet_p) * g.n).sum() / max(w.sum(), 1e-9))
    out = df.groupby(['turbine', 'month']).apply(dev, include_groups=False).rename('dev').reset_index()
    rpm = df.groupby(['turbine', 'month'])[['rpm_ws', 'ws_mid']].first().reset_index()
    return out.merge(rpm, on=['turbine', 'month'], how='left')

def _ts_slope(y):
    """Theil-Sen (%/月); 仅形态描述."""
    y = [v for v in y if v == v]
    n = len(y)
    if n < 4: return np.nan
    sl = [(y[j] - y[i]) / (j - i) for i in range(n) for j in range(i + 1, n)]
    return float(np.median(sl)) * 100

def registry(cfg=None):
    cfg = cfg or farm()
    md = monthly_dev(cfg)
    piv = md.pivot_table(index='month', columns='turbine', values='dev')
    m25 = [m for m in piv.index if m < BREAK]
    m25h2 = [m for m in piv.index if '2025-07' <= m < BREAK]
    m26 = [m for m in piv.index if m >= BREAK]
    rows = []
    for t in piv.columns:
        s = piv[t]
        sl25, sl26 = _ts_slope(s.loc[m25]), _ts_slope(s.loc[m26])
        d_break = float(s.loc[m26].mean() - s.loc[m25h2].mean()) * 100 if len(m26) and len(m25h2) else np.nan
        seg26 = s.loc[m26].dropna()
        persist = bool(len(seg26) >= 4 and (seg26.diff().dropna() < 0).mean() >= 0.7)
        verdict = '平稳'
        if sl26 == sl26 and sl26 <= -0.8:
            verdict = '2026段下行(形态描述级)'
        if d_break == d_break and d_break < -3:
            verdict = '断点下移(先判工况: 2026起限电常态化, 逐月占比55-73% vs 2025间歇)'
        if sl26 == sl26 and sl26 < -0.8 and persist:
            verdict = '渐进退化候选(筛查级, 2026段内持续下行)'
        rows.append(dict(turbine=t, slope25=round(sl25, 2) if sl25 == sl25 else None,
                         slope26=round(sl26, 2) if sl26 == sl26 else None,
                         break_delta=round(d_break, 2) if d_break == d_break else None,
                         dev_last=round(float(s.dropna().iloc[-1]) * 100, 2) if len(s.dropna()) else None,
                         判=verdict))
    out = pd.DataFrame(rows)
    return out.sort_values('dev_last').reset_index(drop=True)

def attribution_router(cfg=None):
    """欠发候选逐轴排除 (SOP 1.3c; 机制候选级封顶)."""
    cfg = cfg or farm()
    st = pathlib.Path(cfg['store'])
    pcd = pd.read_parquet(st / 'powercurve_dev.parquet').set_index('turbine')
    att = pd.read_parquet(st / 'underperf_pitchmid.parquet').set_index('turbine')
    md = monthly_dev(cfg)
    rpiv = md.pivot_table(index='month', columns='turbine', values='rpm_ws')
    rdev = (rpiv - rpiv.median(axis=1).values[:, None]).mean()
    from ..subsys import temp_nbm
    treg = temp_nbm.registry()
    tdf = treg[0] if isinstance(treg, tuple) else treg
    out = []
    for t, r in pcd.iterrows():
        if r['判别'] == '—': continue
        axes = {}
        axes['风速计'] = f"ws偏置 {r.ws_bias:+.2f} m/s" + (' → A类先校风' if 'A类' in r['判别'] else ' (小, 非主因)')
        pd_dev = float(att.loc[t, 'pitch_dev']) if t in att.index else np.nan
        axes['桨角'] = f"1-2MW档 {pd_dev:+.2f}° vs 全场 → {'挂桨候选' if pd_dev == pd_dev and pd_dev > 0.5 else '排除(全场σ=0.2°)'}"
        rv = float(rdev.get(t, np.nan))
        axes['转速跟踪'] = f"中载转速偏 {rv:+.1f} rpm vs 全场 → {'跟踪偏低候选' if rv == rv and rv < -4 else '排除'}"
        # 第五轴 扇区形态 (R3, DS claim5 复现指令落库): 30°扇区×1m/s档 vs fleet同格,
        # ★绝对量因本台机舱风自指配档而放大(风速计偏置台尤甚), 只看**形态**不引绝对量
        sec_note = ''
        try:
            sp = pd.read_parquet(st / 'sector_power.parquet')
            me = sp[sp.turbine == t]
            rest = sp[sp.turbine != t].groupby(['sec', 'wsb']).p.median().rename('fleet')
            j = me.join(rest, on=['sec', 'wsb']).dropna(subset=['fleet'])
            if len(j) >= 60:
                j = j.assign(sdev=(j.p - j.fleet) / j.fleet)
                ps = j.groupby('sec').apply(lambda g: np.average(g.sdev, weights=g.n), include_groups=False)
                spread = float(ps.max() - ps.min()); best = float(ps.max())
                if spread >= 0.10 and best > -0.03:
                    axes['扇区形态'] = (f"强方向依赖(极差{spread:.0%}, 最好扇区{best:+.0%}≥0) → 尾流/布局效应候选(需机位图核), "
                                       f"机组自身欠发证据弱化")
                elif spread < 0.08:
                    axes['扇区形态'] = f"全向均匀(极差{spread:.0%}) → 与仪器/机组自身因一致, 非尾流形态"
                else:
                    axes['扇区形态'] = f"全向基底+个别扇区加深(极差{spread:.0%}, 最差{float(ps.min()):+.0%}@{int(ps.idxmin())}°) → 扇区分量存在, 布局核对可判"
        except Exception:
            pass
        brk = tdf[(tdf.turbine == t) & tdf.channel.str.contains('BrkTmp')]
        hot = brk[brk['dev_K'] >= 4]
        axes['刹车拖曳'] = ('候选: ' + '; '.join(f"{x.channel.split('_')[2]} {x.dev_K:+.1f}K" for _, x in hot.iterrows())) if len(hot) else '排除(刹车温正常)'
        open_axes = [k for k, v in axes.items() if '候选' in v and '排除' not in v and k != '扇区形态']
        wake = '尾流/布局' in axes.get('扇区形态', '')
        # 2026-08-26 外审(Gemini)逮逻辑错并采纳: "A类→不计欠发"混淆了数据事实与根因诊断 —
        # 欠发 -14.7% 是 SCADA 面客观事实, "风速计偏高"是它的根因; 损失仍计入, 归因传感器, 修复路径=校风
        # v3 (双模审两票收敛): dev = SCADA面自指口径的 fleet 相对偏差 — 是测量面事实;
        # 但机舱风自指 ⇒ **真实(来流面)性能欠发与否 INSUFFICIENT**, 不下绝对欠发结论 (Gemini v2 错误3 采纳)
        conclusion = ('SCADA面相对偏差成立(自指口径); 真实性能欠发与否 INSUFFICIENT(需独立风源) → 先校风复测, 校后重判'
                      if 'A类' in r['判别'] else \
            ('扇区形态=强方向依赖 → 主体或为尾流/布局效应(需机位图定谳), 机组自身欠发证据不足' if wake else
             (f"机制候选: {'/'.join(open_axes)}" if open_axes else
              '五轴皆无机组自身指向 → 余因(气动/自由来流) INSUFFICIENT, 需 mast/激光雷达')))
        out.append(dict(turbine=t, dev=f"{r.dev_w:+.1%}", 判别=r['判别'], 轴=axes, 结论=conclusion))
    return out