Compressor/backend/scan_pks_clean_evolution.py

#!/usr/bin/env python3
"""Simple PKS evolution report.

This scanner intentionally avoids the previous state-machine implementation.
It reports smoothed hourly features first, then derives oil-pressure stages and
vibration segments from those features.  It is read-only and uses the actual
site_point alarm limits for each machine.
"""
from __future__ import annotations

import argparse
import csv
import os
import sys
from collections import defaultdict
from datetime import datetime, timedelta
from pathlib import Path
from statistics import median

import pymysql

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from app.config import settings  # noqa: E402

BATCH_TO_UNIT = {30: "7号机", 31: "8号机", 32: "9号机"}
EXPECTED_PER_HOUR = 720
MIN_COVERAGE = 0.80
SMOOTH_HOURS = 6
CONFIRM_HOURS = 3
MERGE_GAP_HOURS = 6

SQL = """
SELECT
  FROM_UNIXTIME((UNIX_TIMESTAMP(sample_time) DIV 3600) * 3600) AS hour_start,
  COUNT(*) AS samples,
  SUM(YSJ_41 > 0) AS running_samples,
  AVG(CASE WHEN YSJ_41 > 0 THEN YSJ_5 END) AS oil_avg,
  MIN(CASE WHEN YSJ_41 > 0 THEN YSJ_5 END) AS oil_min,
  MAX(CASE WHEN YSJ_41 > 0 THEN YSJ_5 END) AS oil_max,
  AVG(CASE WHEN YSJ_41 > 0 THEN YSJ_10 END) AS coupling_avg,
  MAX(CASE WHEN YSJ_41 > 0 THEN YSJ_10 END) AS coupling_max,
  AVG(CASE WHEN YSJ_41 > 0 THEN YSJ_11 END) AS chain_avg,
  MAX(CASE WHEN YSJ_41 > 0 THEN YSJ_11 END) AS chain_max,
  SUM(YSJ_41 > 0 AND YSJ_5 <= %s) AS oil_low_count,
  SUM(YSJ_41 > 0 AND YSJ_5 <= %s) AS oil_low_low_count,
  SUM(YSJ_41 > 0 AND YSJ_10 >= %s) AS coupling_high_count,
  SUM(YSJ_41 > 0 AND YSJ_10 >= %s) AS coupling_high_high_count,
  SUM(YSJ_41 > 0 AND YSJ_11 >= %s) AS chain_high_count,
  SUM(YSJ_41 > 0 AND YSJ_11 >= %s) AS chain_high_high_count
FROM pks_long_sample
WHERE import_batch_id = %s AND sample_time >= %s AND sample_time < %s
GROUP BY FROM_UNIXTIME((UNIX_TIMESTAMP(sample_time) DIV 3600) * 3600)
ORDER BY hour_start
"""


def connect():
    return pymysql.connect(
        host=settings.db_host, port=settings.db_port, user=settings.db_user,
        password=settings.db_password, database=settings.db_name,
        charset="utf8mb4", cursorclass=pymysql.cursors.DictCursor,
        connect_timeout=settings.db_connect_timeout, read_timeout=600,
        write_timeout=60, autocommit=True,
    )


def alarm_value(row, alarm_type):
    for i in range(1, 5):
        if str(row.get(f"AlarmType{i}") or "").strip() == alarm_type:
            return float(row[f"AlarmLimit{i}"])
    raise RuntimeError(f"site_point 缺少 {alarm_type}")


def load_limits():
    names = [f"YSJ{unit}_{point}" for unit in (7, 8, 9) for point in (5, 10, 11)]
    conn = connect()
    try:
        with conn.cursor() as cur:
            marks = ",".join(["%s"] * len(names))
            cur.execute(
                "SELECT ItemName, ItemDescription, Units, AlarmType1, AlarmType2, AlarmType3, AlarmType4, "
                "AlarmLimit1, AlarmLimit2, AlarmLimit3, AlarmLimit4 "
                f"FROM site_point WHERE ItemName IN ({marks})", names,
            )
            rows = {r["ItemName"]: r for r in cur.fetchall()}
    finally:
        conn.close()
    limits = {}
    for batch, unit in BATCH_TO_UNIT.items():
        n = unit[0]
        oil = rows[f"YSJ{n}_5"]
        coupling = rows[f"YSJ{n}_10"]
        chain = rows[f"YSJ{n}_11"]
        limits[batch] = {
            "unit": unit,
            "oil_low": alarm_value(oil, "PVLow"),
            "oil_low_low": alarm_value(oil, "PVLowLow"),
            "coupling_high": alarm_value(coupling, "PVHigh"),
            "coupling_high_high": alarm_value(coupling, "PVHighHigh"),
            "chain_high": alarm_value(chain, "PVHigh"),
            "chain_high_high": alarm_value(chain, "PVHighHigh"),
        }
    return limits


