|
@@ -2,15 +2,17 @@
|
|
|
"""Phase-level features from one compressor indicator diagram cycle.
|
|
"""Phase-level features from one compressor indicator diagram cycle.
|
|
|
|
|
|
|
|
One stored cycle (wave_sample_one) is the pressure of a single cylinder chamber
|
|
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.
|
|
|
|
|
|
|
+over one crankshaft revolution (~60 ms at ~996 RPM).
|
|
|
|
|
+
|
|
|
|
|
+Design goals:
|
|
|
|
|
+ * robust two-plateau (suction P_s / discharge P_d) detection via Otsu split;
|
|
|
|
|
+ * physical phase alignment (search the crank offset that maximises the PV loop
|
|
|
|
|
+ area) so indicated work and a pressure-conversion efficiency are meaningful;
|
|
|
|
|
+ * a normalised phase profile for peer / double-chamber / self-baseline review.
|
|
|
|
|
|
|
|
Compression staging (confirmed with the site):
|
|
Compression staging (confirmed with the site):
|
|
|
一级(1st) = 2/4/6 缸, 二级(2nd) = 3/5 缸, 三级(3rd) = 1 缸.
|
|
一级(1st) = 2/4/6 缸, 二级(2nd) = 3/5 缸, 三级(3rd) = 1 缸.
|
|
|
-No model is trained here; the goal is to expose interpretable features.
|
|
|
|
|
|
|
+Absolute levels depend on the process, so features are meant to be compared.
|
|
|
"""
|
|
"""
|
|
|
from __future__ import annotations
|
|
from __future__ import annotations
|
|
|
|
|
|
|
@@ -21,6 +23,8 @@ CYLINDER_STAGE = {"一": "三级", "二": "一级", "三": "二级", "四": "一
|
|
|
PISTON_STROKE_MM = 148.0
|
|
PISTON_STROKE_MM = 148.0
|
|
|
CONNECTING_ROD_LENGTH_MM = 460.0
|
|
CONNECTING_ROD_LENGTH_MM = 460.0
|
|
|
CLEARANCE_VOLUME_L_BY_BORE = {490.0: 1.51, 390.0: 0.74, 360.0: 0.62}
|
|
CLEARANCE_VOLUME_L_BY_BORE = {490.0: 1.51, 390.0: 0.74, 360.0: 0.62}
|
|
|
|
|
+PROFILE_BINS = 12
|
|
|
|
|
+K_ISENTROPIC = 1.30
|
|
|
|
|
|
|
|
|
|
|
|
|
def smooth(values: np.ndarray, window: int) -> np.ndarray:
|
|
def smooth(values: np.ndarray, window: int) -> np.ndarray:
|
|
@@ -31,12 +35,37 @@ def smooth(values: np.ndarray, window: int) -> np.ndarray:
|
|
|
return np.convolve(values, kernel, mode="same")
|
|
return np.convolve(values, kernel, mode="same")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
+def otsu_threshold(values: np.ndarray, bins: int = 96) -> float:
|
|
|
|
|
+ hist, edges = np.histogram(values, bins=bins)
|
|
|
|
|
+ centers = (edges[:-1] + edges[1:]) / 2.0
|
|
|
|
|
+ total = float(values.size)
|
|
|
|
|
+ if total <= 0:
|
|
|
|
|
+ return float(np.mean(values))
|
|
|
|
|
+ sum_total = float(np.sum(centers * hist))
|
|
|
|
|
+ sum_before = 0.0
|
|
|
|
|
+ weight_before = 0.0
|
|
|
|
|
+ best_variance = -1.0
|
|
|
|
|
+ threshold = float(centers[0])
|
|
|
|
|
+ for index in range(len(hist)):
|
|
|
|
|
+ weight_before += float(hist[index])
|
|
|
|
|
+ if weight_before == 0.0:
|
|
|
|
|
+ continue
|
|
|
|
|
+ weight_after = total - weight_before
|
|
|
|
|
+ if weight_after == 0.0:
|
|
|
|
|
+ break
|
|
|
|
|
+ sum_before += centers[index] * hist[index]
|
|
|
|
|
+ mean_before = sum_before / weight_before
|
|
|
|
|
+ mean_after = (sum_total - sum_before) / weight_after
|
|
|
|
|
+ variance = weight_before * weight_after * (mean_before - mean_after) ** 2
|
|
|
|
|
+ if variance > best_variance:
|
|
|
|
|
+ best_variance = variance
|
|
|
|
|
+ threshold = float(centers[index])
|
|
|
|
|
+ return threshold
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
def volume_from_angle(count: int, bore_mm: float) -> np.ndarray:
|
|
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
|
|
crank_radius = PISTON_STROKE_MM / 2.0
|
|
|
- clearance = CLEARANCE_VOLUME_L_BY_BORE.get(bore_mm)
|
|
|
|
|
- if clearance is None:
|
|
|
|
|
- clearance = 1.0
|
|
|
|
|
|
|
+ clearance = CLEARANCE_VOLUME_L_BY_BORE.get(bore_mm, 1.0)
|
|
|
theta = np.linspace(0.0, 2.0 * np.pi, count, endpoint=False)
|
|
theta = np.linspace(0.0, 2.0 * np.pi, count, endpoint=False)
|
|
|
travel = (
|
|
travel = (
|
|
|
crank_radius * (1.0 - np.cos(theta))
|
|
crank_radius * (1.0 - np.cos(theta))
|
|
@@ -47,50 +76,85 @@ def volume_from_angle(count: int, bore_mm: float) -> np.ndarray:
|
|
|
return clearance + piston_area * travel / 1_000_000.0
|
|
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 _plateau_levels(p: np.ndarray) -> tuple[float, float, float]:
|
|
|
|
|
+ threshold = otsu_threshold(p)
|
|
|
|
|
+ low = p[p <= threshold]
|
|
|
|
|
+ high = p[p > threshold]
|
|
|
|
|
+ suction = float(np.median(low)) if low.size else float(np.min(p))
|
|
|
|
|
+ discharge = float(np.median(high)) if high.size else float(np.max(p))
|
|
|
|
|
+ return suction, discharge, threshold
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+def _align_offset(p: np.ndarray, dvolume: np.ndarray) -> tuple[int, float]:
|
|
|
|
|
+ """Offset (in samples) that maximises the indicated work -sum(p dV)."""
|
|
|
|
|
+ count = len(p)
|
|
|
|
|
+ cross = np.real(np.fft.ifft(np.conj(np.fft.fft(p)) * np.fft.fft(dvolume)))
|
|
|
|
|
+ work = -cross
|
|
|
|
|
+ offset = int(np.argmax(work))
|
|
|
|
|
+ return offset, float(work[offset])
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+def _polytropic_index(volume: np.ndarray, pressure: np.ndarray, low: float, high: float, rising: bool) -> float | None:
|
|
|
|
|
+ dp = np.diff(pressure, prepend=pressure[0])
|
|
|
|
|
+ band = (pressure > low + 0.15 * (high - low)) & (pressure < low + 0.85 * (high - low))
|
|
|
|
|
+ direction = dp > 0 if rising else dp < 0
|
|
|
|
|
+ mask = band & direction
|
|
|
|
|
+ v = volume[mask]
|
|
|
|
|
+ p = pressure[mask]
|
|
|
|
|
+ valid = (v > 0) & (p > 0)
|
|
|
|
|
+ if np.count_nonzero(valid) < 8:
|
|
|
|
|
+ return None
|
|
|
|
|
+ slope, _ = np.polyfit(np.log(v[valid]), np.log(p[valid]), 1)
|
|
|
|
|
+ return float(-slope)
|
|
|
|
|
|
|
|
|
|
|
|
|
def cycle_features(signal: np.ndarray, cylinder: str, smooth_window: int = 21) -> dict:
|
|
def cycle_features(signal: np.ndarray, cylinder: str, smooth_window: int = 21) -> dict:
|
|
|
signal = np.asarray(signal, dtype=float)
|
|
signal = np.asarray(signal, dtype=float)
|
|
|
signal = signal[np.isfinite(signal)]
|
|
signal = signal[np.isfinite(signal)]
|
|
|
count = len(signal)
|
|
count = len(signal)
|
|
|
- if count < 50:
|
|
|
|
|
|
|
+ if count < 60:
|
|
|
return {}
|
|
return {}
|
|
|
p = smooth(signal, smooth_window)
|
|
p = smooth(signal, smooth_window)
|
|
|
p = np.where(np.isfinite(p), p, signal)
|
|
p = np.where(np.isfinite(p), p, signal)
|
|
|
- suction, discharge = _plateau_levels(p)
|
|
|
|
|
|
|
+ suction, discharge, threshold = _plateau_levels(p)
|
|
|
span = discharge - suction
|
|
span = discharge - suction
|
|
|
- threshold = suction + 0.5 * span
|
|
|
|
|
- high = p >= threshold
|
|
|
|
|
|
|
+ volume = volume_from_angle(count, CYLINDER_BORE_MM.get(cylinder, 490.0))
|
|
|
|
|
+ dvolume = np.diff(volume, append=volume[0])
|
|
|
|
|
+
|
|
|
|
|
+ offset, indicated_work = _align_offset(p, dvolume)
|
|
|
|
|
+ aligned_volume = np.roll(volume, -offset)
|
|
|
|
|
+ aligned_dvolume = np.roll(dvolume, -offset)
|
|
|
|
|
+
|
|
|
|
|
+ # Isentropic ideal work (MPa*L = kJ) and pressure-conversion efficiency.
|
|
|
|
|
+ v_suction = float(np.max(aligned_volume))
|
|
|
|
|
+ ratio = (discharge / suction) if suction > 0 else None
|
|
|
|
|
+ if ratio and ratio > 1:
|
|
|
|
|
+ ideal = (K_ISENTROPIC / (K_ISENTROPIC - 1.0)) * suction * v_suction * (
|
|
|
|
|
+ ratio ** ((K_ISENTROPIC - 1.0) / K_ISENTROPIC) - 1.0
|
|
|
|
|
+ )
|
|
|
|
|
+ else:
|
|
|
|
|
+ ideal = None
|
|
|
|
|
+ efficiency = (indicated_work / ideal) if ideal and ideal > 1e-9 else None
|
|
|
|
|
+
|
|
|
|
|
+ # Phase boundaries by pressure crossing (aligned).
|
|
|
|
|
+ high_mask = p >= threshold
|
|
|
|
|
+ suction_frac = float(np.mean(~high_mask))
|
|
|
|
|
+ discharge_frac = float(np.mean(high_mask))
|
|
|
dp = np.diff(p)
|
|
dp = np.diff(p)
|
|
|
rising = dp[dp > 0]
|
|
rising = dp[dp > 0]
|
|
|
falling = 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
|
|
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
|
|
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)
|
|
d2 = np.diff(p, 2)
|
|
|
band = np.abs(p[1:-1] - threshold) <= 0.1 * span
|
|
band = np.abs(p[1:-1] - threshold) <= 0.1 * span
|
|
|
corner = float(np.mean(np.abs(d2[band]))) if band.any() else 0.0
|
|
corner = float(np.mean(np.abs(d2[band]))) if band.any() else 0.0
|
|
|
- return {
|
|
|
|
|
|
|
+
|
|
|
|
|
+ # Normalised phase profile (0..1) on the aligned cycle.
|
|
|
|
|
+ idx = np.linspace(0, count - 1, PROFILE_BINS).round().astype(int)
|
|
|
|
|
+ prof = (p[idx] - suction) / span if span else np.zeros(PROFILE_BINS)
|
|
|
|
|
+ profile = {f"prof_{i:02d}": float(prof[i]) for i in range(PROFILE_BINS)}
|
|
|
|
|
+
|
|
|
|
|
+ features = {
|
|
|
"n_samples": count,
|
|
"n_samples": count,
|
|
|
"p_min": float(np.min(p)),
|
|
"p_min": float(np.min(p)),
|
|
|
"p_max": float(np.max(p)),
|
|
"p_max": float(np.max(p)),
|
|
@@ -98,13 +162,18 @@ def cycle_features(signal: np.ndarray, cylinder: str, smooth_window: int = 21) -
|
|
|
"p_suction": suction,
|
|
"p_suction": suction,
|
|
|
"p_discharge": discharge,
|
|
"p_discharge": discharge,
|
|
|
"p_span": span,
|
|
"p_span": span,
|
|
|
- "pressure_ratio": (discharge / suction) if suction else None,
|
|
|
|
|
- "pressure_rise": span,
|
|
|
|
|
|
|
+ "pressure_ratio": ratio,
|
|
|
"suction_frac": suction_frac,
|
|
"suction_frac": suction_frac,
|
|
|
"discharge_frac": discharge_frac,
|
|
"discharge_frac": discharge_frac,
|
|
|
"compression_slope": compression_slope,
|
|
"compression_slope": compression_slope,
|
|
|
"expansion_slope": expansion_slope,
|
|
"expansion_slope": expansion_slope,
|
|
|
- "indicated_work": work,
|
|
|
|
|
- "imep": imep,
|
|
|
|
|
"corner_sharpness": corner,
|
|
"corner_sharpness": corner,
|
|
|
|
|
+ "n_comp": _polytropic_index(aligned_volume, p, suction, discharge, rising=True),
|
|
|
|
|
+ "n_exp": _polytropic_index(aligned_volume, p, suction, discharge, rising=False),
|
|
|
|
|
+ "indicated_work": indicated_work,
|
|
|
|
|
+ "ideal_work_isentropic": ideal,
|
|
|
|
|
+ "efficiency": efficiency,
|
|
|
|
|
+ "phase_offset": offset,
|
|
|
}
|
|
}
|
|
|
|
|
+ features.update(profile)
|
|
|
|
|
+ return features
|