guanlan-src/app_algorithmModel/app_algorithmModel_guanlan/subsys/pitch.py

# -*- coding: utf-8 -*-
"""变桨系统四维度判级 (v1, 2026-08-24). 判据规范出处: 变桨分册 V2.0 (四维度/§3.2归并) + 零位专项 20260821 口径.
全部确定性; 模型无判级通道. 阈值: 分册场标定优先, 无者 fleet 相对 [暂行标注]. 数据止 2026-06-30 (边界).

★2026-09-19 (用户令「所有的计算均要形成观澜的源代码」): 本面的输入 `pitch/pitch_daily.parquet` 与
`pitch/pitch_zero_monthly.parquet` 过去是**随包快照**(无生成端) ⇒ 换台机器重算后必然缺件、本面恒"不可判"。
现在由 `scripts/pitch_face_build.py`(重算链第 ④c 步)从 `data/raw` 的 `scada_10min` + `scada_1min` 算出,
口径与证据写在同目录 `pitch_face_manifest.json`。零位轴按用户令**分列**: `zero_dev`(停机顺桨段, 机械止挡位)、
`zero_dev_full`(满发段, 并列参考)、`zero_dev_run`(运行段同工况分档 = 零位偏差的行业口径, 实测能把 19# 拎出来);
D3 判据按"停机段与运行段**取更不利**"(§3.2 同一纪律), 依据行三个数都写。
★★**用户令 2026-09-19 就"判据轴怎么定"的裁决 = 维持本口径**(三口径都存、依据行都写、判据取更不利的那一个):
既保留用户选定的"停机段中位"这一列, 又不会漏掉 19# 这类**运行区**零位偏差 —— 本机实测三个口径对 19# 分别
是 +0.11°(停机顺桨段, 机械止挡抹平差异) / +0.61°(满发段) / −0.46°(运行段同工况分档), 只有第三个把它拎出来。
改动口沿时**必须同步**: `scripts/pitch_face_build.py` 的口径注释与本文件这段说明。
"""
from app_common.app_common_guanlan.api import paths as P
import numpy as np, pandas as pd, pathlib, json
from src.windscada.config import ROOT

PITCH_DIR = P.pitch()
DIMS = ['蓄能与压力调节', '变桨轴承与轮毂润滑', '执行与位置反馈', '油路与密封']
# 场标定常数 (分册§5.1/5.2): 正常锯齿带宽≈52bar; 润滑泵全场中位171s/日, 310s+为高位; 零位专项: 精度门0.25°, 19#确诊-0.87°
NORM_BAND, HUBLUB_FLEET_MED = 52.0, 171.0
ZERO_GATE_STD, ZERO_ALARM, ZERO_WATCH = 0.25, 0.6, 0.3   # [暂行] 报警阈=3×2025离散度上界0.2°

def _cur_win(d, days=60, span=None):
    """判级窗：默认 = 数据末端往前 `days` 天（原口径）；`span=(a, b)` 给了就按**所选时间窗**（含两端）。

    ★2026-09-21 用户令「判级也按所选时间窗重算」：本面输入 `pitch_daily.parquet` 是**日粒度**序列
    （每台每日一行），所以按窗重算是**精确**的——只是换一个日期切片，判据不动。
    """
    if span:
        a, b = span
        lo, hi = pd.Timestamp(a).date(), pd.Timestamp(b).date()
        dd = pd.to_datetime(d['date']).dt.date
        return d[(dd >= lo) & (dd <= hi)]
    end = d['date'].max()
    return d[d['date'] >= (pd.Timestamp(end) - pd.Timedelta(days=days)).date()]


def _rate(v) -> str:
    """越线率排版: 缺 1min 计数来源时是 NaN ⇒ 写「—」, 不写成 0.000 (那是"没有越线"的意思)。"""
    return f"{v:.3f}" if v == v else '—'

def load():
    """读变桨面产物 —— **缺件时返回空表而不是抛异常** (2026-09-17, 见下)。

    生成端: `scripts/pitch_face_build.py` (重算链 ④c, 输入 = `data/raw/<场站>/{scada_10min,scada_1min}`)。
    为什么仍要容缺: ① 目标机上"放了原始件但还没重算"是常态; ② 现场没给 1min 导出件时, 零位两口径
    与越线计数会按缺 —— 缺件就该显示缺件, 不能造数。原先这里直接 `pd.read_parquet` ⇒ `FileNotFoundError`
    一路上抛到本体层 `taxonomy.system_matrix()`, 把"整条重算链"卡在第 ⑦ 步
    (2026-09-17 实测: 清产物后重算, 2267 s 后 rc=1 断在这里)。
    现在缺件返回空表, 下游按"该面无数据"降级(页面显示缺件, 不静默造数)。
    """
    p = PITCH_DIR / 'pitch_daily.parquet'
    if p.is_file():
        daily = pd.read_parquet(p)
    else:
        daily = pd.DataFrame(columns=['date', 'turbine', 'n', 'hydlevel_timeon', 'hydfilt_timeon'])
    zp = PITCH_DIR / 'pitch_zero_monthly.parquet'
    zero = pd.read_parquet(zp) if zp.is_file() else pd.DataFrame()
    return daily, zero

