Compressor/models/vibration/vibration_algorithm.py

#!/usr/bin/env python3
"""PKS vibration anomaly-segment algorithm (vibration half of the scanner).

Vibration-only reconstruction of ``backend/scan_pks_fault_evolution.py``.  The
ported functions are byte-for-byte equivalent to the original so the exported
``振动异常段.csv`` stays identical.  Shared query/IO/statistics helpers live in
:mod:`models.common`.

CLI
---
::

    python models/vibration/vibration_algorithm.py \
        --hourly-csv models/datas/hourly/pks_hourly_full.csv \
        --output-dir models/results/vibration
"""
from __future__ import annotations

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

if __package__ in (None, ""):
    sys.path.insert(0, str(Path(__file__).resolve().parents[2]))

from models.common import db, devices, io, paths, queries
from models.common.constants import (
    MIN_HOUR_COVERAGE,
    VIBRATION_BASELINE_DAYS,
    VIBRATION_MERGE_GAP_HOURS,
    VIBRATION_RESTART_GAP_HOURS,
)
from models.common.stats import median_value, robust_baseline


def vibration_hour_signals(rows: list[dict], configs: dict[int, dict]):
    signals = []
    grouped = defaultdict(list)
    for row in rows:
        if row["运行覆盖率"] >= MIN_HOUR_COVERAGE:
            grouped[row["批次"]].append(row)
    for batch, series in grouped.items():
        series.sort(key=lambda row: row["_dt"])
        config = configs[batch]
        valid = []
        last_valid_dt = None
        for row in series:
            if row["运行覆盖率"] < MIN_HOUR_COVERAGE:
                continue
            restart = last_valid_dt is None or row["_dt"] - last_valid_dt > timedelta(hours=VIBRATION_RESTART_GAP_HOURS)
            if restart:
                # The first complete running hour after a restart is a
                # transition hour. Keep older valid hours for the long-term
                # robust reference, but do not evaluate this transition hour.
                last_valid_dt = row["_dt"]
                continue
            baseline_start = row["_dt"] - timedelta(days=VIBRATION_BASELINE_DAYS)
            prior = [item for item in valid if item["_dt"] >= baseline_start]
            side_signals = []
            for side in ("coupling", "chain"):
                avg_center, avg_scale = robust_baseline([item[f"{side}_avg"] for item in prior])
                max_center, max_scale = robust_baseline([item[f"{side}_max"] for item in prior])
                current_avg = row[f"{side}_avg"]
                current_max = row[f"{side}_max"]
                if current_max is None:
                    continue
                avg_z = ((current_avg - avg_center) / avg_scale) if avg_center is not None and current_avg is not None else 0.0
                peak_z = (current_max - max_center) / max_scale if max_center is not None and max_scale else 0.0
                running = max(row["运行样本数"], 1)
                high_ratio = (row[f"{side}_high_count"] or 0) / running
                high_high_ratio = (row[f"{side}_high_high_count"] or 0) / running
                limits = config[side]
                local_prior = valid[:6] if len(valid) >= 6 else []
                local_avg = median_value([item[f"{side}_avg"] for item in local_prior])
                local_max = median_value([item[f"{side}_max"] for item in local_prior])
                local_rise = False
                severe = (
                    high_high_ratio >= 0.01
                    or current_max >= limits["high_high"]
                    or peak_z >= 8
                )
                abnormal = (
                    high_ratio >= 0.01
                    or current_max >= limits["high"]
                    or peak_z >= 5
                    or avg_z >= 4
                    or local_rise
                )
                mild = peak_z >= 4 or avg_z >= 3 or local_rise
                if severe or abnormal or mild:
                    side_signals.append({
                        "side": side,
                        "avg_z": avg_z,
                        "peak_z": peak_z,
                        "high_ratio": high_ratio,
                        "high_high_ratio": high_high_ratio,
                        "severe": severe,
                        "abnormal": abnormal,
                        "value": current_max,
                        "avg": current_avg,
                        "prior_max_center": max_center,
                        "prior_max_scale": max_scale,
                        "prior_avg_center": avg_center,
                    })
            if side_signals:
                strongest = max(side_signals, key=lambda item: (item["severe"], item["peak_z"], item["avg_z"]))
                signal_level = 3 if any(item["severe"] for item in side_signals) else 2 if any(item["abnormal"] for item in side_signals) else 1
                reasons = []
                for item in side_signals:
                    label = "联轴器端" if item["side"] == "coupling" else "链轮端"
                    reasons.append(
                        f"{label}均值偏离{item['avg_z']:.1f}倍、峰值偏离{item['peak_z']:.1f}倍、"
                        f"高报警比例{item['high_ratio']:.1%}"
                    )
                signals.append({
                    "批次": batch,
                    "机组": config["unit"],
                    "_dt": row["_dt"],
                    "_level": signal_level,
                    "_severe": signal_level == 3,
                    "_side_signals": side_signals,
                    "_strongest": strongest,
                    "_reason": "；".join(reasons),
                })
            valid.append(row)
            last_valid_dt = row["_dt"]
    return signals


