#!/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, OIL_PRECURSOR_DECLINE_HOURS, OIL_PRECURSOR_SLOPE, OIL_STARTUP_CANDIDATE_HOURS, OIL_STARTUP_OBSERVE_HOURS, OIL_STARTUP_RECOVERY_RATIO, OIL_STARTUP_STABLE_HOURS, RESET_GAP_HOURS, ) from models.common.stats import consecutive_decline_hours, linear_slope, median_value def _zero_ratio(row: dict) -> float | None: ratio = row.get("oil_zero_ratio") if ratio is not None: return float(ratio) zero_count = row.get("oil_zero_count") running = row.get("运行样本数") if zero_count is None or not running: return None return float(zero_count) / float(running) def _dropout_invalid_hours(rows: list[dict]) -> set[tuple[int, datetime]]: """Hours with >=5% zero oil samples next to a nearly all-zero hour.""" full_zero: set[tuple[int, datetime]] = set() for row in rows: ratio = _zero_ratio(row) if ratio is not None and ratio >= 0.5: full_zero.add((row["批次"], row["_dt"])) invalid: set[tuple[int, datetime]] = set() for row in rows: ratio = _zero_ratio(row) if ratio is None or ratio < 0.05: continue moment = row["_dt"] if ( (row["批次"], moment - timedelta(hours=1)) in full_zero or (row["批次"], moment + timedelta(hours=1)) in full_zero ): invalid.add((row["批次"], moment)) return invalid def _contiguous_slope(history: list[dict], current: dict, hours: int = 6): sequence = [current] newer = current for older in reversed(history): if newer["_dt"] - older["_dt"] != timedelta(hours=1): break if older["oil_avg"] is None: break sequence.append(older) newer = older if len(sequence) >= hours: break sequence.reverse() if len(sequence) < 3: return None return linear_slope([row["oil_avg"] for row in sequence]) def assign_hourly_cycles(rows: list[dict]): """Create oil-pressure cycles directly from valid running hours.""" invalid = _dropout_invalid_hours(rows) grouped = defaultdict(list) for row in rows: if ( row["运行覆盖率"] >= MIN_HOUR_COVERAGE and row["oil_avg"] is not None and row["oil_avg"] > 0 and (row["批次"], row["_dt"]) not in invalid ): 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 = _contiguous_slope(history, row) 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 extend_oil_segments_with_precursor(hourly_levels: list[dict], segments: list[dict]) -> list[dict]: """Move a mild stage start earlier when a confirmed slow decline precedes it.""" by_cycle = defaultdict(list) for row in hourly_levels: if row["压力等级"] != "数据不足": by_cycle[(row["批次"], row["周期编号"])].append(row) for rows in by_cycle.values(): rows.sort(key=lambda row: row["_dt"]) result = [] for segment in segments: if segment["压力阶段"] != "轻微": result.append(segment) continue rows = by_cycle.get((segment["批次"], segment["周期编号"]), []) start_index = next( (index for index, row in enumerate(rows) if row["小时"] == segment["开始小时"]), None, ) end_index = next( (index for index, row in enumerate(rows) if row["小时"] == segment["结束小时"]), None, ) if start_index is None or end_index is None or start_index == 0: result.append(segment) continue extension = [] anchor = rows[start_index] index = start_index - 1 while index >= 0: row = rows[index] deviation = row["相对周期基线变化"] decline = row["连续下降有效小时数"] slope = row["近6有效小时油压斜率"] if (anchor["_dt"] - row["_dt"]) != timedelta(hours=1): break if row["压力等级"] != "正常": break if deviation is None or not (0 <= -deviation < 0.05): break if decline is None or decline < OIL_PRECURSOR_DECLINE_HOURS: break if slope is None or slope >= OIL_PRECURSOR_SLOPE: break extension.insert(0, row) anchor = row index -= 1 if not extension: result.append(segment) continue base_series = [ row for row in rows[start_index:end_index + 1] if row["压力等级"] == segment["压力阶段"] ] extended = _build_oil_hour_stage_segment(extension + base_series, segment["压力阶段"]) first = extension[0] extended["阶段触发依据"] = ( f"前置下降：连续下降{first['连续下降有效小时数']}个有效小时，" f"6小时斜率{first['近6有效小时油压斜率']:.5f}，相对周期基线下降不足5%；" + extended["阶段触发依据"] ) result.append(extended) return sorted(result, key=lambda row: (row["批次"], row["开始小时"])) def filter_startup_recovery_segments(hourly_levels: list[dict], segments: list[dict]): """Remove early-cycle labels only when pressure later stabilises near baseline. A fixed post-startup exclusion window is deliberately not used: the label is removed only if at least ``OIL_STARTUP_STABLE_HOURS`` consecutive valid hours recover to ``OIL_STARTUP_RECOVERY_RATIO`` of the cycle baseline after the candidate segment. The transient hours are also returned so training can exclude them from the normal class. """ by_cycle = defaultdict(list) for row in hourly_levels: if row["压力等级"] != "数据不足": by_cycle[(row["批次"], row["周期编号"])].append(row) for rows in by_cycle.values(): rows.sort(key=lambda row: row["_dt"]) removed = [] excluded_hours: set[tuple[int, datetime]] = set() kept = [] for segment in segments: key = (segment["批次"], segment["周期编号"]) rows = by_cycle.get(key, []) if not rows: kept.append(segment) continue cycle_start = rows[0]["_dt"] start_dt = queries.parse_dt(segment["开始小时"]) if (start_dt - cycle_start).total_seconds() / 3600 > OIL_STARTUP_CANDIDATE_HOURS: kept.append(segment) continue end_index = next( (index for index, row in enumerate(rows) if row["小时"] == segment["结束小时"]), None, ) if end_index is None: kept.append(segment) continue baseline = segment.get("周期基线油压") baseline = float(baseline) if baseline not in (None, "") else None if not baseline or baseline <= 0: kept.append(segment) continue stable_start = None streak = 0 stable_run_start = None previous_dt = rows[end_index]["_dt"] for row in rows[end_index + 1:]: if (row["_dt"] - cycle_start).total_seconds() / 3600 > OIL_STARTUP_OBSERVE_HOURS: break if row["_dt"] - previous_dt > timedelta(hours=1): streak = 0 stable_run_start = None previous_dt = row["_dt"] pressure = row.get("小时油压均值") is_stable = ( row.get("压力等级") == "正常" and pressure is not None and pressure / baseline >= OIL_STARTUP_RECOVERY_RATIO and row.get("是否低报警小时") != "是" ) if is_stable: if streak == 0: stable_run_start = row["_dt"] streak += 1 if streak >= OIL_STARTUP_STABLE_HOURS: stable_start = stable_run_start break else: streak = 0 stable_run_start = None if stable_start is None: kept.append(segment) continue removed.append({ **segment, "稳定恢复开始小时": stable_start.strftime("%Y-%m-%d %H:%M:%S"), "连续稳定小时数": streak, "开机恢复判定": ( f"开机后{OIL_STARTUP_CANDIDATE_HOURS}小时内出现；" f"段后连续{streak}个有效小时恢复至周期基线" f"{OIL_STARTUP_RECOVERY_RATIO:.0%}以上且无低报警；" "判定为开机恢复瞬态，不纳入故障标签" ), }) for row in rows: if row["_dt"] >= stable_start: break excluded_hours.add((row["批次"], row["_dt"])) return kept, removed, excluded_hours 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_startup=False): """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) hour_segments = extend_oil_segments_with_precursor(hourly_levels, hour_segments) hour_segments, removed, excluded = filter_startup_recovery_segments(hourly_levels, hour_segments) cycles = summarize_hourly_oil_cycles(hourly_levels, hour_segments) if return_startup: return hourly_levels, hour_segments, cycles, removed, excluded return hourly_levels, hour_segments, cycles def write_results(output_dir: Path, hourly_levels, hour_segments, cycles, extra_fields=()): output_dir = Path(output_dir) stage_fields = [field for field in HOUR_STAGE_FIELDS if field not in extra_fields] + list(extra_fields) 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), stage_fields, {field: field for field in stage_fields}) io.write_csv(output_dir / "润滑油压力演化周期.csv", io.clean_rows(cycles), CYCLE_FIELDS, {field: field for field in CYCLE_FIELDS}) def write_startup_exclusions(output_dir: Path, hourly_levels, removed, excluded): output_dir = Path(output_dir) unit_by_batch = {row["批次"]: row["机组"] for row in hourly_levels} excluded_rows = [ { "机组": unit_by_batch.get(batch, ""), "批次": batch, "小时": moment.strftime("%Y-%m-%d %H:%M:%S"), "排除原因": "开机恢复瞬态", } for batch, moment in sorted(excluded, key=lambda item: (item[0], item[1])) ] io.write_csv( output_dir / "压力开机恢复剔除小时.csv", excluded_rows, ["机组", "批次", "小时", "排除原因"], {field: field for field in ("机组", "批次", "小时", "排除原因")}, ) removed_fields = [ "机组", "批次", "周期编号", "压力阶段", "开始小时", "结束小时", "周期基线油压", "稳定恢复开始小时", "连续稳定小时数", "开机恢复判定", ] io.write_csv( output_dir / "压力开机恢复剔除标签.csv", removed, removed_fields, {field: field for field in removed_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, removed, excluded = run(rows, configs, return_startup=True) output_dir = Path(args.output_dir) write_results(output_dir, hourly_levels, hour_segments, cycles) write_startup_exclusions(output_dir, hourly_levels, removed, excluded) print(f"设备：{devices_selected}") print(f"小时特征：{len(rows):,} 行") print(f"润滑油压力小时等级：{len(hourly_levels):,} 行") print(f"润滑油压力小时阶段段：{len(hour_segments):,} 段") print(f"开机恢复瞬态剔除：{len(removed)} 段 / {len(excluded)} 小时") print(f"润滑油压力周期：{len(cycles):,} 个") print(f"输出目录：{args.output_dir}") if __name__ == "__main__": main()