ALARM_FAM = {  # M4b 报警轴 (windscada alarms.parquet; 分册监测量对齐)
    '执行反馈': r'跟踪|pawl|反馈|停止位|制转杆|电磁阀|安全阀',
    '蓄能压力': r'高压|蓄能器|压力开关',
    '润滑': r'变桨润滑|无润滑|润滑警告',
    '油路密封': r'液压油位|液压油温|变桨液压',
}


def _alarm_counts(win_start, win_end, store=None):
    import re as _re
    # ★store 显式传入 (2026-09-16 修): 原写 `P.store()` —— 它跟的是环境变量 WINDSCADA_FARM,
    #   而调用方 registry(cfg) 拿的是**显式选定的场**; 多场部署下这里会读到 rudong 的 alarms.parquet
    #   而其余输入都来自当前场 ⇒ 跨场串数据, 且因为"文件存在、列名兼容"而完全不报错。
    ap = (pathlib.Path(store) if store else P.store()) / 'alarms.parquet'
    if not ap.exists():
        return None
    al = pd.read_parquet(ap)
    al = al[(al.t_on >= str(win_start)) & (al.t_on <= str(win_end)) & al.text.str.contains('桨|液压|润滑|高压|蓄能', na=False, regex=True)]
    out = {}
    for fam, pat in ALARM_FAM.items():
        sub = al[al.text.str.contains(pat, na=False, regex=True)]
        out[fam] = sub.groupby('turbine').size()
    return out


