Compressor/models/generalized/pressure_algorithm.py

#!/usr/bin/env python3
"""PKS lubricating-oil-pressure evolution algorithm (pressure half of the scanner).

This is the pressure-only reconstruction of
``backend/scan_pks_fault_evolution.py``.  Only the code path that actually
produces the exported pressure results is retained, copied verbatim so the
output stays identical to the original scanner.  Database/query/IO/statistics
helpers live in :mod:`models.common`.

CLI
---
::

    python models/generalized/pressure_algorithm.py \
        --start 2025-01-01 --end 2026-09-01 \
        --hourly-csv models/datas/hourly/pks_hourly_full.csv \
        --output-dir models/results/pressure
"""
from __future__ import annotations

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

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 (
    HOURLY_BASELINE_HOURS,
    MIN_HOUR_COVERAGE,
    HOUR_ALARM_SAMPLE_RATIO,
    OIL_CONFIRM_HOURS,
    OIL_STAGE_MERGE_GAP_HOURS,
    RESET_GAP_HOURS,
)
from models.common.stats import consecutive_decline_hours, linear_slope, median_value


def assign_hourly_cycles(rows: list[dict]):
    """Create oil-pressure cycles directly from valid running hours."""
    grouped = defaultdict(list)
    for row in rows:
        if row["运行覆盖率"] >= MIN_HOUR_COVERAGE and row["oil_avg"] is not None:
            grouped[row["批次"]].append(row)
    for batch, series in grouped.items():
        series.sort(key=lambda row: row["_dt"])
        cycle_number = 0
        previous = None
        for row in series:
            gap_hours = None if previous is None else (row["_dt"] - previous["_dt"]).total_seconds() / 3600
            # A one-hour pressure recovery is not a restart.  It can be a
            # sensor dropout or part of the same fault episode.  Start a new
            # cycle only after a long gap in valid running hours.
            if previous is None or (gap_hours is not None and gap_hours > RESET_GAP_HOURS):
                cycle_number += 1
            row["_cycle_id"] = f"{row['机组']}-P{cycle_number:02d}"
            previous = row
    return grouped