def fetch(batch, start, end, chunk_days, limit):
    result = []
    conn = connect()
    try:
        with conn.cursor() as cur:
            cursor = start
            while cursor < end:
                nxt = min(cursor + timedelta(days=chunk_days), end)
                params = (limit["oil_low"], limit["oil_low_low"], limit["coupling_high"],
                          limit["coupling_high_high"], limit["chain_high"], limit["chain_high_high"],
                          batch, cursor, nxt)
                cur.execute("EXPLAIN " + SQL, params)
                plan = cur.fetchone()
                if not plan or plan.get("key") != "PRIMARY":
                    raise RuntimeError(f"EXPLAIN 未使用 PRIMARY: batch={batch}, chunk={cursor}")
                cur.execute(SQL, params)
                for raw in cur.fetchall():
                    row = {"机组": limit["unit"], "批次": batch,
                           "小时": raw["hour_start"].strftime("%Y-%m-%d %H:%M:%S"),
                           "_dt": raw["hour_start"], "样本数": int(raw["samples"] or 0),
                           "运行样本数": int(raw["running_samples"] or 0)}
                    row["运行覆盖率"] = row["运行样本数"] / EXPECTED_PER_HOUR
                    for key in ("oil_avg", "oil_min", "oil_max", "coupling_avg", "coupling_max",
                                "chain_avg", "chain_max", "oil_low_count", "oil_low_low_count",
                                "coupling_high_count", "coupling_high_high_count",
                                "chain_high_count", "chain_high_high_count"):
                        row[key] = float(raw[key]) if raw[key] is not None else None
                    result.append(row)
                cursor = nxt
    finally:
        conn.close()
    return result


def rolling_median(values):
    values = [v for v in values if v is not None]
    return median(values) if values else None


def slope(values):
    pairs = [(i, v) for i, v in enumerate(values) if v is not None]
    if len(pairs) < 3:
        return None
    xm = sum(i for i, _ in pairs) / len(pairs)
    ym = sum(v for _, v in pairs) / len(pairs)
    den = sum((i - xm) ** 2 for i, _ in pairs)
    return sum((i - xm) * (v - ym) for i, v in pairs) / den if den else None


def prepare(rows, limits):
    by_batch = defaultdict(list)
    for row in rows:
        if row["运行覆盖率"] >= MIN_COVERAGE and row["oil_avg"] is not None:
            by_batch[row["批次"]].append(row)
    output = []
    for batch, series in by_batch.items():
        series.sort(key=lambda r: r["_dt"])
        limit = limits[batch]
        history = []
        cycle = 0
        previous = None
        for row in series:
            if previous is None or row["_dt"] - previous > timedelta(hours=12):
                cycle += 1
                cycle_rows = []
            cycle_rows.append(row)
            baseline = rolling_median([r["oil_avg"] for r in cycle_rows[:24]])
            smooth = rolling_median([r["oil_avg"] for r in history[-(SMOOTH_HOURS - 1):]] + [row["oil_avg"]])
            trend = slope([r["oil_avg"] for r in history[-5:]] + [row["oil_avg"]])
            stable_values = [r["oil_avg"] for r in cycle_rows[:24]]
            stable_mad = median(abs(value - baseline) for value in stable_values) if baseline is not None else 0
            fluctuation_band = max(3 * stable_mad, 0.0015)
            low = smooth is not None and smooth <= limit["oil_low"]
            low_low = smooth is not None and smooth <= limit["oil_low_low"]
            mild_zone = baseline is not None and smooth <= baseline - fluctuation_band
            abnormal_zone = baseline is not None and smooth <= baseline - 2 * fluctuation_band
            current = (
                "严重异常" if low_low
                else "异常" if low or (abnormal_zone and trend is not None and trend < 0)
                else "轻微" if mild_zone and trend is not None and trend < 0
                else "正常"
            )
            row.update({"周期编号": f"{limit['unit']}-P{cycle:02d}", "周期基线油压": baseline,
                        "6小时油压中位数": smooth, "6小时油压趋势": trend,
                        "正常波动带": fluctuation_band,
                        "低报警阈值": limit["oil_low"], "低低报警阈值": limit["oil_low_low"],
                        "原始压力等级": current,
                        "低报警样本比例": (row["oil_low_count"] or 0) / max(row["运行样本数"], 1),
                        "低低报警样本比例": (row["oil_low_low_count"] or 0) / max(row["运行样本数"], 1)})
            output.append(row)
            history.append(row)
            previous = row["_dt"]
    return sorted(output, key=lambda r: (r["批次"], r["_dt"]))


