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

# -*- coding: utf-8 -*-
"""A翼 特性曲线多镜头件 (windscada M11; SOP §4.6c 正本 + skill curve-diagnostic).
镜头 (X 轴优先功率=直测干净量, 机舱风只作形状/相对):
  L1 风速-功率      功率曲线/欠发/限电平台
  L2 功率-桨距      控制律/削峰/per机偏离 (+三叶互比: 集距/不平衡)
  L3 功率-发电机转速 转速饱和点/降功率型
  L4 转矩-转速      T=P/ω_gen, Kω² 跟踪
  L5 风速-风轮转速   TSR λ=ω_r·R/v, Region-2 跟踪
  L6 功率-Cp        气动效率; ★Betz 0.593 硬闸: 超限=测风读低不是气动神机
  L7 转速比 gen/rot  联轴器/传动链/编码器 (铭牌齿比 119.752)
铁律 (§4.6c): ①列 liveness 先验(nuniq>1 ∧ 有方差 ∧ stuck_frac<0.5 ∧ 物理域) ②X轴选干净参考量
  ③机舱风 self-ref → 只判形状/相对, 绝对达成率须现场测风 ④per机异常必算同型机群残差(禁眼估)
前置: 限电剥离 (只用正常发电态行, SOP §4.3); 窗=2025-H2 (与曲线判别窗一致)."""
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

WIN = ('2025-07-01', '2026-01-01')
R_ROTOR = 65.0                    # SWT-4.0-130 叶轮半径 m
A_ROTOR = np.pi * R_ROTOR ** 2    # 13273 m²
GEAR_I = 119.752                  # 铭牌齿比 (cms-tcm-swt4000)
BETZ = 0.593
WSB = np.arange(3.0, 15.5, 0.5)   # 风速档
PWB = np.arange(0, 4250, 250)     # 功率档
GRB = np.arange(600, 1760, 40)    # 发电机转速档
# 七镜头真正取用的列（slim10min 窄仓里必须有这些；缺列 ⇒ 该列 liveness 判否，不造数）
CURVE_COLS = ['grd_wtc_ActPower_mean', 'tur_wtc_PowerRef_endvalue', 'tur_wtc_SecAnemo_mean',
              'tur_wtc_GenRpm_mean', 'tur_wtc_MainSRpm_mean', 'tur_wtc_PitcPosA_mean',
              'tur_wtc_PitcPosB_mean', 'tur_wtc_PitcPosC_mean', 'tmp_wtc_AmbieTmp_mean']

def _liveness(sr, lo=None, hi=None):
    s = sr.dropna()
    if len(s) < 100: return False, '样本<100'
    nu = s.nunique()
    if nu <= 2: return False, f'恒值(nuniq={nu})'
    stuck = float(s.value_counts(normalize=True).iloc[0])
    if stuck > 0.5: return False, f'卡值{stuck:.0%}'
    if lo is not None and (s.median() < lo or s.median() > hi): return False, f'物理域外(中位{s.median():.1f})'
    return True, f'活(nuniq={nu}, 卡值{stuck:.0%})'

