Compressor/backend/indicator_features/indicator_features.py

#!/usr/bin/env python3
"""Phase-level features from one compressor indicator diagram cycle.

One stored cycle (wave_sample_one) is the pressure of a single cylinder chamber
over one crankshaft revolution (~60 ms at ~996 RPM).  Absolute pressure levels
depend on the process, so these features are meant to be compared:
  * against the same cylinder's own recent baseline (trend), and
  * against the peer cylinders of the same compression stage at the same time,
  * and between the cover/shaft chambers of the same cylinder.

Compression staging (confirmed with the site):
  一级(1st) = 2/4/6 缸,  二级(2nd) = 3/5 缸,  三级(3rd) = 1 缸.
No model is trained here; the goal is to expose interpretable features.
"""
from __future__ import annotations

import numpy as np

CYLINDER_BORE_MM = {"一": 360.0, "二": 490.0, "三": 390.0, "四": 490.0, "五": 390.0, "六": 490.0}
CYLINDER_STAGE = {"一": "三级", "二": "一级", "三": "二级", "四": "一级", "五": "二级", "六": "一级"}
PISTON_STROKE_MM = 148.0
CONNECTING_ROD_LENGTH_MM = 460.0
CLEARANCE_VOLUME_L_BY_BORE = {490.0: 1.51, 390.0: 0.74, 360.0: 0.62}


def smooth(values: np.ndarray, window: int) -> np.ndarray:
    window = max(1, int(window))
    if window <= 1 or len(values) < window:
        return values
    kernel = np.ones(window, dtype=float) / window
    return np.convolve(values, kernel, mode="same")


def volume_from_angle(count: int, bore_mm: float) -> np.ndarray:
    """Chamber volume (L) over one revolution, from crank-slider geometry."""
    crank_radius = PISTON_STROKE_MM / 2.0
    clearance = CLEARANCE_VOLUME_L_BY_BORE.get(bore_mm)
    if clearance is None:
        clearance = 1.0
    theta = np.linspace(0.0, 2.0 * np.pi, count, endpoint=False)
    travel = (
        crank_radius * (1.0 - np.cos(theta))
        + CONNECTING_ROD_LENGTH_MM
        - np.sqrt(CONNECTING_ROD_LENGTH_MM**2 - (crank_radius * np.sin(theta)) ** 2)
    )
    piston_area = np.pi * (bore_mm / 2.0) ** 2
    return clearance + piston_area * travel / 1_000_000.0


def _plateau_levels(p: np.ndarray) -> tuple[float, float]:
    """Suction (low) and discharge (high) plateau levels via robust quantiles."""
    low = np.quantile(p, 0.12)
    high = np.quantile(p, 0.88)
    near_low = p[p <= low + 0.15 * (high - low)]
    near_high = p[p >= high - 0.15 * (high - low)]
    suction = float(np.median(near_low)) if near_low.size else float(np.min(p))
    discharge = float(np.median(near_high)) if near_high.size else float(np.max(p))
    return suction, discharge


def cycle_features(signal: np.ndarray, cylinder: str, smooth_window: int = 21) -> dict:
    signal = np.asarray(signal, dtype=float)
    signal = signal[np.isfinite(signal)]
    count = len(signal)
    if count < 50:
        return {}
    p = smooth(signal, smooth_window)
    p = np.where(np.isfinite(p), p, signal)
    suction, discharge = _plateau_levels(p)
    span = discharge - suction
    threshold = suction + 0.5 * span
    high = p >= threshold
    dp = np.diff(p)
    rising = dp[dp > 0]
    falling = dp[dp < 0]
    # Compression/expansion line steepness (per sample).
    compression_slope = float(np.quantile(rising, 0.9)) if rising.size else 0.0
    expansion_slope = float(-np.quantile(falling, 0.1)) if falling.size else 0.0
    # Plateau occupancy: fraction of the cycle spent near high/low plateau.
    near = 0.08 * span + 1e-9
    suction_frac = float(np.mean(np.abs(p - suction) <= near))
    discharge_frac = float(np.mean(np.abs(p - discharge) <= near))
    # Indicated work over the revolution (net area) and mean effective pressure.
    volume = volume_from_angle(count, CYLINDER_BORE_MM.get(cylinder, 490.0))
    dvolume = np.diff(volume, append=volume[0])
    work = float(-np.sum(p * dvolume))  # net work done on the gas (positive)
    sweep = float(np.max(volume) - np.min(volume))
    imep = work / sweep if sweep else 0.0
    # Corner sharpness: mean |curvature| of pressure where it crosses the mid band.
    d2 = np.diff(p, 2)
    band = np.abs(p[1:-1] - threshold) <= 0.1 * span
    corner = float(np.mean(np.abs(d2[band]))) if band.any() else 0.0
    return {
        "n_samples": count,
        "p_min": float(np.min(p)),
        "p_max": float(np.max(p)),
        "p_mean": float(np.mean(p)),
        "p_suction": suction,
        "p_discharge": discharge,
        "p_span": span,
        "pressure_ratio": (discharge / suction) if suction else None,
        "pressure_rise": span,
        "suction_frac": suction_frac,
        "discharge_frac": discharge_frac,
        "compression_slope": compression_slope,
        "expansion_slope": expansion_slope,
        "indicated_work": work,
        "imep": imep,
        "corner_sharpness": corner,
    }