def classify_oil_hourly_segments(rows: list[dict], configs: dict[int, dict]):
    """Classify pressure by amplitude zone first, then confirm by duration.

    The configured low/low-low limits define the engineering zones.  A short
    excursion is retained as mild; only a run of three or more consecutive
    valid hours in the low zone is confirmed abnormal, and the low-low zone is
    confirmed severe.  No fixed 5/10/20 percent thresholds are used.
    """
    cycles = assign_hourly_cycles(rows)
    output = []
    for batch, batch_series in cycles.items():
        config = configs[batch]
        by_cycle = defaultdict(list)
        for row in batch_series:
            by_cycle[row["_cycle_id"]].append(row)
        for cycle_id, series in sorted(by_cycle.items(), key=lambda item: item[1][0]["_dt"]):
            series.sort(key=lambda item: item["_dt"])
            baseline_values = [row["oil_avg"] for row in series[:HOURLY_BASELINE_HOURS]]
            cycle_base = median_value(baseline_values)
            low = config["oil"]["low"]
            low_low = config["oil"]["low_low"]
            candidates = []
            history = []
            for row in series:
                current = row["oil_avg"]
                trend = linear_slope([item["oil_avg"] for item in history[-5:]] + [current])
                decline = consecutive_decline_hours(history, row)
                deviation = (cycle_base - current) / cycle_base if cycle_base else 0.0
                low_ratio = (row["oil_low_count"] or 0) / max(row["运行样本数"], 1)
                low_low_ratio = (row["oil_low_low_count"] or 0) / max(row["运行样本数"], 1)
                low_zone = current <= low
                low_low_zone = current <= low_low
                downward = (
                    decline >= 3
                    or (trend is not None and trend < -0.0005)
                )
                # Relative decline is a data-derived evolution signal.  The
                # engineering alarm limits remain higher priority, while the
                # 5/10/20% bands describe progression before the alarm limit.
                severe_candidate = low_low_zone or (deviation >= 0.20 and downward)
                abnormal_candidate = low_zone or (deviation >= 0.10 and downward)
                mild_candidate = deviation >= 0.05 and downward
                if severe_candidate:
                    raw = "严重异常"
                elif abnormal_candidate:
                    raw = "异常"
                elif mild_candidate:
                    raw = "轻微"
                else:
                    raw = "正常"
                candidates.append({
                    "row": row, "raw": raw, "trend": trend, "decline": decline,
                    "deviation": deviation, "low_ratio": low_ratio, "low_low_ratio": low_low_ratio,
                    "low_zone": low_zone, "low_low_zone": low_low_zone,
                })
                history.append(row)

            index = 0
            while index < len(candidates):
                raw = candidates[index]["raw"]
                if raw == "正常":
                    final = "正常"
                    end = index
                else:
                    end = index
                    while end + 1 < len(candidates) and candidates[end + 1]["raw"] == raw:
                        end += 1
                    run_length = end - index + 1
                    # Alarm limits have highest priority for the hourly
                    # candidate and immediate alarm fields.  Long-term stage
                    # labels still require persistence: a one/two-hour
                    # low-low excursion is a short severe warning, not a
                    # confirmed long-term severe stage.
                    if run_length >= OIL_CONFIRM_HOURS:
                        final = raw
                    elif raw == "严重异常":
                        final = "正常"
                    else:
                        final = "轻微"
                for item in candidates[index:end + 1]:
                    row = item["row"]
                    level_reason = {
                        "轻微": "相对周期基线下降至少5%并持续下降，或普通异常持续不足3小时",
                        "异常": "相对周期基线下降至少10%并持续下降，或低报警区间连续达到3小时",
                        "严重异常": "相对周期基线下降至少20%并持续下降，或低低报警区间连续达到3小时",
                        "正常": "未形成长期确认等级" if raw != "正常" else "未达到压力异常条件",
                    }[final]
                    if item["low_low_zone"]:
                        level_reason += f"；低低报警样本比例{item['low_low_ratio']:.1%}"
                    elif item["low_zone"]:
                        level_reason += f"；低报警样本比例{item['low_ratio']:.1%}"
                    output.append({
                        "机组": config["unit"], "批次": batch, "周期编号": cycle_id,
                        "小时": row["小时"], "运行样本数": row["运行样本数"], "运行覆盖率": row["运行覆盖率"],
                        "小时油压均值": row["oil_avg"], "小时油压最小值": row["oil_min"], "小时油压最大值": row["oil_max"],
                        "周期基线油压": cycle_base, "近期24有效小时基线油压": "",
                        "相对周期基线变化": (row["oil_avg"] / cycle_base - 1) if cycle_base else None,
                        "相对近期基线变化": "", "近6有效小时油压斜率": item["trend"],
                        "连续下降有效小时数": item["decline"], "历史有效运行小时数": index,
                        "低报警样本比例": item["low_ratio"], "低低报警样本比例": item["low_low_ratio"],
                        "低报警阈值": low, "低低报警阈值": low_low,
                        "是否低报警小时": "是" if item["low_zone"] else "否",
                        "是否低低报警小时": "是" if item["low_low_zone"] else "否",
                        "低低报警连续小时数": 0,
                        "压力等级": final, "原始小时等级": raw,
                        "立即报警": "严重预警" if item["low_low_zone"] else "低压预警" if item["low_zone"] else "",
                        "正常波动标记": "是" if raw != "正常" and final == "轻微" else "否",
                        "是否形成确认等级": "是" if final in ("异常", "严重异常") else "否",
                        "小时判定依据": f"油压{row['oil_avg']:.3f}MPa，周期基线{cycle_base:.3f}MPa，"
                        f"相对周期基线变化{((row['oil_avg'] / cycle_base - 1) if cycle_base else 0):+.1%}，"
                        f"{level_reason}",
                        "_dt": row["_dt"],
                    })
                index = end + 1
    return sorted(output, key=lambda row: (row["批次"], row["_dt"]))