def build_store(cfg=None, span=None, write=True):
    """七镜头分箱件。`span=(起, 止)` 给了就按**所选时间窗**（含两端）重算（用户令 2026-09-21）。

    span 模式走 `slim10min` 窄仓（同一份 `load_10min` 的列子集 ⇒ 同列同值），并且 **不写盘**：
    按窗重算是服务时算给页面看的，不能覆盖正式产物 `curve_lenses.parquet`
    （它的窗是"限电前干净判别窗 2025-07~2026-01"，另有消费者 `physics_check`）。
    """
    cfg = cfg or farm()
    rows, live_log = [], []
    for t in cfg['turbines']:
        if span:
            from app_ETL.app_ETL_guanlan.api import slim
            d = slim.load(t, cfg, span=span, columns=CURVE_COLS)
        else:
            d = load_10min(t, cfg, groups=['A.功率', 'A.风况', 'A.转速', 'B.变桨', 'B.温度NBM'])
            d = d[(d.ts >= WIN[0]) & (d.ts < WIN[1])]
        d['state'] = classify(d, cfg['rated_kw'])
        d = d[d.state == '正常发电']          # ★限电剥离前置
        # ① 列 liveness 先验
        L = {}
        for c, dom in (('grd_wtc_ActPower_mean', (0, 4200)), ('tur_wtc_SecAnemo_mean', (2, 20)),
                       ('tur_wtc_GenRpm_mean', (500, 1800)), ('tur_wtc_MainSRpm_mean', (2, 20)),
                       ('tur_wtc_PitcPosA_mean', (-5, 95)), ('tur_wtc_PitcPosB_mean', (-5, 95)),
                       ('tur_wtc_PitcPosC_mean', (-5, 95)), ('tmp_wtc_AmbieTmp_mean', (-20, 45))):
            ok, why = _liveness(d[c], *dom) if c in d.columns else (False, '无此列')
            L[c] = ok; live_log.append(dict(turbine=t, col=c, ok=ok, why=why))
        P = d['grd_wtc_ActPower_mean']; V = d['tur_wtc_SecAnemo_mean']
        # 空气密度 (环温 + 标准气压; 海上无气压列 → 口径声明)
        rho = 101325.0 / (287.05 * (d['tmp_wtc_AmbieTmp_mean'] + 273.15)) if L['tmp_wtc_AmbieTmp_mean'] else pd.Series(1.225, index=d.index)
        d = d.assign(rho=rho,
                     cp=(P * 1000) / (0.5 * rho * A_ROTOR * V.clip(lower=0.5) ** 3),
                     tq=(P * 1000) / (2 * np.pi * d['tur_wtc_GenRpm_mean'].clip(lower=1) / 60) / 1000,   # kNm 发电机侧
                     lam=(2 * np.pi * d['tur_wtc_MainSRpm_mean'] / 60) * R_ROTOR / V.clip(lower=0.5),
                     ratio=d['tur_wtc_GenRpm_mean'] / d['tur_wtc_MainSRpm_mean'].clip(lower=0.1),
                     p3=d['tur_wtc_PitcPosA_mean'].combine(d['tur_wtc_PitcPosB_mean'], max).combine(d['tur_wtc_PitcPosC_mean'], max)
                        - d['tur_wtc_PitcPosA_mean'].combine(d['tur_wtc_PitcPosB_mean'], min).combine(d['tur_wtc_PitcPosC_mean'], min),
                     wsb=pd.cut(V, WSB, labels=WSB[:-1]).astype(float),
                     pwb=pd.cut(P, PWB, labels=PWB[:-1]).astype(float),
                     grb=pd.cut(d['tur_wtc_GenRpm_mean'], GRB, labels=GRB[:-1]).astype(float))
        for xcol, ycols in (('wsb', ['grd_wtc_ActPower_mean', 'cp', 'lam', 'tur_wtc_MainSRpm_mean']),
                            ('pwb', ['tur_wtc_PitcPosA_mean', 'p3', 'tur_wtc_GenRpm_mean', 'ratio', 'cp']),
                            ('grb', ['tq'])):
            g = d.groupby(xcol, observed=True)[ycols].median()
            n = d.groupby(xcol, observed=True).size()
            for x, r in g.iterrows():
                if n[x] < 12: continue
                for y in ycols:
                    v = r[y]
                    if v == v: rows.append(dict(turbine=t, lens=xcol, x=float(x), y=y, v=float(v), n=int(n[x])))
        if write:
            print(t, flush=True)
    df = pd.DataFrame(rows); st = pathlib.Path(cfg['store'])
    if write:
        df.to_parquet(st / 'curve_lenses.parquet')
        pd.DataFrame(live_log).to_parquet(st / 'curve_liveness.parquet')
        print('→ curve_lenses.parquet', df.shape)
    return df

# 显著门 = 物理绝对锚 (z 只说"相对同类偏离", 定级必须过绝对量; MAD≈0 的镜头会造 z=200 的伪离群)
ANCHOR = {'grd_wtc_ActPower_mean': (40, 'kW', '=1%额定'), 'tur_wtc_PitcPosA_mean': (0.3, '°', '=变桨分册报警门半值'),
          'p3': (0.2, '°', '=三叶不平衡可辨门'), 'tur_wtc_GenRpm_mean': (5, 'rpm', '≈0.3%额定转速'),
          'tq': (0.25, 'kNm', '≈1%额定转矩'), 'lam': (0.2, '—', '≈2%λ_opt'),
          'cp': (0.01, '—', '≈3%Cp峰'), 'ratio': (0.6, '—', '=0.5%齿比')}