def merge_vibration_signals(signals: list[dict], configs: dict[int, dict]):
    grouped = defaultdict(list)
    for signal in signals:
        grouped[signal["批次"]].append(signal)
    events = []
    for batch, series in grouped.items():
        series.sort(key=lambda item: item["_dt"])
        current = []
        for signal in series:
            if not current or signal["_dt"] - current[-1]["_dt"] <= timedelta(hours=VIBRATION_MERGE_GAP_HOURS):
                current.append(signal)
            else:
                events.append(build_vibration_event(current, configs[batch]))
                current = [signal]
        if current:
            events.append(build_vibration_event(current, configs[batch]))

    events.sort(key=lambda row: (queries.unit_sort_key(row["机组"]), row["开始小时"], row["批次"]))
    unit_indexes = defaultdict(int)
    for event in events:
        unit_indexes[event["机组"]] += 1
        event["周期编号"] = f'{event["机组"]}-V{unit_indexes[event["机组"]]:02d}'

    filtered = []
    for event in events:
        # A single post-start relative-rise signal is not sufficient evidence;
        # absolute severe signals remain eligible on their own.
        if event["_严重信号小时数"] > 0 or event["_触发小时数"] >= 2:
            filtered.append(event)
    return sorted(filtered, key=io.output_sort_key)


def vibration_score(duration_hours: int, valid_hours: int, avg_change: float, max_value: float, peak_z: float, high_hours: int, high_high_hours: int, config: dict, side: str):
    duration_score = min(40.0, max(0.0, (duration_hours - 1) / 12 * 40))
    avg_score = min(20.0, max(0.0, (avg_change - 0.10) / 0.40 * 20)) if avg_change >= 0.10 else 0.0
    high = config[side]["high"]
    high_high = config[side]["high_high"]
    if max_value < high:
        max_score = min(10.0, max(0.0, (max_value / max(high, 1e-6) - 0.8) / 0.2 * 10))
    elif max_value < high_high:
        max_score = 10.0 + (max_value - high) / max(high_high - high, 1e-6) * 10.0
    else:
        max_score = 20.0
    peak_score = min(20.0, max(0.0, (peak_z - 4.0) / 8.0 * 20.0)) if peak_z >= 4.0 else 0.0
    total = duration_score + avg_score + max_score + peak_score
    hard_severe = max_value >= high_high or high_high_hours > 0 or (peak_z >= 8.0 and valid_hours >= 2)
    if hard_severe or total >= 60.0:
        level = "严重异常"
    elif total >= 30.0:
        level = "异常"
    else:
        level = "轻微"
    return total, level, duration_score, avg_score, max_score, peak_score


