#!/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()