def build_oil_segments(rows):
    result = []
    for batch, series in defaultdict(list).items():
        pass
    grouped = defaultdict(list)
    for row in rows:
        if row["原始压力等级"] != "正常":
            grouped[(row["批次"], row["周期编号"], row["原始压力等级"])].append(row)
    for (_batch, _cycle, level), series in grouped.items():
        series.sort(key=lambda r: r["_dt"])
        current = []
        for row in series:
            if not current or row["_dt"] - current[-1]["_dt"] <= timedelta(hours=MERGE_GAP_HOURS):
                current.append(row)
            else:
                result.append(oil_segment(current, level)); current = [row]
        if current:
            result.append(oil_segment(current, level))
    return sorted(result, key=lambda r: (r["机组"], r["开始小时"]))


def oil_segment(series, level):
    first, last = series[0], series[-1]
    confirmed = len(series) >= CONFIRM_HOURS or level == "严重异常"
    return {"机组": first["机组"], "批次": first["批次"], "周期编号": first["周期编号"],
            "压力阶段": level if confirmed else "轻微", "开始小时": first["小时"], "结束小时": last["小时"],
            "持续小时数": len(series), "周期基线油压": first["周期基线油压"],
            "阶段开始油压": first["6小时油压中位数"], "阶段结束油压": last["6小时油压中位数"],
            "阶段最低油压": min(r["6小时油压中位数"] for r in series),
            "低报警小时数": sum(r["6小时油压中位数"] <= r["低报警阈值"] for r in series),
            "低低报警小时数": sum(r["6小时油压中位数"] <= r["低低报警阈值"] for r in series),
            "是否形成确认等级": "是" if confirmed else "否",
            "阶段判定依据": f"6小时中位数{first['6小时油压中位数']:.3f}～{last['6小时油压中位数']:.3f}MPa，"
            f"连续{len(series)}个有效运行小时；基线波动带{first.get('正常波动带', 0):.3f}MPa；"
            f"{'达到低低报警阈值' if level == '严重异常' else '超过正常波动带并持续下降'}"}


def build_vibration_segments(rows, limits):
    signals = []
    by_batch = defaultdict(list)
    for row in rows:
        by_batch[row["批次"]].append(row)
    for batch, series in by_batch.items():
        series.sort(key=lambda r: r["_dt"])
        history = []
        limit = limits[batch]
        for row in series:
            for side in ("coupling", "chain"):
                avg = rolling_median([r[f"{side}_avg"] for r in history[-5:]] + [row[f"{side}_avg"]])
                peak = row[f"{side}_max"]
                base = rolling_median([r[f"{side}_avg"] for r in history[-24:]])
                if avg is None or peak is None or base is None:
                    continue
                if peak >= limit[f"{side}_high_high"] or avg >= base * 1.20:
                    level = "严重异常"
                elif peak >= limit[f"{side}_high"] or avg >= base * 1.10:
                    level = "异常"
                elif avg > base * 1.05:
                    level = "轻微"
                else:
                    continue
                signals.append({"机组": row["机组"], "批次": batch, "_dt": row["_dt"], "侧": side,
                                "等级": level, "均值": avg, "峰值": peak})
            history.append(row)
    result = []
    grouped = defaultdict(list)
    for signal in signals: grouped[(signal["批次"], signal["侧"])].append(signal)
    for (batch, side), series in grouped.items():
        series.sort(key=lambda r: r["_dt"]); current=[]
        for signal in series:
            if not current or signal["_dt"] - current[-1]["_dt"] <= timedelta(hours=MERGE_GAP_HOURS): current.append(signal)
            else: result.append(vibration_segment(current, limits[batch])); current=[signal]
        if current: result.append(vibration_segment(current, limits[batch]))
    return sorted(result, key=lambda r: (r["机组"], r["开始时间"]))


