|
|
@@ -0,0 +1,137 @@
|
|
|
+#!/usr/bin/env python3
|
|
|
+"""Scan every repair event for a single-cylinder slow drift that leads it.
|
|
|
+
|
|
|
+For each event it looks back --days before the repair, builds the daily
|
|
|
+peer-relative series (pressure_ratio / p_discharge / efficiency within the same
|
|
|
+stage and side) and the self series, fits a linear drift, and ranks by strength.
|
|
|
+Read-only; uses feature CSVs produced by extract_from_db.py.
|
|
|
+"""
|
|
|
+from __future__ import annotations
|
|
|
+
|
|
|
+import argparse
|
|
|
+import csv
|
|
|
+from collections import defaultdict
|
|
|
+from datetime import datetime, timedelta
|
|
|
+from pathlib import Path
|
|
|
+
|
|
|
+import numpy as np
|
|
|
+
|
|
|
+EVENTS = [
|
|
|
+ ("7", "2025-08-08", "气阀"),
|
|
|
+ ("7", "2025-10-07", "气阀"),
|
|
|
+ ("7", "2026-05-21", "气阀"),
|
|
|
+ ("8", "2025-02-17", "气阀"),
|
|
|
+ ("8", "2025-12-06", "活塞"),
|
|
|
+ ("8", "2026-03-10", "十字头"),
|
|
|
+ ("8", "2026-03-29", "气阀"),
|
|
|
+ ("8", "2026-04-01", "气阀"),
|
|
|
+ ("9", "2026-01-18", "活塞"),
|
|
|
+]
|
|
|
+STAGE_MEMBERS = {"一级": ("二", "四", "六"), "二级": ("三", "五"), "三级": ("一",)}
|
|
|
+CYL_STAGE = {"一": "三级", "二": "一级", "三": "二级", "四": "一级", "五": "二级", "六": "一级"}
|
|
|
+FEATURES = ("pressure_ratio", "p_discharge", "efficiency")
|
|
|
+
|
|
|
+
|
|
|
+def load(path):
|
|
|
+ rows = []
|
|
|
+ for row in csv.DictReader(open(path, encoding="utf-8-sig")):
|
|
|
+ row["day"] = row["sample_time"][:10]
|
|
|
+ for name in FEATURES:
|
|
|
+ try:
|
|
|
+ row[name] = float(row[name])
|
|
|
+ except (TypeError, ValueError):
|
|
|
+ row[name] = None
|
|
|
+ rows.append(row)
|
|
|
+ return rows
|
|
|
+
|
|
|
+
|
|
|
+def daily_median(rows, cylinder, side, feature):
|
|
|
+ groups = defaultdict(list)
|
|
|
+ for row in rows:
|
|
|
+ if row["cylinder"] == cylinder and row["side"] == side and row[feature] is not None:
|
|
|
+ groups[row["day"]].append(row[feature])
|
|
|
+ return {day: float(np.median(v)) for day, v in groups.items()}
|
|
|
+
|
|
|
+
|
|
|
+def peer_relative(rows, feature):
|
|
|
+ by_key = defaultdict(dict)
|
|
|
+ for row in rows:
|
|
|
+ if row[feature] is not None:
|
|
|
+ by_key[(row["day"], row["side"])][row["cylinder"]] = row[feature]
|
|
|
+ series = defaultdict(dict)
|
|
|
+ for (day, side), members in by_key.items():
|
|
|
+ clean = {c: v for c, v in members.items() if v is not None}
|
|
|
+ if len(clean) < 2:
|
|
|
+ continue
|
|
|
+ med = float(np.median(list(clean.values())))
|
|
|
+ if abs(med) < 1e-9:
|
|
|
+ continue
|
|
|
+ for cylinder, value in clean.items():
|
|
|
+ series[(cylinder, side)][day] = (value - med) / med
|
|
|
+ return series
|
|
|
+
|
|
|
+
|
|
|
+def fit(daily):
|
|
|
+ days = sorted(daily)
|
|
|
+ if len(days) < 10:
|
|
|
+ return None
|
|
|
+ dates = [datetime.strptime(d, "%Y-%m-%d") for d in days]
|
|
|
+ span = (dates[-1] - dates[0]).days
|
|
|
+ if span <= 0 or span / (len(dates) - 1) > 15:
|
|
|
+ return None
|
|
|
+ xd = np.array([(d - dates[0]).days for d in dates], dtype=float)
|
|
|
+ y = np.array([daily[d] for d in days], dtype=float)
|
|
|
+ if xd[-1] <= 0:
|
|
|
+ return None
|
|
|
+ slope, intercept = np.polyfit(xd, y, 1)
|
|
|
+ pred = slope * xd + intercept
|
|
|
+ ss_res = float(np.sum((y - pred) ** 2))
|
|
|
+ ss_tot = float(np.sum((y - y.mean()) ** 2))
|
|
|
+ r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 0.0
|
|
|
+ scale = float(np.median(np.abs(y - np.median(y)))) * 1.4826 or float(np.std(y))
|
|
|
+ if scale <= 1e-9:
|
|
|
+ return None
|
|
|
+ total_mads = abs(slope) * xd[-1] / scale # total drift over the window in MADs
|
|
|
+ return total_mads, r2, len(days), float(y[0]), float(y[-1])
|
|
|
+
|
|
|
+
|
|
|
+def main():
|
|
|
+ parser = argparse.ArgumentParser(description="扫描所有维修事件的单缸慢漂移")
|
|
|
+ parser.add_argument("--data-dir", default=str(Path(__file__).resolve().parents[1] / "cache" / "indicator"))
|
|
|
+ parser.add_argument("--days", type=int, default=180, help="维修前回看天数")
|
|
|
+ args = parser.parse_args()
|
|
|
+
|
|
|
+ cache = {}
|
|
|
+ for unit in ("7", "8", "9"):
|
|
|
+ path = Path(args.data_dir) / f"all_{unit}_2025_2026.csv"
|
|
|
+ cache[unit] = load(path) if path.is_file() else []
|
|
|
+
|
|
|
+ for unit, date, kind in EVENTS:
|
|
|
+ event_day = datetime.strptime(date, "%Y-%m-%d")
|
|
|
+ start = (event_day - timedelta(days=args.days)).strftime("%Y-%m-%d")
|
|
|
+ rows = [r for r in cache[unit] if start <= r["day"] <= date]
|
|
|
+ days = sorted({r["day"] for r in rows})
|
|
|
+ print(f"\n== {unit}号机 {date} {kind} 回看{args.days}天 数据日 {len(days)} ({days[0] if days else '-'}~{days[-1] if days else '-'})")
|
|
|
+ if len(days) < 8:
|
|
|
+ print(" 数据不足,跳过")
|
|
|
+ continue
|
|
|
+ ranked = []
|
|
|
+ for feature in FEATURES:
|
|
|
+ for (cylinder, side), daily in peer_relative(rows, feature).items():
|
|
|
+ result = fit(daily)
|
|
|
+ if result and result[0] >= 2.0 and result[1] >= 0.4:
|
|
|
+ ranked.append((result[0], f"peer {cylinder}缸{side} {feature}", result))
|
|
|
+ for cylinder in ("一", "二", "三", "四", "五", "六"):
|
|
|
+ for side in ("盖侧", "轴侧"):
|
|
|
+ result = fit(daily_median(rows, cylinder, side, feature))
|
|
|
+ if result and result[0] >= 3.0 and result[1] >= 0.5:
|
|
|
+ ranked.append((result[0], f"self {cylinder}缸{side} {feature}", result))
|
|
|
+ ranked.sort(reverse=True)
|
|
|
+ if not ranked:
|
|
|
+ print(" 无显著慢漂移")
|
|
|
+ for _, label, (mads, r2, n, first, last) in ranked[:5]:
|
|
|
+ print(f" {label:<34} 漂移{mads:4.1f}MAD R2={r2:.2f} (n={n}, {first:.3f}->{last:.3f})")
|
|
|
+
|
|
|
+
|
|
|
+if __name__ == "__main__":
|
|
|
+ main()
|