def build_oil_hour_stage_segments(hourly_levels: list[dict]):
    """Merge confirmed levels and short unconfirmed triggers into segments."""
    grouped = defaultdict(list)
    for row in hourly_levels:
        if row["压力等级"] in ("轻微", "异常", "严重异常"):
            level = row["压力等级"]
        elif row["压力等级"] == "轻微" and row.get("是否形成确认等级") == "否":
            level = "轻微"
        else:
            continue
        grouped[(row["批次"], row["周期编号"], level)].append(row)
    segments = []
    for (batch, cycle_id, level), series in grouped.items():
        series.sort(key=lambda row: row["_dt"])
        current = []
        for row in series:
            if not current or (
                row["_dt"] - current[-1]["_dt"] <= timedelta(hours=OIL_STAGE_MERGE_GAP_HOURS)
                and row["压力等级"] == current[-1]["压力等级"]
            ):
                current.append(row)
            else:
                segments.append(_build_oil_hour_stage_segment(current, level))
                current = [row]
        if current:
            segments.append(_build_oil_hour_stage_segment(current, level))
    return sorted(segments, key=lambda row: (row["批次"], row["开始小时"]))


def _build_oil_hour_stage_segment(series: list[dict], level: str):
    first = series[0]
    last = series[-1]
    first_reason = series[0]["小时判定依据"]
    return {
        "机组": first["机组"],
        "批次": first["批次"],
        "周期编号": first["周期编号"],
        "压力阶段": level,
        "开始小时": first["小时"],
        "结束小时": last["小时"],
        "持续自然小时数": int((last["_dt"] - first["_dt"]).total_seconds() / 3600) + 1,
        "有效运行小时数": len(series),
        "周期基线油压": first["周期基线油压"],
        "正常波动带": first.get("正常波动带", ""),
        "阶段开始小时油压": first["小时油压均值"],
        "阶段结束小时油压": last["小时油压均值"],
        "阶段最低小时油压": min(row["小时油压均值"] for row in series),
        "阶段开始相对周期基线变化": first["相对周期基线变化"],
        "阶段结束相对周期基线变化": last["相对周期基线变化"],
        "阶段最大连续下降有效小时数": max(row["连续下降有效小时数"] for row in series),
        "阶段低报警小时数": sum(row["是否低报警小时"] == "是" for row in series),
        "阶段低低报警小时数": sum(row["是否低低报警小时"] == "是" for row in series),
        "阶段触发依据": (
            f"油压{first['小时油压均值']:.3f}～{last['小时油压均值']:.3f}MPa，"
            f"最低{min(row['小时油压均值'] for row in series):.3f}MPa，"
            f"连续{len(series)}个有效运行小时；{first_reason}"
        ),
        "是否形成确认等级": "是" if level != "轻微" or any(row.get("是否形成确认等级") == "是" for row in series) else "否",
    }


def summarize_hourly_oil_cycles(hourly_levels: list[dict], hourly_segments: list[dict]):
    by_cycle = defaultdict(list)
    for row in hourly_levels:
        if row["压力等级"] != "数据不足":
            by_cycle[row["周期编号"]].append(row)
    segments_by_cycle = defaultdict(list)
    for segment in hourly_segments:
        segments_by_cycle[segment["周期编号"]].append(segment)
    result = []
    for cycle_id, series in sorted(by_cycle.items(), key=lambda item: item[1][0]["_dt"]):
        series.sort(key=lambda row: row["_dt"])
        segments = segments_by_cycle[cycle_id]
        levels = ["正常波动", "轻微", "异常", "严重异常"]
        fields = {
            "机组": series[0]["机组"], "批次": series[0]["批次"], "周期编号": cycle_id,
            "周期开始小时": series[0]["小时"], "周期结束小时": series[-1]["小时"],
            "周期有效运行小时数": len(series), "周期基线油压": series[0]["周期基线油压"],
            "正常波动带": series[0].get("正常波动带", ""),
            "阶段演化顺序": "→".join(segment["压力阶段"] for segment in segments),
            "周期最低小时油压": min(row["小时油压均值"] for row in series),
            "是否达到轻微": "是" if any(row["压力等级"] == "轻微" for row in series) else "否",
            "是否达到异常": "是" if any(row["压力等级"] == "异常" for row in series) else "否",
            "是否达到严重异常": "是" if any(row["压力等级"] == "严重异常" for row in series) else "否",
            "周期短时低压小时数": sum(row["立即报警"] != "" and row["压力等级"] == "正常" for row in series),
            "正常波动段数": sum(segment["压力阶段"] == "正常波动" for segment in segments),
            "正常波动小时数": sum(segment["有效运行小时数"] for segment in segments if segment["压力阶段"] == "正常波动"),
        }
        for level in levels:
            selected = [segment for segment in segments if segment["压力阶段"] == level]
            fields[f"{level}阶段开始小时"] = selected[0]["开始小时"] if selected else ""
            fields[f"{level}阶段结束小时"] = selected[-1]["结束小时"] if selected else ""
            fields[f"{level}阶段运行小时数"] = sum(segment["有效运行小时数"] for segment in selected)
        result.append(fields)
    return result