def vibration_segment(series, limit):
    level_order={"轻微":1,"异常":2,"严重异常":3}; first=series[0]; last=series[-1]
    strongest=max(series,key=lambda r:(level_order[r["等级"]],r["峰值"]))
    side_name="联轴器端" if first["侧"]=="coupling" else "链轮端"
    return {"机组":first["机组"],"批次":first["批次"],"开始时间":first["_dt"].strftime("%Y-%m-%d %H:%M:%S"),"结束时间":last["_dt"].strftime("%Y-%m-%d %H:%M:%S"),"峰值时间":strongest["_dt"].strftime("%Y-%m-%d %H:%M:%S"),"振动等级":max(series,key=lambda r:level_order[r["等级"]])["等级"],"主要振动侧":side_name,"持续小时数":int((last["_dt"]-first["_dt"]).total_seconds()/3600)+1,"触发小时数":len(series),"最大振动值":max(r["峰值"] for r in series),"高报警阈值":limit[f"{first['侧']}_high"],"高高报警阈值":limit[f"{first['侧']}_high_high"],"阶段判定依据":f"{side_name}6小时中位数/小时峰值触发，连续{len(series)}个有效运行小时"}


def write(path, rows, fields):
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", encoding="utf-8-sig", newline="") as f:
        w=csv.DictWriter(f,fieldnames=fields,extrasaction="ignore");w.writeheader()
        for row in rows:
            w.writerow({k: "" if row.get(k) is None else round(row[k],6) if isinstance(row.get(k),float) else row.get(k) for k in fields})


def main():
    p=argparse.ArgumentParser();p.add_argument("--batch",nargs="+",type=int,choices=sorted(BATCH_TO_UNIT),default=list(BATCH_TO_UNIT));p.add_argument("--start",default="2025-04-01 00:00:00");p.add_argument("--end",default="2026-05-01 00:00:00");p.add_argument("--chunk-days",type=int,default=14);p.add_argument("--output-dir",default="cache/pks_scanner/clean_evolution");a=p.parse_args();start=datetime.strptime(a.start,"%Y-%m-%d %H:%M:%S");end=datetime.strptime(a.end,"%Y-%m-%d %H:%M:%S")
    limits=load_limits(); rows=[]
    for batch in a.batch: rows.extend(fetch(batch,start,end,a.chunk_days,limits[batch]))
    rows=prepare(rows,limits); oil=build_oil_segments(rows); vib=build_vibration_segments(rows,limits); out=Path(a.output_dir)
    write(out/"润滑油压力小时特征.csv",rows,["机组","批次","周期编号","小时","运行样本数","运行覆盖率","油压均值","油压最小值","油压最大值","6小时油压中位数","6小时油压趋势","周期基线油压","正常波动带","低报警阈值","低低报警阈值","原始压力等级"])
    for r in rows:r["油压均值"]=r.pop("oil_avg");r["油压最小值"]=r.pop("oil_min");r["油压最大值"]=r.pop("oil_max")
    write(out/"润滑油压力小时特征.csv",rows,["机组","批次","周期编号","小时","运行样本数","运行覆盖率","油压均值","油压最小值","油压最大值","6小时油压中位数","6小时油压趋势","周期基线油压","低报警阈值","低低报警阈值","原始压力等级"])
    write(out/"润滑油压力小时阶段段.csv",oil,list(oil[0]) if oil else ["机组","压力阶段","开始小时","结束小时"]);write(out/"振动异常段.csv",vib,list(vib[0]) if vib else ["机组","振动等级","开始时间","结束时间"])
    print(f"小时特征：{len(rows):,} 行");print(f"油压阶段：{len(oil):,} 段");print(f"振动阶段：{len(vib):,} 段");print(f"输出：{out}")


if __name__ == "__main__": main()