From ad9aeeb4e9dae4685fccc92ad3b033d33c54b9c3 Mon Sep 17 00:00:00 2001 From: dmitei Date: Tue, 10 Nov 2026 23:02:57 +0300 Subject: [PATCH] add physics degradation model with open-meteo and residual Co-Authored-By: Claude Opus 5.5 --- .env.example | 2 + .github/workflows/ci.yml | 10 +- README.md | 3 + backend/app/api/equipment.py | 15 +- backend/app/schemas/equipment.py | 7 + backend/app/services/equipment_service.py | 21 +++ backend/pytest.ini | 3 + backend/tests/conftest.py | 76 ++++++++++ backend/tests/test_physics.py | 61 ++++++++ frontend/src/api.ts | 8 ++ frontend/src/format.ts | 9 ++ frontend/src/index.css | 11 -- frontend/src/pages/EquipmentCard.tsx | 30 +++- ml/physics.py | 166 ++++++++++++++++++++++ ml/physics_check.py | 77 ++++++++++ ml/worker.py | 40 +++++- 16 files changed, 517 insertions(+), 22 deletions(-) create mode 100644 backend/pytest.ini create mode 100644 backend/tests/conftest.py create mode 100644 backend/tests/test_physics.py create mode 100644 ml/physics_check.py diff --git a/.env.example b/.env.example index b78e297..ecde501 100644 --- a/.env.example +++ b/.env.example @@ -10,5 +10,7 @@ CORS_ORIGINS= RUL_WARNING=75 RUL_CRITICAL=30 RUL_CLIP=125 +WEATHER_LATITUDE=55.7558 +WEATHER_LONGITUDE=37.6173 BOT_TOKEN= CHAT_ID= diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 066a9bf..1527329 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -90,6 +90,10 @@ jobs: working-directory: ml run: python baseline_model_check.py ../data/test_FD001.txt --rul ../data/RUL_FD001.txt + - name: Run physics check + working-directory: ml + run: python physics_check.py ../data/train_FD001.txt + pipeline: needs: [migrations, model] runs-on: ubuntu-latest @@ -154,5 +158,7 @@ jobs: EQ=$(psql "$PSQL_URL" -tAc "SELECT count(*) FROM equipment") SR=$(psql "$PSQL_URL" -tAc "SELECT count(*) FROM sensor_readings") ML=$(psql "$PSQL_URL" -tAc "SELECT count(*) FROM ml_predictions") - echo "equipment=$EQ sensor_readings=$SR ml_predictions=$ML" - [ "$EQ" -eq 20 ] && [ "$SR" -gt 0 ] && [ "$ML" -eq "$SR" ] + PH=$(psql "$PSQL_URL" -tAc "SELECT count(*) FROM physics_predictions") + RS=$(psql "$PSQL_URL" -tAc "SELECT count(*) FROM ml_predictions WHERE residual <> 0") + echo "equipment=$EQ sensor_readings=$SR ml_predictions=$ML physics_predictions=$PH nonzero_residual=$RS" + [ "$EQ" -eq 20 ] && [ "$SR" -gt 0 ] && [ "$ML" -eq "$SR" ] && [ "$PH" -eq "$SR" ] && [ "$RS" -gt 0 ] diff --git a/README.md b/README.md index de2d9f4..2661cab 100644 --- a/README.md +++ b/README.md @@ -50,6 +50,8 @@ Frontend: React **Статусы оборудования по прогнозу RUL:** `normal` — более 75 циклов, `warning` — от 31 до 75, `critical` — 30 циклов и меньше. +**Физическая модель** (`ml/physics.py`) считается воркером вместе с RUL. Фактическая деградация — проекция нормированных показаний 14 датчиков на первую главную компоненту, отмасштабированная так, что 0 — исправный двигатель, 1 — отказ (по двигателям 1–80), со сглаживанием EWMA. Ожидаемая деградация — экспоненциальный закон накопления повреждений D(τ) = a·(e^(b·τ) − 1), которым NASA моделировала износ в C-MAPSS (Saxena et al., 2008). Эквивалентная наработка τ за цикл растёт на коэффициент ускорения по Аррениусу от температуры воздуха из Open-Meteo: при 15 °C он равен 1, при 0 °C — 0,83, при 30 °C — 1,18. Координаты площадки — `WEATHER_LATITUDE` и `WEATHER_LONGITUDE` (по умолчанию Москва). Ожидаемая деградация пишется в `physics_predictions`, невязка «фактическая минус ожидаемая» — в `ml_predictions.residual`. Проверка модели: `python physics_check.py ../data/train_FD001.txt` из `ml/`. + ## Масштаб системы Целевая аудитория системы - инженеры по эксплуатации и ответственные за производство на промышленных предприятиях. При оценке масштаба в 10 000 уникальных пользователей в сутки предполагается, что система обслуживает порядка N предприятий с M единицами оборудования на каждом, где каждый пользователь в среднем совершает K запросов в день (просмотр дашборда, проверка истории, получение алертов). Основная нагрузка - read-запросы (просмотр данных), запись происходит асинхронно от симулятора/датчиков оборудования, а не от действий пользователей напрямую. @@ -125,6 +127,7 @@ api/v1 (эндпоинты, валидация) → services (бизнес-ло | GET | `/equipment/{id}/status` | Текущее состояние: статус, RUL, вероятность дефекта, невязка | < 150 мс | | GET | `/equipment/{id}/history` | История показаний датчиков (`from_cycle`, `to_cycle`, `limit`) | < 500 мс | | GET | `/equipment/{id}/forecast` | Прогноз деградации: RUL по будущим циклам с доверительным интервалом | < 500 мс | +| GET | `/equipment/{id}/physics` | Физическая модель против фактических показаний: ожидаемая и фактическая деградация, невязка (`limit`) | < 500 мс | ### Алерты diff --git a/backend/app/api/equipment.py b/backend/app/api/equipment.py index d82907e..f7bb676 100644 --- a/backend/app/api/equipment.py +++ b/backend/app/api/equipment.py @@ -2,7 +2,13 @@ from app.core.db import DB from app.core.errors import NOT_FOUND, equipment_not_found -from app.schemas.equipment import EquipmentOut, ForecastOut, ReadingOut, StatusOut +from app.schemas.equipment import ( + EquipmentOut, + ForecastOut, + PhysicsPointOut, + ReadingOut, + StatusOut, +) from app.services import equipment_service router = APIRouter(prefix="/equipment", tags=["Оборудование"]) @@ -42,6 +48,13 @@ def get_history(equipment_id: int, db: DB, return _found(equipment_service.get_history(db, equipment_id, from_cycle, to_cycle, limit), equipment_id) +@router.get("/{equipment_id}/physics", response_model=list[PhysicsPointOut], responses=NOT_FOUND, + summary="Физическая модель против фактических показаний: ожидаемая и фактическая деградация") +def get_physics(equipment_id: int, db: DB, + limit: int = Query(500, ge=1, le=5000, description="Максимум записей (берутся последние)")): + return _found(equipment_service.get_physics(db, equipment_id, limit), equipment_id) + + @router.get("/{equipment_id}/forecast", response_model=ForecastOut, responses=NOT_FOUND, summary="Прогноз деградации: RUL по будущим циклам с доверительным интервалом") def get_forecast(equipment_id: int, db: DB): diff --git a/backend/app/schemas/equipment.py b/backend/app/schemas/equipment.py index 634101d..7986fc4 100644 --- a/backend/app/schemas/equipment.py +++ b/backend/app/schemas/equipment.py @@ -31,6 +31,13 @@ class ReadingOut(BaseModel): sensors: dict[str, float] = Field(description="Показания 14 рабочих датчиков: код → значение") +class PhysicsPointOut(BaseModel): + cycle: int + expected: float = Field(description="Ожидаемая деградация по физической модели: 0 — исправный двигатель, 1 — отказ") + actual: float = Field(description="Фактическая деградация по показаниям датчиков") + residual: float = Field(description="Невязка: фактическая деградация минус ожидаемая") + + class ForecastPoint(BaseModel): cycle: int rul: float = Field(description="Ожидаемый остаточный ресурс на этом цикле") diff --git a/backend/app/services/equipment_service.py b/backend/app/services/equipment_service.py index 717142a..d501aac 100644 --- a/backend/app/services/equipment_service.py +++ b/backend/app/services/equipment_service.py @@ -3,6 +3,7 @@ from app.core.config import RUL_CLIP, RUL_CRITICAL, RUL_WARNING from app.models.equipment import Equipment from app.models.ml_predictions import MlPrediction +from app.models.physics_predictions import PhysicsPrediction from app.models.sensor_readings import SensorReading SENSOR_CODES = ["s2", "s3", "s4", "s7", "s8", "s9", "s11", "s12", "s13", "s14", "s15", "s17", "s20", "s21"] @@ -115,6 +116,26 @@ def get_history(db: Connection, equipment_id: int, from_cycle: int | None, to_cy return rows +def get_physics(db: Connection, equipment_id: int, limit: int) -> list[dict] | None: + if not equipment_exists(db, equipment_id): + return None + + query = ( + select(PhysicsPrediction.cycle, PhysicsPrediction.predicted_degradation, MlPrediction.residual) + .join(MlPrediction, and_(MlPrediction.equipment_id == PhysicsPrediction.equipment_id, + MlPrediction.cycle == PhysicsPrediction.cycle)) + .where(PhysicsPrediction.equipment_id == equipment_id) + .order_by(PhysicsPrediction.cycle.desc()) + .limit(limit) + ) + + points = [{"cycle": row.cycle, "expected": row.predicted_degradation, + "actual": row.predicted_degradation + row.residual, "residual": row.residual} + for row in db.execute(query)] + points.reverse() + return points + + def get_forecast(db: Connection, equipment_id: int) -> dict | None: items = _latest_predictions(db, equipment_id) if not items: diff --git a/backend/pytest.ini b/backend/pytest.ini new file mode 100644 index 0000000..c7b23ec --- /dev/null +++ b/backend/pytest.ini @@ -0,0 +1,3 @@ +[pytest] +pythonpath = . +testpaths = tests diff --git a/backend/tests/conftest.py b/backend/tests/conftest.py new file mode 100644 index 0000000..13285c2 --- /dev/null +++ b/backend/tests/conftest.py @@ -0,0 +1,76 @@ +import os +from datetime import date + +import pytest +from fastapi.testclient import TestClient +from sqlalchemy import insert, text + +# Тесты очищают таблицы, поэтому БД для них задаётся отдельной переменной: +# рабочий DATABASE_URL здесь намеренно не используется. +os.environ["DATABASE_URL"] = os.environ["TEST_DATABASE_URL"] +os.environ["RUL_WARNING"] = "75" +os.environ["RUL_CRITICAL"] = "30" +os.environ["RUL_CLIP"] = "125" + +TABLES = "alerts, ml_predictions, physics_predictions, sensor_readings, equipment" + + +@pytest.fixture +def db(): + from app.core.db import engine + + with engine.begin() as connection: + connection.execute(text(f"TRUNCATE {TABLES} RESTART IDENTITY CASCADE")) + with engine.connect() as connection: + yield connection + + +@pytest.fixture +def client(db): + from main import app + + return TestClient(app) + + +@pytest.fixture +def add_equipment(db): + from app.models import Equipment + + def add(unit_number: int = 81) -> int: + statement = insert(Equipment).values( + name=f"Turbofan №{unit_number}", type="turbofan", dataset="FD001", + unit_number=unit_number, install_date=date(2024, 1, 1), + ) + equipment_id = db.execute(statement.returning(Equipment.id)).scalar_one() + db.commit() + return equipment_id + + return add + + +@pytest.fixture +def add_prediction(db): + from app.models import MlPrediction + + def add(equipment_id: int, cycle: int, rul: float) -> None: + db.execute(insert(MlPrediction).values( + equipment_id=equipment_id, cycle=cycle, predicted_rul=rul, + defect_probability=0, residual=0, model_version="test", + )) + db.commit() + + return add + + +@pytest.fixture +def add_reading(db): + from app.models import SensorReading + from app.services.equipment_service import SENSOR_CODES + + def add(equipment_id: int, cycle: int) -> dict[str, float]: + sensors = {code: float(i + cycle) for i, code in enumerate(SENSOR_CODES)} + db.execute(insert(SensorReading).values(equipment_id=equipment_id, cycle=cycle, **sensors)) + db.commit() + return sensors + + return add diff --git a/backend/tests/test_physics.py b/backend/tests/test_physics.py new file mode 100644 index 0000000..2238a82 --- /dev/null +++ b/backend/tests/test_physics.py @@ -0,0 +1,61 @@ +import pytest +from sqlalchemy import insert + + +@pytest.fixture +def add_physics(db): + from app.models import MlPrediction, PhysicsPrediction + + def add(equipment_id: int, cycle: int, expected: float, residual: float) -> None: + db.execute(insert(PhysicsPrediction).values( + equipment_id=equipment_id, cycle=cycle, predicted_degradation=expected, + )) + db.execute(insert(MlPrediction).values( + equipment_id=equipment_id, cycle=cycle, predicted_rul=100, + defect_probability=0, residual=residual, model_version="test", + )) + db.commit() + + return add + + +def test_physics_returns_expected_actual_and_residual(client, add_equipment, add_physics): + equipment_id = add_equipment() + add_physics(equipment_id, 2, expected=0.2, residual=-0.05) + add_physics(equipment_id, 1, expected=0.1, residual=0.05) + + response = client.get(f"/equipment/{equipment_id}/physics") + + assert response.status_code == 200 + points = response.json() + assert [p["cycle"] for p in points] == [1, 2] + assert points[0]["expected"] == pytest.approx(0.1) + assert points[0]["actual"] == pytest.approx(0.15) + assert points[0]["residual"] == pytest.approx(0.05) + assert points[1]["actual"] == pytest.approx(0.15) + + +def test_physics_limit_keeps_latest_cycles(client, add_equipment, add_physics): + equipment_id = add_equipment() + for cycle in (1, 2, 3): + add_physics(equipment_id, cycle, expected=cycle / 10, residual=0) + + response = client.get(f"/equipment/{equipment_id}/physics", params={"limit": 2}) + + assert [p["cycle"] for p in response.json()] == [2, 3] + + +def test_physics_is_empty_before_worker(client, add_equipment, add_reading): + equipment_id = add_equipment() + add_reading(equipment_id, 1) + + response = client.get(f"/equipment/{equipment_id}/physics") + + assert response.status_code == 200 + assert response.json() == [] + + +def test_physics_unknown_equipment(client, db): + response = client.get("/equipment/999/physics") + + assert response.status_code == 404 diff --git a/frontend/src/api.ts b/frontend/src/api.ts index 5819da1..b1cda21 100644 --- a/frontend/src/api.ts +++ b/frontend/src/api.ts @@ -25,6 +25,13 @@ export interface Reading { sensors: Record; } +export interface PhysicsPoint { + cycle: number; + expected: number; + actual: number; + residual: number; +} + export interface ForecastPoint { cycle: number; rul: number; @@ -69,6 +76,7 @@ export const api = { getEquipment: (id: number) => request(`/equipment/${id}`), getHistory: (id: number, limit = 500) => request(`/equipment/${id}/history?limit=${limit}`), getForecast: (id: number) => request(`/equipment/${id}/forecast`), + getPhysics: (id: number, limit = 500) => request(`/equipment/${id}/physics?limit=${limit}`), listAlerts: (limit = 100) => request(`/alerts?limit=${limit}`), alertsByEquipment: (id: number, limit = 50) => request(`/alerts/${id}?limit=${limit}`), predict: (id: number) => diff --git a/frontend/src/format.ts b/frontend/src/format.ts index 2d65f43..1ded1d9 100644 --- a/frontend/src/format.ts +++ b/frontend/src/format.ts @@ -41,3 +41,12 @@ export function formatDateTime(value: string): string { export function formatRul(value: number): string { return Math.round(value).toString(); } + +export function formatPercent(value: number): string { + return `${Math.round(value * 100)}%`; +} + +export function formatResidual(value: number): string { + const points = Math.round(value * 100); + return `${points > 0 ? "+" : ""}${points} п.п.`; +} diff --git a/frontend/src/index.css b/frontend/src/index.css index 9cca51c..843f5c4 100644 --- a/frontend/src/index.css +++ b/frontend/src/index.css @@ -205,17 +205,6 @@ h2 { .grid > .card:last-child { grid-column: 1 / -1; } -.placeholder { - height: 260px; - display: flex; - align-items: center; - justify-content: center; - text-align: center; - padding: 24px; - color: #6b7280; - border: 2px dashed #e5e7eb; - border-radius: 6px; -} .muted { color: #6b7280; diff --git a/frontend/src/pages/EquipmentCard.tsx b/frontend/src/pages/EquipmentCard.tsx index ef5eaa3..ee4e401 100644 --- a/frontend/src/pages/EquipmentCard.tsx +++ b/frontend/src/pages/EquipmentCard.tsx @@ -4,6 +4,7 @@ import { Area, CartesianGrid, ComposedChart, + Legend, Line, LineChart, ReferenceLine, @@ -18,7 +19,7 @@ import { AlertsTable } from "../components/AlertsTable"; import { LiveIndicator } from "../components/LiveIndicator"; import { Empty, ErrorMessage, Loading } from "../components/State"; import { StatusBadge } from "../components/StatusBadge"; -import { SENSORS, formatDate, formatRul } from "../format"; +import { SENSORS, formatDate, formatPercent, formatResidual, formatRul } from "../format"; import { useLiveFetch } from "../useLiveFetch"; const RUL_WARNING = 75; @@ -84,7 +85,7 @@ export function EquipmentCard() {
- +
@@ -158,15 +159,32 @@ function ForecastChart({ id }: { id: number }) { ); } -function PhysicsChart() { +function PhysicsChart({ id }: { id: number }) { + const { data, error } = useLiveFetch(() => api.getPhysics(id), [id]); + const last = data && data.length > 0 ? data[data.length - 1] : null; + return (

Физическая модель против фактических показаний

+ {last && невязка {formatResidual(last.residual)}}
-
- Заглушка -
+ {error && } + {!data && !error && } + {data && data.length === 0 && } + {data && data.length > 0 && ( + + + + + + [formatPercent(v), name]} labelFormatter={(c) => `Цикл ${c}`} /> + + + + + + )}
); } diff --git a/ml/physics.py b/ml/physics.py index e69de29..1cd4f6a 100644 --- a/ml/physics.py +++ b/ml/physics.py @@ -0,0 +1,166 @@ +"""Физическая модель деградации турбовентиляторного двигателя. + +Фактическая деградация по датчикам: проекция нормированных показаний 14 рабочих датчиков на главное +направление их изменения (первая главная компонента), отмасштабированная так, что 0 соответствует +исправному двигателю (первые 20 циклов), а 1 — отказу (последние 5 циклов) на обучающих двигателях. + +Ожидаемая деградация: экспоненциальный закон накопления повреждений D(τ) = a·(e^(b·τ) − 1), которым +NASA моделировала износ при генерации C-MAPSS (Saxena et al., 2008). τ — эквивалентная наработка: +каждый цикл добавляет к ней коэффициент ускорения по Аррениусу от температуры воздуха из Open-Meteo, +AF = exp(Ea/R · (1/T_ref − 1/T)). T_ref = 288,15 K (15 °C) — стандартная атмосфера, в которой сняты +данные C-MAPSS, поэтому при 15 °C AF = 1. + +Невязка (residual) — фактическая деградация минус ожидаемая: больше нуля, если двигатель изнашивается +быстрее, чем предсказывает физическая модель. +""" + +import json +import logging +import math +import time +import urllib.parse +import urllib.request +from collections.abc import Callable +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import pandas as pd +from prepare_rul import load_train +from scipy.optimize import curve_fit + +logger = logging.getLogger("physics") + +SENSORS = ["s2", "s3", "s4", "s7", "s8", "s9", "s11", "s12", "s13", "s14", "s15", "s17", "s20", "s21"] +TRAIN_UNITS = 80 # двигатели 81–100 проигрывает симулятор, в обучение они не попадают +HEALTHY_CYCLES = 20 +FAILED_CYCLES = 5 + +GAS_CONSTANT = 8.314 # Дж/(моль·К) +T_REF = 288.15 # К +ACTIVATION_ENERGY = 8000.0 # Дж/моль: AF ≈ 0,63 при −20 °C и ≈ 1,24 при +35 °C +SMOOTHING = 0.3 # коэффициент экспоненциального сглаживания фактической деградации + +OPEN_METEO_URL = "https://api.open-meteo.com/v1/forecast" +WEATHER_TTL_SEC = 3600 +WEATHER_RETRY_SEC = 300 + + +def damage(tau, a, b): + return a * (np.exp(b * tau) - 1.0) + + +@dataclass(frozen=True) +class PhysicsModel: + mean: np.ndarray + std: np.ndarray + direction: np.ndarray + healthy: float + failed: float + a: float + b: float + + def observed(self, sensors: np.ndarray) -> np.ndarray: + """Фактическая деградация по показаниям: строки — циклы, колонки — SENSORS.""" + score = ((sensors - self.mean) / self.std) @ self.direction + return (score - self.healthy) / (self.failed - self.healthy) + + def expected(self, tau: float) -> float: + return float(damage(tau, self.a, self.b)) + + def age(self, degradation: float) -> float: + """Эквивалентная наработка по ожидаемой деградации — функция, обратная expected.""" + return math.log(degradation / self.a + 1.0) / self.b + + +def fit(path: Path) -> PhysicsModel: + df = load_train(path) + df = df[df["unit"] <= TRAIN_UNITS] + + x = df[SENSORS].to_numpy(dtype=float) + mean, std = x.mean(axis=0), x.std(axis=0) + z = (x - mean) / std + _, _, vt = np.linalg.svd(z, full_matrices=False) + direction = vt[0] + + score = z @ direction + cycle = df["cycle"].to_numpy(dtype=float) + life = df.groupby("unit")["cycle"].transform("max").to_numpy() + healthy = score[cycle <= HEALTHY_CYCLES].mean() + failed = score[cycle > life - FAILED_CYCLES].mean() + observed = (score - healthy) / (failed - healthy) + + (a, b), _ = curve_fit(damage, cycle, observed, p0=(0.5, 0.004), maxfev=20000) + logger.info("physics model: D(t) = %.4f * (exp(%.5f * t) - 1)", a, b) + return PhysicsModel(mean, std, direction, float(healthy), float(failed), float(a), float(b)) + + +def acceleration(temperature: float | None, activation_energy: float = ACTIVATION_ENERGY) -> float: + """Коэффициент ускорения износа по Аррениусу; без погоды — 1.""" + if temperature is None: + return 1.0 + kelvin = temperature + 273.15 + return math.exp(activation_energy / GAS_CONSTANT * (1.0 / T_REF - 1.0 / kelvin)) + + +def predict(model: PhysicsModel, readings: pd.DataFrame, state: dict[int, tuple[int, float, float]], + factor: float, smoothing: float = SMOOTHING) -> pd.DataFrame: + """Ожидаемая деградация и невязка для новых показаний. + + readings — новые строки sensor_readings, отсортированные по (equipment_id, cycle). + state — последнее посчитанное состояние по оборудованию: цикл, ожидаемая деградация и сглаженная + фактическая. Функция дополняет его, чтобы следующая пачка продолжила накопление. + """ + observed = model.observed(readings[SENSORS].to_numpy(dtype=float)) + rows = [] + for equipment_id, cycle, raw in zip(readings["equipment_id"], readings["cycle"], observed): + equipment_id, cycle, raw = int(equipment_id), int(cycle), float(raw) + previous = state.get(equipment_id) + if previous is None: + tau, actual = cycle * factor, raw + else: + last_cycle, last_expected, last_actual = previous + tau = model.age(last_expected) + (cycle - last_cycle) * factor + actual = smoothing * raw + (1.0 - smoothing) * last_actual + expected = model.expected(tau) + state[equipment_id] = (cycle, expected, actual) + rows.append((equipment_id, cycle, expected, actual - expected)) + return pd.DataFrame(rows, columns=["equipment_id", "cycle", "predicted_degradation", "residual"]) + + +class Weather: + """Температура воздуха из Open-Meteo с кэшем на час. + + Если API недоступен, используется последнее полученное значение; если значения ещё нет, + модель работает без погодного фактора. + """ + + def __init__(self, latitude: float, longitude: float, fetch: Callable[[], float] | None = None, + ttl: float = WEATHER_TTL_SEC, retry: float = WEATHER_RETRY_SEC): + self.latitude = latitude + self.longitude = longitude + self.ttl = ttl + self.retry = retry + self._fetch = fetch or self._fetch_open_meteo + self._value: float | None = None + self._next_fetch = 0.0 + + def _fetch_open_meteo(self) -> float: + query = urllib.parse.urlencode({"latitude": self.latitude, "longitude": self.longitude, + "current": "temperature_2m"}) + with urllib.request.urlopen(f"{OPEN_METEO_URL}?{query}", timeout=5) as response: + return float(json.load(response)["current"]["temperature_2m"]) + + def temperature(self) -> float | None: + now = time.monotonic() + if now < self._next_fetch: + return self._value + try: + self._value = self._fetch() + self._next_fetch = now + self.ttl + logger.info("open-meteo: %.1f °C, acceleration %.3f", self._value, acceleration(self._value)) + except (OSError, ValueError, KeyError, TypeError) as error: + self._next_fetch = now + self.retry + fallback = "last cached value" if self._value is not None else "no weather factor" + logger.warning("open-meteo unavailable (%s), using %s", error, fallback) + return self._value diff --git a/ml/physics_check.py b/ml/physics_check.py new file mode 100644 index 0000000..a2f3f51 --- /dev/null +++ b/ml/physics_check.py @@ -0,0 +1,77 @@ +"""Проверка физмодели: обучение на двигателях 1–80, проверка на 81–100, которые проигрывает симулятор.""" + +import argparse +import logging +import sys +from pathlib import Path + +import numpy as np +import pandas as pd +import physics +from prepare_rul import load_train + +logger = logging.getLogger("physics_check") + + +def check_model(path: Path) -> None: + model = physics.fit(path) + assert model.a > 0 and model.b > 0, "ожидаемая деградация должна расти с наработкой" + + df = load_train(path) + df = df[df["unit"] > physics.TRAIN_UNITS].sort_values(["unit", "cycle"]).rename(columns={"unit": "equipment_id"}) + result = physics.predict(model, df, {}, factor=1.0) + result["actual"] = result["predicted_degradation"] + result["residual"] + + life = result.groupby("equipment_id")["cycle"].transform("max") + spearman = result["actual"].corr(result["cycle"] / life, method="spearman") + logger.info("spearman(actual degradation, life fraction): %.3f", spearman) + assert spearman > 0.6, "фактическая деградация должна расти к отказу" + + at_failure = result[result["cycle"] == life] + corr = float(np.corrcoef(at_failure["cycle"], at_failure["residual"])[0, 1]) + logger.info("corr(life, residual at failure): %.3f", corr) + assert corr < -0.5, "у быстро изнашивающихся двигателей невязка к отказу должна быть больше" + + for _, row in at_failure.iterrows(): + logger.info("unit %d: life %d, actual %.3f, expected %.3f, residual %+.3f", row["equipment_id"], + row["cycle"], row["actual"], row["predicted_degradation"], row["residual"]) + + +def check_weather() -> None: + assert physics.acceleration(None) == 1.0 + assert abs(physics.acceleration(15.0) - 1.0) < 1e-9 + assert physics.acceleration(0.0) < 1.0 < physics.acceleration(30.0) + + responses = [] + + def flaky() -> float: + if responses: + return responses.pop() + raise OSError("network is unreachable") + + assert physics.Weather(0, 0, fetch=flaky).temperature() is None, "без кэша — без погодного фактора" + + responses.append(20.0) + weather = physics.Weather(0, 0, fetch=flaky, ttl=0, retry=0) + assert weather.temperature() == 20.0 + assert weather.temperature() == 20.0, "при недоступности API используется последнее значение" + logger.info("weather fallback: ok") + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("path", type=Path) + a = ap.parse_args() + + pd.set_option("display.width", 200) + check_model(a.path) + check_weather() + + +if __name__ == "__main__": + logging.basicConfig( + level=logging.INFO, + format="%(asctime)s %(levelname)s %(name)s: %(message)s", + stream=sys.stdout, + ) + main() diff --git a/ml/worker.py b/ml/worker.py index a02a5f2..d0ef904 100644 --- a/ml/worker.py +++ b/ml/worker.py @@ -5,6 +5,7 @@ import numpy as np import pandas as pd +import physics import xgboost as xgb from dotenv import load_dotenv from sqlalchemy import create_engine, text @@ -16,6 +17,10 @@ warning_rul = int(os.getenv("RUL_WARNING")) critical_rul = int(os.getenv("RUL_CRITICAL")) +TRAIN_PATH = os.getenv("TRAIN_PATH", "../data/train_FD001.txt") +WEATHER_LATITUDE = float(os.getenv("WEATHER_LATITUDE", "55.7558")) +WEATHER_LONGITUDE = float(os.getenv("WEATHER_LONGITUDE", "37.6173")) + engine = create_engine(os.getenv("DATABASE_URL")) def get_model() -> xgb.XGBRegressor: @@ -33,6 +38,7 @@ def get_data_from_simulator() -> tuple[pd.DataFrame, pd.DataFrame] | tuple[None, ON sensor_readings.equipment_id = ml_predictions.equipment_id AND sensor_readings.cycle = ml_predictions.cycle WHERE ml_predictions.id IS NULL + ORDER BY sensor_readings.equipment_id, sensor_readings.cycle """) result = connection.execute(query) @@ -50,14 +56,39 @@ def model_predict(df: pd.DataFrame, model: xgb.XGBRegressor) -> np.ndarray: return arr +def get_physics_state() -> dict[int, tuple[int, float, float]]: + """Последний посчитанный цикл по каждому оборудованию: ожидаемая и сглаженная фактическая деградация.""" + query = text(""" + SELECT DISTINCT ON (physics_predictions.equipment_id) + physics_predictions.equipment_id, physics_predictions.cycle, + physics_predictions.predicted_degradation, ml_predictions.residual + FROM physics_predictions + JOIN ml_predictions + ON physics_predictions.equipment_id = ml_predictions.equipment_id + AND physics_predictions.cycle = ml_predictions.cycle + ORDER BY physics_predictions.equipment_id, physics_predictions.cycle DESC + """) + + with engine.connect() as connection: + return {eq_id: (cycle, expected, expected + residual) + for eq_id, cycle, expected, residual in connection.execute(query)} + def write_to_tables(): model = get_model() + physics_model = physics.fit(TRAIN_PATH) + weather = physics.Weather(WEATHER_LATITUDE, WEATHER_LONGITUDE) last_id = 0 + query_first = text(""" + INSERT INTO physics_predictions (equipment_id, cycle, predicted_degradation) + VALUES (:equipment_id, :cycle, :predicted_degradation) + ON CONFLICT DO NOTHING + """) + query_second = text(""" INSERT INTO ml_predictions (equipment_id, cycle, predicted_rul, defect_probability, residual, model_version) - VALUES (:eq_id, :cycle, :rul, 0, 0, :model) + VALUES (:eq_id, :cycle, :rul, 0, :residual, :model) ON CONFLICT DO NOTHING """) @@ -71,11 +102,16 @@ def write_to_tables(): else: t0 = time.perf_counter() arr = model_predict(df_pred, model) + factor = physics.acceleration(weather.temperature()) + df_physics = physics.predict(physics_model, df, get_physics_state(), factor) + with engine.connect() as connection: + connection.execute(query_first, df_physics[['equipment_id', 'cycle', 'predicted_degradation']].to_dict('records')) for i in range (arr.size): df_row = df.iloc[i] - connection.execute(query_second, {'eq_id': df_row['equipment_id'], 'cycle': df_row['cycle'], 'rul':arr[i], 'model': 'first_xgboost_model'}) + connection.execute(query_second, {'eq_id': df_row['equipment_id'], 'cycle': df_row['cycle'], 'rul':arr[i], + 'residual': df_physics['residual'].iat[i], 'model': 'first_xgboost_model'}) connection.commit()