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

# -*- coding: utf-8 -*-
"""A翼 故障分析件 (windscada M8; SOP 5.1/5.2/5.3 锚定).
轴:
  F1 停机事件化: 10min 五态序列 → 连续停机段 (≥1h) → stop_events.parquet; 调度令段另标;
  F2 停机归因: 段起点 ±30min 窗内首发报警码 = 归因码 (无码=无记录停机); 帕累托按"停机时长"加权 (SOP 5.2 禁 raw 计数单口径);
  F3 MTBF/MDT (GB/Z 178 口径族): SCADA 全停机口径 (无工单不可拆计划/故障 → 【参考级】+口径声明, SOP 5.1 定论须台账);
  F4 季节性: 报警族×月 (2025-01→2026-07, 19个月跨2冬1夏; 显著性未做χ²前=参考).
边界: 事件时长=状态段口径 (非报警闩锁口径); 停机含计划检修不可分; 调度令停机单列不入 MTBF 分子."""
import numpy as np, pandas as pd, pathlib
from src.windscada.config import farm
from app_ETL.app_ETL_guanlan.api import load_10min
from .curtail import classify

def build_stop_events(cfg=None):
    cfg = cfg or farm()
    rows = []
    al = pd.read_parquet(pathlib.Path(cfg['store']) / 'alarms.parquet')
    for t in cfg['turbines']:
        try:
            d = load_10min(t, cfg, groups=['A.功率', 'A.风况', 'A.转速', 'A.状态']).sort_values('ts').reset_index(drop=True)
            d['state'] = classify(d, cfg['rated_kw'])
            stop = d['state'].eq('停机').astype(int)
            seg = (stop.diff() != 0).cumsum()
            alt = al[al.turbine == t]
            for _, g in d[stop == 1].groupby(seg[stop == 1]):
                if len(g) < 6: continue          # ≥1h
                t0, t1 = g.ts.iloc[0], g.ts.iloc[-1] + pd.Timedelta(minutes=10)
                dur_h = (t1 - t0).total_seconds() / 3600
                # 两级归因 (2026-08-26 外审逮"无记录38%"→成因调查后改):
                #   T1 紧窗 ±30min 首发码 = 跳机型停机的直接归因 (高置信)
                #   T2 宽窗 前6h 首发码 = 长停型 (待料/待船窗口: 报警首发后进入等待, 期间无新报警)
                #      实测: 紧窗未命中的段里 53% 在前6h内有报警; 中位段长2.5h, 最长174h
                win = alt[(alt.t_on >= t0 - pd.Timedelta(minutes=30)) & (alt.t_on <= t0 + pd.Timedelta(minutes=30))]
                cause_code, cause_n, tier = '', 0, ''
                if len(win):
                    first = win.sort_values('t_on').iloc[0]
                    cause_code = str(first.code); cause_n = len(win); tier = 'T1紧窗'
                else:
                    w2 = alt[(alt.t_on >= t0 - pd.Timedelta(hours=6)) & (alt.t_on <= t0 + pd.Timedelta(minutes=30))]
                    if len(w2):
                        last = w2.sort_values('t_on').iloc[-1]      # 取最近一条 (最可能是触发)
                        cause_code = str(last.code); cause_n = len(w2); tier = 'T2宽窗(前6h)'
                # 调度令判别: 窗内含 1007/1008
                disp = bool(len(win) and win.code.isin(['1007', '1008']).any())
                rows.append(dict(turbine=t, t_start=t0, t_end=t1, dur_h=round(dur_h, 2),
                                 cause_code=cause_code, n_alarm_win=cause_n, tier=tier, dispatch=disp))
            print(t, flush=True)
        except Exception as e:
            print(t, 'ERR', str(e)[:60], flush=True)
    df = pd.DataFrame(rows)
    df.to_parquet(pathlib.Path(cfg['store']) / 'stop_events.parquet')
    print('→ stop_events.parquet', df.shape)

def _win_filter(df, months):
    m = df['t_start'].dt.to_period('M').astype(str)
    return df[m.isin(months)]

