212 KiB
212 KiB
🎯 [MNT04] 고장실적·수리부속 소모 기반 PBL 품목관리
Google Colab 실습 환경 (장비정비·수리부속 / Level 3)
📌 1. 프로젝트 개요 및 배경
| 항목 | 내용 |
|---|---|
| 과정 분류 | 장비정비·수리부속 (Level 3) |
| 학습 소요 시간 | 10개 미션 (5일 / 총 28.5시간) |
| 환경 사양 | Python 3.10+, numpy, pandas, matplotlib, scipy, scikit-learn |
📖 프로젝트 목표
항공·기갑 장비의 엔진과 구동계는 고장이 나기 전에 진동, 온도, 압력에 미세한 변화가 먼저 나타나지만, 고정 임계값만으로는 고부하 운용일을 이상으로 오인하거나 초기 열화를 놓쳐 수리부속을 제때 확보하지 못한다. 본 프로젝트는 고장 실적과 수리부속 소모를 연결해 정비와 부품 관리의 연계 문제를 학습하는 것을 목표로, 교육용 가상 기동체 48대의 센서 시계열과 부품 소모를 코드로 만들어 분석한다. 센서 데이터셋 정의, 잡음 제거와 FFT 특징량 추출, 임계값 Baseline과 부품 소모 파레토 분석, 오토인코더 기반 미세 변화 감지와 파생 피처 모델, 잔존수명(RUL) 예측과 이상 신뢰도 점수, 설명가능한 경보 화면, 검사관 승인(HITL) 기반 조기 입고 시뮬레이터, 대용량 스트림 스트레스 시험, AAR 오경보 분석과 임계 재보정, 재학습을 포함한 통합 도구까지 단계별로 구현한다. 최종 목표는 고장 전에 이상 징후를 일찍 잡아 정비 시점과 수리부속 소요를 예고하되, 조기 입고 결정은 검사관이 근거를 보고 내리도록 지원하는 것이다.
⚠️ [교육용 가상 데이터셋 원칙]
- 본 노트북의 모든 데이터(제원, 로그, 시나리오 등)는 보안 규정을 준수하여 코드 내 난수 기반 가상 데이터로 생성됩니다.
- 실제 부대 위치, 작전 인가, 실제 전장 정보, 군사 기밀 데이터는 일체 포함되지 않습니다.
In [ ]:
# ============================================================================
# [Google Colab 환경 설정] 필수 라이브러리 설치 및 임포트
# ============================================================================
!pip install -q numpy pandas matplotlib scipy scikit-learn
import os
import sys
import time
import math
import json
import re
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
# 그래프 영문 폰트 기본 설정 (환경 독립성 보장)
plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
print(f"실행 환경 준비 완료 (Python {sys.version.split()[0]})")
In [ ]:
"""헬기/기갑 엔진 및 구동계의 진동·온도·압력 센서 시계열, 고장 실적, 수리부속 소모 기록의 데이터셋 스키마·데이터 딕셔너리(피처 20개 이상)와 예지정비·PBL 품목관리 업무 프로세스 맵을 정의하고, 가상 센서 자료로 검증한다."""
import json
import re
import numpy as np
import pandas as pd
from scipy.stats import kurtosis
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# 1) 업무 프로세스 맵: 센서 수집부터 수리부속 소요 환류까지 (각 단계에서 생기는 데이터)
PROCESS = [("센서 수집", "개체 ID, 측정일, 진동 파형, 온도, 압력, 부하"), ("전처리", "잡음 제거, 단위·표기 통일, 결측 보정"),
("특징 추출", "RMS, 첨도, 파고율, FFT 대역 에너지"), ("이상 감지", "이상 점수, 경보 단계"),
("잔존수명(RUL) 예측", "고장 임박 시점, 신뢰 구간"), ("검사관 승인", "조기 입고 제안, 승인·반려 사유"),
("정비·부품 소모", "교체 부품, 수량, 정비 소요시간"), ("고장 실적 연결", "정비 건 ID로 센서 이력과 고장 실적 연결"),
("부품 소요·가동률 집계", "품목별 소모, 가용도(PBL 지표)"), ("사후검토(AAR)", "오경보 원인, 임계 파라미터 조정")]
# 2) 스키마(데이터 딕셔너리): 그룹 / 필드 / 자료형 / 규칙 / 설명
S = {
"sensor": {
"unit_id": (str, {"pattern": r"^U\d{2}$"}, "개체 ID(가상 기동체)"),
"day": (int, {"min": 0, "max": 260}, "운용 경과일(측정일)"),
"vib_rms": (float, {"min": 0, "max": 10}, "진동 RMS(가상 단위)"),
"vib_kurt": (float, {"min": 1, "max": 30}, "진동 첨도(충격성 지표)"),
"temp": (float, {"min": -30, "max": 150}, "엔진 온도(섭씨)"),
"press": (float, {"min": 0, "max": 700}, "오일 압력(kPa)"),
"load": (float, {"min": 0, "max": 1.2}, "운용 부하율"),
},
"feature": {
"band_50": (float, {"min": 0}, "50Hz 부근 FFT 대역 에너지"),
"band_165": (float, {"min": 0}, "165Hz 부근(베어링 결함 대역) 에너지"),
"crest": (float, {"min": 1, "max": 30}, "파고율(최대값/RMS)"),
"hf_ratio": (float, {"min": 0}, "고주파/기본주파 에너지 비"),
"temp_res": (float, {"min": -60, "max": 80}, "부하·외기 대비 온도 잔차"),
},
"failure": {
"event_id": (str, {"pattern": r"^F\d{3}$"}, "고장·정비 건 ID(센서 이력과 연결하는 키)"),
"mode": (str, {"enum": list(MODES)}, "열화(고장) 유형"),
"life": (int, {"min": 60, "max": 260}, "고장까지 운용일수"),
"unplanned": (bool, {}, "비계획 정비 여부"),
},
"parts": {
"part_name": (str, {"enum": [v[2] for v in MODES.values()]}, "소모 수리부속"),
"qty": (int, {"min": 1, "max": 10}, "소모 수량"),
"unit_cost": (float, {"min": 0, "max": 50}, "가상 단가지수"),
},
"target": {
"rul": (float, {"min": 0, "max": 260}, "잔존수명(일) - 예측대상 1"),
"anomaly_score": (float, {"min": 0}, "이상 점수(재구성 오차) - 예측대상 2"),
"availability": (float, {"min": 0, "max": 1}, "가용도(가동일/전체일) - PBL 성과 지표"),
},
}
def check(name, v, typ, rule):
if not isinstance(v, typ) or (typ in (int, float) and isinstance(v, bool)):
return [f"{name}: 자료형 오류({type(v).__name__})"]
e = []
if "enum" in rule and v not in rule["enum"]:
e.append(f"{name}: 허용되지 않은 값 '{v}'")
if "min" in rule and v < rule["min"]:
e.append(f"{name}: 최소값 {rule['min']} 미만")
if "max" in rule and v > rule["max"]:
e.append(f"{name}: 최대값 {rule['max']} 초과")
if "pattern" in rule and not re.match(rule["pattern"], v):
e.append(f"{name}: 형식 불일치 '{v}'")
return e
def validate(rec):
errs = []
for g, fields in S.items():
for k, (typ, rule, _) in fields.items():
if k not in rec.get(g, {}):
errs.append(f"{g}.{k}: 필수 필드 누락")
else:
errs += check(f"{g}.{k}", rec[g][k], typ, rule)
return errs
# 3) [교육용가상자료] 생성 자료에서 20건을 레코드로 만들어 검증
units, snap, W = make_runs()
F = extract(W)
snap = pd.concat([snap, F], axis=1)
snap["temp_res"] = snap["temp"] - (55 + 22 * snap["load"] + 0.5 * snap["ambient"])
snap["hf_ratio"] = (snap["b150"] + snap["b165"] + snap["b200"]) / (snap["b50"] + 1e-6)
ui = units.set_index("unit")
recs = []
for i, r in enumerate(snap.sample(20, random_state=0).itertuples()):
m = r.mode
recs.append({
"sensor": {"unit_id": r.unit, "day": int(r.day), "vib_rms": float(round(r.rms, 3)), "vib_kurt": float(round(r.kurt, 2)),
"temp": float(round(r.temp, 1)), "press": float(round(r.press, 1)), "load": float(round(r.load, 2))},
"feature": {"band_50": float(round(r.b50, 4)), "band_165": float(round(r.b165, 4)), "crest": float(round(r.crest, 2)),
"hf_ratio": float(round(r.hf_ratio, 3)), "temp_res": float(round(r.temp_res, 2))},
"failure": {"event_id": f"F{int(r.unit[1:]):03d}", "mode": m, "life": int(ui.loc[r.unit, "life"]), "unplanned": bool(r.rul < 5)},
"parts": {"part_name": MODES[m][2], "qty": 1, "unit_cost": float(MODES[m][3])},
"target": {"rul": float(r.rul), "anomaly_score": float(round(r.deg, 3)), "availability": 0.92}})
ok = sum(not validate(r) for r in recs)
print(f"가상 센서 자료 {len(snap)}건(개체 {len(units)}대), 표본 {len(recs)}건 중 검증 통과 {ok}건, 오류 {len(recs) - ok}건")
print(f"전체 파형 {W.shape[0]}개 x {W.shape[1]}표본, 온도·압력 결측 {int(snap[['temp', 'press']].isna().sum().sum())}건")
# 4) 오류 주입 시험
bad = json.loads(json.dumps(recs[0]))
bad["sensor"]["unit_id"] = "ENGINE-3" # ID 형식 오류
bad["sensor"]["temp"] = 480.0 # 센서 범위 초과(단위 혼동 의심)
bad["failure"]["mode"] = "turbine" # 정의되지 않은 열화 유형
del bad["parts"]["qty"] # 필수 필드 누락
bad["target"]["rul"] = "30일" # 자료형 오류
print("\n[오류 주입 결과]")
for e in validate(bad):
print(" -", e)
# 5) 프로세스 맵과 데이터 딕셔너리
print("\n[업무 프로세스 맵]")
print(" -> ".join(f"{i}.{p}" for i, (p, _) in enumerate(PROCESS, 1)))
for p, d in PROCESS:
print(f" {p:12s}: {d}")
rows = [(g, k, t.__name__, d) for g, fl in S.items() for k, (t, r, d) in fl.items()]
dd = pd.DataFrame(rows, columns=["그룹", "필드", "자료형", "설명"])
print(f"\n[데이터 딕셔너리] 피처 {len(dd)}개 (목표 20개 이상: {'달성' if len(dd) >= 20 else '미달'})")
print(dd.to_string(index=False))
In [ ]:
"""표기 불일치·중복·단위 혼재·불가능값·결측·단발 스파이크가 섞인 가상 센서 원천 자료에 정제 규칙 10종(파형 스파이크 제거, 저역통과 잡음 제거, KNN 결측 보정 포함)을 적용하고, FFT 변환으로 RMS·첨도 등 주요 특징량을 추출한다. 고장 임박 시점의 긴급 점검 측정은 보존하며 정합성 검증 통과율과 정제 전후 통계를 확인한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import kurtosis
from scipy.signal import medfilt, butter, sosfiltfilt
from sklearn.impute import KNNImputer
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# 1) [교육용가상자료] 원천 자료: 생성 자료에 일부러 오류를 섞는다
units, snap, W0 = make_runs()
rng = np.random.default_rng(2)
S = len(snap)
W = W0.copy()
raw = pd.DataFrame({"unit": snap["unit"].values, "date": [(pd.Timestamp("2026-01-01") + pd.Timedelta(days=int(c))).strftime("%Y-%m-%d") for c in snap["cal"]],
"temp": snap["temp"].round(1).values, "temp_unit": "C", "press": snap["press"].round(1).values, "press_unit": "kPa",
"load": snap["load"].round(2).values, "emergency": snap["rul"].values <= 3})
idx = rng.choice(S, 1500, replace=False)
a, b, c, d_, e, f, g, h_ = np.split(idx, [150, 350, 600, 750, 800, 850, 1050])
raw.loc[a, "unit"] = rng.choice(["u", " U", "unit-"], len(a)) + raw.loc[a, "unit"].str[1:] # ID 표기 불일치
raw.loc[b, "date"] = pd.to_datetime(raw.loc[b, "date"]).dt.strftime("%Y/%m/%d") # 날짜 형식 혼재
raw.loc[c, "temp"] = (raw.loc[c, "temp"] * 9 / 5 + 32).round(1); raw.loc[c, "temp_unit"] = "F" # 온도 단위 혼재
raw.loc[d_, "press"] = (raw.loc[d_, "press"] / 6.895).round(1); raw.loc[d_, "press_unit"] = "psi" # 압력 단위 혼재
raw.loc[e, "temp"] = -999.0; raw.loc[f, "press"] = -5.0 # 불가능 값
raw.loc[g, ["temp", "press"]] = np.nan # 결측
sp = h_[:200] # 파형 단발 스파이크(센서 글리치)
for i in sp:
W[i, rng.integers(0, NS, 2)] += rng.choice([-9.0, 9.0], 2)
nz = h_[200:] # 전자기 간섭: 파형 전체 잡음 증가
W[nz] += rng.normal(0, 0.9, (len(nz), NS))
raw = pd.concat([raw, raw.sample(80, random_state=1)], ignore_index=True) # 중복 80건
print(f"원천 자료 {len(raw)}행 (정상 {S}행 + 중복 80행, 오류 주입 {len(idx)}행, 파형 스파이크 {len(sp)}개·간섭 {len(nz)}개 포함)")
# 2) 정제 규칙 10종 (표 규칙 R1~R6, 결측 R9, 예외 R10 / 파형 규칙 R7~R8)
df, applied = raw.iloc[:S].copy(), []
def rule(name, n):
applied.append((name, int(n)))
rule("R1 중복 행 제거", raw.duplicated().sum())
uid = "U" + df["unit"].str.extract(r"(\d+)")[0]; rule("R2 개체 ID 표기 통일", (uid != df["unit"]).sum()); df["unit"] = uid
t1 = pd.to_datetime(df["date"], format="%Y-%m-%d", errors="coerce"); t2 = pd.to_datetime(df["date"], format="%Y/%m/%d", errors="coerce")
rule("R3 측정일 형식 통일", t1.isna().sum()); df["date"] = t1.fillna(t2)
fm = df["temp_unit"] == "F"; df.loc[fm, "temp"] = ((df.loc[fm, "temp"] - 32) * 5 / 9).round(1); df["temp_unit"] = "C"; rule("R4 온도 단위 통일(F -> C)", fm.sum())
pm = df["press_unit"] == "psi"; df.loc[pm, "press"] = (df.loc[pm, "press"] * 6.895).round(1); df["press_unit"] = "kPa"; rule("R5 압력 단위 통일(psi -> kPa)", pm.sum())
bt, bp = ~df["temp"].between(-30, 150) & df["temp"].notna(), (df["press"] < 0) & df["press"].notna()
df.loc[bt, "temp"] = np.nan; df.loc[bp, "press"] = np.nan; rule("R6 불가능 값 결측 처리", bt.sum() + bp.sum())
Wc = W.copy()
glitch = np.abs(Wc - np.median(Wc, axis=1, keepdims=True)).max(axis=1) > 8 # 단발 스파이크 창
for i in np.where(glitch)[0]:
Wc[i] = medfilt(Wc[i], 3)
rule("R7 파형 단발 스파이크 제거(중앙값 필터)", glitch.sum())
sos = butter(4, 250, fs=FS, output="sos")
hf = (Wc - sosfiltfilt(sos, Wc, axis=1)).std(axis=1) # 250Hz 이상 성분의 크기
noisy = hf > 0.4 # 간섭 잡음 창
Wc[noisy] = sosfiltfilt(sos, Wc[noisy], axis=1); rule("R8 간섭 잡음 창 저역통과(250Hz) 필터", noisy.sum())
miss = df[["temp", "press"]].isna().sum().sum()
X = df.assign(un=df["unit"].str[1:].astype(int) / 50)[["temp", "press", "load", "un"]]
df[["temp", "press"]] = KNNImputer(n_neighbors=5).fit_transform(X)[:, :2]; rule("R9 센서 결측 KNN 보정", miss)
df["rms"] = np.sqrt(((Wc - Wc.mean(axis=1, keepdims=True)) ** 2).mean(axis=1)); lim = df["rms"].quantile(0.99)
df["outlier_flag"] = (df["rms"] > lim) & ~df["emergency"]; rule("R10 긴급 점검 측정 보존(일반 이상값은 표시만)", (df["rms"] > lim).sum() - df["outlier_flag"].sum())
print("\n[적용 규칙과 처리 건수]")
for nme, n in applied:
print(f" {nme:36s} {n:5d}건")
# 3) 정합성 검증
checks = {"중복 0건": not df.duplicated(["unit", "date"]).any(), "ID 형식": df["unit"].str.match(r"^U\d{2}$").all(),
"날짜 변환": df["date"].notna().all(), "온도 단위 C": (df["temp_unit"] == "C").all(), "압력 단위 kPa": (df["press_unit"] == "kPa").all(),
"온도 범위": df["temp"].between(-30, 150).all(), "압력 0 이상": (df["press"] >= 0).all(), "결측 0건": df[["temp", "press"]].notna().all().all(),
"파형 스파이크 잔존 0개": not (np.abs(Wc - np.median(Wc, axis=1, keepdims=True)).max(axis=1) > 8).any(),
f"긴급 점검 {int(df['emergency'].sum())}건 보존": int(df["emergency"].sum()) == int(snap["rul"].le(3).sum())}
rate = np.mean(list(checks.values()))
print(f"\n정합성 검증 {sum(checks.values())}/{len(checks)} 통과, 통과율 {rate:.0%} (목표 100%: {'달성' if rate == 1 else '미달'})")
for k, v in checks.items():
if not v:
print(" 실패:", k)
# 4) FFT 특징량: 정제 효과(깨끗한 원본 파형 대비 오차)
def feats(Wx):
x = Wx - Wx.mean(axis=1, keepdims=True)
return np.sqrt((x ** 2).mean(axis=1)), kurtosis(x, axis=1, fisher=False)
r0, k0 = feats(W0); r1, k1 = feats(W); r2, k2 = feats(Wc)
bad = glitch | noisy
print("\n[특징량 정제 효과] 스파이크·간섭이 있던 창 " + f"{int(bad.sum())}개의 깨끗한 원본 대비 평균 오차")
print(f" RMS 정제 전 {np.abs(r1 - r0)[bad].mean():.3f} -> 정제 후 {np.abs(r2 - r0)[bad].mean():.3f}")
print(f" 첨도 정제 전 {np.abs(k1 - k0)[bad].mean():.2f} -> 정제 후 {np.abs(k2 - k0)[bad].mean():.2f}")
cmp = pd.DataFrame({"정제 전": [len(raw), pd.to_numeric(raw["temp"], errors="coerce").mean(), raw["temp"].min(), raw["press"].min(), raw["press"].isna().sum()],
"정제 후": [len(df), df["temp"].mean(), df["temp"].min(), df["press"].min(), df["press"].isna().sum()]},
index=["행 수", "온도 평균", "온도 최소", "압력 최소", "압력 결측"]).round(1)
print("\n[정제 전후 기술통계]"); print(cmp.to_string())
i = int(np.where(glitch)[0][0])
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5))
ax[0].plot(W[i], label="raw", alpha=0.6); ax[0].plot(Wc[i], label="cleaned"); ax[0].set_title("Waveform with glitch"); ax[0].legend()
f = np.fft.rfftfreq(NS, 1 / FS); ax[1].plot(f, np.abs(np.fft.rfft(Wc[i])) / NS); ax[1].set_title("FFT spectrum (cleaned)"); ax[1].set_xlabel("Hz")
plt.tight_layout(); plt.show()
In [ ]:
"""임계값 초과(진동 RMS·첨도·온도·압력) 기반 1차 센서 이상 징후 감지 Baseline 알고리즘을 구현해 탐지율·오경보율·고장 전 확보 시간을 확인하고, 탐색적 분석(EDA 차트 8종, 피처와 잔존수명의 상관계수·p-value, 수리부속 소모 파레토 ABC 분류)으로 핵심 관리 대상을 식별한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import kurtosis, spearmanr
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# ---------------------------------------------------------------------------
# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)
# ---------------------------------------------------------------------------
BASE = ["rms", "kurt", "crest", "p2p", "b50", "b100", "b150", "b165", "b200", "temp", "press"]
ROLL = ["rms", "kurt", "b150", "b165", "b200", "temp", "press"]
def build_features(snap, W):
D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)
D["temp_res"] = D["temp"] - (55 + 22 * D["load"] + 0.5 * D["ambient"]) # 부하·외기 대비 온도 잔차
D["press_res"] = D["press"] - (450 - 12 * D["load"]) # 부하 대비 압력 잔차
D["rms_n"] = D["rms"] / (0.5 + D["load"]) # 부하 정규화 RMS
D["hf_ratio"] = (D["b150"] + D["b165"] + D["b200"]) / (D["b50"] + 1e-6) # 고주파/기본주파 비
g = D.groupby("unit")
for c in ROLL:
D[c + "_m3"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())
for c in ("rms_n", "temp_res", "press_res"):
D[c + "_d"] = g[c].diff().fillna(0)
return D
DERIVED = ["temp_res", "press_res", "rms_n", "hf_ratio"] + [c + "_m3" for c in ROLL] + ["rms_n_d", "temp_res_d", "press_res_d"]
FEATS = BASE + DERIVED
AE_COLS = ["rms_n", "kurt", "b150", "b165", "b200", "hf_ratio", "temp_res", "press_res"]
def split_units(units, n_train=24):
return list(units["unit"][:n_train]), list(units["unit"][n_train:])
class AE:
"""오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다."""
def __init__(self, cols=AE_COLS, q=0.99, seed=0):
self.cols, self.q, self.seed = cols, q, seed
def fit(self, healthy):
self.sc = StandardScaler().fit(healthy[self.cols])
Z = self.sc.transform(healthy[self.cols])
self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)
self.thr = float(np.quantile(self.err(healthy), self.q))
return self
def contrib(self, d):
Z = self.sc.transform(d[self.cols])
return (Z - self.net.predict(Z)) ** 2 # 입력 항목별 재구성 오차(이상 기여도)
def err(self, d):
return self.contrib(d).mean(axis=1)
def latch(d, e, thr, k=2):
"""같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다."""
out, e = np.zeros(len(d), bool), np.asarray(e)
for idx in d.reset_index(drop=True).groupby("unit").indices.values():
run, on = 0, False
for i in idx:
run = run + 1 if e[i] > thr else 0
on = on or run >= k
out[i] = on
return out
units, snap, W = make_runs()
D = build_features(snap, W)
# 1) Baseline: 단순 임계값(하나라도 넘으면 이상). 임계값은 정상(열화 전) 자료의 분포로 정한다.
healthy = D[D["deg"] == 0]
TH = {"rms": 1.25, "kurt": 3.0, "temp": 92.0, "press_low": 405.0}
D["alarm"] = (D["rms"] > TH["rms"]) | (D["kurt"] > TH["kurt"]) | (D["temp"] > TH["temp"]) | (D["press"] < TH["press_low"])
D["anomaly"] = D["deg"] >= 0.1 # 정답: 열화가 일정 수준 이상 진행된 측정
tp, fp = (D["alarm"] & D["anomaly"]).sum(), (D["alarm"] & ~D["anomaly"]).sum()
fn = (~D["alarm"] & D["anomaly"]).sum()
print(f"[교육용가상자료] 측정 {len(D)}건, 개체 {len(units)}대, 열화 구간 {int(D['anomaly'].sum())}건")
print(f"Baseline 임계값 {TH}")
print(f" 탐지율 {tp / (tp + fn):.1%}, 정밀도 {tp / (tp + fp):.1%}, 정상 측정 오경보율 {fp / (~D['anomaly']).sum():.1%} (고부하 운용일 오경보율 {D.loc[D['hi_load'] & ~D['anomaly'], 'alarm'].mean():.1%})")
lead = []
for u, g in D.groupby("unit"):
a = g[g["alarm"] & g["anomaly"]]
lead.append(a["rul"].max() if len(a) else 0)
units["lead"] = lead
print(f" 고장 전 첫 경보 확보 시간 중앙값 {np.median(lead):.0f}일, 경보 못 낸 개체 {int((units['lead'] == 0).sum())}대/{len(units)}대")
print(" 열화 유형별 탐지율: " + ", ".join(f"{m} {D.loc[D['anomaly'] & (D['mode'] == m), 'alarm'].mean():.0%}" for m in MODES))
# 2) 상관분석: 피처와 잔존수명(RUL) (스피어만)
rows = []
for c in ["rms", "kurt", "crest", "b150", "b165", "b200", "temp", "press", "hf_ratio", "temp_res", "press_res", "load", "ambient"]:
r, p = spearmanr(D[c], D["rul"]); rows.append((c, r, p))
cor = pd.DataFrame(rows, columns=["피처", "rho", "p"]); cor["유의"] = (cor["p"] < 0.05) & (cor["rho"].abs() >= 0.2)
print("\n[상관분석] 피처 | 스피어만 rho | p-value | 유의(p<0.05, |rho|>=0.2)")
for x in cor.sort_values("rho", key=abs, ascending=False).head(8).itertuples():
print(f" {x.피처:10s} {x.rho:+.2f} {x.p:.4f} {'유의' if x.유의 else '-'}")
nsig = int(cor["유의"].sum())
# 3) 파레토(ABC): 고장별 수리부속 소모 비용(가상 단가지수)
cons = {}
for m in units["mode"]:
cons[MODES[m][2]] = cons.get(MODES[m][2], 0) + MODES[m][3]
cons["씰·필터 키트"] = cons.get("씰·필터 키트", 0) + KIT_COST
par = pd.Series(cons).sort_values(ascending=False).to_frame("소모 비용지수")
par["누적비율"] = par["소모 비용지수"].cumsum() / par["소모 비용지수"].sum()
par["등급"] = np.where(par["누적비율"] <= 0.8, "A", np.where(par["누적비율"] <= 0.95, "B", "C")); par.iloc[0, 2] = "A"
print("\n[파레토] 수리부속 소모 비용(고장 1회 교체 기준)")
print(par.assign(누적비율=lambda d: d["누적비율"].map("{:.0%}".format)).round(0).to_string())
top20 = par.head(int(np.ceil(len(par) * 0.2)))
print(f"상위 20% 품목({', '.join(top20.index)})이 전체 소모 비용의 {top20['소모 비용지수'].sum() / par['소모 비용지수'].sum():.0%}")
# 4) EDA 차트 8종
fig, ax = plt.subplots(2, 4, figsize=(16, 7))
ax[0, 0].scatter(D["rul"], D["rms"], s=3); ax[0, 0].set_title("1 RMS vs RUL"); ax[0, 0].invert_xaxis()
for m in MODES:
g = D[D["mode"] == m]; ax[0, 1].scatter(g["deg"], g["kurt"], s=2, label=m)
ax[0, 1].set_title("2 Kurtosis vs degradation"); ax[0, 1].legend(fontsize=6)
for m in MODES:
g = D[D["mode"] == m]; ax[0, 2].scatter(g["deg"], g["temp_res"], s=2)
ax[0, 2].set_title("3 Temp residual vs degradation")
for m in MODES:
g = D[D["mode"] == m]; ax[0, 3].scatter(g["deg"], g["press_res"], s=2)
ax[0, 3].set_title("4 Pressure residual vs degradation")
f = np.fft.rfftfreq(NS, 1 / FS); bm = D["mode"] == "bearing"
ax[1, 0].plot(f, np.abs(np.fft.rfft(W[(bm & (D["deg"] == 0)).values], axis=1)).mean(axis=0) / NS, label="healthy")
ax[1, 0].plot(f, np.abs(np.fft.rfft(W[(bm & (D["deg"] > 0.6)).values], axis=1)).mean(axis=0) / NS, label="degraded"); ax[1, 0].legend(); ax[1, 0].set_title("5 Mean spectrum (bearing)")
ax[1, 1].hist([units.loc[units["mode"] == m, "life"] for m in MODES], bins=8, stacked=True, label=list(MODES)); ax[1, 1].set_title("6 Life by failure mode")
ax[1, 2].scatter(healthy["load"], healthy["rms"], s=2); ax[1, 2].set_title("7 Load vs RMS (healthy)")
ax[1, 3].bar(range(len(par)), par["소모 비용지수"]); ax[1, 3].set_title("8 Spare-part consumption (Pareto)")
plt.tight_layout(); plt.show()
print(f"\nEDA 차트 8종 생성(목표 8종 이상: 달성), 유의한 피처 {nsig}개 (목표 5개 이상: {'달성' if nsig >= 5 else '미달'})")
print("핵심 발견: 1) 고정 임계값은 고부하일 오경보와 온도·압력 유형 놓침이 있음 2) 부하·외기를 보정한 잔차 피처가 열화와 더 뚜렷함 3) 베어링·기어 부품이 소모 비용 대부분")
In [ ]:
"""정상 자료만으로 학습한 오토인코더의 재구성 오차로 미세 진동·온도 변화를 감지(2일 연속 시 경보)해 임계값 Baseline과 비교하고, 파생 피처 14개로 잔존수명(RUL) 기준 모델을 만들어 개체 단위 시계열 교차검증 MAPE를 평가하며 순열 중요도로 피처를 줄이는 예지정비(PdM) 모델을 구현한다."""
import numpy as np
import pandas as pd
from scipy.stats import kurtosis
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.inspection import permutation_importance
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# ---------------------------------------------------------------------------
# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)
# ---------------------------------------------------------------------------
BASE = ["rms", "kurt", "crest", "p2p", "b50", "b100", "b150", "b165", "b200", "temp", "press"]
ROLL = ["rms", "kurt", "b150", "b165", "b200", "temp", "press"]
def build_features(snap, W):
D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)
D["temp_res"] = D["temp"] - (55 + 22 * D["load"] + 0.5 * D["ambient"]) # 부하·외기 대비 온도 잔차
D["press_res"] = D["press"] - (450 - 12 * D["load"]) # 부하 대비 압력 잔차
D["rms_n"] = D["rms"] / (0.5 + D["load"]) # 부하 정규화 RMS
D["hf_ratio"] = (D["b150"] + D["b165"] + D["b200"]) / (D["b50"] + 1e-6) # 고주파/기본주파 비
g = D.groupby("unit")
for c in ROLL:
D[c + "_m3"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())
for c in ("rms_n", "temp_res", "press_res"):
D[c + "_d"] = g[c].diff().fillna(0)
return D
DERIVED = ["temp_res", "press_res", "rms_n", "hf_ratio"] + [c + "_m3" for c in ROLL] + ["rms_n_d", "temp_res_d", "press_res_d"]
FEATS = BASE + DERIVED
AE_COLS = ["rms_n", "kurt", "b150", "b165", "b200", "hf_ratio", "temp_res", "press_res"]
def split_units(units, n_train=24):
return list(units["unit"][:n_train]), list(units["unit"][n_train:])
class AE:
"""오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다."""
def __init__(self, cols=AE_COLS, q=0.99, seed=0):
self.cols, self.q, self.seed = cols, q, seed
def fit(self, healthy):
self.sc = StandardScaler().fit(healthy[self.cols])
Z = self.sc.transform(healthy[self.cols])
self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)
self.thr = float(np.quantile(self.err(healthy), self.q))
return self
def contrib(self, d):
Z = self.sc.transform(d[self.cols])
return (Z - self.net.predict(Z)) ** 2 # 입력 항목별 재구성 오차(이상 기여도)
def err(self, d):
return self.contrib(d).mean(axis=1)
def latch(d, e, thr, k=2):
"""같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다."""
out, e = np.zeros(len(d), bool), np.asarray(e)
for idx in d.reset_index(drop=True).groupby("unit").indices.values():
run, on = 0, False
for i in idx:
run = run + 1 if e[i] > thr else 0
on = on or run >= k
out[i] = on
return out
units, snap, W = make_runs()
D = build_features(snap, W)
tr_u, va_u = split_units(units)
tr, va = D[D["unit"].isin(tr_u)].reset_index(drop=True), D[D["unit"].isin(va_u)].reset_index(drop=True)
print(f"[교육용가상자료] 학습 {len(tr_u)}대({len(tr)}건), 검증 {len(va_u)}대({len(va)}건), 파생 피처 {len(DERIVED)}개 (목표 10개 이상: {'달성' if len(DERIVED) >= 10 else '미달'})")
print(" 파생 피처:", ", ".join(DERIVED))
# 1) 오토인코더 이상 감지 vs 임계값 Baseline (검증 개체)
ae = AE().fit(tr[tr["deg"] == 0])
e = ae.err(va)
va["anomaly"] = va["deg"] >= 0.1
va["ae_alarm"] = latch(va, e, ae.thr, k=2)
va["th_alarm"] = (va["rms"] > 1.25) | (va["kurt"] > 3.0) | (va["temp"] > 92.0) | (va["press"] < 405.0)
va["ae_raw"] = e > ae.thr
def report(name, col):
a = va[col]
tpr, far = (a & va["anomaly"]).sum() / va["anomaly"].sum(), (a & ~va["anomaly"]).sum() / (~va["anomaly"]).sum()
hi = a[va["hi_load"] & ~va["anomaly"]].mean()
lead = [g.loc[g[col] & g["anomaly"], "rul"].max() if (g[col] & g["anomaly"]).any() else 0 for _, g in va.groupby("unit")]
print(f" {name:22s} 탐지율 {tpr:5.1%} 오경보율 {far:5.1%} 고부하일 오경보 {hi:5.1%} 고장 전 확보 중앙값 {np.median(lead):.0f}일")
return tpr, far
print(f"\n[오토인코더 입력 {len(AE_COLS)}개, 정상 자료 재구성 오차 99분위 {ae.thr:.2f}를 기준으로 사용]")
report("임계값 Baseline", "th_alarm")
report("오토인코더(1회 초과)", "ae_raw")
tpr, far = report("오토인코더(2일 연속)", "ae_alarm")
print(f" 미세 변화 구간(열화 0.1~0.3) 탐지율: Baseline {va.loc[va['deg'].between(0.1, 0.3), 'th_alarm'].mean():.1%} -> 오토인코더 {va.loc[va['deg'].between(0.1, 0.3), 'ae_alarm'].mean():.1%}")
# 2) RUL 기준 모델: 학습 개체 일부로 학습, 개체 단위 시계열 교차검증(앞선 개체로 학습 -> 뒤 개체 검증)
CAP = 120 # RUL 상한(일): 열화 초기에는 120일 이상은 구분하지 않는다
def fit_rul(d, cols):
return GradientBoostingRegressor(n_estimators=100, max_depth=3, learning_rate=0.1, subsample=0.8, random_state=0).fit(d[cols], np.minimum(d["rul"], CAP))
def mape(d, p):
m = (d["deg"] >= 0.1) & (d["rul"] >= 10) # 열화가 시작된 구간, 고장 직전 10일 미만은 분모 문제로 제외
y = np.minimum(d["rul"], CAP)[m]
return float(np.mean(np.abs(y - p[m]) / y))
order = list(units["unit"])
folds, cut = [], [12, 24, 36]
for i, c in enumerate(cut):
trn, val = D[D["unit"].isin(order[:c])], D[D["unit"].isin(order[c:c + 12])]
if len(val):
folds.append(mape(val, fit_rul(trn, FEATS).predict(val[FEATS])))
base_only = fit_rul(tr, BASE)
full = fit_rul(tr, FEATS)
m_base, m_full = mape(va, base_only.predict(va[BASE])), mape(va, full.predict(va[FEATS]))
print(f"\n[RUL 모델] 개체 단위 시계열 교차검증 폴드별 MAPE: {', '.join(f'{f:.1%}' for f in folds)} -> 평균 {np.mean(folds):.1%} (목표 25% 이하: {'달성' if np.mean(folds) <= 0.25 else '미달'})")
print(f" 검증 24대 MAPE: 원본 센서 {len(BASE)}개 {m_base:.1%} -> 파생 포함 {len(FEATS)}개 {m_full:.1%} (열화 구간, 잔존 10일 이상)")
# 3) 순열 중요도로 차원 축소
pi = permutation_importance(full, va[FEATS], np.minimum(va["rul"], CAP), n_repeats=3, random_state=0, scoring="neg_mean_absolute_error")
imp = pd.Series(pi.importances_mean, FEATS).sort_values(ascending=False)
print("피처 중요도 상위 5: " + ", ".join(f"{k}({v:.1f})" for k, v in imp.head(5).items()))
top = list(imp.index[:8])
m_top = mape(va, fit_rul(tr, top).predict(va[top]))
print(f"차원 축소 {len(FEATS)}개 -> {len(top)}개: 검증 MAPE {m_full:.1%} -> {m_top:.1%}")
In [ ]:
"""진동·온도·압력 특징으로 부품 고장 임박 시점(RUL)을 예측하는 앙상블(그래디언트 부스팅+랜덤 포레스트+릿지) 모델을 하이퍼파라미터 탐색으로 만들어 단순 선형 기준 모델과 MAE·RMSE·R2를 비교하고, 오토인코더 이상 점수와 예측 구간으로 상태 이상 신뢰도 점수를 산출해 신뢰도 구간별 실제 정확도와 점검 카드를 출력한다."""
import numpy as np
import pandas as pd
from scipy.stats import kurtosis
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# ---------------------------------------------------------------------------
# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)
# ---------------------------------------------------------------------------
BASE = ["rms", "kurt", "crest", "p2p", "b50", "b100", "b150", "b165", "b200", "temp", "press"]
ROLL = ["rms", "kurt", "b150", "b165", "b200", "temp", "press"]
def build_features(snap, W):
D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)
D["temp_res"] = D["temp"] - (55 + 22 * D["load"] + 0.5 * D["ambient"]) # 부하·외기 대비 온도 잔차
D["press_res"] = D["press"] - (450 - 12 * D["load"]) # 부하 대비 압력 잔차
D["rms_n"] = D["rms"] / (0.5 + D["load"]) # 부하 정규화 RMS
D["hf_ratio"] = (D["b150"] + D["b165"] + D["b200"]) / (D["b50"] + 1e-6) # 고주파/기본주파 비
g = D.groupby("unit")
for c in ROLL:
D[c + "_m3"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())
for c in ("rms_n", "temp_res", "press_res"):
D[c + "_d"] = g[c].diff().fillna(0)
return D
DERIVED = ["temp_res", "press_res", "rms_n", "hf_ratio"] + [c + "_m3" for c in ROLL] + ["rms_n_d", "temp_res_d", "press_res_d"]
FEATS = BASE + DERIVED
AE_COLS = ["rms_n", "kurt", "b150", "b165", "b200", "hf_ratio", "temp_res", "press_res"]
def split_units(units, n_train=24):
return list(units["unit"][:n_train]), list(units["unit"][n_train:])
class AE:
"""오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다."""
def __init__(self, cols=AE_COLS, q=0.99, seed=0):
self.cols, self.q, self.seed = cols, q, seed
def fit(self, healthy):
self.sc = StandardScaler().fit(healthy[self.cols])
Z = self.sc.transform(healthy[self.cols])
self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)
self.thr = float(np.quantile(self.err(healthy), self.q))
return self
def contrib(self, d):
Z = self.sc.transform(d[self.cols])
return (Z - self.net.predict(Z)) ** 2 # 입력 항목별 재구성 오차(이상 기여도)
def err(self, d):
return self.contrib(d).mean(axis=1)
def latch(d, e, thr, k=2):
"""같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다."""
out, e = np.zeros(len(d), bool), np.asarray(e)
for idx in d.reset_index(drop=True).groupby("unit").indices.values():
run, on = 0, False
for i in idx:
run = run + 1 if e[i] > thr else 0
on = on or run >= k
out[i] = on
return out
CAP = 120
units, snap, W = make_runs()
D = build_features(snap, W)
D["hf_max"] = D.groupby("unit")["hf_ratio"].cummax() # 지금까지 가장 높았던 고주파 비(열화는 되돌아가지 않는다)
D["temp_res_max"] = D.groupby("unit")["temp_res"].cummax()
FEATS = FEATS + ["day", "hf_max", "temp_res_max"] # 운용일수(day)는 정비 이력에서 아는 값
tr_u, va_u = split_units(units)
tr, va = D[D["unit"].isin(tr_u)].reset_index(drop=True), D[D["unit"].isin(va_u)].reset_index(drop=True)
ytr, yva = np.minimum(tr["rul"], CAP), np.minimum(va["rul"], CAP)
print(f"[교육용가상자료] 학습 {len(tr_u)}대 {len(tr)}건, Hold-out 검증 {len(va_u)}대 {len(va)}건 (RUL은 {CAP}일 상한)")
# 1) 기준 모델: 진동 RMS와 온도만 쓰는 선형 회귀
base = LinearRegression().fit(tr[["rms", "temp"]], ytr)
pb = base.predict(va[["rms", "temp"]])
# 2) 고도화: 하이퍼파라미터 탐색(학습 개체의 뒤쪽 6대를 내부 검증) 후 앙상블
inner_tr, inner_va = tr[tr["unit"].isin(tr_u[:18])], tr[tr["unit"].isin(tr_u[18:])]
def ens(d_tr, y, d_va, depth, alpha):
ms = [GradientBoostingRegressor(n_estimators=100, max_depth=depth, learning_rate=0.1, subsample=0.8, random_state=0),
RandomForestRegressor(n_estimators=40, max_depth=10, min_samples_leaf=5, n_jobs=-1, random_state=0),
make_pipeline(StandardScaler(), Ridge(alpha=alpha))]
for m in ms:
m.fit(d_tr[FEATS], y)
return ms, np.mean([m.predict(d_va[FEATS]) for m in ms], axis=0)
score = {(dp, al): mean_absolute_error(np.minimum(inner_va["rul"], CAP), ens(inner_tr, np.minimum(inner_tr["rul"], CAP), inner_va, dp, al)[1])
for dp in (3, 4) for al in (1, 10)}
dp, al = min(score, key=score.get)
models, pe = ens(tr, ytr, va, dp, al)
names = ["GBR", "RF", "Ridge"]
print(f"선택 하이퍼파라미터: GBR depth={dp}, ridge alpha={al} (내부 검증 MAE {score[(dp, al)]:.1f}일)")
def metrics(p):
return mean_absolute_error(yva, p), mean_squared_error(yva, p) ** 0.5, r2_score(yva, p)
mb, me = metrics(pb), metrics(pe)
print("\n[RUL 예측 성능] 모델 | MAE(일) | RMSE(일) | R2")
print(f" 기준(RMS+온도 선형) {mb[0]:6.2f} {mb[1]:6.2f} {mb[2]:.3f}")
for n, m in zip(names, models):
x = metrics(m.predict(va[FEATS]))
print(f" 단일 {n:12s} {x[0]:6.2f} {x[1]:6.2f} {x[2]:.3f}")
print(f" 고도화(앙상블) {me[0]:6.2f} {me[1]:6.2f} {me[2]:.3f}")
imp_mae, imp_rmse = 1 - me[0] / mb[0], 1 - me[1] / mb[1]
print(f"개선율 MAE {imp_mae:.1%}, RMSE {imp_rmse:.1%} (목표 20% 이상: {'달성' if min(imp_mae, imp_rmse) >= 0.2 else '미달'}), R2 {me[2]:.3f} (목표 0.85 이상: {'달성' if me[2] >= 0.85 else '미달'})")
for lo, hi, nm in [(0, 30, "고장 30일 이내"), (30, 90, "30~90일"), (90, 9999, "90일 이상")]:
m = (va["rul"] >= lo) & (va["rul"] < hi)
print(f" 구간 {nm:10s} MAE {mean_absolute_error(yva[m], pe[m]):5.1f}일 (기준 {mean_absolute_error(yva[m], pb[m]):5.1f}일)")
# 3) 상태 이상 신뢰도 점수 = 오토인코더 이상 점수(정상 대비 배수)와 앙상블 구성 모델 간 의견 일치도를 결합
ae = AE().fit(tr[tr["deg"] == 0])
e = ae.err(va)
z = np.log(e / ae.thr)
conf = 1 / (1 + np.exp(-2.5 * z)) # 재구성 오차가 기준과 같으면 0.5
spread = np.std([m.predict(va[FEATS]) for m in models], axis=0)
agree = np.clip(1 - spread / 40, 0, 1) # 모델 간 RUL 의견 차이가 작을수록 높음
score_c = conf * (0.7 + 0.3 * agree)
truth = va["deg"] >= 0.1
print("\n[상태 이상 신뢰도] 구간 | 건수 | 평균 점수 | 실제 열화 비율")
bins = [0, 0.2, 0.5, 0.8, 1.01]
for a, b in zip(bins[:-1], bins[1:]):
m = (score_c >= a) & (score_c < b)
print(f" {a:.1f}~{min(b, 1):.1f} {m.sum():5d} {score_c[m].mean():.2f} {truth[m].mean():.2f}")
ece = sum(((score_c >= a) & (score_c < b)).mean() * abs(score_c[(score_c >= a) & (score_c < b)].mean() - truth[(score_c >= a) & (score_c < b)].mean())
for a, b in zip(bins[:-1], bins[1:]) if ((score_c >= a) & (score_c < b)).any())
print(f"기대 보정 오차(ECE) {ece:.3f}")
# 4) 점검 카드 예시(개체 하나의 최근 상태)
print("\n[점검 카드 예시] 개체 | 측정일 | 이상 신뢰도 | 예측 RUL(일) | 실제 RUL(일)")
va["conf"], va["rul_hat"] = score_c, pe
for u in va_u[:3]:
g = va[va["unit"] == u]
r = g[g["rul"] == 40].iloc[0] if (g["rul"] == 40).any() else g.iloc[-1]
print(f" {u} {int(r['day']):3d}일째 {r['conf']:.2f} {r['rul_hat']:5.1f} {int(r['rul'])} -> {'정비 계획 수립 권고' if r['conf'] >= 0.8 else '추가 관찰'}(검사관 확인)")
In [ ]:
"""장비 센서 상태를 개체별로 모니터링해 이상 신호를 경보로 띄우는 정비병용 화면(상태표·점검 카드)을 구성하고, 오토인코더 재구성 오차에서 경보 건마다 기여도 상위 3개 인자를 자동 산출하며(근거 유효성 검증), 의심 부품과 자연어 요약 문장을 생성하는 설명가능성(XAI) 모듈을 구현한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import kurtosis
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# ---------------------------------------------------------------------------
# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)
# ---------------------------------------------------------------------------
BASE = ["rms", "kurt", "crest", "p2p", "b50", "b100", "b150", "b165", "b200", "temp", "press"]
ROLL = ["rms", "kurt", "b150", "b165", "b200", "temp", "press"]
def build_features(snap, W):
D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)
D["temp_res"] = D["temp"] - (55 + 22 * D["load"] + 0.5 * D["ambient"]) # 부하·외기 대비 온도 잔차
D["press_res"] = D["press"] - (450 - 12 * D["load"]) # 부하 대비 압력 잔차
D["rms_n"] = D["rms"] / (0.5 + D["load"]) # 부하 정규화 RMS
D["hf_ratio"] = (D["b150"] + D["b165"] + D["b200"]) / (D["b50"] + 1e-6) # 고주파/기본주파 비
g = D.groupby("unit")
for c in ROLL:
D[c + "_m3"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())
for c in ("rms_n", "temp_res", "press_res"):
D[c + "_d"] = g[c].diff().fillna(0)
return D
DERIVED = ["temp_res", "press_res", "rms_n", "hf_ratio"] + [c + "_m3" for c in ROLL] + ["rms_n_d", "temp_res_d", "press_res_d"]
FEATS = BASE + DERIVED
AE_COLS = ["rms_n", "kurt", "b150", "b165", "b200", "hf_ratio", "temp_res", "press_res"]
def split_units(units, n_train=24):
return list(units["unit"][:n_train]), list(units["unit"][n_train:])
class AE:
"""오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다."""
def __init__(self, cols=AE_COLS, q=0.99, seed=0):
self.cols, self.q, self.seed = cols, q, seed
def fit(self, healthy):
self.sc = StandardScaler().fit(healthy[self.cols])
Z = self.sc.transform(healthy[self.cols])
self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)
self.thr = float(np.quantile(self.err(healthy), self.q))
return self
def contrib(self, d):
Z = self.sc.transform(d[self.cols])
return (Z - self.net.predict(Z)) ** 2 # 입력 항목별 재구성 오차(이상 기여도)
def err(self, d):
return self.contrib(d).mean(axis=1)
def latch(d, e, thr, k=2):
"""같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다."""
out, e = np.zeros(len(d), bool), np.asarray(e)
for idx in d.reset_index(drop=True).groupby("unit").indices.values():
run, on = 0, False
for i in idx:
run = run + 1 if e[i] > thr else 0
on = on or run >= k
out[i] = on
return out
units, snap, W = make_runs()
D = build_features(snap, W)
tr_u, va_u = split_units(units)
tr, va = D[D["unit"].isin(tr_u)].reset_index(drop=True), D[D["unit"].isin(va_u)].reset_index(drop=True)
ae = AE().fit(tr[tr["deg"] == 0])
C = ae.contrib(va) # 건 x 입력항목 기여도
e = C.mean(axis=1)
va["alarm"] = latch(va, e, ae.thr, k=2)
va["score"] = e / ae.thr
LABEL = {"rms_n": "부하 보정 진동 크기(RMS)", "kurt": "진동 충격성(첨도)", "b150": "150Hz 부근 진동(기어 맞물림 대역)", "b165": "165Hz 부근 진동(베어링 결함 대역)",
"b200": "200Hz 부근 진동(펌프 대역)", "hf_ratio": "고주파/기본주파 에너지 비", "temp_res": "부하 대비 온도 상승", "press_res": "부하 대비 압력 저하"}
PART = {"b165": "주 베어링 세트", "kurt": "주 베어링 세트", "b150": "구동 기어 세트", "b200": "오일 펌프", "press_res": "오일 펌프", "temp_res": "냉각 팬 모터"}
TRUE_PART = {m: MODES[m][2] for m in MODES}
cols = ae.cols
# 1) 경보 건마다 기여도 상위 3개 인자 자동 산출
idx = np.where(va["alarm"].values)[0]
top3 = np.argsort(-C[idx], axis=1)[:, :3]
ok = sum(len(set(t)) == 3 for t in top3)
print(f"[교육용가상자료] 검증 {len(va_u)}대 {len(va)}건 중 경보 {len(idx)}건: 기여도 상위 3개 인자 산출 {ok}건 (산출률 {ok / len(idx):.0%}, 목표 100%: {'달성' if ok == len(idx) else '미달'})")
# 2) 근거 유효성: 1순위 인자를 정상 평균(표준화 0)으로 되돌리면 이상 점수가 낮아지는가
Z = ae.sc.transform(va[cols].iloc[idx])
Z2 = Z.copy()
Z2[np.arange(len(idx)), top3[:, 0]] = 0.0
e1, e2 = ((Z - ae.net.predict(Z)) ** 2).mean(axis=1), ((Z2 - ae.net.predict(Z2)) ** 2).mean(axis=1)
print(f"근거 유효성: 1순위 인자를 정상값으로 되돌리면 이상 점수가 낮아지는 비율 {(e2 < e1).mean():.1%} (평균 {e1.mean():.2f} -> {e2.mean():.2f})")
# 3) 의심 부품: 1순위 인자 -> 부품 매핑 (열화가 뚜렷한 경보에서 정답 비교)
sub = va.iloc[idx].copy()
sub["top1"] = [cols[t[0]] for t in top3]
sub["hint"] = sub["top1"].map(PART).fillna("미분류(점검 필요)")
sub["true_part"] = sub["mode"].map(TRUE_PART)
strong = sub[sub["deg"] >= 0.3]
print(f"의심 부품 추정 정확도(열화 0.3 이상 경보 {len(strong)}건): {(strong['hint'] == strong['true_part']).mean():.1%}")
print(" 1순위 인자 분포: " + ", ".join(f"{LABEL[k].split('(')[0]} {v}" for k, v in sub["top1"].value_counts().head(4).items()))
def summarize(r, t):
f = [LABEL[cols[i]] for i in t]
part = PART.get(cols[t[0]], "원인 불명")
return f"{part} 이상 가능성(이상 점수 {r['score']:.1f}배). 주요 근거: {f[0]}, {f[1]}, {f[2]}."
# 4) 정비병용 화면: 점검 카드와 상태표
print("\n[점검 카드] 개체 | 측정일 | 실제 열화 유형")
for m in MODES:
k = sub[(sub["mode"] == m) & (sub["deg"] > 0.5)].index[0]
j = list(sub.index).index(k)
r = sub.loc[k]
print(f" {r['unit']} {int(r['day'])}일째 ({m})")
print(" " + summarize(r, top3[j]))
print(" 근거 기여도: " + " / ".join(f"{LABEL[cols[i]].split('(')[0]}({C[idx[j], i]:.1f})" for i in top3[j]))
CAL = int(va.groupby('cal')['alarm'].sum().idxmax()) # 경보가 가장 많은 날을 가상 기준일로 쓴다
act = va[(va["cal"] == CAL)]
reg = {"정상": 0, "주의": 0, "경보": 0}
rows = []
for r in act.itertuples():
st = "경보" if r.alarm else ("주의" if r.score > 0.7 else "정상")
reg[st] += 1
for st in ("경보", "주의"):
for r in act[(act["alarm"] if st == "경보" else (~act["alarm"] & (act["score"] > 0.7)))].itertuples():
rows.append((r.unit, st, f"{r.score:.1f}", int(r.rul)))
print(f"\n[운용 상태표] 가상 기준일 {CAL}: 가동 중 {len(act)}대 -> " + ", ".join(f"{k} {v}대" for k, v in reg.items()))
for r in rows[:5]:
print(f" {r[0]} {r[1]} 이상 점수 {r[2]}배 (실제 잔존 {r[3]}일)")
plt.figure(figsize=(8, 4))
for u in sub["unit"].unique()[:3]:
g = va[va["unit"] == u]
plt.plot(g["day"], g["score"], label=u)
plt.axhline(1.0, ls="--", color="k"); plt.yscale("log"); plt.xlabel("day"); plt.ylabel("anomaly score (x threshold)"); plt.legend(); plt.title("Monitoring view: anomaly score")
plt.tight_layout(); plt.show()
In [ ]:
"""센서 이상 경보가 나면 예방 정비 조기 입고일을 제안하되 자동 확정하지 않고 항공/기갑 검사관 승인 대기열에 올려 승인·반려·무응답(기본 계획 유지)을 사유와 함께 해시 체인 감사 로그에 남기는 HITL 승인 모듈을 시뮬레이션하고, 정기 교체 방식 대비 정비·고장 비용 절감률과 승인 지연·전시 전환(가동 강도 급증) 민감도를 분석한다."""
import hashlib
import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import kurtosis
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# ---------------------------------------------------------------------------
# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)
# ---------------------------------------------------------------------------
BASE = ["rms", "kurt", "crest", "p2p", "b50", "b100", "b150", "b165", "b200", "temp", "press"]
ROLL = ["rms", "kurt", "b150", "b165", "b200", "temp", "press"]
def build_features(snap, W):
D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)
D["temp_res"] = D["temp"] - (55 + 22 * D["load"] + 0.5 * D["ambient"]) # 부하·외기 대비 온도 잔차
D["press_res"] = D["press"] - (450 - 12 * D["load"]) # 부하 대비 압력 잔차
D["rms_n"] = D["rms"] / (0.5 + D["load"]) # 부하 정규화 RMS
D["hf_ratio"] = (D["b150"] + D["b165"] + D["b200"]) / (D["b50"] + 1e-6) # 고주파/기본주파 비
g = D.groupby("unit")
for c in ROLL:
D[c + "_m3"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())
for c in ("rms_n", "temp_res", "press_res"):
D[c + "_d"] = g[c].diff().fillna(0)
return D
DERIVED = ["temp_res", "press_res", "rms_n", "hf_ratio"] + [c + "_m3" for c in ROLL] + ["rms_n_d", "temp_res_d", "press_res_d"]
FEATS = BASE + DERIVED
AE_COLS = ["rms_n", "kurt", "b150", "b165", "b200", "hf_ratio", "temp_res", "press_res"]
def split_units(units, n_train=24):
return list(units["unit"][:n_train]), list(units["unit"][n_train:])
class AE:
"""오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다."""
def __init__(self, cols=AE_COLS, q=0.99, seed=0):
self.cols, self.q, self.seed = cols, q, seed
def fit(self, healthy):
self.sc = StandardScaler().fit(healthy[self.cols])
Z = self.sc.transform(healthy[self.cols])
self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)
self.thr = float(np.quantile(self.err(healthy), self.q))
return self
def contrib(self, d):
Z = self.sc.transform(d[self.cols])
return (Z - self.net.predict(Z)) ** 2 # 입력 항목별 재구성 오차(이상 기여도)
def err(self, d):
return self.contrib(d).mean(axis=1)
def latch(d, e, thr, k=2):
"""같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다."""
out, e = np.zeros(len(d), bool), np.asarray(e)
for idx in d.reset_index(drop=True).groupby("unit").indices.values():
run, on = 0, False
for i in idx:
run = run + 1 if e[i] > thr else 0
on = on or run >= k
out[i] = on
return out
PERIOD = 90 # 기본 계획: 정기 교체 주기(운용일). 가장 빠른 고장 사례보다 짧게 정해 사후 고장을 막는 보수적 주기
UNPLANNED = (2.0, 12.0) # [교육용가상자료] 고장 후 정비 비용 = 부품 단가 x 2.0 + 연쇄 손상·가동중단 12.0
PLANNED_FIXED = 1.0 # 입고 정비 작업 고정 비용(지수)
INSPECT = {"base": 0.3, "pdm": 0.05} # 30일당 점검 비용: 정기 방식은 분해 점검, 센서 감시 방식은 경보 시 확인 점검 위주(가정)
class Audit:
def __init__(self):
self.log = []
def add(self, **body):
prev = self.log[-1]["hash"] if self.log else "GENESIS"
body = {"seq": len(self.log), **body}
self.log.append({**body, "prev": prev, "hash": hashlib.sha256((prev + json.dumps(body, sort_keys=True, ensure_ascii=False)).encode()).hexdigest()})
def verify(self):
prev = "GENESIS"
for e in self.log:
b = {k: v for k, v in e.items() if k not in ("prev", "hash")}
if e["prev"] != prev or hashlib.sha256((prev + json.dumps(b, sort_keys=True, ensure_ascii=False)).encode()).hexdigest() != e["hash"]:
return e["seq"]
prev = e["hash"]
return None
def add_feat(D):
D["hf_max"] = D.groupby("unit")["hf_ratio"].cummax(); D["temp_res_max"] = D.groupby("unit")["temp_res"].cummax()
return D
F7 = FEATS + ["day", "hf_max", "temp_res_max"]
u1, s1, w1 = make_runs(seed=1)
D1 = add_feat(build_features(s1, w1))
tr_u, _ = split_units(u1)
tr = D1[D1["unit"].isin(tr_u)]
ae = AE().fit(tr[tr["deg"] == 0])
gbr = GradientBoostingRegressor(n_estimators=100, max_depth=3, learning_rate=0.1, subsample=0.8, random_state=0).fit(tr[F7], np.minimum(tr["rul"], 120))
# 시험 개체: 다른 시드로 만든 가상 개체 3묶음(144대)에서 개체별 첫 경보일과 그때의 RUL 예측을 모은다
cases = []
for sd in (2, 3, 4):
un, sn, ww = make_runs(seed=sd)
D = add_feat(build_features(sn, ww))
D["alarm"] = latch(D, ae.err(D), ae.thr, k=2)
D["rul_hat"] = gbr.predict(D[F7])
for u in un.itertuples():
g = D[D["unit"] == u.unit]
a = g[g["alarm"]]
cases.append({"id": f"S{sd}-{u.unit}", "mode": u.mode, "life": int(u.life), "a": int(a["day"].iloc[0]) if len(a) else None,
"rul_hat": float(a["rul_hat"].iloc[0]) if len(a) else None, "load": float(a["load"].iloc[0]) if len(a) else None})
C = pd.DataFrame(cases)
print(f"[교육용가상자료] 시험 개체 {len(C)}대(수명 {C['life'].min()}~{C['life'].max()}일), 경보 발생 {int(C['a'].notna().sum())}대")
def cost(mode, life, maint, kind="pdm"):
p = MODES[mode][3] + KIT_COST + INSPECT[kind] * min(maint if maint is not None else life, life) / 30
if maint is None or maint >= life: # 정비 전에 고장
return MODES[mode][3] * UNPLANNED[0] + KIT_COST + UNPLANNED[1] + INSPECT[kind] * life / 30, True
return p + PLANNED_FIXED + MODES[mode][3] * (life - maint) / life, False # 남은 수명 낭비분 포함
def commander(row, rng):
"""[모의 검사관] 잔존 예측이 길면 반려, 10%는 무응답."""
if rng.random() < 0.10:
return "무응답", "기한 내 응답 없음 -> 기본 계획(정기 교체) 유지"
if row["rul_hat"] > 90:
return "반려", f"예측 잔존 {row['rul_hat']:.0f}일로 충분 -> 현 계획 유지"
if row["load"] > 1.0:
return "반려", "고부하 운용일 신호로 판단 -> 추가 관찰"
return "승인", f"이상 지속, 예측 잔존 {row['rul_hat']:.0f}일 -> 조기 입고"
def run(mode, delay=2, scale=1.0, seed=0, audit=None):
"""mode: base(정기 교체) / hitl(제안 -> 검사관 결정) / auto(모든 제안 적용: 비교용 상한, 실제 운용 금지)"""
rng = np.random.default_rng(seed)
tot, fails, dec = 0.0, 0, {"승인": 0, "반려": 0, "무응답": 0}
for _, r in C.iterrows():
life = r["life"] / scale
maint, kind = PERIOD, "base"
if mode != "base" and r["a"] is not None:
a, rh = r["a"] / scale, r["rul_hat"] / scale
inbound = a + max(int(0.4 * rh), 3)
row = {"rul_hat": rh, "load": r["load"]}
if mode == "auto":
d, why = "승인", "비교용 전부 적용"
else:
d, why = commander(row, rng)
audit.add(event="PROPOSE", id=r["id"], day=round(a), inbound=int(inbound), status="PENDING")
audit.add(event="DECIDE", id=r["id"], decision=d, reason=why, by="검사관(모의)")
dec[d] += 1
if d == "승인":
maint, kind = inbound + (delay if mode == "hitl" else 0), "pdm"
c, f = cost(r["mode"], life, maint, kind)
tot += c; fails += f
return tot, fails, dec
b_cost, b_fail, _ = run("base")
audit = Audit()
h_cost, h_fail, dec = run("hitl", seed=3, audit=audit)
a_cost, a_fail, _ = run("auto")
print("\n[비용 비교] 방식 | 비용지수 | 사후 고장 건수")
for nm, c, f in [("정기 교체(기본 계획)", b_cost, b_fail), ("경보 제안 + 검사관 승인", h_cost, h_fail), ("제안 전부 적용(비교용)", a_cost, a_fail)]:
print(f" {nm:22s} {c:8.0f} {f:3d}")
cut = 1 - h_cost / b_cost
print(f"HITL 적용 비용 절감률 {cut:.1%} (목표 15% 이상: {'달성' if cut >= 0.15 else '미달'}), 고장 {b_fail} -> {h_fail}건")
print(f"결정 {sum(dec.values())}건: " + ", ".join(f"{k} {v}건" for k, v in dec.items()))
n_prop = sum(e["event"] == "PROPOSE" for e in audit.log); n_dec = sum(e["event"] == "DECIDE" for e in audit.log)
print(f"제안 {n_prop}건 중 결정 기록 {n_dec}건(로깅률 {n_dec / n_prop:.0%}), 자동 승인 0건, 무응답 시 기본 계획 유지")
rej = [e["reason"] for e in audit.log if e.get("decision") == "반려"]
print("반려 사유 예:", rej[0] if rej else "-")
print(f"감사 로그 무결성: {'정상' if audit.verify() is None else '이상'}", end=", ")
audit.log[5]["reason"] = "사후 수정"
print(f"사후 변조 시 seq {audit.verify()}에서 탐지")
print("\n[민감도 1] 승인 후 입고까지 지연(일) -> 비용 절감률")
for d in (0, 2, 7, 14, 21):
c, f, _ = run("hitl", delay=d, seed=3, audit=Audit())
print(f" {d:2d}일: {1 - c / b_cost:.1%} (고장 {f}건)")
print("[민감도 2] 가동 강도 급증(전시 전환 가정, 열화가 빨라져 수명이 1/배수로 단축) -> 절감률")
for s in (1.0, 1.3, 1.6, 2.0):
bc, bf, _ = run("base", scale=s); hc, hf, _ = run("hitl", scale=s, seed=3, audit=Audit())
print(f" x{s}: 정기 교체 비용 {bc:.0f}(고장 {bf}) -> HITL {hc:.0f}(고장 {hf}), 절감 {1 - hc / bc:.1%}")
plt.figure(figsize=(6, 4))
plt.bar(["periodic", "HITL", "auto(ref)"], [b_cost, h_cost, a_cost]); plt.ylabel("cost index"); plt.title("HITL early-inbound simulator")
plt.tight_layout(); plt.show()
In [ ]:
"""장시간 기동 센서 데이터를 1,000~100,000건 규모로 한꺼번에 주입해 검증·특징 추출·오토인코더 점수·경보 판정을 배치로 처리하는 스트림 시험을 수행하고, 처리 시간과 상태 화면 구성 시간, 경보 정확도, 불량 패킷(결측·고정값·포화·중복) 격리와 유실 0건 여부를 확인한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import kurtosis
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# ---------------------------------------------------------------------------
# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)
# ---------------------------------------------------------------------------
BASE = ["rms", "kurt", "crest", "p2p", "b50", "b100", "b150", "b165", "b200", "temp", "press"]
ROLL = ["rms", "kurt", "b150", "b165", "b200", "temp", "press"]
def build_features(snap, W):
D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)
D["temp_res"] = D["temp"] - (55 + 22 * D["load"] + 0.5 * D["ambient"]) # 부하·외기 대비 온도 잔차
D["press_res"] = D["press"] - (450 - 12 * D["load"]) # 부하 대비 압력 잔차
D["rms_n"] = D["rms"] / (0.5 + D["load"]) # 부하 정규화 RMS
D["hf_ratio"] = (D["b150"] + D["b165"] + D["b200"]) / (D["b50"] + 1e-6) # 고주파/기본주파 비
g = D.groupby("unit")
for c in ROLL:
D[c + "_m3"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())
for c in ("rms_n", "temp_res", "press_res"):
D[c + "_d"] = g[c].diff().fillna(0)
return D
DERIVED = ["temp_res", "press_res", "rms_n", "hf_ratio"] + [c + "_m3" for c in ROLL] + ["rms_n_d", "temp_res_d", "press_res_d"]
FEATS = BASE + DERIVED
AE_COLS = ["rms_n", "kurt", "b150", "b165", "b200", "hf_ratio", "temp_res", "press_res"]
def split_units(units, n_train=24):
return list(units["unit"][:n_train]), list(units["unit"][n_train:])
class AE:
"""오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다."""
def __init__(self, cols=AE_COLS, q=0.99, seed=0):
self.cols, self.q, self.seed = cols, q, seed
def fit(self, healthy):
self.sc = StandardScaler().fit(healthy[self.cols])
Z = self.sc.transform(healthy[self.cols])
self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)
self.thr = float(np.quantile(self.err(healthy), self.q))
return self
def contrib(self, d):
Z = self.sc.transform(d[self.cols])
return (Z - self.net.predict(Z)) ** 2 # 입력 항목별 재구성 오차(이상 기여도)
def err(self, d):
return self.contrib(d).mean(axis=1)
def latch(d, e, thr, k=2):
"""같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다."""
out, e = np.zeros(len(d), bool), np.asarray(e)
for idx in d.reset_index(drop=True).groupby("unit").indices.values():
run, on = 0, False
for i in idx:
run = run + 1 if e[i] > thr else 0
on = on or run >= k
out[i] = on
return out
units, snap, W = make_runs()
D = build_features(snap, W)
tr_u, va_u = split_units(units)
tr, va = D[D["unit"].isin(tr_u)].reset_index(drop=True), D[D["unit"].isin(va_u)].reset_index(drop=True)
Wva = W[D["unit"].isin(va_u).values]
ae = AE().fit(tr[tr["deg"] == 0])
KEEP = ["unit", "day", "load", "ambient", "temp", "press"]
def score_batch(w, meta):
"""스트림 한 묶음: 파형 -> 특징 -> 이상 점수(배수). 점수는 개체 연속성(래치)과 무관한 단건 값."""
f = extract(w)
d = pd.DataFrame({"rms_n": f["rms"] / (0.5 + meta["load"]), "kurt": f["kurt"], "b150": f["b150"], "b165": f["b165"], "b200": f["b200"],
"hf_ratio": (f["b150"] + f["b165"] + f["b200"]) / (f["b50"] + 1e-6),
"temp_res": meta["temp"] - (55 + 22 * meta["load"] + 0.5 * meta["ambient"]), "press_res": meta["press"] - (450 - 12 * meta["load"])})
return ae.err(d) / ae.thr
def valid(w, meta):
"""불량 패킷 판정: 결측(NaN), 센서 고정(분산 0), 포화(|값|>12), 값 범위 이탈."""
bad = np.isnan(w).any(axis=1) | (np.nanstd(w, axis=1) < 1e-6) | (np.abs(np.nan_to_num(w)).max(axis=1) > 12)
bad |= ~meta["temp"].between(20, 160).values | ~meta["press"].between(100, 600).values
return ~bad
def make_stream(n, bad_rate, seed=0):
rng = np.random.default_rng(seed)
pick = rng.integers(0, len(va), n)
w = Wva[pick] + rng.normal(0, 0.02, (n, NS))
meta = va.iloc[pick][KEEP + ["deg"]].reset_index(drop=True)
kind = np.array([""] * n, dtype=object)
bad = rng.random(n) < bad_rate
k = rng.integers(0, 4, n)
for t, nm in enumerate(["결측 패킷", "센서 고정", "포화", "범위 이탈"]):
m = bad & (k == t)
kind[m] = nm
if t == 0: w[m, 10:40] = np.nan
if t == 1: w[m] = 0.5
if t == 2: w[m, ::7] = 50.0
if t == 3: meta.loc[m, "temp"] = 400.0
return w, meta, kind
# 1) 1,000건 동시 주입
t0 = time.perf_counter()
w, meta, kind = make_stream(1000, 0.06)
ok = valid(w, meta)
s = np.full(1000, np.nan)
s[ok] = score_batch(w[ok], meta[ok].reset_index(drop=True))
sec = time.perf_counter() - t0
inj = int((kind != "").sum())
print(f"[교육용가상자료] 1,000건 동시 주입: 정상 처리 {int(ok.sum())}건, 격리 {int((~ok).sum())}건, 유실 {1000 - int(ok.sum()) - int((~ok).sum())}건, 처리 시간 {sec:.2f}초 (목표 2초 이내: {'달성' if sec <= 2 else '미달'})")
print(f" 불량 주입 {inj}건 중 격리 {int(((kind != '') & ~ok).sum())}건 (" + ", ".join(f"{k} {int((kind == k).sum())}" for k in ["결측 패킷", "센서 고정", "포화", "범위 이탈"]) + ")")
# 2) 경보 정확도: 정상 구간과 열화가 뚜렷한 구간(deg>=0.3)을 판정 대상으로 하고, 열화 초기 회색 구간은 정확도에서 제외하고 따로 표시
al = s > 1.0
y = (meta["deg"] >= 0.3).values
eval_m = ok & ((meta["deg"] == 0).values | y)
acc = (al[eval_m] == y[eval_m]).mean()
tp, fp, fn = (al & y & ok).sum(), (al & ~y & eval_m).sum(), (~al & y & ok).sum()
print(f"경보 정확도 {acc:.1%} (목표 95% 이상: {'달성' if acc >= 0.95 else '미달'}), 재현율 {tp / (tp + fn):.1%}, 정밀도 {tp / (tp + fp):.1%}, 경보 {int((al & ok).sum())}건")
gray = ok & (meta["deg"] > 0).values & (meta["deg"] < 0.3).values
print(f" 참고: 열화 초기(0~0.3) 구간 {int(gray.sum())}건 중 경보 {al[gray].mean():.0%}(회색 구간은 정확도에서 제외)")
# 3) 상태 화면 구성: 최근 30일 x 24대 요약(개체별 최댓값·경보 여부) 집계 시간
t0 = time.perf_counter()
last = va[va["day"] >= va.groupby("unit")["day"].transform("max") - 30]
sc = score_batch(Wva[last.index.values], last.reset_index(drop=True))
view = last.assign(score=sc).groupby("unit").agg(max_score=("score", "max"), last_score=("score", "last"))
view["상태"] = np.where(view["last_score"] > 1, "경보", np.where(view["max_score"] > 1, "주의", "정상"))
tv = time.perf_counter() - t0
print(f"상태 화면 구성(24대 x 최근 30일 {len(last)}건) {tv:.2f}초 (목표 2초 이내: {'달성' if tv <= 2 else '미달'}), " + ", ".join(f"{k} {v}대" for k, v in view["상태"].value_counts().items()))
# 4) 규모 확대: 배치 처리량과 메모리 한계 대응(배치 크기 5,000)
print("\n[주입 규모별 처리] 건수 | 시간 | 처리량 | 격리 | 유실")
sizes = (1000, 10000, 50000, 100000)
ts = []
for n in sizes:
w, meta, kind = make_stream(n, 0.06, seed=n)
t0 = time.perf_counter(); done = quarantine = 0
for i in range(0, n, 5000):
wb, mb = w[i:i + 5000], meta.iloc[i:i + 5000].reset_index(drop=True)
okb = valid(wb, mb)
if okb.any():
score_batch(wb[okb], mb[okb].reset_index(drop=True))
done += int(okb.sum()); quarantine += int((~okb).sum())
t = time.perf_counter() - t0; ts.append(t)
print(f" {n:7d} {t:5.2f}초 {n / t:8.0f}건/초 {quarantine:5d} {n - done - quarantine}")
print("운용 기준: 전체 격리율이 20%를 넘으면 센서 계통 점검 알림, 처리 지연 시 FFT 대역 5개만 계산하는 축소 모드로 전환")
plt.figure(figsize=(6, 4))
plt.plot(sizes, ts, marker="o"); plt.xscale("log"); plt.xlabel("stream size"); plt.ylabel("seconds"); plt.title("Stream stress test")
plt.tight_layout(); plt.show()
In [ ]:
"""AAR 기간의 예지정비 오경보를 원인(고부하 운용일·외기 급변·센서 스파이크·기타)별로 분석하고, 정상 분위수·연속 일수·고부하일 임계 완화 파라미터를 AAR 기간 자료로만 격자 탐색해 재보정한 뒤 검증 기간에서 오경보·고장 전 확보 시간 변화와 비용 절감·점검 시간 단축·확보 시간 증가의 3대 정량 효과 지표를 가정 단가로 산출한다."""
import itertools
import numpy as np
import pandas as pd
from scipy.stats import kurtosis
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# ---------------------------------------------------------------------------
# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)
# ---------------------------------------------------------------------------
BASE = ["rms", "kurt", "crest", "p2p", "b50", "b100", "b150", "b165", "b200", "temp", "press"]
ROLL = ["rms", "kurt", "b150", "b165", "b200", "temp", "press"]
def build_features(snap, W):
D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)
D["temp_res"] = D["temp"] - (55 + 22 * D["load"] + 0.5 * D["ambient"]) # 부하·외기 대비 온도 잔차
D["press_res"] = D["press"] - (450 - 12 * D["load"]) # 부하 대비 압력 잔차
D["rms_n"] = D["rms"] / (0.5 + D["load"]) # 부하 정규화 RMS
D["hf_ratio"] = (D["b150"] + D["b165"] + D["b200"]) / (D["b50"] + 1e-6) # 고주파/기본주파 비
g = D.groupby("unit")
for c in ROLL:
D[c + "_m3"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())
for c in ("rms_n", "temp_res", "press_res"):
D[c + "_d"] = g[c].diff().fillna(0)
return D
DERIVED = ["temp_res", "press_res", "rms_n", "hf_ratio"] + [c + "_m3" for c in ROLL] + ["rms_n_d", "temp_res_d", "press_res_d"]
FEATS = BASE + DERIVED
AE_COLS = ["rms_n", "kurt", "b150", "b165", "b200", "hf_ratio", "temp_res", "press_res"]
def split_units(units, n_train=24):
return list(units["unit"][:n_train]), list(units["unit"][n_train:])
class AE:
"""오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다."""
def __init__(self, cols=AE_COLS, q=0.99, seed=0):
self.cols, self.q, self.seed = cols, q, seed
def fit(self, healthy):
self.sc = StandardScaler().fit(healthy[self.cols])
Z = self.sc.transform(healthy[self.cols])
self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)
self.thr = float(np.quantile(self.err(healthy), self.q))
return self
def contrib(self, d):
Z = self.sc.transform(d[self.cols])
return (Z - self.net.predict(Z)) ** 2 # 입력 항목별 재구성 오차(이상 기여도)
def err(self, d):
return self.contrib(d).mean(axis=1)
def latch(d, e, thr, k=2):
"""같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다."""
out, e = np.zeros(len(d), bool), np.asarray(e)
for idx in d.reset_index(drop=True).groupby("unit").indices.values():
run, on = 0, False
for i in idx:
run = run + 1 if e[i] > thr else 0
on = on or run >= k
out[i] = on
return out
FA_COST, FA_HOURS, LATE_COST = 1.0, 3.0, 0.15 # [교육용가상자료] 오경보 1회 점검 비용지수, 점검 소요시간(h), 확보 시간 1일 손실 비용지수
units, snap, W = make_runs(seed=5)
D = build_features(snap, W)
order = list(units["unit"])
tr_u, aar_u, val_u = order[:24], order[24:36], order[36:]
tr = D[D["unit"].isin(tr_u)]
ae = AE().fit(tr[tr["deg"] == 0])
healthy_err = ae.err(tr[tr["deg"] == 0])
E = {u: (D[D["unit"] == u].reset_index(drop=True), ae.err(D[D["unit"] == u])) for u in order}
def evaluate(us, q, k, relax):
"""반환: 오경보 에피소드 수, 원인 분류 건수, 고장 전 확보 시간(일) 목록, 미탐 개체 수, 손실 확보 일수"""
base_thr = float(np.quantile(healthy_err, q))
fa, causes, leads, miss, lost = 0, {"고부하 운용일": 0, "외기 급변·고온": 0, "센서 스파이크": 0, "열화 초기 경계": 0, "기타": 0}, [], 0, 0.0
for u in us:
g, e = E[u]
thr = base_thr * np.where(g["load"].values > 0.95, 1 + relax, 1.0)
over = e > thr
run = np.zeros(len(g), int)
for i in range(len(g)):
run[i] = run[i - 1] + 1 if over[i] and i else int(over[i])
al = run >= k
an = (g["deg"] >= 0.1).values
start = al & ~np.r_[False, al[:-1]]
for i in np.where(start & ~an)[0]:
fa += 1
r = g.iloc[i]
c = "고부하 운용일" if r["hi_load"] else ("외기 급변·고온" if r["ambient"] > 24 else ("센서 스파이크" if r["kurt"] > 3.2 else ("열화 초기 경계" if r["deg"] > 0 else "기타")))
causes[c] += 1
first_an = np.argmax(an)
hit = np.where(al & an)[0]
ideal = g["rul"].iloc[first_an]
if len(hit):
leads.append(int(g["rul"].iloc[hit[0]])); lost += ideal - g["rul"].iloc[hit[0]]
else:
miss += 1; lost += ideal
return fa, causes, leads, miss, lost
def J(res):
return FA_COST * res[0] + LATE_COST * res[4]
cur = dict(q=0.99, k=1, relax=0.0) # 현행(AAR 이전) 파라미터
r_aar = evaluate(aar_u, **cur)
print(f"[교육용가상자료] AAR 기간 {len(aar_u)}대: 현행 파라미터(정상 {cur['q']:.0%} 분위, {cur['k']}일, 고부하 완화 {cur['relax']:.0%}) 오경보 {r_aar[0]}건")
tot = sum(r_aar[1].values())
for c, v in sorted(r_aar[1].items(), key=lambda x: -x[1]):
print(f" {c:10s} {v:3d}건 ({v / max(tot, 1):.0%})")
# 재보정: AAR 기간 자료만으로 격자 탐색 (미탐 개체가 없는 조합 중 비용 최소)
grid = []
for q, k, rl in itertools.product((0.98, 0.99, 0.995, 0.999), (1, 2, 3), (0.0, 0.3, 0.6)):
r = evaluate(aar_u, q, k, rl)
if r[3] == 0:
grid.append((J(r), q, k, rl))
grid.sort()
best = dict(q=grid[0][1], k=grid[0][2], relax=grid[0][3])
print(f"\n재보정 결과(격자 {4 * 3 * 3}개 중 미탐 0개 조합 {len(grid)}개 탐색): 정상 {best['q']:.1%} 분위, 연속 {best['k']}일, 고부하일 임계 {best['relax']:.0%} 완화")
rows = []
for nm, us in (("AAR 기간(튜닝)", aar_u), ("검증 기간(미사용)", val_u)):
a, b = evaluate(us, **cur), evaluate(us, **best)
rows.append((nm, a, b))
print(f" {nm:14s} 오경보 {a[0]:3d} -> {b[0]:3d}건, 고장 전 확보 중앙값 {np.median(a[2]):.0f} -> {np.median(b[2]):.0f}일, 미탐 {a[3]} -> {b[3]}대")
# 3대 정량 효과 지표(검증 3개월 대신 검증 12대 기준, 연환산은 개체 수 비례)
a, b = rows[1][1], rows[1][2]
jc, jn = J(a), J(b)
print("\n[3대 정량 효과 지표] (검증 12대, 가정 단가)")
print(f" 1) 비용 절감: 비용지수 {jc:.1f} -> {jn:.1f} ({1 - jn / jc:.1%} 절감), 48대 연환산 약 {(jc - jn) * 4:.0f}")
print(f" 2) 시간 단축: 불필요 점검 {a[0] * FA_HOURS:.0f}h -> {b[0] * FA_HOURS:.0f}h ({1 - b[0] / max(a[0], 1):.1%} 단축)")
print(f" 3) 확보 시간: 고장 전 첫 경보 확보 중앙값 {np.median(a[2]):.0f}일 -> {np.median(b[2]):.0f}일, 정상 오경보 {a[0]}건 -> {b[0]}건")
print("[도입 로드맵] 단기(1~3월): 재보정 임계 적용 + 검사관 확인 병행 / 중기(3~12월): 분기 AAR로 임계 재보정, 부하 보정 항 추가 / 장기(1년~): 군수 정보체계 정비 이력과 연동해 임계 자동 갱신")
In [ ]:
"""센서 파형 검증·정제, 오토인코더 이상 감지 모델의 자동 재학습(챔피언·도전자 비교 후 교체), 재보정 임계 적용 경보, 잔존수명(RUL) 예측, 검사관 승인 대기 입고 제안, 수리부속 소요 예고까지 예지정비(PdM) 이상 감지 시스템을 한 파일로 통합 실행하고 자체 점검·운영 문서 완비 점검·정비 담당 장교 브리핑 요약을 출력한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import kurtosis
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)
# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.
# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)
# ---------------------------------------------------------------------------
FS, NS = 1024, 256 # 가상 표본화 주파수(Hz), 파형 1개당 표본 수
TT = np.arange(NS) / FS
MODES = { # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)
"bearing": (0.40, 160, "주 베어링 세트", 6.0), "gear": (0.25, 175, "구동 기어 세트", 9.0),
"oil_pump": (0.20, 150, "오일 펌프", 4.0), "cooling": (0.15, 140, "냉각 팬 모터", 3.0)}
KIT_COST = 1.0 # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수
BANDS = {"b50": (40, 60), "b100": (90, 110), "b150": (140, 156), "b165": (156, 176), "b200": (190, 215)}
def make_runs(n=48, seed=1):
"""반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])"""
rng = np.random.default_rng(seed)
names = list(MODES)
mode = rng.choice(names, n, p=[MODES[m][0] for m in names])
life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)
start = np.sort(rng.integers(0, 240, n))
units = pd.DataFrame({"unit": [f"U{i:02d}" for i in range(n)], "mode": mode, "life": life, "start": start, "base_load": rng.uniform(0.4, 0.75, n)})
rows = []
for u in units.itertuples():
d = np.arange(u.life + 1)
onset = 0.45 * u.life
deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5
load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)
hi = rng.random(len(d)) < 0.07 # 고부하 운용일(정상이지만 진동·온도 상승)
load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)
rows.append(pd.DataFrame({"unit": u.unit, "day": d, "cal": u.start + d, "rul": u.life - d, "deg": deg, "load": load, "hi_load": hi, "mode": u.mode}))
snap = pd.concat(rows, ignore_index=True)
S = len(snap)
snap["ambient"] = 15 + 12 * np.sin(2 * np.pi * snap["cal"] / 365) + rng.normal(0, 2, S)
dg, ld, md = snap["deg"].values, snap["load"].values, snap["mode"].values
isb, isg, iso, isc = (md == "bearing"), (md == "gear"), (md == "oil_pump"), (md == "cooling")
snap["temp"] = 55 + 22 * ld + 0.5 * snap["ambient"] + 9 * dg * iso + 26 * dg * isc + 3 * dg * (isb | isg) + rng.normal(0, 1.5, S)
snap["press"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)
ph = rng.uniform(0, 2 * np.pi, (S, 3))
a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)
W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])
+ rng.normal(0, 0.35, (S, NS)))
W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])
W += (0.8 * dg * isg)[:, None] * (np.sin(2 * np.pi * 150 * TT + ph[:, 2:3]) + 0.8 * np.sin(2 * np.pi * 100 * TT))
W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])
k = rng.poisson(10 * dg * isb) # 베어링 열화: 충격성 펄스(첨도 상승)
pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]
amp = amp * (np.arange(30)[None, :] < k[:, None])
np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())
return units, snap, W
def extract(W):
"""진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)"""
x = W - W.mean(axis=1, keepdims=True)
rms = np.sqrt((x ** 2).mean(axis=1))
out = {"rms": rms, "kurt": kurtosis(x, axis=1, fisher=False), "crest": np.abs(x).max(axis=1) / rms, "p2p": x.max(axis=1) - x.min(axis=1)}
A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS
f = np.fft.rfftfreq(NS, 1 / FS)
for nm, (lo, hi) in BANDS.items():
out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))
return pd.DataFrame(out)
# ---------------------------------------------------------------------------
# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)
# ---------------------------------------------------------------------------
BASE = ["rms", "kurt", "crest", "p2p", "b50", "b100", "b150", "b165", "b200", "temp", "press"]
ROLL = ["rms", "kurt", "b150", "b165", "b200", "temp", "press"]
def build_features(snap, W):
D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)
D["temp_res"] = D["temp"] - (55 + 22 * D["load"] + 0.5 * D["ambient"]) # 부하·외기 대비 온도 잔차
D["press_res"] = D["press"] - (450 - 12 * D["load"]) # 부하 대비 압력 잔차
D["rms_n"] = D["rms"] / (0.5 + D["load"]) # 부하 정규화 RMS
D["hf_ratio"] = (D["b150"] + D["b165"] + D["b200"]) / (D["b50"] + 1e-6) # 고주파/기본주파 비
g = D.groupby("unit")
for c in ROLL:
D[c + "_m3"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())
for c in ("rms_n", "temp_res", "press_res"):
D[c + "_d"] = g[c].diff().fillna(0)
return D
DERIVED = ["temp_res", "press_res", "rms_n", "hf_ratio"] + [c + "_m3" for c in ROLL] + ["rms_n_d", "temp_res_d", "press_res_d"]
FEATS = BASE + DERIVED
AE_COLS = ["rms_n", "kurt", "b150", "b165", "b200", "hf_ratio", "temp_res", "press_res"]
def split_units(units, n_train=24):
return list(units["unit"][:n_train]), list(units["unit"][n_train:])
class AE:
"""오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다."""
def __init__(self, cols=AE_COLS, q=0.99, seed=0):
self.cols, self.q, self.seed = cols, q, seed
def fit(self, healthy):
self.sc = StandardScaler().fit(healthy[self.cols])
Z = self.sc.transform(healthy[self.cols])
self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)
self.thr = float(np.quantile(self.err(healthy), self.q))
return self
def contrib(self, d):
Z = self.sc.transform(d[self.cols])
return (Z - self.net.predict(Z)) ** 2 # 입력 항목별 재구성 오차(이상 기여도)
def err(self, d):
return self.contrib(d).mean(axis=1)
def latch(d, e, thr, k=2):
"""같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다."""
out, e = np.zeros(len(d), bool), np.asarray(e)
for idx in d.reset_index(drop=True).groupby("unit").indices.values():
run, on = 0, False
for i in idx:
run = run + 1 if e[i] > thr else 0
on = on or run >= k
out[i] = on
return out
CONFIG = {"q": 0.995, "k": 1, "relax": 0.3, "retrain_tolerance": 0.02, "approve_rul": 90, "spare_horizon": 60} # 미션 9 재보정 결과
OPS_GUIDE = { # docs/mnt04_operations_guide.md 요약
"1. 목적과 범위": "엔진·구동계 센서로 이상 징후와 잔존수명을 보조 판단. 입고·정비 결정은 검사관.",
"2. 입력 데이터": "일 1회 진동 파형·온도·압력·부하. 실제 장비 자료 대신 가상 식별코드 사용.",
"3. 일일 운용 절차": "수집 -> 검증·격리 -> 특징 추출 -> 이상 점수 -> 경보 -> 입고 제안(대기) -> 검사관 결정 -> 결과 기록.",
"4. 승인(HITL) 기준": "조기 입고 제안은 검사관 승인 후 확정, 무응답 시 기본 계획(정기 교체) 유지, 사유 기록.",
"5. 재학습 절차": "분기 1회 신규 정상·고장 실적으로 도전자 학습, 검증 F1이 나빠지지 않을 때만 교체.",
"6. 장애 대응": "모델 오류 시 정기 교체 주기 기준으로 대체, 손상 센서 자료는 수동 확인 목록으로.",
"7. 연동 인터페이스": "국방 군수 정보체계와 정비 이력·수리부속 소모 표준 필드(장비코드·분류군·일시)를 CSV/API로 교환하는 안을 차기 과제로 둔다.",
}
F10 = FEATS + ["day", "hf_max", "temp_res_max"]
def prep(units, snap, W):
"""검증: 파형 결측·고정·포화, 값 범위 이탈 측정은 격리하고 나머지로 특징표를 만든다."""
bad = np.isnan(W).any(axis=1) | (np.nanstd(W, axis=1) < 1e-6) | (np.abs(np.nan_to_num(W)).max(axis=1) > 12) | ~snap["temp"].between(20, 160).values
ok = ~bad
D = build_features(snap[ok].reset_index(drop=True), W[ok])
D["hf_max"] = D.groupby("unit")["hf_ratio"].cummax(); D["temp_res_max"] = D.groupby("unit")["temp_res"].cummax()
return D, int(bad.sum())
def alarms(D, e, thr, cfg=CONFIG):
over = np.asarray(e) > thr * np.where(D["load"].values > 0.95, 1 + cfg["relax"], 1.0)
out = np.zeros(len(D), bool)
for idx in D.reset_index(drop=True).groupby("unit").indices.values():
run = 0
for i in idx:
run = run + 1 if over[i] else 0
out[i] = run >= cfg["k"]
return out
def f1(D, al):
y = (D["deg"] >= 0.1).values
tp, fp, fn = (al & y).sum(), (al & ~y).sum(), (~al & y).sum()
return float(2 * tp / (2 * tp + fp + fn))
def retrain(champion, D_old, D_new):
"""챔피언 AE(지난달 자료) vs 도전자 AE(지난달+이번 달 학습 개체) -> 이번 달 검증 개체 F1 비교."""
new_tr, new_va = D_new[D_new["unit"].isin(list(D_new["unit"].unique())[:24])], D_new[D_new["unit"].isin(list(D_new["unit"].unique())[24:])]
ch = AE(seed=1).fit(pd.concat([D_old[D_old["deg"] == 0], new_tr[new_tr["deg"] == 0]]))
thr_c = float(np.quantile(champion.err(D_old[D_old["deg"] == 0]), CONFIG["q"]))
thr_n = float(np.quantile(ch.err(pd.concat([D_old[D_old["deg"] == 0], new_tr[new_tr["deg"] == 0]])), CONFIG["q"]))
f_old, f_new = f1(new_va, alarms(new_va, champion.err(new_va), thr_c)), f1(new_va, alarms(new_va, ch.err(new_va), thr_n))
swap = f_new >= f_old - CONFIG["retrain_tolerance"]
return (ch if swap else champion), (thr_n if swap else thr_c), {"챔피언 F1": round(f_old, 3), "도전자 F1": round(f_new, 3), "교체": bool(swap)}
def predict_safe(model, D):
try:
return np.maximum(model.predict(D[F10]), 0), "모델"
except Exception:
return np.maximum(90 - D["day"].values, 0).astype(float), "정기 주기 기준 대체"
# ---------------- 1. 통합 실행 ----------------
t_start = time.perf_counter()
Dold, _ = prep(*make_runs(seed=11)) # 지난달 실적
Dnew, _ = prep(*make_runs(seed=1)) # 이번 달 실적
u3, s3, w3 = make_runs(seed=21) # 다음 달 운용(신규 측정)
w3[[100, 2000], 20:50] = np.nan # 손상 측정 2건 주입(결측 패킷)
Dop, n_bad = prep(u3, s3, w3)
champion = AE().fit(Dold[Dold["deg"] == 0])
champ_thr = float(np.quantile(champion.err(Dold[Dold["deg"] == 0]), CONFIG["q"]))
model_ae, thr, rt = retrain(champion, Dold, Dnew)
print(f"[교육용가상자료] 재학습: {rt} -> {'도전자 교체' if rt['교체'] else '챔피언 유지'}")
gbr = GradientBoostingRegressor(n_estimators=100, max_depth=3, learning_rate=0.1, subsample=0.8, random_state=0).fit(Dnew[F10], np.minimum(Dnew["rul"], 120))
e = model_ae.err(Dop)
Dop["alarm"] = alarms(Dop, e, thr)
Dop["rul_hat"], src = predict_safe(gbr, Dop)
Dop["score"] = e / thr
sec = time.perf_counter() - t_start
print(f"운용 측정 {len(s3)}건: 처리 {len(Dop)}건, 격리 {n_bad}건, 경보 측정 {int(Dop['alarm'].sum())}건, RUL 출처 '{src}', 전체 실행 {sec:.1f}초")
# 입고 제안(대기)과 수리부속 소요 예고: 경보 개체의 마지막 측정 기준
PART = {"b165": "주 베어링 세트", "kurt": "주 베어링 세트", "b150": "구동 기어 세트", "b200": "오일 펌프", "press_res": "오일 펌프", "temp_res": "냉각 팬 모터"}
first = Dop[Dop["alarm"]].groupby("unit").first()
cols = model_ae.cols
C = model_ae.contrib(Dop)
prop = []
for u, r in first.iterrows():
i = Dop.index[(Dop["unit"] == u) & (Dop["day"] == r["day"])][0]
prop.append({"unit": u, "day": int(r["day"]), "rul_hat": float(r["rul_hat"]), "rul": int(r["rul"]), "part": PART.get(cols[int(np.argmax(C[i]))], "미분류"),
"mode": r["mode"], "status": "PENDING" if r["rul_hat"] <= CONFIG["approve_rul"] else "관찰"})
P = pd.DataFrame(prop)
P["true_part"] = P["mode"].map({m: MODES[m][2] for m in MODES})
# 수리부속 소요 예고: 예측 잔존 30일 이내로 내려온 개체의 의심 부품
need = P[P["rul_hat"] <= CONFIG["spare_horizon"]]
# ---------------- 2. 자체 점검 ----------------
_, bad_src = predict_safe(None, Dop.head(3))
y = (Dop["deg"] >= 0.1).values
checks = {
"측정 유실 0건": len(Dop) + n_bad == len(s3),
"손상 측정 2건 격리": n_bad == 2,
"재학습 파이프라인 실행": "교체" in rt,
"RUL 예측값 음수 없음": bool((Dop["rul_hat"] >= 0).all()),
"모델 장애 시 정기 주기 대체": bad_src == "정기 주기 기준 대체",
"열화 개체 미탐 0대": bool(Dop.groupby("unit").apply(lambda g: (g["alarm"] & (g["deg"] >= 0.1)).any() or not (g["deg"] >= 0.1).any()).all()),
"경보 개체 모두 제안 또는 관찰 분류": len(P) == int(Dop[Dop["alarm"]]["unit"].nunique()),
"자동 승인 0건(승인 대상은 대기)": bool(P["status"].isin(["PENDING", "관찰"]).all()),
}
print(f"\n자체 점검 {sum(checks.values())}/{len(checks)} 통과, 자동 재학습 파이프라인 {'정상 실행' if checks['재학습 파이프라인 실행'] else '실패'}")
for k, v in checks.items():
if not v:
print(" 실패:", k)
doc = sum(bool(v.strip()) for v in OPS_GUIDE.values()) / len(OPS_GUIDE)
print(f"운영 문서 {len(OPS_GUIDE)}개 장 작성률 {doc:.0%} (목표 100%: {'달성' if doc == 1 else '미달'}) - " + ", ".join(k.split(". ")[1] for k in OPS_GUIDE))
# ---------------- 3. 화면과 브리핑 ----------------
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
for u in list(P["unit"])[:3]:
g = Dop[Dop["unit"] == u]
ax[0].plot(g["day"], g["score"], label=u)
ax[0].axhline(1.0, ls="--", color="k"); ax[0].set_yscale("log"); ax[0].set_title("Anomaly score (x threshold)"); ax[0].set_xlabel("day"); ax[0].legend()
ax[1].scatter(P["rul"], P["rul_hat"]); ax[1].plot([0, 120], [0, 120], "k--"); ax[1].set_title("RUL at first alarm: actual vs predicted"); ax[1].set_xlabel("actual (days)")
plt.tight_layout(); plt.show()
lead = P["rul"]
print("\n[정비 담당 장교 브리핑 요약] [교육용가상자료]")
print(f"1. 운용 개체 {len(u3)}대 측정 {len(s3)}건 중 손상 {n_bad}건 격리, 열화 개체 {len(P)}대에서 경보 발령, 고장 전 확보 시간 중앙값 {lead.median():.0f}일(최소 {lead.min()}일)")
print(f"2. 이상 감지 모델: {'신규 실적으로 교체' if rt['교체'] else '기존 유지'}(검증 F1 {max(rt['챔피언 F1'], rt['도전자 F1']):.2f}), 재보정 임계 적용으로 정상 개체 오경보 {int((Dop['alarm'] & (Dop['deg'] == 0)).sum())}건")
print(f"3. 조기 입고 제안 {int((P['status'] == 'PENDING').sum())}건은 검사관 승인 대기(자동 승인 0건), 예측 잔존 {CONFIG['approve_rul']}일 초과 {int((P['status'] == '관찰').sum())}건은 관찰 유지")
print("4. 수리부속 소요 예고(예측 잔존 " + str(CONFIG["spare_horizon"]) + "일 이내 개체의 의심 부품): " + ", ".join(f"{k} {v}건" for k, v in need["part"].value_counts().items())
+ f" (첫 경보 시점 의심 부품 정확도 {(P['part'] == P['true_part']).mean():.0%})")
print(f"5. RUL 예측: 경보 시점 평균 절대오차 {np.mean(np.abs(P['rul_hat'] - np.minimum(P['rul'], 120))):.1f}일")
print("6. 건의: 구동 기어·주 베어링 세트 우선 확보(소모 비용 상위), 국방 군수 정보체계와 정비 이력·수리부속 소모 표준 필드 연동 인터페이스 정의")