def registry(cfg=None, span=None):
    """变桨面登记表。`cfg` 给定时, 其 `store` 决定读哪个场的 alarms (见 _alarm_counts 的注释)。

    `span=(起, 止)` 给了就按**所选时间窗**（含两端）重算判级（用户令 2026-09-21）；不给 = 原口径
    （数据末端往前 60 天）。日粒度件 ⇒ 按窗重算是精确切片，判据/阈值一字未动。
    """
    daily, zero = load()
    if not len(daily):
        # 缺"无生成端"的 pitch_daily (交付包不随产物 ⇒ 目标机上必然缺): 该面整体**不可判**,
        # 但**结构要在**(列齐、零行), 否则下游 `reg.set_index('机组')` 会 KeyError, 又把整条链卡住。
        # 页面据此显示"变桨面缺件: outputs/<场>/pitch/pitch_daily.parquet (docs §7)", 不静默造数。
        return pd.DataFrame(columns=['机组', '系统级', *DIMS]), {}
    daily['date'] = pd.to_datetime(daily['date']).dt.date
    cur = _cur_win(daily, span=span).copy()
    # 语义修正 (2026-08-24 数据实逮): HydLevel/HydrFilt 开关 1=正常 → 报警秒 = 当日应测秒(n×600) − timeon;
    # 润滑泵/柱塞用窗内日均 (日中位恒0); PitchPum 使能位恒开=无判据力, 废弃
    cur['lvl_alarm'] = (cur['n'] * 600 - cur['hydlevel_timeon']).clip(lower=0)
    cur['filt_alarm'] = (cur['n'] * 600 - cur['hydfilt_timeon']).clip(lower=0)
    # 双通道同时近整日置零 = 断链/停机伪影日 (24# 实逮 05-08), 剔出油路报警秒
    _cm = (cur['lvl_alarm'] > 600) & (cur['lvl_alarm'] == cur['filt_alarm'])   # 两独立开关逐秒相等=共因伪影(断链/停机)
    cur.loc[_cm, ['lvl_alarm', 'filt_alarm']] = 0.0
    fleet = cur.groupby('turbine').agg(band=('hyd_band', 'median'), over=('hyd_over_relief', 'sum'), under=('hyd_under_pump', 'sum'),
                                       lub=('hublub_timeon', 'mean'), pump=('pitchpum_timeon', 'median'),
                                       lvl=('lvl_alarm', 'sum'), filt=('filt_alarm', 'sum'),
                                       pa=('pistA', 'mean'), pb=('pistB', 'mean'), pc=('pistC', 'mean'),
                                       nd=('n', 'count'), nop=('n_op', 'sum'))
    fleet['over_rate'] = fleet['over'] / (fleet['nop'] * 10).clip(lower=1)   # 越线**分钟**/运行分钟 (工况归一)
    fleet['under_rate'] = fleet['under'] / (fleet['nop'] * 10).clip(lower=1)
    med = fleet.median()
    # 零位: 精度门逐月 (全场逐台偏差 std ≤0.25° 才可判)
    # ★2026-09-19: 输入件现由 scripts/pitch_face_build.py 自算, 零位有**三个口径**(见该脚本 docstring):
    #   停机顺桨段 `zero_dev`(机械止挡位, 量不出运行区零位差) / 满发段 `zero_dev_full`(混风况, 参考) /
    #   运行段分档 `zero_dev_run`(同工况读数差 = 行业口径的零位, 分档归一后能把 19# 这类台单独拎出来)。
    #   判据按"**取更不利**"(§3.2 同一纪律): 停机段与运行段谁偏得多用谁, 依据行三个数都写出来。
    zsum = {}
    if len(zero):
        mstd = zero.groupby('month')['zero_dev'].std()
        ok_months = mstd[mstd <= ZERO_GATE_STD].index
        zok = zero[zero['month'].isin(ok_months)]
        for t, g in zok.groupby('turbine'):
            g = g.sort_values('month')
            _full = float(g['zero_dev_full'].median()) if 'zero_dev_full' in g.columns else float('nan')
            _run = float(g['zero_dev_run'].median()) if 'zero_dev_run' in g.columns else float('nan')
            _n_alarm = int((g['zero_dev'].abs() >= ZERO_ALARM).sum())
            if _run == _run:
                _n_alarm = int(((g['zero_dev'].abs() >= ZERO_ALARM)
                                | (g['zero_dev_run'].abs() >= ZERO_ALARM)).sum())
            zsum[t] = dict(months=int(len(g)), dev_med=float(g['zero_dev'].median()), last=float(g['zero_dev'].iloc[-1]),
                           dev_full_med=_full, dev_run_med=_run, n_alarm=_n_alarm,
                           last_month=str(g['month'].iloc[-1]))
    # 报警轴同窗（用户令 2026-09-21：按所选窗重算）：span 给了就用 span 的两端，
    # 否则仍是"数据末端往前 60 天"（原口径）。
    if span:
        win_start, win_end = pd.Timestamp(span[0]), pd.Timestamp(span[1])
    else:
        win_end = daily['date'].max(); win_start = pd.Timestamp(win_end) - pd.Timedelta(days=60)
    ac = _alarm_counts(win_start, win_end, store=(cfg or {}).get('store') if cfg else None)
    rows, detail = [], {}
    for t, r in fleet.iterrows():
        d = dict(维度={}, 依据={})
        blind = (r['nd'] < 10) or (r['nop'] < 144)   # 现窗运行不足1天等效 → 不可判
        # D1 蓄能与压力调节
        band_x = r['band'] / max(med['band'], 1e-9)
        st1, why1 = '优秀', f"带宽中位 {r['band']:.0f}bar ({band_x:.2f}×fleet, 正常参考≈{NORM_BAND:.0f}bar 1min口径)"
        if r['band'] >= 2.0 * med['band']:
            st1 = '报警'; why1 += (f"; 越卸荷 {_rate(r['over_rate'])}/跌启泵 {_rate(r['under_rate'])}"
                                  f"(运行分钟占比)作佐证 [暂行: 带宽≥2×单证报警; 蓄能失效佐证=锯齿失稳]")
        elif r['band'] >= 1.4 * med['band']:
            st1 = '良好'
        # D2 润滑
        lub_x = r['lub'] / max(med['lub'], 1e-9)
        st2, why2 = '优秀', f"轮毂润滑泵日均 {r['lub']:.0f}s ({lub_x:.2f}×fleet中位{med['lub']:.0f}s; 分册: 全场中位171s/高位310s)"
        pist = [r['pa'], r['pb'], r['pc']]
        if all(x == x for x in pist) and min(pist) > 0 and max(pist) / min(pist) > 1.5:
            why2 += f"; 三叶柱塞不对称 {max(pist)/min(pist):.2f}×"
        if r['lub'] >= 1.8 * med['lub']:
            st2 = '报警'
        elif r['lub'] >= 1.4 * med['lub'] or (all(x == x for x in pist) and min(pist) > 0 and max(pist) / min(pist) > 1.5):
            st2 = '良好'
        # D3 执行与位置反馈 (零位轴并入 + 带宽骤变签名)
        st3, why3 = '优秀', ''
        z = zsum.get(t)
        if z:
            _stop, _full, _run = z['dev_med'], z.get('dev_full_med'), z.get('dev_run_med')
            _mag = abs(_stop)
            _lead = '停机顺桨段'
            if _run == _run and abs(_run) > _mag:
                _mag, _lead = abs(_run), '运行段'
            why3 = f"零位偏差 停机顺桨段 {_stop:+.2f}°"
            if _full == _full:
                why3 += f" / 满发段 {_full:+.2f}°"
            if _run == _run:
                why3 += f" / 运行段(同工况) {_run:+.2f}°"
            why3 += f" (判据轴取更不利的 {_lead}, 可判月 {z['months']}, 精度门0.25°)"
            if z['n_alarm'] >= 3 and _mag >= ZERO_ALARM:
                st3 = '报警'; why3 += f" ≥{ZERO_ALARM}°持续{z['n_alarm']}月 → 偏差确诊·成因待诊 [专项口径]"
            elif _mag >= ZERO_WATCH:
                st3 = '良好'
            if _mag >= ZERO_ALARM and abs(z['last']) < ZERO_WATCH:
                st3 = '良好'; why3 += ' (末月已归零→已闭环观察)'
        else:
            why3 = '零位: 可判月不足 (精度门/样本)'
        prev = daily[(daily['turbine'] == t) & (~daily['date'].isin(cur['date']))]
        if len(prev) > 20:
            pb_ = prev['hyd_band'].median()
            if pb_ > 0 and (r['band'] / pb_ < 0.5):
                st3 = st3 if st3 == '报警' else '良好'
                why3 += f"; 带宽骤降 {pb_:.0f}→{r['band']:.0f}bar (整定改动签名)"
        # D4 油路与密封
        st4, why4 = '优秀', f"油位报警 {r['lvl']:.0f}s/窗, 滤网报警 {r['filt']:.0f}s/窗 (开关1=正常, 报警秒=应测−timeon; 变桨泵使能位恒开无判据力已弃)"
        if r['lvl'] > 3600 or r['filt'] > 3600:
            st4 = '报警' if max(r['lvl'], r['filt']) > 6 * 3600 else '良好'
        # M4b 报警轴佐证: 族内计数 ≥10 且 ≥3×fleet中位 → 该维度升报警 (分册监测量的报警面, 7#/38# 类)
        if ac is not None:
            for fam, st_ref in (('蓄能压力', 'st1'), ('润滑', 'st2'), ('执行反馈', 'st3'), ('油路密封', 'st4')):
                cnt = int(ac[fam].get(t, 0)); fam_med = float(ac[fam].median()) if len(ac[fam]) else 0.0
                if cnt >= 10 and cnt >= 3 * max(fam_med, 1):
                    if fam == '蓄能压力': st1 = '报警'; why1 += f'; 报警轴 {cnt}条 (fleet中位{fam_med:.0f})'
                    if fam == '润滑': st2 = '报警'; why2 += f'; 报警轴 {cnt}条 (fleet中位{fam_med:.0f})'
                    if fam == '执行反馈': st3 = '报警'; why3 += f'; 报警轴 {cnt}条 (fleet中位{fam_med:.0f})'
                    if fam == '油路密封': st4 = '报警'; why4 += f'; 报警轴 {cnt}条 (fleet中位{fam_med:.0f})'
        for dim, st, why in zip(DIMS, (st1, st2, st3, st4), (why1, why2, why3, why4)):
            d['维度'][dim] = '不可判' if blind else st
            d['依据'][dim] = why
        order = {'报警': 0, '良好': 1, '优秀': 2, '不可判': 3}
        judged = [v for v in d['维度'].values() if v != '不可判']
        sys_st = sorted(judged, key=lambda x: order[x])[0] if judged else '不可判'   # §3.2: 取最不利; 不可判不传染
        rows.append(dict(机组=t, 系统级=sys_st, **d['维度']))
        detail[t] = d
    reg = pd.DataFrame(rows).sort_values('机组').reset_index(drop=True)
    return reg, detail

def reconcile_delivered():
    ref = json.load(open(ROOT / 'reference/rudong/pitch_delivered_grades.json', encoding='utf-8'))
    reg, _ = registry()
    ix = reg.set_index('机组')
    m = {'蓄能压力': '蓄能与压力调节', '润滑': '变桨轴承与轮毂润滑', '执行反馈': '执行与位置反馈', '油路密封': '油路与密封'}
    rows = []
    for t, v in ref.items():
        if t not in ix.index: continue
        for k, dim in m.items():
            a, b = v[k], ix.loc[t, dim]
            rows.append(dict(台=t, 维度=k, 交付=a, 现算=b, 一致=(a == '报警') == (b == '报警')))
    df = pd.DataFrame(rows)
    return df