HOURLY_LEVEL_FIELDS = [
    "机组", "批次", "周期编号", "小时", "运行样本数", "运行覆盖率", "小时油压均值",
    "小时油压最小值", "小时油压最大值", "周期基线油压", "近期24有效小时基线油压",
    "相对周期基线变化", "相对近期基线变化", "近6有效小时油压斜率", "连续下降有效小时数",
    "历史有效运行小时数", "低报警样本比例", "低低报警样本比例", "低报警阈值", "低低报警阈值",
    "是否低报警小时", "是否低低报警小时", "压力等级", "原始小时等级", "立即报警",
    "正常波动标记", "低低报警连续小时数", "是否纳入等级判定", "小时判定依据",
]

HOUR_STAGE_FIELDS = [
    "机组", "批次", "周期编号", "压力阶段", "开始小时", "结束小时", "持续自然小时数",
    "有效运行小时数", "周期基线油压", "正常波动带", "阶段开始小时油压", "阶段结束小时油压", "阶段最低小时油压",
    "阶段开始相对周期基线变化", "阶段结束相对周期基线变化", "阶段最大连续下降有效小时数",
    "阶段低报警小时数", "阶段低低报警小时数", "是否形成确认等级", "阶段触发依据",
]

CYCLE_FIELDS = [
    "机组", "批次", "周期编号", "周期开始小时", "周期结束小时", "周期有效运行小时数",
    "周期基线油压", "正常波动带", "阶段演化顺序", "周期最低小时油压", "是否达到轻微", "是否达到异常",
    "是否达到严重异常", "周期短时低压小时数", "正常波动段数", "正常波动小时数", "轻微阶段开始小时", "轻微阶段结束小时",
    "轻微阶段运行小时数", "异常阶段开始小时", "异常阶段结束小时", "异常阶段运行小时数",
    "严重异常阶段开始小时", "严重异常阶段结束小时", "严重异常阶段运行小时数",
]


def run(rows: list[dict], configs: dict[int, dict]):
    """Return ``(hourly_levels, hour_segments, cycles)`` for pressure."""
    hourly_levels = classify_oil_hourly_segments(rows, configs)
    hour_segments = build_oil_hour_stage_segments(hourly_levels)
    cycles = summarize_hourly_oil_cycles(hourly_levels, hour_segments)
    return hourly_levels, hour_segments, cycles


def write_results(output_dir: Path, hourly_levels, hour_segments, cycles):
    output_dir = Path(output_dir)
    io.write_csv(output_dir / "润滑油压力小时等级.csv", io.clean_rows(hourly_levels), HOURLY_LEVEL_FIELDS, {field: field for field in HOURLY_LEVEL_FIELDS})
    io.write_csv(output_dir / "润滑油压力小时阶段段.csv", io.clean_rows(hour_segments), HOUR_STAGE_FIELDS, {field: field for field in HOUR_STAGE_FIELDS})
    io.write_csv(output_dir / "润滑油压力演化周期.csv", io.clean_rows(cycles), CYCLE_FIELDS, {field: field for field in CYCLE_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.PRESSURE_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)
    hourly_levels, hour_segments, cycles = run(rows, configs)
    write_results(Path(args.output_dir), hourly_levels, hour_segments, cycles)
    print(f"设备：{devices_selected}")
    print(f"小时特征：{len(rows):,} 行")
    print(f"润滑油压力小时等级：{len(hourly_levels):,} 行")
    print(f"润滑油压力小时阶段段：{len(hour_segments):,} 段")
    print(f"润滑油压力周期：{len(cycles):,} 个")
    print(f"输出目录：{args.output_dir}")


if __name__ == "__main__":
    main()