#!/usr/bin/env python3 """Backtest the cylinder-temperature scan against labelled 4#/5#/6# faults. Reads ``_tmp_temperature/out/温度小时特征.csv`` and ``温度异常段.csv`` and reports, for every temperature-related fault: * the maximum same-stage spread change / peer deviation in the 7 days before; * whether a mild/abnormal/severe segment hit the pre-fault window; * the lead time. It also prints normal-period percentiles and a small threshold sweep so the "how many degC is a fault" question can be answered from data instead of a hand-picked 3 degC rule. """ from __future__ import annotations import argparse import sys from datetime import timedelta from pathlib import Path import pandas as pd ROOT = Path(__file__).resolve().parents[1] if str(ROOT) not in sys.path: sys.path.insert(0, str(ROOT)) from models.common.db import connect OUT = Path(__file__).resolve().parent / "out" LOOKBACK_DAYS = 7 def load_events() -> pd.DataFrame: connection = connect() try: with connection.cursor() as cursor: cursor.execute( "SELECT id, unit_name, fault_date, fault_category, fault_description " "FROM compressor_fault " "WHERE unit_name IN ('4号机组','5号机组','6号机组') " "AND (fault_description LIKE %s OR fault_description LIKE %s " " OR fault_description LIKE %s OR fault_description LIKE %s) " "ORDER BY unit_name, fault_date", ("%缸温度%", "%排气温度%", "%排温%", "%温度异常%"), ) rows = cursor.fetchall() finally: connection.close() frame = pd.DataFrame(rows) frame["event_date"] = pd.to_datetime(frame["fault_date"]) return frame def main() -> None: parser = argparse.ArgumentParser() parser.add_argument("--dir", default=str(OUT)) args = parser.parse_args() directory = Path(args.dir) features = pd.read_csv(directory / "温度小时特征.csv", low_memory=False) features["dt"] = pd.to_datetime(features["小时"]) features["unit"] = features["机组"].str.replace("号机", "号机组", regex=False) for column in ("级差变化", "最大同级偏离", "最大同级偏低", "最大稳健z", "级差", "最高温度"): features[column] = pd.to_numeric(features[column], errors="coerce") segments = pd.read_csv(directory / "温度异常段.csv") segments = segments[segments["特征类型"] != "高频波动"].copy() segments["start"] = pd.to_datetime(segments["开始小时"]) segments["end"] = pd.to_datetime(segments["结束小时"]) segments["unit"] = segments["机组"].str.replace("号机", "号机组", regex=False) channel_path = directory / "温度通道异常段.csv" if channel_path.exists(): channels = pd.read_csv(channel_path) channels["start"] = pd.to_datetime(channels["开始小时"]) channels["end"] = pd.to_datetime(channels["结束小时"]) channels["unit"] = channels["机组"].str.replace("号机", "号机组", regex=False) else: channels = pd.DataFrame(columns=["unit", "start", "end", "异常类型", "点位描述"]) events = load_events() print(f"温度类故障事件：{len(events)} 起\n") def _closest_start(frame, moment): if frame.empty: return "" values = frame["start"].tolist() before = [value for value in values if value <= moment] chosen = max(before) if before else min(values) return chosen.strftime("%Y-%m-%d %H:%M") rows = [] for _, event in events.iterrows(): unit = event["unit_name"] date = event["event_date"] start = date - timedelta(days=LOOKBACK_DAYS) end = date + timedelta(days=1) window = features[(features["unit"] == unit) & (features["dt"] >= start) & (features["dt"] <= end)] hits = segments[(segments["unit"] == unit) & (segments["start"] <= end) & (segments["end"] >= start)] channel_hits = channels[(channels["unit"] == unit) & (channels["start"] <= end) & (channels["end"] >= start)] starts = [] if not hits.empty: starts.extend(hits["start"].tolist()) if not channel_hits.empty: starts.extend(channel_hits["start"].tolist()) if starts: before = [value for value in starts if value <= end] first = max(before) if before else min(starts) else: first = pd.NaT lead = (date - first).total_seconds() / 3600 if starts else None def peak(stage: str, column: str): part = window[window["级"] == stage][column].dropna() return round(float(part.max()), 2) if not part.empty else None def peak_row(stage: str, column: str): part = window[(window["级"] == stage) & window[column].notna()] if part.empty: return None row = part.loc[part[column].idxmax()] return row stage = "一级" if "一级" in str(event["fault_description"]) else "二级" if "二级" in str(event["fault_description"]) else "" stages = [stage] if stage else ["一级", "二级"] best = None for candidate in stages: peak_peer = peak(candidate, "最大同级偏离") peak_spread = peak(candidate, "级差变化") score = max(value for value in (peak_peer, peak_spread) if value is not None) if (peak_peer is not None or peak_spread is not None) else -1 if best is None or score > best[0]: best = (score, candidate, peak_peer, peak_spread) hit = not hits.empty or not channel_hits.empty status = "数据不足" if not window.empty and window["最高温度"].dropna().empty else ("命中" if hit else "未命中") rows.append({ "机组": unit, "故障日期": date.date().isoformat(), "描述": event["fault_description"], "关注级": best[1] if best else "", "最大同级偏离℃": best[2] if best else None, "最大级差变化℃": best[3] if best else None, "命中": status, "首个事件开始": first.strftime("%Y-%m-%d %H:%M") if not hits.empty else "", "提前小时": round(lead, 1) if lead is not None else "", "命中等级": "、".join(sorted(set(pd.concat([hits["温度阶段"]]).astype(str)))) if not hits.empty else ("通道异常" if not channel_hits.empty else ""), "通道异常命中": "是" if not channel_hits.empty else "否", "通道异常类型": "、".join(sorted(set(channel_hits["异常类型"].astype(str)))) if not channel_hits.empty else "", "通道异常开始": _closest_start(channel_hits, end), "同日命中": "是" if starts and first.date() == date.date() else "否", }) report = pd.DataFrame(rows) print(report.to_string(index=False)) report.to_csv(directory / "回测_温度事件.csv", index=False, encoding="utf-8-sig") # ---- normal-period distribution ---- event_mask = pd.Series(False, index=features.index) for _, event in events.iterrows(): start = event["event_date"] - timedelta(days=3) end = event["event_date"] + timedelta(days=3) event_mask |= (features["unit"] == event["unit_name"]) & (features["dt"] >= start) & (features["dt"] <= end) normal = features[~event_mask] print("\n正常运行期分位数（排除故障日±3天）：") for column in ("级差变化", "最大同级偏离", "最大同级偏低", "最大稳健z", "最高温度"): series = normal[column].dropna() if series.empty: continue print(f" {column:>6}: " f"P50={series.quantile(0.5):6.2f} P90={series.quantile(0.9):6.2f} " f"P99={series.quantile(0.99):6.2f} P99.5={series.quantile(0.995):6.2f} " f"P99.9={series.quantile(0.999):6.2f} max={series.max():6.2f}") # ---- threshold sweep ---- print("\n阈值扫描（事件召回 vs 正常期小时触发率）：") pre = [] for _, event in events.iterrows(): start = event["event_date"] - timedelta(days=LOOKBACK_DAYS) end = event["event_date"] + timedelta(days=1) part = features[(features["unit"] == event["unit_name"]) & (features["dt"] >= start) & (features["dt"] <= end)] pre.append(part) print(" peer\\spread | " + " ".join(f"{spread:>6}" for spread in (2, 3, 4, 5, 6, 8))) for peer in (2, 3, 4, 5, 6, 8, 10): cells = [] for spread in (2, 3, 4, 5, 6, 8): hit_count = 0 for part in pre: mask = (part["最大同级偏离"] >= peer) | (part["级差变化"] >= spread) if mask.any(): hit_count += 1 recall = hit_count / len(pre) if pre else 0 normal_mask = (normal["最大同级偏离"] >= peer) | (normal["级差变化"] >= spread) normal_rate = normal_mask.mean() if not normal.empty else 0 cells.append(f"{recall:.2f}/{normal_rate*100:.2f}%") print(f" {peer:>4} | " + " ".join(f"{cell:>7}" for cell in cells)) print(" 单元格 = 事件召回 / 正常小时触发率（百分比）") if __name__ == "__main__": main()