def build_vibration_event(signals: list[dict], config: dict):
    start = signals[0]["_dt"]
    end = signals[-1]["_dt"]
    all_side = [item for signal in signals for item in signal["_side_signals"]]
    strongest = max(all_side, key=lambda item: (item["severe"], item["peak_z"], item["value"]))
    main_side = "联轴器端" if strongest["side"] == "coupling" else "链轮端"
    max_value = max(item["value"] for item in all_side if item["value"] is not None)
    max_peak_z = max(item["peak_z"] for item in all_side)
    high_hours = sum(
        1 for signal in signals
        if any(item["high_ratio"] > 0 for item in signal["_side_signals"])
    )
    high_high_hours = sum(
        1 for signal in signals
        if any(item["high_high_ratio"] > 0 for item in signal["_side_signals"])
    )
    severe_hours = sum(1 for signal in signals if signal["_severe"])
    duration_hours = int((end - start).total_seconds() / 3600) + 1
    running_days = len({signal["_dt"].date() for signal in signals})
    side_avg_changes = []
    for side in ("coupling", "chain"):
        values = [item["avg"] for item in all_side if item["side"] == side and item["avg"] is not None]
        prior_values = [item["prior_avg_center"] for item in all_side if item["side"] == side and item["prior_avg_center"] is not None]
        if values and prior_values and median(prior_values) != 0:
            side_avg_changes.append((median(values) / median(prior_values) - 1, side))
    avg_change, _ = max(side_avg_changes, key=lambda item: item[0], default=(0.0, strongest["side"]))
    peak_ratio = max_value / max(strongest["prior_max_center"], 1e-6)

    if max_value >= config[strongest["side"]]["high_high"] or severe_hours > 0 and (high_high_hours > 0 or max_peak_z >= 8):
        level = "严重异常"
    elif max_value >= config[strongest["side"]]["high"] and high_hours >= 2:
        level = "严重异常"
    elif max_value >= config[strongest["side"]]["high"] or max_peak_z >= 5 or avg_change >= 0.10:
        level = "异常"
    else:
        level = "轻微"

    if peak_ratio >= 1.5 and avg_change < 0.25:
        shape = "间歇冲击型"
    elif avg_change >= 0.10 and len(signals) >= 3:
        shape = "持续抬升型"
    else:
        shape = "混合型"

    peak_signal = max(signals, key=lambda signal: signal["_strongest"]["peak_z"])
    first_signal = signals[0]
    last_signal = signals[-1]
    first_value = first_signal["_strongest"]["value"]
    last_value = last_signal["_strongest"]["value"]
    baseline = strongest["prior_max_center"]
    scale = max(strongest.get("prior_max_scale", 0.0), 0.0)
    normal_band = f"{baseline - 2 * scale:.3f}～{baseline + 2 * scale:.3f} mm/s" if baseline is not None else ""
    consecutive = 0
    longest_consecutive = 0
    for signal in signals:
        if signal["_level"] >= 2:
            consecutive += 1
            longest_consecutive = max(longest_consecutive, consecutive)
        else:
            consecutive = 0
    score, level, duration_score, avg_score, max_score, peak_score = vibration_score(
        duration_hours, len(signals), avg_change, max_value, max_peak_z,
        high_hours, high_high_hours, config, strongest["side"]
    )
    confirmed = "是" if level != "轻微" else "否"
    if shape == "持续抬升型":
        basis = f"{main_side}振动均值持续抬升，阶段内{longest_consecutive or len(signals)}个有效运行小时触发"
    elif shape == "间歇冲击型":
        basis = f"{main_side}重复出现振动冲击，阶段内{len(signals)}个有效运行小时触发"
    else:
        basis = f"{main_side}振动峰值或均值偏离近期基线，阶段内{len(signals)}个有效运行小时触发"
    if high_high_hours:
        basis += f"；高高报警小时{high_high_hours}个"
    elif high_hours:
        basis += f"；高报警小时{high_hours}个"
    basis += f"；最大振动值{max_value:.3f}mm/s，峰值高于近期基线{max_peak_z:.2f}个稳健尺度（峰值/基线{peak_ratio:.2f}倍）；综合评分{score:.1f}分"
    return {
        "机组": config["unit"],
        "批次": config["batch"],
        "周期编号": "",
        "振动阶段": level,
        "开始小时": start.strftime("%Y-%m-%d %H:%M:%S"),
        "结束小时": end.strftime("%Y-%m-%d %H:%M:%S"),
        "持续自然小时数": duration_hours,
        "有效运行小时数": len(signals),
        "近期基线振动": baseline,
        "正常波动带": normal_band,
        "阶段开始小时振动": first_value,
        "阶段结束小时振动": last_value,
        "阶段最大振动值": max_value,
        "阶段开始相对近期基线变化": (first_value / baseline - 1) if baseline else None,
        "阶段结束相对近期基线变化": (last_value / baseline - 1) if baseline else None,
        "阶段最大连续异常有效小时数": longest_consecutive,
        "阶段综合评分": score,
        "持续时间得分": duration_score,
        "平均振动得分": avg_score,
        "最大振动值得分": max_score,
        "峰值偏离得分": peak_score,
        "最强异常侧": main_side,
        "最强侧近期峰值基线": baseline,
        "最强侧近期稳健尺度": scale,
        "阶段最大峰值/基线比例": peak_ratio,
        "阶段高报警小时数": high_hours,
        "阶段高高报警小时数": high_high_hours,
        "_严重信号小时数": severe_hours,
        "_触发小时数": len(signals),
        "是否形成确认等级": confirmed,
        "阶段触发依据": basis,
    }