def mtbf_summary(months, cfg=None):
    """GB/Z 178 口径族, SCADA 全停机口径【参考级】: 无工单不可拆计划/故障停机."""
    cfg = cfg or farm()
    ev = pd.read_parquet(pathlib.Path(cfg['store']) / 'stop_events.parquet')
    ev = _win_filter(ev, months)
    fault = ev[~ev.dispatch]
    n_t = len(cfg['turbines'])
    # ★分母必须用实测台时, 不能用理想日历 (2026-08-28 自审逮)。
    # 原式 len(months)*30.4*24*n_t 把残月当整月、把数据缺口当有数据:
    # 2026-01~07 理想日历 194,050 台时 vs 实测 165,475 台时 (差 17%) — 2026-07 只有 6 天,
    # 另有几个月覆盖 90~96%。分子虚高而分母(事件数)是实测的 ⇒ MTBO 虚高 18% (239h 实为 203h)。
    # 与故障统计的残月归一是同一类问题, 那里已修, 这里遗漏。
    cal_h = None
    try:
        lm = pd.read_parquet(pathlib.Path(cfg['store']) / 'loss_monthly.parquet')
        lm = lm[lm.month.astype(str).isin([str(m) for m in months])]
        if len(lm):
            cal_h = float(lm['rows_'].sum()) / 6      # 10min 行 → 台时
    except Exception:
        pass
    cal_src = '实测台时'
    if not cal_h:
        cal_h = len(months) * 30.4 * 24 * n_t
        cal_src = '理想日历(实测台时不可得)'
    stop_h = float(fault.dur_h.sum()); disp_h = float(ev[ev.dispatch].dur_h.sum())
    n_ev = int(len(fault))
    mtbf = (cal_h - stop_h - disp_h) / max(n_ev, 1)
    mdt = stop_h / max(n_ev, 1)
    # 剔除口径: 去 1020就地模式(检修代理) 与 1161切入边界段(判据风SecAnemo vs 控制器风的边界争议, 非故障)
    core = fault[~fault.cause_code.isin(['1020', '1161'])]
    core_h = float(core.dur_h.sum()); n_core = int(len(core))
    mtbf2 = (cal_h - stop_h - disp_h) / max(n_core, 1)
    per_t = fault.groupby('turbine').agg(次=('dur_h', 'size'), 时长h=('dur_h', 'sum')).sort_values('时长h', ascending=False)
    n_drop = n_ev - n_core
    return dict(口径=('平均停机间隔(MTBO口径): 停机=≥1h停机状态段(不含调度令/低风待机), 含计划检修不可拆 → '
                     '**不可与行业MTBF(纯故障口径)比**; 分母=窗内实测台时(非理想日历, 残月不按整月计); 真MTBF须2025-26工单定性故障【参考级】'),
                停机事件数=n_ev, 停机总时长h=round(stop_h, 0), 调度令时长h=round(disp_h, 0),
                统计台时=round(cal_h, 0), 台时口径=cal_src,
                MTBF_h=round(mtbf, 0), MDT_h=round(mdt, 1),
                剔除口径=dict(说明='剔1020就地(检修代理)+1161切入边界段', 事件数=n_core, 时长h=round(core_h, 0), MTBF_h=round(mtbf2, 0),
                              次数变化=f'{n_ev}→{n_core} (剔{n_drop}起, {n_drop/max(n_ev,1):.0%}) — 间隔升高主因=分母减小, 非可靠性变好'),
                台级top=[dict(t=i, 次=int(r.次), 时长h=round(float(r.时长h), 0)) for i, r in per_t.head(8).iterrows()])

def stop_pareto(months, cfg=None, topn=12):
    """停机归因帕累托 (按停机时长加权, 双口径列次数)."""
    cfg = cfg or farm()
    from src.windscada import i18n
    ev = pd.read_parquet(pathlib.Path(cfg['store']) / 'stop_events.parquet')
    ev = _win_filter(ev, months)
    ev = ev[~ev.dispatch]
    ev['cause'] = ev.cause_code.replace('', '无归因停机(前6h无报警; 实测98%段前72h内有报警=长停远窗, 非链路静默)')
    g = ev.groupby('cause').agg(时长h=('dur_h', 'sum'), 次=('dur_h', 'size')).sort_values('时长h', ascending=False)
    out = []
    for code, r in g.head(topn).iterrows():
        label = code if '无记录' in code else f"{code} {i18n.alarm_label(code)}"
        out.append(dict(k=label, h=round(float(r.时长h), 0), n=int(r.次)))
    return out

def seasonal(cfg=None, fams=None):
    """报警族×月 季节矩阵 (全程 19 个月; 参考级)."""
    cfg = cfg or farm()
    from src.windscada import i18n
    al = pd.read_parquet(pathlib.Path(cfg['store']) / 'alarms.parquet')
    al['month'] = al.t_on.dt.to_period('M').astype(str)
    top_codes = al.groupby('code').size().sort_values(ascending=False).head(8).index
    ms = sorted(al.month.unique())
    out = []
    for c in top_codes:
        g = al[al.code == c].groupby('month').size().reindex(ms).fillna(0)
        out.append(dict(name=f"{c} {i18n.alarm_label(c)[:14]}", months=ms, vals=[int(v) for v in g]))
    return out