def lens(xcol, ycol, cfg=None, store_df=None, win_label=None):
    """回 fleet 中位曲线 + 逐台残差 (§4.6c④: per机必算残差, 禁眼估).

    `store_df` 给了就用它（= 按**所选时间窗**重算出来的分箱件），否则读正式产物
    `curve_lenses.parquet`（固定判别窗）。`win_label` 只进图注，说明这张图的窗是什么。
    """
    cfg = cfg or farm()
    if store_df is not None:
        df = store_df
    else:
        df = pd.read_parquet(pathlib.Path(cfg['store']) / 'curve_lenses.parquet')
    g = df[(df.lens == xcol) & (df.y == ycol)]
    # ★工况段闸: 剔非发电段档 (0kW 档=切入/空转混合态, 转速散度极大; 低转速档同理)
    #   实逮: 未剔时 19#/29#/21#/31# 因单个 0kW 档 -66~-81rpm 被误判"转速跟踪偏低族",
    #   而真发电段 500~3500kW 差仅 ±5rpm ⇒ 单档极值经均值放大 = 伪信号 (残差改用逐档中位)
    if xcol == 'pwb': g = g[g.x >= 250]
    if xcol == 'grb': g = g[g.x >= 900]
    if not len(g): return None
    piv = g.pivot_table(index='x', columns='turbine', values='v')
    med = piv.median(axis=1)
    resid = piv.sub(med, axis=0)
    w = piv.notna().sum(axis=1)
    keep = med.index[w >= max(10, int(0.5 * piv.shape[1]))]
    piv, med, resid = piv.loc[keep], med.loc[keep], resid.loc[keep]
    rmean = resid.median()          # 逐档中位: 单档极值不许主导
    mad = float((rmean - rmean.median()).abs().median()) or 1e-9
    z = (rmean - rmean.median()) / (1.4826 * mad)
    thr, unit, thr_note = ANCHOR.get(ycol, (0, '', ''))
    out_t = {t: dict(resid=round(float(rmean[t]), 4), z=round(float(z[t]), 2))
             for t in piv.columns if abs(z[t]) >= 3 and abs(rmean[t]) >= thr}
    # 样本量必须随曲线出来 (2026-08-28 审核逮): parquet 本来就有 n 列, lens() 没传出去
    # → 19 张图全部无样本量声明, 违反第一性原理②"Sample size declared"。
    npv = g.pivot_table(index='x', columns='turbine', values='n', aggfunc='sum').reindex(keep)
    n_by_x = [int(v) for v in npv.sum(axis=1).fillna(0)]
    return dict(x=[float(v) for v in med.index], fleet=[float(v) for v in med],
                n_total=int(sum(n_by_x)), n_by_x=n_by_x, n_min_bin=int(min(n_by_x) if n_by_x else 0),
                n_turbines=int(piv.shape[1]), win=(win_label or f"{WIN[0]}~{WIN[1]}"),
                anchor=dict(门=thr, 单位=unit, 说明=thr_note), 离群=out_t,
                q1=[float(v) for v in piv.quantile(0.25, axis=1)], q3=[float(v) for v in piv.quantile(0.75, axis=1)],
                per_t={t: [None if v != v else float(v) for v in piv[t]] for t in piv.columns},
                resid={t: round(float(rmean[t]), 4) for t in piv.columns},
                z={t: round(float(z[t]), 2) for t in piv.columns})

def physics_check(cfg=None):
    """第一性原理硬闸: Betz + 齿比恒等 + λ 合理域."""
    cfg = cfg or farm()
    df = pd.read_parquet(pathlib.Path(cfg['store']) / 'curve_lenses.parquet')
    out = {}
    cp = df[(df.lens == 'wsb') & (df.y == 'cp')]
    cpmax = cp.groupby('turbine').v.max()
    over = cpmax[cpmax > BETZ]
    out['Betz'] = dict(判据=f'Cp ≤ {BETZ} (贝兹极限)', 全场最大Cp=round(float(cpmax.max()), 3),
                       中位峰值Cp=round(float(cpmax.median()), 3),
                       超限台=[dict(t=k, cp=round(float(v), 3)) for k, v in over.items()],
                       判=('超限=机舱风读低(NTF)或密度口径, 非气动超越' if len(over) else '未见荒谬输出 (sanity check: 本闸只防极端口径错误, 不构成功率曲线绝对准确性背书 — 机舱风自指, Cp绝对值不可信)'))
    ra = df[(df.lens == 'pwb') & (df.y == 'ratio')]
    rmed = ra.groupby('turbine').v.median()
    dev = (rmed - GEAR_I).abs()
    out['齿比恒等'] = dict(铭牌=GEAR_I, 全场中位=round(float(rmed.median()), 3),
                          离群台=[dict(t=k, r=round(float(rmed[k]), 2)) for k in dev[dev > 1.5].index],
                          判='破恒等=编码器/传感, 非联轴器打滑本身')
    lam = df[(df.lens == 'wsb') & (df.y == 'lam')]
    r2 = lam[(lam.x >= 5) & (lam.x <= 9)].groupby('turbine').v.median()          # Region-2 平台
    floor = lam[lam.x <= 3.5].groupby('turbine').v.median()
    plateau = float(r2.median())
    out['λ叶尖速比'] = dict(Region2平台_中位=round(plateau, 2), 台间范围=[round(float(r2.min()), 2), round(float(r2.max()), 2)],
                          切入段_中位=round(float(floor.median()), 2),
                          判=('Region-2 平台 λ≈%.1f ∈ 合理域 6~10 ✔; 切入段 λ 被最小转速地板顶高 (风轮 5.5rpm 恒定) 属设计非异常, '
                             '故 λ 判据取平台不取峰值' % plateau))
    return out