VIBRATION_FIELDS = [
    "机组", "批次", "周期编号", "振动阶段", "开始小时", "结束小时", "持续自然小时数", "有效运行小时数",
    "近期基线振动", "正常波动带", "阶段开始小时振动", "阶段结束小时振动", "阶段最大振动值",
    "阶段开始相对近期基线变化", "阶段结束相对近期基线变化", "阶段最大连续异常有效小时数",
    "阶段高报警小时数", "阶段高高报警小时数", "是否形成确认等级", "阶段触发依据",
    "阶段综合评分", "持续时间得分", "平均振动得分", "最大振动值得分", "峰值偏离得分",
    "最强异常侧", "最强侧近期峰值基线", "最强侧近期稳健尺度", "阶段最大峰值/基线比例",
]


def run(rows: list[dict], configs: dict[int, dict]):
    """Return vibration anomaly events."""
    signals = vibration_hour_signals(rows, configs)
    return merge_vibration_signals(signals, configs)


def write_results(output_dir: Path, events):
    output_dir = Path(output_dir)
    io.write_csv(output_dir / "振动异常段.csv", io.clean_rows(events), VIBRATION_FIELDS, {field: field for field in VIBRATION_FIELDS})


def load_hourly(args, selected_batches, unit_by_batch, configs):
    if args.hourly_csv:
        selected_set = set(selected_batches)
        return [row for row in queries.load_hourly_csv(Path(args.hourly_csv)) if row["批次"] in selected_set]
    start, end = queries.parse_dt(args.start), queries.parse_dt(args.end)
    return queries.fetch_all_hourly(start, end, args.chunk_days, unit_by_batch, configs, selected_batches)


def main() -> None:
    parser = argparse.ArgumentParser(description="PKS 振动异常段算法（独立于压力）")
    parser.add_argument("--device", nargs="+", default=list(devices.DEVICE_CODES), help="压缩机编号，如 7# 8# 9#")
    parser.add_argument("--start", default="2025-01-01 00:00:00")
    parser.add_argument("--end", default="2026-09-01 00:00:00")
    parser.add_argument("--chunk-days", type=int, default=14)
    parser.add_argument("--hourly-csv", default=None)
    parser.add_argument("--output-dir", default=str(paths.VIBRATION_RESULTS))
    args = parser.parse_args()
    devices_selected = devices.normalize_devices(args.device)
    device_batches = devices.load_device_batches(devices_selected)
    unit_by_batch = devices.batch_to_unit(device_batches)
    configs = db.load_alarm_config(unit_by_batch)
    selected_batches = sorted(unit_by_batch)
    rows = load_hourly(args, selected_batches, unit_by_batch, configs)
    events = run(rows, configs)
    write_results(Path(args.output_dir), events)
    print(f"设备：{devices_selected}")
    print(f"小时特征：{len(rows):,} 行")
    print(f"振动异常段：{len(events):,} 段")
    print(f"输出目录：{args.output_dir}")


if __name__ == "__main__":
    main()