{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# 🎯 [MNT04] 고장실적·수리부속 소모 기반 PBL 품목관리\n",
        "### Google Colab 실습 환경 (장비정비·수리부속 / Level 3)\n",
        "\n",
        "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/edunlab/milicontent2/blob/main/uploads/MNT04/MNT04_solution.ipynb)\n",
        "\n",
        "---\n",
        "\n",
        "## 📌 1. 프로젝트 개요 및 배경\n",
        "\n",
        "| 항목 | 내용 |\n",
        "| :--- | :--- |\n",
        "| **과정 분류** | 장비정비·수리부속 (Level 3) |\n",
        "| **학습 소요 시간** | 10개 미션 (5일 / 총 28.5시간) |\n",
        "| **환경 사양** | Python 3.10+, numpy, pandas, matplotlib, scipy, scikit-learn |\n",
        "\n",
        "### 📖 프로젝트 목표\n",
        "항공·기갑 장비의 엔진과 구동계는 고장이 나기 전에 진동, 온도, 압력에 미세한 변화가 먼저 나타나지만, 고정 임계값만으로는 고부하 운용일을 이상으로 오인하거나 초기 열화를 놓쳐 수리부속을 제때 확보하지 못한다. 본 프로젝트는 고장 실적과 수리부속 소모를 연결해 정비와 부품 관리의 연계 문제를 학습하는 것을 목표로, 교육용 가상 기동체 48대의 센서 시계열과 부품 소모를 코드로 만들어 분석한다. 센서 데이터셋 정의, 잡음 제거와 FFT 특징량 추출, 임계값 Baseline과 부품 소모 파레토 분석, 오토인코더 기반 미세 변화 감지와 파생 피처 모델, 잔존수명(RUL) 예측과 이상 신뢰도 점수, 설명가능한 경보 화면, 검사관 승인(HITL) 기반 조기 입고 시뮬레이터, 대용량 스트림 스트레스 시험, AAR 오경보 분석과 임계 재보정, 재학습을 포함한 통합 도구까지 단계별로 구현한다. 최종 목표는 고장 전에 이상 징후를 일찍 잡아 정비 시점과 수리부속 소요를 예고하되, 조기 입고 결정은 검사관이 근거를 보고 내리도록 지원하는 것이다.\n",
        "\n",
        "> **⚠️ [교육용 가상 데이터셋 원칙]**\n",
        "> - 본 노트북의 모든 데이터(제원, 로그, 시나리오 등)는 보안 규정을 준수하여 **코드 내 난수 기반 가상 데이터로 생성**됩니다.\n",
        "> - 실제 부대 위치, 작전 인가, 실제 전장 정보, 군사 기밀 데이터는 일체 포함되지 않습니다.\n",
        "\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# ============================================================================\n",
        "# [Google Colab 환경 설정] 필수 라이브러리 설치 및 임포트\n",
        "# ============================================================================\n",
        "!pip install -q numpy pandas matplotlib scipy scikit-learn\n",
        "\n",
        "import os\n",
        "import sys\n",
        "import time\n",
        "import math\n",
        "import json\n",
        "import re\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "import warnings\n",
        "warnings.filterwarnings('ignore')\n",
        "\n",
        "# 그래프 영문 폰트 기본 설정 (환경 독립성 보장)\n",
        "plt.rcParams['font.sans-serif'] = 'DejaVu Sans'\n",
        "plt.rcParams['axes.unicode_minus'] = False\n",
        "\n",
        "print(f\"실행 환경 준비 완료 (Python {sys.version.split()[0]})\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 01. 헬기/기갑 엔진 및 구동계 진동, 온도, 압력 센서 시계열 데이터셋 정의\n",
        "**소요시간**: 약 1.5시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "헬기/기갑 엔진 진동, 온도, 압력 센서 시계열 데이터셋을 정의한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **헬기/기갑 엔진 및 구동계의 진동, 온도, 압력 센서 시계열과 고장·수리부속 소모 실적의 데이터셋 스키마, 데이터 딕셔너리, 업무 프로세스 맵을 정의한다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **센서 수집부터 전처리, 특징 추출, 이상 감지, RUL 예측, 검사관 승인, 정비·부품 소모, 사후검토까지 업무 흐름을 단계별로 정리하고, 센서·특징·고장·부품 묶음으로 피처 20개 이상 정의하고 검증하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 4개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `docs/mnt04_logistics_domain_spec.md 및 데이터 딕셔너리` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 군수 실무 업무 프로세스와의 일치성 및 분석 목표의 명확성\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: 업무 프로세스 맵 작성 완료, 필요 피처(Feature) 목록 20개 이상 정의\n",
        "\n",
        "\n",
        "- 진행 기준: 도메인 분석서 및 데이터 요구명세서 승인\n",
        "\n",
        "\n",
        "- 검증 방법: 군수 담당 실무자 인터뷰 및 기존 전산체계 데이터 구조 대조\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"헬기/기갑 엔진 및 구동계의 진동·온도·압력 센서 시계열, 고장 실적, 수리부속 소모 기록의 데이터셋 스키마·데이터 딕셔너리(피처 20개 이상)와 예지정비·PBL 품목관리 업무 프로세스 맵을 정의하고, 가상 센서 자료로 검증한다.\"\"\"\n",
        "import json\n",
        "import re\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "from scipy.stats import kurtosis\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# 1) 업무 프로세스 맵: 센서 수집부터 수리부속 소요 환류까지 (각 단계에서 생기는 데이터)\n",
        "PROCESS = [(\"센서 수집\", \"개체 ID, 측정일, 진동 파형, 온도, 압력, 부하\"), (\"전처리\", \"잡음 제거, 단위·표기 통일, 결측 보정\"),\n",
        "           (\"특징 추출\", \"RMS, 첨도, 파고율, FFT 대역 에너지\"), (\"이상 감지\", \"이상 점수, 경보 단계\"),\n",
        "           (\"잔존수명(RUL) 예측\", \"고장 임박 시점, 신뢰 구간\"), (\"검사관 승인\", \"조기 입고 제안, 승인·반려 사유\"),\n",
        "           (\"정비·부품 소모\", \"교체 부품, 수량, 정비 소요시간\"), (\"고장 실적 연결\", \"정비 건 ID로 센서 이력과 고장 실적 연결\"),\n",
        "           (\"부품 소요·가동률 집계\", \"품목별 소모, 가용도(PBL 지표)\"), (\"사후검토(AAR)\", \"오경보 원인, 임계 파라미터 조정\")]\n",
        "\n",
        "# 2) 스키마(데이터 딕셔너리): 그룹 / 필드 / 자료형 / 규칙 / 설명\n",
        "S = {\n",
        "    \"sensor\": {\n",
        "        \"unit_id\": (str, {\"pattern\": r\"^U\\d{2}$\"}, \"개체 ID(가상 기동체)\"),\n",
        "        \"day\": (int, {\"min\": 0, \"max\": 260}, \"운용 경과일(측정일)\"),\n",
        "        \"vib_rms\": (float, {\"min\": 0, \"max\": 10}, \"진동 RMS(가상 단위)\"),\n",
        "        \"vib_kurt\": (float, {\"min\": 1, \"max\": 30}, \"진동 첨도(충격성 지표)\"),\n",
        "        \"temp\": (float, {\"min\": -30, \"max\": 150}, \"엔진 온도(섭씨)\"),\n",
        "        \"press\": (float, {\"min\": 0, \"max\": 700}, \"오일 압력(kPa)\"),\n",
        "        \"load\": (float, {\"min\": 0, \"max\": 1.2}, \"운용 부하율\"),\n",
        "    },\n",
        "    \"feature\": {\n",
        "        \"band_50\": (float, {\"min\": 0}, \"50Hz 부근 FFT 대역 에너지\"),\n",
        "        \"band_165\": (float, {\"min\": 0}, \"165Hz 부근(베어링 결함 대역) 에너지\"),\n",
        "        \"crest\": (float, {\"min\": 1, \"max\": 30}, \"파고율(최대값/RMS)\"),\n",
        "        \"hf_ratio\": (float, {\"min\": 0}, \"고주파/기본주파 에너지 비\"),\n",
        "        \"temp_res\": (float, {\"min\": -60, \"max\": 80}, \"부하·외기 대비 온도 잔차\"),\n",
        "    },\n",
        "    \"failure\": {\n",
        "        \"event_id\": (str, {\"pattern\": r\"^F\\d{3}$\"}, \"고장·정비 건 ID(센서 이력과 연결하는 키)\"),\n",
        "        \"mode\": (str, {\"enum\": list(MODES)}, \"열화(고장) 유형\"),\n",
        "        \"life\": (int, {\"min\": 60, \"max\": 260}, \"고장까지 운용일수\"),\n",
        "        \"unplanned\": (bool, {}, \"비계획 정비 여부\"),\n",
        "    },\n",
        "    \"parts\": {\n",
        "        \"part_name\": (str, {\"enum\": [v[2] for v in MODES.values()]}, \"소모 수리부속\"),\n",
        "        \"qty\": (int, {\"min\": 1, \"max\": 10}, \"소모 수량\"),\n",
        "        \"unit_cost\": (float, {\"min\": 0, \"max\": 50}, \"가상 단가지수\"),\n",
        "    },\n",
        "    \"target\": {\n",
        "        \"rul\": (float, {\"min\": 0, \"max\": 260}, \"잔존수명(일) - 예측대상 1\"),\n",
        "        \"anomaly_score\": (float, {\"min\": 0}, \"이상 점수(재구성 오차) - 예측대상 2\"),\n",
        "        \"availability\": (float, {\"min\": 0, \"max\": 1}, \"가용도(가동일/전체일) - PBL 성과 지표\"),\n",
        "    },\n",
        "}\n",
        "\n",
        "\n",
        "def check(name, v, typ, rule):\n",
        "    if not isinstance(v, typ) or (typ in (int, float) and isinstance(v, bool)):\n",
        "        return [f\"{name}: 자료형 오류({type(v).__name__})\"]\n",
        "    e = []\n",
        "    if \"enum\" in rule and v not in rule[\"enum\"]:\n",
        "        e.append(f\"{name}: 허용되지 않은 값 '{v}'\")\n",
        "    if \"min\" in rule and v < rule[\"min\"]:\n",
        "        e.append(f\"{name}: 최소값 {rule['min']} 미만\")\n",
        "    if \"max\" in rule and v > rule[\"max\"]:\n",
        "        e.append(f\"{name}: 최대값 {rule['max']} 초과\")\n",
        "    if \"pattern\" in rule and not re.match(rule[\"pattern\"], v):\n",
        "        e.append(f\"{name}: 형식 불일치 '{v}'\")\n",
        "    return e\n",
        "\n",
        "\n",
        "def validate(rec):\n",
        "    errs = []\n",
        "    for g, fields in S.items():\n",
        "        for k, (typ, rule, _) in fields.items():\n",
        "            if k not in rec.get(g, {}):\n",
        "                errs.append(f\"{g}.{k}: 필수 필드 누락\")\n",
        "            else:\n",
        "                errs += check(f\"{g}.{k}\", rec[g][k], typ, rule)\n",
        "    return errs\n",
        "\n",
        "\n",
        "# 3) [교육용가상자료] 생성 자료에서 20건을 레코드로 만들어 검증\n",
        "units, snap, W = make_runs()\n",
        "F = extract(W)\n",
        "snap = pd.concat([snap, F], axis=1)\n",
        "snap[\"temp_res\"] = snap[\"temp\"] - (55 + 22 * snap[\"load\"] + 0.5 * snap[\"ambient\"])\n",
        "snap[\"hf_ratio\"] = (snap[\"b150\"] + snap[\"b165\"] + snap[\"b200\"]) / (snap[\"b50\"] + 1e-6)\n",
        "ui = units.set_index(\"unit\")\n",
        "recs = []\n",
        "for i, r in enumerate(snap.sample(20, random_state=0).itertuples()):\n",
        "    m = r.mode\n",
        "    recs.append({\n",
        "        \"sensor\": {\"unit_id\": r.unit, \"day\": int(r.day), \"vib_rms\": float(round(r.rms, 3)), \"vib_kurt\": float(round(r.kurt, 2)),\n",
        "                   \"temp\": float(round(r.temp, 1)), \"press\": float(round(r.press, 1)), \"load\": float(round(r.load, 2))},\n",
        "        \"feature\": {\"band_50\": float(round(r.b50, 4)), \"band_165\": float(round(r.b165, 4)), \"crest\": float(round(r.crest, 2)),\n",
        "                    \"hf_ratio\": float(round(r.hf_ratio, 3)), \"temp_res\": float(round(r.temp_res, 2))},\n",
        "        \"failure\": {\"event_id\": f\"F{int(r.unit[1:]):03d}\", \"mode\": m, \"life\": int(ui.loc[r.unit, \"life\"]), \"unplanned\": bool(r.rul < 5)},\n",
        "        \"parts\": {\"part_name\": MODES[m][2], \"qty\": 1, \"unit_cost\": float(MODES[m][3])},\n",
        "        \"target\": {\"rul\": float(r.rul), \"anomaly_score\": float(round(r.deg, 3)), \"availability\": 0.92}})\n",
        "ok = sum(not validate(r) for r in recs)\n",
        "print(f\"가상 센서 자료 {len(snap)}건(개체 {len(units)}대), 표본 {len(recs)}건 중 검증 통과 {ok}건, 오류 {len(recs) - ok}건\")\n",
        "print(f\"전체 파형 {W.shape[0]}개 x {W.shape[1]}표본, 온도·압력 결측 {int(snap[['temp', 'press']].isna().sum().sum())}건\")\n",
        "\n",
        "# 4) 오류 주입 시험\n",
        "bad = json.loads(json.dumps(recs[0]))\n",
        "bad[\"sensor\"][\"unit_id\"] = \"ENGINE-3\"           # ID 형식 오류\n",
        "bad[\"sensor\"][\"temp\"] = 480.0                  # 센서 범위 초과(단위 혼동 의심)\n",
        "bad[\"failure\"][\"mode\"] = \"turbine\"             # 정의되지 않은 열화 유형\n",
        "del bad[\"parts\"][\"qty\"]                        # 필수 필드 누락\n",
        "bad[\"target\"][\"rul\"] = \"30일\"                  # 자료형 오류\n",
        "print(\"\\n[오류 주입 결과]\")\n",
        "for e in validate(bad):\n",
        "    print(\" -\", e)\n",
        "\n",
        "# 5) 프로세스 맵과 데이터 딕셔너리\n",
        "print(\"\\n[업무 프로세스 맵]\")\n",
        "print(\" -> \".join(f\"{i}.{p}\" for i, (p, _) in enumerate(PROCESS, 1)))\n",
        "for p, d in PROCESS:\n",
        "    print(f\"   {p:12s}: {d}\")\n",
        "rows = [(g, k, t.__name__, d) for g, fl in S.items() for k, (t, r, d) in fl.items()]\n",
        "dd = pd.DataFrame(rows, columns=[\"그룹\", \"필드\", \"자료형\", \"설명\"])\n",
        "print(f\"\\n[데이터 딕셔너리] 피처 {len(dd)}개 (목표 20개 이상: {'달성' if len(dd) >= 20 else '미달'})\")\n",
        "print(dd.to_string(index=False))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 01 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "가상 센서 자료 7733건(개체 48대), 표본 20건 중 검증 통과 20건, 오류 0건\n",
        "전체 파형 7733개 x 256표본, 온도·압력 결측 0건\n",
        "\n",
        "[오류 주입 결과]\n",
        " - sensor.unit_id: 형식 불일치 'ENGINE-3'\n",
        " - sensor.temp: 최대값 150 초과\n",
        " - failure.mode: 허용되지 않은 값 'turbine'\n",
        " - parts.qty: 필수 필드 누락\n",
        " - target.rul: 자료형 오류(str)\n",
        "\n",
        "[업무 프로세스 맵]\n",
        "1.센서 수집 -> 2.전처리 -> 3.특징 추출 -> 4.이상 감지 -> 5.잔존수명(RUL) 예측 -> 6.검사관 승인 -> 7.정비·부품 소모 -> 8.고장 실적 연결 -> 9.부품 소요·가동률 집계 -> 10.사후검토(AAR)\n",
        "   센서 수집       : 개체 ID, 측정일, 진동 파형, 온도, 압력, 부하\n",
        "   전처리         : 잡음 제거, 단위·표기 통일, 결측 보정\n",
        "   특징 추출       : RMS, 첨도, 파고율, FFT 대역 에너지\n",
        "   이상 감지       : 이상 점수, 경보 단계\n",
        "   잔존수명(RUL) 예측: 고장 임박 시점, 신뢰 구간\n",
        "   검사관 승인      : 조기 입고 제안, 승인·반려 사유\n",
        "   정비·부품 소모    : 교체 부품, 수량, 정비 소요시간\n",
        "   고장 실적 연결    : 정비 건 ID로 센서 이력과 고장 실적 연결\n",
        "   부품 소요·가동률 집계: 품목별 소모, 가용도(PBL 지표)\n",
        "   사후검토(AAR)   : 오경보 원인, 임계 파라미터 조정\n",
        "\n",
        "[데이터 딕셔너리] 피처 22개 (목표 20개 이상: 달성)\n",
        "     그룹            필드   자료형                        설명\n",
        " sensor       unit_id   str             개체 ID(가상 기동체)\n",
        " sensor           day   int               운용 경과일(측정일)\n",
        " sensor       vib_rms float             진동 RMS(가상 단위)\n",
        " sensor      vib_kurt float             진동 첨도(충격성 지표)\n",
        " sensor          temp float                 엔진 온도(섭씨)\n",
        " sensor         press float                오일 압력(kPa)\n",
        " sensor          load float                    운용 부하율\n",
        "feature       band_50 float        50Hz 부근 FFT 대역 에너지\n",
        "feature      band_165 float   165Hz 부근(베어링 결함 대역) 에너지\n",
        "feature         crest float              파고율(최대값/RMS)\n",
        "feature      hf_ratio float            고주파/기본주파 에너지 비\n",
        "feature      temp_res float            부하·외기 대비 온도 잔차\n",
        "failure      event_id   str 고장·정비 건 ID(센서 이력과 연결하는 키)\n",
        "failure          mode   str                 열화(고장) 유형\n",
        "failure          life   int                 고장까지 운용일수\n",
        "... (이하 생략)\n",
        "```\n",
        "</details>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 02. 고주파 센서 신호 잡음 제거, FFT 변환 및 주요 특징량(RMS, Kurtosis) 전처리 구현\n",
        "**소요시간**: 약 2.0시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "센서 잡음 제거, FFT 변환 및 주요 특징량(RMS 등) 추출 전처리 모듈을 구현한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **고주파 센서 신호의 잡음을 제거하고 FFT 변환으로 RMS, 첨도 등 주요 특징량을 추출하는 전처리 모듈을 구현한다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **표기·단위(온도 F/C, 압력 psi/kPa) 통일, 중복 제거, 불가능값 결측 처리, 파형 단발 스파이크의 중앙값 필터 제거, 간섭 잡음 창의 저역통과 필터, KNN 결측 보정 등 정제 규칙 10종을 적용 수행하고 검증하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 4개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `src/etl/mnt04_data_cleansing.py 및 정제 품질 보고서` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 군수 특수 예외(비계획 정비, 긴급 불출 등) 데이터 보존의 타당성\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: 데이터 정제 규칙 10종 적용, 데이터 정합성 검증 통과율 100%\n",
        "\n",
        "\n",
        "- 진행 기준: ETL 파이프라인 구현 및 정제된 마스터 데이터셋 구축\n",
        "\n",
        "\n",
        "- 검증 방법: 전후 데이터 통계량(기술통계) 비교 및 무결성 검증\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"표기 불일치·중복·단위 혼재·불가능값·결측·단발 스파이크가 섞인 가상 센서 원천 자료에 정제 규칙 10종(파형 스파이크 제거, 저역통과 잡음 제거, KNN 결측 보정 포함)을 적용하고, FFT 변환으로 RMS·첨도 등 주요 특징량을 추출한다. 고장 임박 시점의 긴급 점검 측정은 보존하며 정합성 검증 통과율과 정제 전후 통계를 확인한다.\"\"\"\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.stats import kurtosis\n",
        "from scipy.signal import medfilt, butter, sosfiltfilt\n",
        "from sklearn.impute import KNNImputer\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# 1) [교육용가상자료] 원천 자료: 생성 자료에 일부러 오류를 섞는다\n",
        "units, snap, W0 = make_runs()\n",
        "rng = np.random.default_rng(2)\n",
        "S = len(snap)\n",
        "W = W0.copy()\n",
        "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\"]],\n",
        "                    \"temp\": snap[\"temp\"].round(1).values, \"temp_unit\": \"C\", \"press\": snap[\"press\"].round(1).values, \"press_unit\": \"kPa\",\n",
        "                    \"load\": snap[\"load\"].round(2).values, \"emergency\": snap[\"rul\"].values <= 3})\n",
        "idx = rng.choice(S, 1500, replace=False)\n",
        "a, b, c, d_, e, f, g, h_ = np.split(idx, [150, 350, 600, 750, 800, 850, 1050])\n",
        "raw.loc[a, \"unit\"] = rng.choice([\"u\", \" U\", \"unit-\"], len(a)) + raw.loc[a, \"unit\"].str[1:]                 # ID 표기 불일치\n",
        "raw.loc[b, \"date\"] = pd.to_datetime(raw.loc[b, \"date\"]).dt.strftime(\"%Y/%m/%d\")                              # 날짜 형식 혼재\n",
        "raw.loc[c, \"temp\"] = (raw.loc[c, \"temp\"] * 9 / 5 + 32).round(1); raw.loc[c, \"temp_unit\"] = \"F\"             # 온도 단위 혼재\n",
        "raw.loc[d_, \"press\"] = (raw.loc[d_, \"press\"] / 6.895).round(1); raw.loc[d_, \"press_unit\"] = \"psi\"          # 압력 단위 혼재\n",
        "raw.loc[e, \"temp\"] = -999.0; raw.loc[f, \"press\"] = -5.0                                                    # 불가능 값\n",
        "raw.loc[g, [\"temp\", \"press\"]] = np.nan                                                                     # 결측\n",
        "sp = h_[:200]                                                                                              # 파형 단발 스파이크(센서 글리치)\n",
        "for i in sp:\n",
        "    W[i, rng.integers(0, NS, 2)] += rng.choice([-9.0, 9.0], 2)\n",
        "nz = h_[200:]                                                                                              # 전자기 간섭: 파형 전체 잡음 증가\n",
        "W[nz] += rng.normal(0, 0.9, (len(nz), NS))\n",
        "raw = pd.concat([raw, raw.sample(80, random_state=1)], ignore_index=True)                                   # 중복 80건\n",
        "print(f\"원천 자료 {len(raw)}행 (정상 {S}행 + 중복 80행, 오류 주입 {len(idx)}행, 파형 스파이크 {len(sp)}개·간섭 {len(nz)}개 포함)\")\n",
        "\n",
        "# 2) 정제 규칙 10종 (표 규칙 R1~R6, 결측 R9, 예외 R10 / 파형 규칙 R7~R8)\n",
        "df, applied = raw.iloc[:S].copy(), []\n",
        "def rule(name, n):\n",
        "    applied.append((name, int(n)))\n",
        "\n",
        "rule(\"R1 중복 행 제거\", raw.duplicated().sum())\n",
        "uid = \"U\" + df[\"unit\"].str.extract(r\"(\\d+)\")[0]; rule(\"R2 개체 ID 표기 통일\", (uid != df[\"unit\"]).sum()); df[\"unit\"] = uid\n",
        "t1 = pd.to_datetime(df[\"date\"], format=\"%Y-%m-%d\", errors=\"coerce\"); t2 = pd.to_datetime(df[\"date\"], format=\"%Y/%m/%d\", errors=\"coerce\")\n",
        "rule(\"R3 측정일 형식 통일\", t1.isna().sum()); df[\"date\"] = t1.fillna(t2)\n",
        "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())\n",
        "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())\n",
        "bt, bp = ~df[\"temp\"].between(-30, 150) & df[\"temp\"].notna(), (df[\"press\"] < 0) & df[\"press\"].notna()\n",
        "df.loc[bt, \"temp\"] = np.nan; df.loc[bp, \"press\"] = np.nan; rule(\"R6 불가능 값 결측 처리\", bt.sum() + bp.sum())\n",
        "Wc = W.copy()\n",
        "glitch = np.abs(Wc - np.median(Wc, axis=1, keepdims=True)).max(axis=1) > 8                               # 단발 스파이크 창\n",
        "for i in np.where(glitch)[0]:\n",
        "    Wc[i] = medfilt(Wc[i], 3)\n",
        "rule(\"R7 파형 단발 스파이크 제거(중앙값 필터)\", glitch.sum())\n",
        "sos = butter(4, 250, fs=FS, output=\"sos\")\n",
        "hf = (Wc - sosfiltfilt(sos, Wc, axis=1)).std(axis=1)                                                       # 250Hz 이상 성분의 크기\n",
        "noisy = hf > 0.4                                                                                          # 간섭 잡음 창\n",
        "Wc[noisy] = sosfiltfilt(sos, Wc[noisy], axis=1); rule(\"R8 간섭 잡음 창 저역통과(250Hz) 필터\", noisy.sum())\n",
        "miss = df[[\"temp\", \"press\"]].isna().sum().sum()\n",
        "X = df.assign(un=df[\"unit\"].str[1:].astype(int) / 50)[[\"temp\", \"press\", \"load\", \"un\"]]\n",
        "df[[\"temp\", \"press\"]] = KNNImputer(n_neighbors=5).fit_transform(X)[:, :2]; rule(\"R9 센서 결측 KNN 보정\", miss)\n",
        "df[\"rms\"] = np.sqrt(((Wc - Wc.mean(axis=1, keepdims=True)) ** 2).mean(axis=1)); lim = df[\"rms\"].quantile(0.99)\n",
        "df[\"outlier_flag\"] = (df[\"rms\"] > lim) & ~df[\"emergency\"]; rule(\"R10 긴급 점검 측정 보존(일반 이상값은 표시만)\", (df[\"rms\"] > lim).sum() - df[\"outlier_flag\"].sum())\n",
        "print(\"\\n[적용 규칙과 처리 건수]\")\n",
        "for nme, n in applied:\n",
        "    print(f\" {nme:36s} {n:5d}건\")\n",
        "\n",
        "# 3) 정합성 검증\n",
        "checks = {\"중복 0건\": not df.duplicated([\"unit\", \"date\"]).any(), \"ID 형식\": df[\"unit\"].str.match(r\"^U\\d{2}$\").all(),\n",
        "          \"날짜 변환\": df[\"date\"].notna().all(), \"온도 단위 C\": (df[\"temp_unit\"] == \"C\").all(), \"압력 단위 kPa\": (df[\"press_unit\"] == \"kPa\").all(),\n",
        "          \"온도 범위\": df[\"temp\"].between(-30, 150).all(), \"압력 0 이상\": (df[\"press\"] >= 0).all(), \"결측 0건\": df[[\"temp\", \"press\"]].notna().all().all(),\n",
        "          \"파형 스파이크 잔존 0개\": not (np.abs(Wc - np.median(Wc, axis=1, keepdims=True)).max(axis=1) > 8).any(),\n",
        "          f\"긴급 점검 {int(df['emergency'].sum())}건 보존\": int(df[\"emergency\"].sum()) == int(snap[\"rul\"].le(3).sum())}\n",
        "rate = np.mean(list(checks.values()))\n",
        "print(f\"\\n정합성 검증 {sum(checks.values())}/{len(checks)} 통과, 통과율 {rate:.0%} (목표 100%: {'달성' if rate == 1 else '미달'})\")\n",
        "for k, v in checks.items():\n",
        "    if not v:\n",
        "        print(\" 실패:\", k)\n",
        "\n",
        "# 4) FFT 특징량: 정제 효과(깨끗한 원본 파형 대비 오차)\n",
        "def feats(Wx):\n",
        "    x = Wx - Wx.mean(axis=1, keepdims=True)\n",
        "    return np.sqrt((x ** 2).mean(axis=1)), kurtosis(x, axis=1, fisher=False)\n",
        "r0, k0 = feats(W0); r1, k1 = feats(W); r2, k2 = feats(Wc)\n",
        "bad = glitch | noisy\n",
        "print(\"\\n[특징량 정제 효과] 스파이크·간섭이 있던 창 \" + f\"{int(bad.sum())}개의 깨끗한 원본 대비 평균 오차\")\n",
        "print(f\" RMS   정제 전 {np.abs(r1 - r0)[bad].mean():.3f} -> 정제 후 {np.abs(r2 - r0)[bad].mean():.3f}\")\n",
        "print(f\" 첨도  정제 전 {np.abs(k1 - k0)[bad].mean():.2f} -> 정제 후 {np.abs(k2 - k0)[bad].mean():.2f}\")\n",
        "cmp = pd.DataFrame({\"정제 전\": [len(raw), pd.to_numeric(raw[\"temp\"], errors=\"coerce\").mean(), raw[\"temp\"].min(), raw[\"press\"].min(), raw[\"press\"].isna().sum()],\n",
        "                    \"정제 후\": [len(df), df[\"temp\"].mean(), df[\"temp\"].min(), df[\"press\"].min(), df[\"press\"].isna().sum()]},\n",
        "                   index=[\"행 수\", \"온도 평균\", \"온도 최소\", \"압력 최소\", \"압력 결측\"]).round(1)\n",
        "print(\"\\n[정제 전후 기술통계]\"); print(cmp.to_string())\n",
        "\n",
        "i = int(np.where(glitch)[0][0])\n",
        "fig, ax = plt.subplots(1, 2, figsize=(10, 3.5))\n",
        "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()\n",
        "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\")\n",
        "plt.tight_layout(); plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 02 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "원천 자료 7813행 (정상 7733행 + 중복 80행, 오류 주입 1500행, 파형 스파이크 200개·간섭 250개 포함)\n",
        "\n",
        "[적용 규칙과 처리 건수]\n",
        " R1 중복 행 제거                              80건\n",
        " R2 개체 ID 표기 통일                         150건\n",
        " R3 측정일 형식 통일                           200건\n",
        " R4 온도 단위 통일(F -> C)                    250건\n",
        " R5 압력 단위 통일(psi -> kPa)                150건\n",
        " R6 불가능 값 결측 처리                         100건\n",
        " R7 파형 단발 스파이크 제거(중앙값 필터)               207건\n",
        " R8 간섭 잡음 창 저역통과(250Hz) 필터              347건\n",
        " R9 센서 결측 KNN 보정                        500건\n",
        " R10 긴급 점검 측정 보존(일반 이상값은 표시만)            11건\n",
        "\n",
        "정합성 검증 10/10 통과, 통과율 100% (목표 100%: 달성)\n",
        "\n",
        "[특징량 정제 효과] 스파이크·간섭이 있던 창 554개의 깨끗한 원본 대비 평균 오차\n",
        " RMS   정제 전 0.267 -> 정제 후 0.116\n",
        " 첨도  정제 전 8.31 -> 정제 후 0.86\n",
        "\n",
        "[정제 전후 기술통계]\n",
        "         정제 전    정제 후\n",
        "행 수    7813.0  7733.0\n",
        "온도 평균    73.2    77.3\n",
        "온도 최소  -999.0    61.5\n",
        "압력 최소    -5.0   342.8\n",
        "압력 결측   201.0     0.0\n",
        "```\n",
        "</details>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 03. 임계값 초과 기반 1차 센서 이상 징후 감지 Baseline 알고리즘 구현\n",
        "**소요시간**: 약 2.5시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "임계값 초과 기반 1차 센서 이상 징후 감지 Baseline 알고리즘을 구현한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **임계값 초과 기반 1차 센서 이상 징후 감지 Baseline을 구현하고, 탐색적 분석으로 고정 임계값의 한계와 수리부속 소모의 핵심 관리 대상을 찾는다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **정상 자료 분포로 RMS, 첨도, 온도, 압력 임계값을 정하고 하나라도 넘으면 이상으로 보는 Baseline의 탐지율, 정밀도, 고부하 운용일 오경보율, 고장 전 확보 시간을 측정 수행하고 검증하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 4개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `notebooks/mnt04_eda_analysis.ipynb 및 EDA 요약 리포트` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 군수 의사결정에 직관적인 인사이트 제공 수준\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: EDA 시각화 차트 8종 이상 생성, 유의미한 상관관계 피처 5개 이상 도출\n",
        "\n",
        "\n",
        "- 진행 기준: EDA 노트북 실행 완료 및 핵심 발견점(Key Findings) 문서화\n",
        "\n",
        "\n",
        "- 검증 방법: 통계적 가설 검정(p-value, 상관계수)을 통한 유의성 검증\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"임계값 초과(진동 RMS·첨도·온도·압력) 기반 1차 센서 이상 징후 감지 Baseline 알고리즘을 구현해 탐지율·오경보율·고장 전 확보 시간을 확인하고, 탐색적 분석(EDA 차트 8종, 피처와 잔존수명의 상관계수·p-value, 수리부속 소모 파레토 ABC 분류)으로 핵심 관리 대상을 식별한다.\"\"\"\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.stats import kurtosis, spearmanr\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)\n",
        "# ---------------------------------------------------------------------------\n",
        "BASE = [\"rms\", \"kurt\", \"crest\", \"p2p\", \"b50\", \"b100\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "ROLL = [\"rms\", \"kurt\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "\n",
        "\n",
        "def build_features(snap, W):\n",
        "    D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)\n",
        "    D[\"temp_res\"] = D[\"temp\"] - (55 + 22 * D[\"load\"] + 0.5 * D[\"ambient\"])        # 부하·외기 대비 온도 잔차\n",
        "    D[\"press_res\"] = D[\"press\"] - (450 - 12 * D[\"load\"])                           # 부하 대비 압력 잔차\n",
        "    D[\"rms_n\"] = D[\"rms\"] / (0.5 + D[\"load\"])                                      # 부하 정규화 RMS\n",
        "    D[\"hf_ratio\"] = (D[\"b150\"] + D[\"b165\"] + D[\"b200\"]) / (D[\"b50\"] + 1e-6)         # 고주파/기본주파 비\n",
        "    g = D.groupby(\"unit\")\n",
        "    for c in ROLL:\n",
        "        D[c + \"_m3\"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())\n",
        "    for c in (\"rms_n\", \"temp_res\", \"press_res\"):\n",
        "        D[c + \"_d\"] = g[c].diff().fillna(0)\n",
        "    return D\n",
        "\n",
        "\n",
        "DERIVED = [\"temp_res\", \"press_res\", \"rms_n\", \"hf_ratio\"] + [c + \"_m3\" for c in ROLL] + [\"rms_n_d\", \"temp_res_d\", \"press_res_d\"]\n",
        "FEATS = BASE + DERIVED\n",
        "AE_COLS = [\"rms_n\", \"kurt\", \"b150\", \"b165\", \"b200\", \"hf_ratio\", \"temp_res\", \"press_res\"]\n",
        "\n",
        "\n",
        "def split_units(units, n_train=24):\n",
        "    return list(units[\"unit\"][:n_train]), list(units[\"unit\"][n_train:])\n",
        "\n",
        "\n",
        "class AE:\n",
        "    \"\"\"오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다.\"\"\"\n",
        "\n",
        "    def __init__(self, cols=AE_COLS, q=0.99, seed=0):\n",
        "        self.cols, self.q, self.seed = cols, q, seed\n",
        "\n",
        "    def fit(self, healthy):\n",
        "        self.sc = StandardScaler().fit(healthy[self.cols])\n",
        "        Z = self.sc.transform(healthy[self.cols])\n",
        "        self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)\n",
        "        self.thr = float(np.quantile(self.err(healthy), self.q))\n",
        "        return self\n",
        "\n",
        "    def contrib(self, d):\n",
        "        Z = self.sc.transform(d[self.cols])\n",
        "        return (Z - self.net.predict(Z)) ** 2                 # 입력 항목별 재구성 오차(이상 기여도)\n",
        "\n",
        "    def err(self, d):\n",
        "        return self.contrib(d).mean(axis=1)\n",
        "\n",
        "\n",
        "def latch(d, e, thr, k=2):\n",
        "    \"\"\"같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다.\"\"\"\n",
        "    out, e = np.zeros(len(d), bool), np.asarray(e)\n",
        "    for idx in d.reset_index(drop=True).groupby(\"unit\").indices.values():\n",
        "        run, on = 0, False\n",
        "        for i in idx:\n",
        "            run = run + 1 if e[i] > thr else 0\n",
        "            on = on or run >= k\n",
        "            out[i] = on\n",
        "    return out\n",
        "\n",
        "units, snap, W = make_runs()\n",
        "D = build_features(snap, W)\n",
        "\n",
        "# 1) Baseline: 단순 임계값(하나라도 넘으면 이상). 임계값은 정상(열화 전) 자료의 분포로 정한다.\n",
        "healthy = D[D[\"deg\"] == 0]\n",
        "TH = {\"rms\": 1.25, \"kurt\": 3.0, \"temp\": 92.0, \"press_low\": 405.0}\n",
        "D[\"alarm\"] = (D[\"rms\"] > TH[\"rms\"]) | (D[\"kurt\"] > TH[\"kurt\"]) | (D[\"temp\"] > TH[\"temp\"]) | (D[\"press\"] < TH[\"press_low\"])\n",
        "D[\"anomaly\"] = D[\"deg\"] >= 0.1                                              # 정답: 열화가 일정 수준 이상 진행된 측정\n",
        "tp, fp = (D[\"alarm\"] & D[\"anomaly\"]).sum(), (D[\"alarm\"] & ~D[\"anomaly\"]).sum()\n",
        "fn = (~D[\"alarm\"] & D[\"anomaly\"]).sum()\n",
        "print(f\"[교육용가상자료] 측정 {len(D)}건, 개체 {len(units)}대, 열화 구간 {int(D['anomaly'].sum())}건\")\n",
        "print(f\"Baseline 임계값 {TH}\")\n",
        "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%})\")\n",
        "lead = []\n",
        "for u, g in D.groupby(\"unit\"):\n",
        "    a = g[g[\"alarm\"] & g[\"anomaly\"]]\n",
        "    lead.append(a[\"rul\"].max() if len(a) else 0)\n",
        "units[\"lead\"] = lead\n",
        "print(f\"  고장 전 첫 경보 확보 시간 중앙값 {np.median(lead):.0f}일, 경보 못 낸 개체 {int((units['lead'] == 0).sum())}대/{len(units)}대\")\n",
        "print(\"  열화 유형별 탐지율: \" + \", \".join(f\"{m} {D.loc[D['anomaly'] & (D['mode'] == m), 'alarm'].mean():.0%}\" for m in MODES))\n",
        "\n",
        "# 2) 상관분석: 피처와 잔존수명(RUL) (스피어만)\n",
        "rows = []\n",
        "for c in [\"rms\", \"kurt\", \"crest\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\", \"hf_ratio\", \"temp_res\", \"press_res\", \"load\", \"ambient\"]:\n",
        "    r, p = spearmanr(D[c], D[\"rul\"]); rows.append((c, r, p))\n",
        "cor = pd.DataFrame(rows, columns=[\"피처\", \"rho\", \"p\"]); cor[\"유의\"] = (cor[\"p\"] < 0.05) & (cor[\"rho\"].abs() >= 0.2)\n",
        "print(\"\\n[상관분석] 피처 | 스피어만 rho | p-value | 유의(p<0.05, |rho|>=0.2)\")\n",
        "for x in cor.sort_values(\"rho\", key=abs, ascending=False).head(8).itertuples():\n",
        "    print(f\" {x.피처:10s} {x.rho:+.2f}  {x.p:.4f}  {'유의' if x.유의 else '-'}\")\n",
        "nsig = int(cor[\"유의\"].sum())\n",
        "\n",
        "# 3) 파레토(ABC): 고장별 수리부속 소모 비용(가상 단가지수)\n",
        "cons = {}\n",
        "for m in units[\"mode\"]:\n",
        "    cons[MODES[m][2]] = cons.get(MODES[m][2], 0) + MODES[m][3]\n",
        "    cons[\"씰·필터 키트\"] = cons.get(\"씰·필터 키트\", 0) + KIT_COST\n",
        "par = pd.Series(cons).sort_values(ascending=False).to_frame(\"소모 비용지수\")\n",
        "par[\"누적비율\"] = par[\"소모 비용지수\"].cumsum() / par[\"소모 비용지수\"].sum()\n",
        "par[\"등급\"] = np.where(par[\"누적비율\"] <= 0.8, \"A\", np.where(par[\"누적비율\"] <= 0.95, \"B\", \"C\")); par.iloc[0, 2] = \"A\"\n",
        "print(\"\\n[파레토] 수리부속 소모 비용(고장 1회 교체 기준)\")\n",
        "print(par.assign(누적비율=lambda d: d[\"누적비율\"].map(\"{:.0%}\".format)).round(0).to_string())\n",
        "top20 = par.head(int(np.ceil(len(par) * 0.2)))\n",
        "print(f\"상위 20% 품목({', '.join(top20.index)})이 전체 소모 비용의 {top20['소모 비용지수'].sum() / par['소모 비용지수'].sum():.0%}\")\n",
        "\n",
        "# 4) EDA 차트 8종\n",
        "fig, ax = plt.subplots(2, 4, figsize=(16, 7))\n",
        "ax[0, 0].scatter(D[\"rul\"], D[\"rms\"], s=3); ax[0, 0].set_title(\"1 RMS vs RUL\"); ax[0, 0].invert_xaxis()\n",
        "for m in MODES:\n",
        "    g = D[D[\"mode\"] == m]; ax[0, 1].scatter(g[\"deg\"], g[\"kurt\"], s=2, label=m)\n",
        "ax[0, 1].set_title(\"2 Kurtosis vs degradation\"); ax[0, 1].legend(fontsize=6)\n",
        "for m in MODES:\n",
        "    g = D[D[\"mode\"] == m]; ax[0, 2].scatter(g[\"deg\"], g[\"temp_res\"], s=2)\n",
        "ax[0, 2].set_title(\"3 Temp residual vs degradation\")\n",
        "for m in MODES:\n",
        "    g = D[D[\"mode\"] == m]; ax[0, 3].scatter(g[\"deg\"], g[\"press_res\"], s=2)\n",
        "ax[0, 3].set_title(\"4 Pressure residual vs degradation\")\n",
        "f = np.fft.rfftfreq(NS, 1 / FS); bm = D[\"mode\"] == \"bearing\"\n",
        "ax[1, 0].plot(f, np.abs(np.fft.rfft(W[(bm & (D[\"deg\"] == 0)).values], axis=1)).mean(axis=0) / NS, label=\"healthy\")\n",
        "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)\")\n",
        "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\")\n",
        "ax[1, 2].scatter(healthy[\"load\"], healthy[\"rms\"], s=2); ax[1, 2].set_title(\"7 Load vs RMS (healthy)\")\n",
        "ax[1, 3].bar(range(len(par)), par[\"소모 비용지수\"]); ax[1, 3].set_title(\"8 Spare-part consumption (Pareto)\")\n",
        "plt.tight_layout(); plt.show()\n",
        "\n",
        "print(f\"\\nEDA 차트 8종 생성(목표 8종 이상: 달성), 유의한 피처 {nsig}개 (목표 5개 이상: {'달성' if nsig >= 5 else '미달'})\")\n",
        "print(\"핵심 발견: 1) 고정 임계값은 고부하일 오경보와 온도·압력 유형 놓침이 있음 2) 부하·외기를 보정한 잔차 피처가 열화와 더 뚜렷함 3) 베어링·기어 부품이 소모 비용 대부분\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 03 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "[교육용가상자료] 측정 7733건, 개체 48대, 열화 구간 3336건\n",
        "Baseline 임계값 {'rms': 1.25, 'kurt': 3.0, 'temp': 92.0, 'press_low': 405.0}\n",
        "  탐지율 25.3%, 정밀도 90.3%, 정상 측정 오경보율 2.1% (고부하 운용일 오경보율 28.6%)\n",
        "  고장 전 첫 경보 확보 시간 중앙값 40일, 경보 못 낸 개체 0대/48대\n",
        "  열화 유형별 탐지율: bearing 29%, gear 9%, oil_pump 54%, cooling 31%\n",
        "\n",
        "[상관분석] 피처 | 스피어만 rho | p-value | 유의(p<0.05, |rho|>=0.2)\n",
        " hf_ratio   -0.63  0.0000  유의\n",
        " temp_res   -0.53  0.0000  유의\n",
        " b165       -0.38  0.0000  유의\n",
        " b150       -0.37  0.0000  유의\n",
        " rms        -0.32  0.0000  유의\n",
        " ambient    +0.31  0.0000  유의\n",
        " kurt       -0.29  0.0000  유의\n",
        " b200       -0.21  0.0000  유의\n",
        "\n",
        "[파레토] 수리부속 소모 비용(고장 1회 교체 기준)\n",
        "          소모 비용지수  누적비율 등급\n",
        "구동 기어 세트    162.0   46%  A\n",
        "주 베어링 세트     90.0   71%  A\n",
        "씰·필터 키트      48.0   85%  B\n",
        "오일 펌프        32.0   94%  B\n",
        "냉각 팬 모터      21.0  100%  C\n",
        "상위 20% 품목(구동 기어 세트)이 전체 소모 비용의 46%\n",
        "\n",
        "EDA 차트 8종 생성(목표 8종 이상: 달성), 유의한 피처 9개 (목표 5개 이상: 달성)\n",
        "핵심 발견: 1) 고정 임계값은 고부하일 오경보와 온도·압력 유형 놓침이 있음 2) 부하·외기를 보정한 잔차 피처가 열화와 더 뚜렷함 3) 베어링·기어 부품이 소모 비용 대부분\n",
        "```\n",
        "</details>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 04. 오토인코더(Autoencoder)/LSTM 기반 미세 진동 변화 감지 예지정비(PdM) 모델 고도화\n",
        "**소요시간**: 약 3.0시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "오토인코더 기반 미세 진동 변화 감지 예지정비(PdM) 모델을 고도화한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **오토인코더 기반 미세 진동·온도 변화 감지 예지정비(PdM) 모델을 고도화하고, 파생 피처로 잔존수명 기준 모델을 만든다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **정상 자료만으로 오토인코더(MLP)를 학습해 재구성 오차를 이상 점수로 쓰고 2일 연속 초과 시 경보를 켠다 과정을 진행하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 5개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `src/features/mnt04_feature_engineering.py 및 기준 모델` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 군수 실무에서 납득 가능한 피처 설계의 논리성\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: 파생 피처 10개 이상 생성, 기준 모델 MAPE 25% 이하 달성\n",
        "\n",
        "\n",
        "- 진행 기준: 피처 파이프라인 및 베이스라인 모델 성능 지표 도출\n",
        "\n",
        "\n",
        "- 검증 방법: 교차 검증(Time-series Cross Validation)을 통한 과적합 방지\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"정상 자료만으로 학습한 오토인코더의 재구성 오차로 미세 진동·온도 변화를 감지(2일 연속 시 경보)해 임계값 Baseline과 비교하고, 파생 피처 14개로 잔존수명(RUL) 기준 모델을 만들어 개체 단위 시계열 교차검증 MAPE를 평가하며 순열 중요도로 피처를 줄이는 예지정비(PdM) 모델을 구현한다.\"\"\"\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "from scipy.stats import kurtosis\n",
        "from sklearn.ensemble import GradientBoostingRegressor\n",
        "from sklearn.inspection import permutation_importance\n",
        "from sklearn.neural_network import MLPRegressor\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)\n",
        "# ---------------------------------------------------------------------------\n",
        "BASE = [\"rms\", \"kurt\", \"crest\", \"p2p\", \"b50\", \"b100\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "ROLL = [\"rms\", \"kurt\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "\n",
        "\n",
        "def build_features(snap, W):\n",
        "    D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)\n",
        "    D[\"temp_res\"] = D[\"temp\"] - (55 + 22 * D[\"load\"] + 0.5 * D[\"ambient\"])        # 부하·외기 대비 온도 잔차\n",
        "    D[\"press_res\"] = D[\"press\"] - (450 - 12 * D[\"load\"])                           # 부하 대비 압력 잔차\n",
        "    D[\"rms_n\"] = D[\"rms\"] / (0.5 + D[\"load\"])                                      # 부하 정규화 RMS\n",
        "    D[\"hf_ratio\"] = (D[\"b150\"] + D[\"b165\"] + D[\"b200\"]) / (D[\"b50\"] + 1e-6)         # 고주파/기본주파 비\n",
        "    g = D.groupby(\"unit\")\n",
        "    for c in ROLL:\n",
        "        D[c + \"_m3\"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())\n",
        "    for c in (\"rms_n\", \"temp_res\", \"press_res\"):\n",
        "        D[c + \"_d\"] = g[c].diff().fillna(0)\n",
        "    return D\n",
        "\n",
        "\n",
        "DERIVED = [\"temp_res\", \"press_res\", \"rms_n\", \"hf_ratio\"] + [c + \"_m3\" for c in ROLL] + [\"rms_n_d\", \"temp_res_d\", \"press_res_d\"]\n",
        "FEATS = BASE + DERIVED\n",
        "AE_COLS = [\"rms_n\", \"kurt\", \"b150\", \"b165\", \"b200\", \"hf_ratio\", \"temp_res\", \"press_res\"]\n",
        "\n",
        "\n",
        "def split_units(units, n_train=24):\n",
        "    return list(units[\"unit\"][:n_train]), list(units[\"unit\"][n_train:])\n",
        "\n",
        "\n",
        "class AE:\n",
        "    \"\"\"오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다.\"\"\"\n",
        "\n",
        "    def __init__(self, cols=AE_COLS, q=0.99, seed=0):\n",
        "        self.cols, self.q, self.seed = cols, q, seed\n",
        "\n",
        "    def fit(self, healthy):\n",
        "        self.sc = StandardScaler().fit(healthy[self.cols])\n",
        "        Z = self.sc.transform(healthy[self.cols])\n",
        "        self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)\n",
        "        self.thr = float(np.quantile(self.err(healthy), self.q))\n",
        "        return self\n",
        "\n",
        "    def contrib(self, d):\n",
        "        Z = self.sc.transform(d[self.cols])\n",
        "        return (Z - self.net.predict(Z)) ** 2                 # 입력 항목별 재구성 오차(이상 기여도)\n",
        "\n",
        "    def err(self, d):\n",
        "        return self.contrib(d).mean(axis=1)\n",
        "\n",
        "\n",
        "def latch(d, e, thr, k=2):\n",
        "    \"\"\"같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다.\"\"\"\n",
        "    out, e = np.zeros(len(d), bool), np.asarray(e)\n",
        "    for idx in d.reset_index(drop=True).groupby(\"unit\").indices.values():\n",
        "        run, on = 0, False\n",
        "        for i in idx:\n",
        "            run = run + 1 if e[i] > thr else 0\n",
        "            on = on or run >= k\n",
        "            out[i] = on\n",
        "    return out\n",
        "\n",
        "units, snap, W = make_runs()\n",
        "D = build_features(snap, W)\n",
        "tr_u, va_u = split_units(units)\n",
        "tr, va = D[D[\"unit\"].isin(tr_u)].reset_index(drop=True), D[D[\"unit\"].isin(va_u)].reset_index(drop=True)\n",
        "print(f\"[교육용가상자료] 학습 {len(tr_u)}대({len(tr)}건), 검증 {len(va_u)}대({len(va)}건), 파생 피처 {len(DERIVED)}개 (목표 10개 이상: {'달성' if len(DERIVED) >= 10 else '미달'})\")\n",
        "print(\"  파생 피처:\", \", \".join(DERIVED))\n",
        "\n",
        "# 1) 오토인코더 이상 감지 vs 임계값 Baseline (검증 개체)\n",
        "ae = AE().fit(tr[tr[\"deg\"] == 0])\n",
        "e = ae.err(va)\n",
        "va[\"anomaly\"] = va[\"deg\"] >= 0.1\n",
        "va[\"ae_alarm\"] = latch(va, e, ae.thr, k=2)\n",
        "va[\"th_alarm\"] = (va[\"rms\"] > 1.25) | (va[\"kurt\"] > 3.0) | (va[\"temp\"] > 92.0) | (va[\"press\"] < 405.0)\n",
        "va[\"ae_raw\"] = e > ae.thr\n",
        "\n",
        "\n",
        "def report(name, col):\n",
        "    a = va[col]\n",
        "    tpr, far = (a & va[\"anomaly\"]).sum() / va[\"anomaly\"].sum(), (a & ~va[\"anomaly\"]).sum() / (~va[\"anomaly\"]).sum()\n",
        "    hi = a[va[\"hi_load\"] & ~va[\"anomaly\"]].mean()\n",
        "    lead = [g.loc[g[col] & g[\"anomaly\"], \"rul\"].max() if (g[col] & g[\"anomaly\"]).any() else 0 for _, g in va.groupby(\"unit\")]\n",
        "    print(f\"  {name:22s} 탐지율 {tpr:5.1%}  오경보율 {far:5.1%}  고부하일 오경보 {hi:5.1%}  고장 전 확보 중앙값 {np.median(lead):.0f}일\")\n",
        "    return tpr, far\n",
        "\n",
        "\n",
        "print(f\"\\n[오토인코더 입력 {len(AE_COLS)}개, 정상 자료 재구성 오차 99분위 {ae.thr:.2f}를 기준으로 사용]\")\n",
        "report(\"임계값 Baseline\", \"th_alarm\")\n",
        "report(\"오토인코더(1회 초과)\", \"ae_raw\")\n",
        "tpr, far = report(\"오토인코더(2일 연속)\", \"ae_alarm\")\n",
        "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%}\")\n",
        "\n",
        "# 2) RUL 기준 모델: 학습 개체 일부로 학습, 개체 단위 시계열 교차검증(앞선 개체로 학습 -> 뒤 개체 검증)\n",
        "CAP = 120                                         # RUL 상한(일): 열화 초기에는 120일 이상은 구분하지 않는다\n",
        "def fit_rul(d, cols):\n",
        "    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))\n",
        "\n",
        "\n",
        "def mape(d, p):\n",
        "    m = (d[\"deg\"] >= 0.1) & (d[\"rul\"] >= 10)      # 열화가 시작된 구간, 고장 직전 10일 미만은 분모 문제로 제외\n",
        "    y = np.minimum(d[\"rul\"], CAP)[m]\n",
        "    return float(np.mean(np.abs(y - p[m]) / y))\n",
        "\n",
        "\n",
        "order = list(units[\"unit\"])\n",
        "folds, cut = [], [12, 24, 36]\n",
        "for i, c in enumerate(cut):\n",
        "    trn, val = D[D[\"unit\"].isin(order[:c])], D[D[\"unit\"].isin(order[c:c + 12])]\n",
        "    if len(val):\n",
        "        folds.append(mape(val, fit_rul(trn, FEATS).predict(val[FEATS])))\n",
        "base_only = fit_rul(tr, BASE)\n",
        "full = fit_rul(tr, FEATS)\n",
        "m_base, m_full = mape(va, base_only.predict(va[BASE])), mape(va, full.predict(va[FEATS]))\n",
        "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 '미달'})\")\n",
        "print(f\"  검증 24대 MAPE: 원본 센서 {len(BASE)}개 {m_base:.1%} -> 파생 포함 {len(FEATS)}개 {m_full:.1%} (열화 구간, 잔존 10일 이상)\")\n",
        "\n",
        "# 3) 순열 중요도로 차원 축소\n",
        "pi = permutation_importance(full, va[FEATS], np.minimum(va[\"rul\"], CAP), n_repeats=3, random_state=0, scoring=\"neg_mean_absolute_error\")\n",
        "imp = pd.Series(pi.importances_mean, FEATS).sort_values(ascending=False)\n",
        "print(\"피처 중요도 상위 5: \" + \", \".join(f\"{k}({v:.1f})\" for k, v in imp.head(5).items()))\n",
        "top = list(imp.index[:8])\n",
        "m_top = mape(va, fit_rul(tr, top).predict(va[top]))\n",
        "print(f\"차원 축소 {len(FEATS)}개 -> {len(top)}개: 검증 MAPE {m_full:.1%} -> {m_top:.1%}\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 04 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "[교육용가상자료] 학습 24대(3847건), 검증 24대(3886건), 파생 피처 14개 (목표 10개 이상: 달성)\n",
        "  파생 피처: temp_res, press_res, rms_n, hf_ratio, rms_m3, kurt_m3, b150_m3, b165_m3, b200_m3, temp_m3, press_m3, rms_n_d, temp_res_d, press_res_d\n",
        "\n",
        "[오토인코더 입력 8개, 정상 자료 재구성 오차 99분위 1.27를 기준으로 사용]\n",
        "  임계값 Baseline           탐지율 22.7%  오경보율  1.7%  고부하일 오경보 22.8%  고장 전 확보 중앙값 40일\n",
        "  오토인코더(1회 초과)           탐지율 80.1%  오경보율  1.2%  고부하일 오경보  0.6%  고장 전 확보 중앙값 66일\n",
        "  오토인코더(2일 연속)           탐지율 79.2%  오경보율  0.0%  고부하일 오경보  0.0%  고장 전 확보 중앙값 56일\n",
        "  미세 변화 구간(열화 0.1~0.3) 탐지율: Baseline 2.0% -> 오토인코더 32.4%\n",
        "\n",
        "[RUL 모델] 개체 단위 시계열 교차검증 폴드별 MAPE: 26.7%, 24.2%, 17.1% -> 평균 22.7% (목표 25% 이하: 달성)\n",
        "  검증 24대 MAPE: 원본 센서 11개 43.3% -> 파생 포함 25개 20.4% (열화 구간, 잔존 10일 이상)\n",
        "피처 중요도 상위 5: hf_ratio(11.8), temp_res(7.0), b150_m3(3.8), b165_m3(2.8), press_m3(2.5)\n",
        "차원 축소 25개 -> 8개: 검증 MAPE 20.4% -> 20.8%\n",
        "```\n",
        "</details>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 05. 부품 고장 임박 시점(RUL) 예측 및 상태 기반 이상 신뢰도 점수 엔진 개발\n",
        "**소요시간**: 약 3.5시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "부품 고장 임박 시점(RUL) 예측 및 상태 이상 신뢰도 점수 엔진을 개발한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **부품 고장 임박 시점(RUL)을 예측하는 앙상블 모델과 상태 이상 신뢰도 점수 엔진을 개발한다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **그래디언트 부스팅, 랜덤 포레스트, 릿지를 평균하는 앙상블을 하이퍼파라미터 탐색으로 만들어 RMS와 온도만 쓰는 선형 기준 모델과 Hold-out 검증 개체에서 MAE, RMSE, R2를 비교 수행하고 검증하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 4개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `src/models/mnt04_advanced_model.py` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 비정형/불규칙 수요에 대한 모델의 적응력 및 예측 안정성\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: 기준 모델 대비 오차율(MAE/RMSE) 20% 이상 개선, R2 Score 0.85 이상\n",
        "\n",
        "\n",
        "- 진행 기준: 모델 하이퍼파라미터 튜닝 완료 및 최적 모델 가중치 저장\n",
        "\n",
        "\n",
        "- 검증 방법: Hold-out 테스트 데이터셋 대상 최종 성능 평가\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"진동·온도·압력 특징으로 부품 고장 임박 시점(RUL)을 예측하는 앙상블(그래디언트 부스팅+랜덤 포레스트+릿지) 모델을 하이퍼파라미터 탐색으로 만들어 단순 선형 기준 모델과 MAE·RMSE·R2를 비교하고, 오토인코더 이상 점수와 예측 구간으로 상태 이상 신뢰도 점수를 산출해 신뢰도 구간별 실제 정확도와 점검 카드를 출력한다.\"\"\"\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "from scipy.stats import kurtosis\n",
        "from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor\n",
        "from sklearn.linear_model import LinearRegression, Ridge\n",
        "from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score\n",
        "from sklearn.neural_network import MLPRegressor\n",
        "from sklearn.pipeline import make_pipeline\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)\n",
        "# ---------------------------------------------------------------------------\n",
        "BASE = [\"rms\", \"kurt\", \"crest\", \"p2p\", \"b50\", \"b100\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "ROLL = [\"rms\", \"kurt\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "\n",
        "\n",
        "def build_features(snap, W):\n",
        "    D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)\n",
        "    D[\"temp_res\"] = D[\"temp\"] - (55 + 22 * D[\"load\"] + 0.5 * D[\"ambient\"])        # 부하·외기 대비 온도 잔차\n",
        "    D[\"press_res\"] = D[\"press\"] - (450 - 12 * D[\"load\"])                           # 부하 대비 압력 잔차\n",
        "    D[\"rms_n\"] = D[\"rms\"] / (0.5 + D[\"load\"])                                      # 부하 정규화 RMS\n",
        "    D[\"hf_ratio\"] = (D[\"b150\"] + D[\"b165\"] + D[\"b200\"]) / (D[\"b50\"] + 1e-6)         # 고주파/기본주파 비\n",
        "    g = D.groupby(\"unit\")\n",
        "    for c in ROLL:\n",
        "        D[c + \"_m3\"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())\n",
        "    for c in (\"rms_n\", \"temp_res\", \"press_res\"):\n",
        "        D[c + \"_d\"] = g[c].diff().fillna(0)\n",
        "    return D\n",
        "\n",
        "\n",
        "DERIVED = [\"temp_res\", \"press_res\", \"rms_n\", \"hf_ratio\"] + [c + \"_m3\" for c in ROLL] + [\"rms_n_d\", \"temp_res_d\", \"press_res_d\"]\n",
        "FEATS = BASE + DERIVED\n",
        "AE_COLS = [\"rms_n\", \"kurt\", \"b150\", \"b165\", \"b200\", \"hf_ratio\", \"temp_res\", \"press_res\"]\n",
        "\n",
        "\n",
        "def split_units(units, n_train=24):\n",
        "    return list(units[\"unit\"][:n_train]), list(units[\"unit\"][n_train:])\n",
        "\n",
        "\n",
        "class AE:\n",
        "    \"\"\"오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다.\"\"\"\n",
        "\n",
        "    def __init__(self, cols=AE_COLS, q=0.99, seed=0):\n",
        "        self.cols, self.q, self.seed = cols, q, seed\n",
        "\n",
        "    def fit(self, healthy):\n",
        "        self.sc = StandardScaler().fit(healthy[self.cols])\n",
        "        Z = self.sc.transform(healthy[self.cols])\n",
        "        self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)\n",
        "        self.thr = float(np.quantile(self.err(healthy), self.q))\n",
        "        return self\n",
        "\n",
        "    def contrib(self, d):\n",
        "        Z = self.sc.transform(d[self.cols])\n",
        "        return (Z - self.net.predict(Z)) ** 2                 # 입력 항목별 재구성 오차(이상 기여도)\n",
        "\n",
        "    def err(self, d):\n",
        "        return self.contrib(d).mean(axis=1)\n",
        "\n",
        "\n",
        "def latch(d, e, thr, k=2):\n",
        "    \"\"\"같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다.\"\"\"\n",
        "    out, e = np.zeros(len(d), bool), np.asarray(e)\n",
        "    for idx in d.reset_index(drop=True).groupby(\"unit\").indices.values():\n",
        "        run, on = 0, False\n",
        "        for i in idx:\n",
        "            run = run + 1 if e[i] > thr else 0\n",
        "            on = on or run >= k\n",
        "            out[i] = on\n",
        "    return out\n",
        "\n",
        "CAP = 120\n",
        "units, snap, W = make_runs()\n",
        "D = build_features(snap, W)\n",
        "D[\"hf_max\"] = D.groupby(\"unit\")[\"hf_ratio\"].cummax()                 # 지금까지 가장 높았던 고주파 비(열화는 되돌아가지 않는다)\n",
        "D[\"temp_res_max\"] = D.groupby(\"unit\")[\"temp_res\"].cummax()\n",
        "FEATS = FEATS + [\"day\", \"hf_max\", \"temp_res_max\"]                    # 운용일수(day)는 정비 이력에서 아는 값\n",
        "tr_u, va_u = split_units(units)\n",
        "tr, va = D[D[\"unit\"].isin(tr_u)].reset_index(drop=True), D[D[\"unit\"].isin(va_u)].reset_index(drop=True)\n",
        "ytr, yva = np.minimum(tr[\"rul\"], CAP), np.minimum(va[\"rul\"], CAP)\n",
        "print(f\"[교육용가상자료] 학습 {len(tr_u)}대 {len(tr)}건, Hold-out 검증 {len(va_u)}대 {len(va)}건 (RUL은 {CAP}일 상한)\")\n",
        "\n",
        "# 1) 기준 모델: 진동 RMS와 온도만 쓰는 선형 회귀\n",
        "base = LinearRegression().fit(tr[[\"rms\", \"temp\"]], ytr)\n",
        "pb = base.predict(va[[\"rms\", \"temp\"]])\n",
        "\n",
        "# 2) 고도화: 하이퍼파라미터 탐색(학습 개체의 뒤쪽 6대를 내부 검증) 후 앙상블\n",
        "inner_tr, inner_va = tr[tr[\"unit\"].isin(tr_u[:18])], tr[tr[\"unit\"].isin(tr_u[18:])]\n",
        "\n",
        "\n",
        "def ens(d_tr, y, d_va, depth, alpha):\n",
        "    ms = [GradientBoostingRegressor(n_estimators=100, max_depth=depth, learning_rate=0.1, subsample=0.8, random_state=0),\n",
        "          RandomForestRegressor(n_estimators=40, max_depth=10, min_samples_leaf=5, n_jobs=-1, random_state=0),\n",
        "          make_pipeline(StandardScaler(), Ridge(alpha=alpha))]\n",
        "    for m in ms:\n",
        "        m.fit(d_tr[FEATS], y)\n",
        "    return ms, np.mean([m.predict(d_va[FEATS]) for m in ms], axis=0)\n",
        "\n",
        "\n",
        "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])\n",
        "         for dp in (3, 4) for al in (1, 10)}\n",
        "dp, al = min(score, key=score.get)\n",
        "models, pe = ens(tr, ytr, va, dp, al)\n",
        "names = [\"GBR\", \"RF\", \"Ridge\"]\n",
        "print(f\"선택 하이퍼파라미터: GBR depth={dp}, ridge alpha={al} (내부 검증 MAE {score[(dp, al)]:.1f}일)\")\n",
        "\n",
        "\n",
        "def metrics(p):\n",
        "    return mean_absolute_error(yva, p), mean_squared_error(yva, p) ** 0.5, r2_score(yva, p)\n",
        "\n",
        "\n",
        "mb, me = metrics(pb), metrics(pe)\n",
        "print(\"\\n[RUL 예측 성능] 모델 | MAE(일) | RMSE(일) | R2\")\n",
        "print(f\" 기준(RMS+온도 선형)  {mb[0]:6.2f}  {mb[1]:6.2f}  {mb[2]:.3f}\")\n",
        "for n, m in zip(names, models):\n",
        "    x = metrics(m.predict(va[FEATS]))\n",
        "    print(f\" 단일 {n:12s}  {x[0]:6.2f}  {x[1]:6.2f}  {x[2]:.3f}\")\n",
        "print(f\" 고도화(앙상블)      {me[0]:6.2f}  {me[1]:6.2f}  {me[2]:.3f}\")\n",
        "imp_mae, imp_rmse = 1 - me[0] / mb[0], 1 - me[1] / mb[1]\n",
        "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 '미달'})\")\n",
        "for lo, hi, nm in [(0, 30, \"고장 30일 이내\"), (30, 90, \"30~90일\"), (90, 9999, \"90일 이상\")]:\n",
        "    m = (va[\"rul\"] >= lo) & (va[\"rul\"] < hi)\n",
        "    print(f\"  구간 {nm:10s} MAE {mean_absolute_error(yva[m], pe[m]):5.1f}일 (기준 {mean_absolute_error(yva[m], pb[m]):5.1f}일)\")\n",
        "\n",
        "# 3) 상태 이상 신뢰도 점수 = 오토인코더 이상 점수(정상 대비 배수)와 앙상블 구성 모델 간 의견 일치도를 결합\n",
        "ae = AE().fit(tr[tr[\"deg\"] == 0])\n",
        "e = ae.err(va)\n",
        "z = np.log(e / ae.thr)\n",
        "conf = 1 / (1 + np.exp(-2.5 * z))                                   # 재구성 오차가 기준과 같으면 0.5\n",
        "spread = np.std([m.predict(va[FEATS]) for m in models], axis=0)\n",
        "agree = np.clip(1 - spread / 40, 0, 1)                              # 모델 간 RUL 의견 차이가 작을수록 높음\n",
        "score_c = conf * (0.7 + 0.3 * agree)\n",
        "truth = va[\"deg\"] >= 0.1\n",
        "print(\"\\n[상태 이상 신뢰도] 구간 | 건수 | 평균 점수 | 실제 열화 비율\")\n",
        "bins = [0, 0.2, 0.5, 0.8, 1.01]\n",
        "for a, b in zip(bins[:-1], bins[1:]):\n",
        "    m = (score_c >= a) & (score_c < b)\n",
        "    print(f\" {a:.1f}~{min(b, 1):.1f}  {m.sum():5d}  {score_c[m].mean():.2f}  {truth[m].mean():.2f}\")\n",
        "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())\n",
        "          for a, b in zip(bins[:-1], bins[1:]) if ((score_c >= a) & (score_c < b)).any())\n",
        "print(f\"기대 보정 오차(ECE) {ece:.3f}\")\n",
        "\n",
        "# 4) 점검 카드 예시(개체 하나의 최근 상태)\n",
        "print(\"\\n[점검 카드 예시] 개체 | 측정일 | 이상 신뢰도 | 예측 RUL(일) | 실제 RUL(일)\")\n",
        "va[\"conf\"], va[\"rul_hat\"] = score_c, pe\n",
        "for u in va_u[:3]:\n",
        "    g = va[va[\"unit\"] == u]\n",
        "    r = g[g[\"rul\"] == 40].iloc[0] if (g[\"rul\"] == 40).any() else g.iloc[-1]\n",
        "    print(f\"  {u}  {int(r['day']):3d}일째  {r['conf']:.2f}  {r['rul_hat']:5.1f}  {int(r['rul'])}  -> {'정비 계획 수립 권고' if r['conf'] >= 0.8 else '추가 관찰'}(검사관 확인)\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 05 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "[교육용가상자료] 학습 24대 3847건, Hold-out 검증 24대 3886건 (RUL은 120일 상한)\n",
        "선택 하이퍼파라미터: GBR depth=4, ridge alpha=1 (내부 검증 MAE 6.8일)\n",
        "\n",
        "[RUL 예측 성능] 모델 | MAE(일) | RMSE(일) | R2\n",
        " 기준(RMS+온도 선형)   35.67   41.43  -0.078\n",
        " 단일 GBR             9.45   13.37  0.888\n",
        " 단일 RF              8.00   12.31  0.905\n",
        " 단일 Ridge           9.98   12.78  0.897\n",
        " 고도화(앙상블)        8.56   12.05  0.909\n",
        "개선율 MAE 76.0%, RMSE 70.9% (목표 20% 이상: 달성), R2 0.909 (목표 0.85 이상: 달성)\n",
        "  구간 고장 30일 이내  MAE   5.7일 (기준  44.8일)\n",
        "  구간 30~90일     MAE  10.2일 (기준  16.9일)\n",
        "  구간 90일 이상     MAE   8.4일 (기준  47.5일)\n",
        "\n",
        "[상태 이상 신뢰도] 구간 | 건수 | 평균 점수 | 실제 열화 비율\n",
        " 0.0~0.2   2242  0.04  0.10\n",
        " 0.2~0.5    282  0.31  0.41\n",
        " 0.5~0.8    161  0.65  0.86\n",
        " 0.8~1.0   1201  0.95  1.00\n",
        "기대 보정 오차(ECE) 0.065\n",
        "\n",
        "[점검 카드 예시] 개체 | 측정일 | 이상 신뢰도 | 예측 RUL(일) | 실제 RUL(일)\n",
        "  U24   94일째  0.01   56.0  40  -> 추가 관찰(검사관 확인)\n",
        "  U25   95일째  0.62   56.6  40  -> 추가 관찰(검사관 확인)\n",
        "  U26  148일째  0.92   28.0  40  -> 정비 계획 수립 권고(검사관 확인)\n",
        "```\n",
        "</details>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 06. 장비 센서 상태 실시간 모니터링 및 결함 징후 경보 대시보드 구현\n",
        "**소요시간**: 약 3.0시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "장비 센서 상태 실시간 모니터링 및 결함 경보 대시보드를 구현한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **장비 센서 상태를 개체별로 모니터링해 결함 경보를 띄우는 정비병용 화면을 구성하고, 경보 건마다 기여도 상위 3개 인자와 자연어 요약을 자동 산출하는 설명가능성(XAI) 모듈을 구현한다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **오토인코더 재구성 오차를 입력 인자별로 나누어 경보 건마다 상위 3개 인자를 뽑고, 1순위 인자를 정상값으로 되돌렸을 때 점수가 내려가는지로 근거 유효성을 점검 수행하고 검증하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 4개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `src/explainability/mnt04_xai_explainer.py 및 XAI 리포트` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 군수 장교/군무원이 신뢰하고 조치할 수 있는 설명의 명확성\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: 주요 예측 건당 기여도 상위 3개 인자 자동 산출률 100%\n",
        "\n",
        "\n",
        "- 진행 기준: SHAP summary plot 및 waterfall plot 생성 파이프라인 완성\n",
        "\n",
        "\n",
        "- 검증 방법: 실무자 대상 설명 근거 유용성 평가 검증\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"장비 센서 상태를 개체별로 모니터링해 이상 신호를 경보로 띄우는 정비병용 화면(상태표·점검 카드)을 구성하고, 오토인코더 재구성 오차에서 경보 건마다 기여도 상위 3개 인자를 자동 산출하며(근거 유효성 검증), 의심 부품과 자연어 요약 문장을 생성하는 설명가능성(XAI) 모듈을 구현한다.\"\"\"\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.stats import kurtosis\n",
        "from sklearn.ensemble import GradientBoostingRegressor\n",
        "from sklearn.neural_network import MLPRegressor\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)\n",
        "# ---------------------------------------------------------------------------\n",
        "BASE = [\"rms\", \"kurt\", \"crest\", \"p2p\", \"b50\", \"b100\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "ROLL = [\"rms\", \"kurt\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "\n",
        "\n",
        "def build_features(snap, W):\n",
        "    D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)\n",
        "    D[\"temp_res\"] = D[\"temp\"] - (55 + 22 * D[\"load\"] + 0.5 * D[\"ambient\"])        # 부하·외기 대비 온도 잔차\n",
        "    D[\"press_res\"] = D[\"press\"] - (450 - 12 * D[\"load\"])                           # 부하 대비 압력 잔차\n",
        "    D[\"rms_n\"] = D[\"rms\"] / (0.5 + D[\"load\"])                                      # 부하 정규화 RMS\n",
        "    D[\"hf_ratio\"] = (D[\"b150\"] + D[\"b165\"] + D[\"b200\"]) / (D[\"b50\"] + 1e-6)         # 고주파/기본주파 비\n",
        "    g = D.groupby(\"unit\")\n",
        "    for c in ROLL:\n",
        "        D[c + \"_m3\"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())\n",
        "    for c in (\"rms_n\", \"temp_res\", \"press_res\"):\n",
        "        D[c + \"_d\"] = g[c].diff().fillna(0)\n",
        "    return D\n",
        "\n",
        "\n",
        "DERIVED = [\"temp_res\", \"press_res\", \"rms_n\", \"hf_ratio\"] + [c + \"_m3\" for c in ROLL] + [\"rms_n_d\", \"temp_res_d\", \"press_res_d\"]\n",
        "FEATS = BASE + DERIVED\n",
        "AE_COLS = [\"rms_n\", \"kurt\", \"b150\", \"b165\", \"b200\", \"hf_ratio\", \"temp_res\", \"press_res\"]\n",
        "\n",
        "\n",
        "def split_units(units, n_train=24):\n",
        "    return list(units[\"unit\"][:n_train]), list(units[\"unit\"][n_train:])\n",
        "\n",
        "\n",
        "class AE:\n",
        "    \"\"\"오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다.\"\"\"\n",
        "\n",
        "    def __init__(self, cols=AE_COLS, q=0.99, seed=0):\n",
        "        self.cols, self.q, self.seed = cols, q, seed\n",
        "\n",
        "    def fit(self, healthy):\n",
        "        self.sc = StandardScaler().fit(healthy[self.cols])\n",
        "        Z = self.sc.transform(healthy[self.cols])\n",
        "        self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)\n",
        "        self.thr = float(np.quantile(self.err(healthy), self.q))\n",
        "        return self\n",
        "\n",
        "    def contrib(self, d):\n",
        "        Z = self.sc.transform(d[self.cols])\n",
        "        return (Z - self.net.predict(Z)) ** 2                 # 입력 항목별 재구성 오차(이상 기여도)\n",
        "\n",
        "    def err(self, d):\n",
        "        return self.contrib(d).mean(axis=1)\n",
        "\n",
        "\n",
        "def latch(d, e, thr, k=2):\n",
        "    \"\"\"같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다.\"\"\"\n",
        "    out, e = np.zeros(len(d), bool), np.asarray(e)\n",
        "    for idx in d.reset_index(drop=True).groupby(\"unit\").indices.values():\n",
        "        run, on = 0, False\n",
        "        for i in idx:\n",
        "            run = run + 1 if e[i] > thr else 0\n",
        "            on = on or run >= k\n",
        "            out[i] = on\n",
        "    return out\n",
        "\n",
        "units, snap, W = make_runs()\n",
        "D = build_features(snap, W)\n",
        "tr_u, va_u = split_units(units)\n",
        "tr, va = D[D[\"unit\"].isin(tr_u)].reset_index(drop=True), D[D[\"unit\"].isin(va_u)].reset_index(drop=True)\n",
        "ae = AE().fit(tr[tr[\"deg\"] == 0])\n",
        "C = ae.contrib(va)                                         # 건 x 입력항목 기여도\n",
        "e = C.mean(axis=1)\n",
        "va[\"alarm\"] = latch(va, e, ae.thr, k=2)\n",
        "va[\"score\"] = e / ae.thr\n",
        "\n",
        "LABEL = {\"rms_n\": \"부하 보정 진동 크기(RMS)\", \"kurt\": \"진동 충격성(첨도)\", \"b150\": \"150Hz 부근 진동(기어 맞물림 대역)\", \"b165\": \"165Hz 부근 진동(베어링 결함 대역)\",\n",
        "         \"b200\": \"200Hz 부근 진동(펌프 대역)\", \"hf_ratio\": \"고주파/기본주파 에너지 비\", \"temp_res\": \"부하 대비 온도 상승\", \"press_res\": \"부하 대비 압력 저하\"}\n",
        "PART = {\"b165\": \"주 베어링 세트\", \"kurt\": \"주 베어링 세트\", \"b150\": \"구동 기어 세트\", \"b200\": \"오일 펌프\", \"press_res\": \"오일 펌프\", \"temp_res\": \"냉각 팬 모터\"}\n",
        "TRUE_PART = {m: MODES[m][2] for m in MODES}\n",
        "cols = ae.cols\n",
        "\n",
        "# 1) 경보 건마다 기여도 상위 3개 인자 자동 산출\n",
        "idx = np.where(va[\"alarm\"].values)[0]\n",
        "top3 = np.argsort(-C[idx], axis=1)[:, :3]\n",
        "ok = sum(len(set(t)) == 3 for t in top3)\n",
        "print(f\"[교육용가상자료] 검증 {len(va_u)}대 {len(va)}건 중 경보 {len(idx)}건: 기여도 상위 3개 인자 산출 {ok}건 (산출률 {ok / len(idx):.0%}, 목표 100%: {'달성' if ok == len(idx) else '미달'})\")\n",
        "\n",
        "# 2) 근거 유효성: 1순위 인자를 정상 평균(표준화 0)으로 되돌리면 이상 점수가 낮아지는가\n",
        "Z = ae.sc.transform(va[cols].iloc[idx])\n",
        "Z2 = Z.copy()\n",
        "Z2[np.arange(len(idx)), top3[:, 0]] = 0.0\n",
        "e1, e2 = ((Z - ae.net.predict(Z)) ** 2).mean(axis=1), ((Z2 - ae.net.predict(Z2)) ** 2).mean(axis=1)\n",
        "print(f\"근거 유효성: 1순위 인자를 정상값으로 되돌리면 이상 점수가 낮아지는 비율 {(e2 < e1).mean():.1%} (평균 {e1.mean():.2f} -> {e2.mean():.2f})\")\n",
        "\n",
        "# 3) 의심 부품: 1순위 인자 -> 부품 매핑 (열화가 뚜렷한 경보에서 정답 비교)\n",
        "sub = va.iloc[idx].copy()\n",
        "sub[\"top1\"] = [cols[t[0]] for t in top3]\n",
        "sub[\"hint\"] = sub[\"top1\"].map(PART).fillna(\"미분류(점검 필요)\")\n",
        "sub[\"true_part\"] = sub[\"mode\"].map(TRUE_PART)\n",
        "strong = sub[sub[\"deg\"] >= 0.3]\n",
        "print(f\"의심 부품 추정 정확도(열화 0.3 이상 경보 {len(strong)}건): {(strong['hint'] == strong['true_part']).mean():.1%}\")\n",
        "print(\"  1순위 인자 분포: \" + \", \".join(f\"{LABEL[k].split('(')[0]} {v}\" for k, v in sub[\"top1\"].value_counts().head(4).items()))\n",
        "\n",
        "\n",
        "def summarize(r, t):\n",
        "    f = [LABEL[cols[i]] for i in t]\n",
        "    part = PART.get(cols[t[0]], \"원인 불명\")\n",
        "    return f\"{part} 이상 가능성(이상 점수 {r['score']:.1f}배). 주요 근거: {f[0]}, {f[1]}, {f[2]}.\"\n",
        "\n",
        "\n",
        "# 4) 정비병용 화면: 점검 카드와 상태표\n",
        "print(\"\\n[점검 카드] 개체 | 측정일 | 실제 열화 유형\")\n",
        "for m in MODES:\n",
        "    k = sub[(sub[\"mode\"] == m) & (sub[\"deg\"] > 0.5)].index[0]\n",
        "    j = list(sub.index).index(k)\n",
        "    r = sub.loc[k]\n",
        "    print(f\" {r['unit']}  {int(r['day'])}일째  ({m})\")\n",
        "    print(\"   \" + summarize(r, top3[j]))\n",
        "    print(\"   근거 기여도: \" + \" / \".join(f\"{LABEL[cols[i]].split('(')[0]}({C[idx[j], i]:.1f})\" for i in top3[j]))\n",
        "\n",
        "CAL = int(va.groupby('cal')['alarm'].sum().idxmax())               # 경보가 가장 많은 날을 가상 기준일로 쓴다\n",
        "act = va[(va[\"cal\"] == CAL)]\n",
        "reg = {\"정상\": 0, \"주의\": 0, \"경보\": 0}\n",
        "rows = []\n",
        "for r in act.itertuples():\n",
        "    st = \"경보\" if r.alarm else (\"주의\" if r.score > 0.7 else \"정상\")\n",
        "    reg[st] += 1\n",
        "for st in (\"경보\", \"주의\"):\n",
        "    for r in act[(act[\"alarm\"] if st == \"경보\" else (~act[\"alarm\"] & (act[\"score\"] > 0.7)))].itertuples():\n",
        "        rows.append((r.unit, st, f\"{r.score:.1f}\", int(r.rul)))\n",
        "print(f\"\\n[운용 상태표] 가상 기준일 {CAL}: 가동 중 {len(act)}대 -> \" + \", \".join(f\"{k} {v}대\" for k, v in reg.items()))\n",
        "for r in rows[:5]:\n",
        "    print(f\"  {r[0]}  {r[1]}  이상 점수 {r[2]}배  (실제 잔존 {r[3]}일)\")\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "for u in sub[\"unit\"].unique()[:3]:\n",
        "    g = va[va[\"unit\"] == u]\n",
        "    plt.plot(g[\"day\"], g[\"score\"], label=u)\n",
        "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\")\n",
        "plt.tight_layout(); plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 06 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "[교육용가상자료] 검증 24대 3886건 중 경보 1327건: 기여도 상위 3개 인자 산출 1327건 (산출률 100%, 목표 100%: 달성)\n",
        "근거 유효성: 1순위 인자를 정상값으로 되돌리면 이상 점수가 낮아지는 비율 98.9% (평균 21.19 -> 6.57)\n",
        "의심 부품 추정 정확도(열화 0.3 이상 경보 1167건): 99.3%\n",
        "  1순위 인자 분포: 150Hz 부근 진동 603, 165Hz 부근 진동 338, 부하 대비 온도 상승 146, 부하 대비 압력 저하 129\n",
        "\n",
        "[점검 카드] 개체 | 측정일 | 실제 열화 유형\n",
        " U27  131일째  (bearing)\n",
        "   주 베어링 세트 이상 가능성(이상 점수 7.1배). 주요 근거: 165Hz 부근 진동(베어링 결함 대역), 150Hz 부근 진동(기어 맞물림 대역), 부하 대비 압력 저하.\n",
        "   근거 기여도: 165Hz 부근 진동(36.6) / 150Hz 부근 진동(22.9) / 부하 대비 압력 저하(5.4)\n",
        " U26  150일째  (gear)\n",
        "   구동 기어 세트 이상 가능성(이상 점수 10.4배). 주요 근거: 150Hz 부근 진동(기어 맞물림 대역), 165Hz 부근 진동(베어링 결함 대역), 진동 충격성(첨도).\n",
        "   근거 기여도: 150Hz 부근 진동(60.3) / 165Hz 부근 진동(19.9) / 진동 충격성(10.9)\n",
        " U25  108일째  (oil_pump)\n",
        "   오일 펌프 이상 가능성(이상 점수 6.0배). 주요 근거: 부하 대비 압력 저하, 200Hz 부근 진동(펌프 대역), 165Hz 부근 진동(베어링 결함 대역).\n",
        "   근거 기여도: 부하 대비 압력 저하(29.3) / 200Hz 부근 진동(22.0) / 165Hz 부근 진동(4.7)\n",
        " U24  107일째  (cooling)\n",
        "   냉각 팬 모터 이상 가능성(이상 점수 4.4배). 주요 근거: 부하 대비 온도 상승, 200Hz 부근 진동(펌프 대역), 진동 충격성(첨도).\n",
        "   근거 기여도: 부하 대비 온도 상승(34.5) / 200Hz 부근 진동(3.9) / 진동 충격성(2.4)\n",
        "\n",
        "[운용 상태표] 가상 기준일 308: 가동 중 18대 -> 정상 6대, 주의 1대, 경보 11대\n",
        "  U26  경보  이상 점수 30.5배  (실제 잔존 1일)\n",
        "  U28  경보  이상 점수 55.0배  (실제 잔존 5일)\n",
        "  U31  경보  이상 점수 85.1배  (실제 잔존 1일)\n",
        "  U32  경보  이상 점수 12.9배  (실제 잔존 14일)\n",
        "  U33  경보  이상 점수 27.4배  (실제 잔존 14일)\n",
        "```\n",
        "</details>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 07. 예방 정비 조기 입고 결정을 위한 항공/기갑 검사관 HITL 승인 모듈 설계\n",
        "**소요시간**: 약 3.5시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "예방 정비 조기 입고 결정을 위한 검사관 HITL 승인 모듈을 설계한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **센서 경보 시 예방 정비 조기 입고를 제안하되 항공/기갑 검사관이 승인해야 확정되는 HITL 승인 모듈을 설계하고 효과를 시뮬레이션으로 검증한다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **경보와 RUL 예측으로 입고일을 제안하고 승인 대기열에 올려 검사관의 승인·반려·무응답(기본 계획인 정기 교체 유지)을 사유와 함께 해시 체인 감사 로그로 남긴다 과정을 진행하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 4개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `src/optimization/mnt04_simulator.py` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 현실적 군수 제약조건 반영도 및 시나리오별 대안 제시 적절성\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: 최적화 적용 시 비용 15% 절감 또는 재고 회전율 20% 향상\n",
        "\n",
        "\n",
        "- 진행 기준: 선형/정수 계획법(LP/MIP) 최적화 해 도출 및 시뮬레이션 완료\n",
        "\n",
        "\n",
        "- 검증 방법: 과거 실적 데이터 대비 최적화 결과 비교 백테스팅\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"센서 이상 경보가 나면 예방 정비 조기 입고일을 제안하되 자동 확정하지 않고 항공/기갑 검사관 승인 대기열에 올려 승인·반려·무응답(기본 계획 유지)을 사유와 함께 해시 체인 감사 로그에 남기는 HITL 승인 모듈을 시뮬레이션하고, 정기 교체 방식 대비 정비·고장 비용 절감률과 승인 지연·전시 전환(가동 강도 급증) 민감도를 분석한다.\"\"\"\n",
        "import hashlib\n",
        "import json\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.stats import kurtosis\n",
        "from sklearn.ensemble import GradientBoostingRegressor\n",
        "from sklearn.neural_network import MLPRegressor\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)\n",
        "# ---------------------------------------------------------------------------\n",
        "BASE = [\"rms\", \"kurt\", \"crest\", \"p2p\", \"b50\", \"b100\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "ROLL = [\"rms\", \"kurt\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "\n",
        "\n",
        "def build_features(snap, W):\n",
        "    D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)\n",
        "    D[\"temp_res\"] = D[\"temp\"] - (55 + 22 * D[\"load\"] + 0.5 * D[\"ambient\"])        # 부하·외기 대비 온도 잔차\n",
        "    D[\"press_res\"] = D[\"press\"] - (450 - 12 * D[\"load\"])                           # 부하 대비 압력 잔차\n",
        "    D[\"rms_n\"] = D[\"rms\"] / (0.5 + D[\"load\"])                                      # 부하 정규화 RMS\n",
        "    D[\"hf_ratio\"] = (D[\"b150\"] + D[\"b165\"] + D[\"b200\"]) / (D[\"b50\"] + 1e-6)         # 고주파/기본주파 비\n",
        "    g = D.groupby(\"unit\")\n",
        "    for c in ROLL:\n",
        "        D[c + \"_m3\"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())\n",
        "    for c in (\"rms_n\", \"temp_res\", \"press_res\"):\n",
        "        D[c + \"_d\"] = g[c].diff().fillna(0)\n",
        "    return D\n",
        "\n",
        "\n",
        "DERIVED = [\"temp_res\", \"press_res\", \"rms_n\", \"hf_ratio\"] + [c + \"_m3\" for c in ROLL] + [\"rms_n_d\", \"temp_res_d\", \"press_res_d\"]\n",
        "FEATS = BASE + DERIVED\n",
        "AE_COLS = [\"rms_n\", \"kurt\", \"b150\", \"b165\", \"b200\", \"hf_ratio\", \"temp_res\", \"press_res\"]\n",
        "\n",
        "\n",
        "def split_units(units, n_train=24):\n",
        "    return list(units[\"unit\"][:n_train]), list(units[\"unit\"][n_train:])\n",
        "\n",
        "\n",
        "class AE:\n",
        "    \"\"\"오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다.\"\"\"\n",
        "\n",
        "    def __init__(self, cols=AE_COLS, q=0.99, seed=0):\n",
        "        self.cols, self.q, self.seed = cols, q, seed\n",
        "\n",
        "    def fit(self, healthy):\n",
        "        self.sc = StandardScaler().fit(healthy[self.cols])\n",
        "        Z = self.sc.transform(healthy[self.cols])\n",
        "        self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)\n",
        "        self.thr = float(np.quantile(self.err(healthy), self.q))\n",
        "        return self\n",
        "\n",
        "    def contrib(self, d):\n",
        "        Z = self.sc.transform(d[self.cols])\n",
        "        return (Z - self.net.predict(Z)) ** 2                 # 입력 항목별 재구성 오차(이상 기여도)\n",
        "\n",
        "    def err(self, d):\n",
        "        return self.contrib(d).mean(axis=1)\n",
        "\n",
        "\n",
        "def latch(d, e, thr, k=2):\n",
        "    \"\"\"같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다.\"\"\"\n",
        "    out, e = np.zeros(len(d), bool), np.asarray(e)\n",
        "    for idx in d.reset_index(drop=True).groupby(\"unit\").indices.values():\n",
        "        run, on = 0, False\n",
        "        for i in idx:\n",
        "            run = run + 1 if e[i] > thr else 0\n",
        "            on = on or run >= k\n",
        "            out[i] = on\n",
        "    return out\n",
        "\n",
        "PERIOD = 90           # 기본 계획: 정기 교체 주기(운용일). 가장 빠른 고장 사례보다 짧게 정해 사후 고장을 막는 보수적 주기\n",
        "UNPLANNED = (2.0, 12.0)   # [교육용가상자료] 고장 후 정비 비용 = 부품 단가 x 2.0 + 연쇄 손상·가동중단 12.0\n",
        "PLANNED_FIXED = 1.0   # 입고 정비 작업 고정 비용(지수)\n",
        "INSPECT = {\"base\": 0.3, \"pdm\": 0.05}   # 30일당 점검 비용: 정기 방식은 분해 점검, 센서 감시 방식은 경보 시 확인 점검 위주(가정)\n",
        "\n",
        "\n",
        "class Audit:\n",
        "    def __init__(self):\n",
        "        self.log = []\n",
        "\n",
        "    def add(self, **body):\n",
        "        prev = self.log[-1][\"hash\"] if self.log else \"GENESIS\"\n",
        "        body = {\"seq\": len(self.log), **body}\n",
        "        self.log.append({**body, \"prev\": prev, \"hash\": hashlib.sha256((prev + json.dumps(body, sort_keys=True, ensure_ascii=False)).encode()).hexdigest()})\n",
        "\n",
        "    def verify(self):\n",
        "        prev = \"GENESIS\"\n",
        "        for e in self.log:\n",
        "            b = {k: v for k, v in e.items() if k not in (\"prev\", \"hash\")}\n",
        "            if e[\"prev\"] != prev or hashlib.sha256((prev + json.dumps(b, sort_keys=True, ensure_ascii=False)).encode()).hexdigest() != e[\"hash\"]:\n",
        "                return e[\"seq\"]\n",
        "            prev = e[\"hash\"]\n",
        "        return None\n",
        "\n",
        "\n",
        "def add_feat(D):\n",
        "    D[\"hf_max\"] = D.groupby(\"unit\")[\"hf_ratio\"].cummax(); D[\"temp_res_max\"] = D.groupby(\"unit\")[\"temp_res\"].cummax()\n",
        "    return D\n",
        "\n",
        "\n",
        "F7 = FEATS + [\"day\", \"hf_max\", \"temp_res_max\"]\n",
        "u1, s1, w1 = make_runs(seed=1)\n",
        "D1 = add_feat(build_features(s1, w1))\n",
        "tr_u, _ = split_units(u1)\n",
        "tr = D1[D1[\"unit\"].isin(tr_u)]\n",
        "ae = AE().fit(tr[tr[\"deg\"] == 0])\n",
        "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))\n",
        "\n",
        "# 시험 개체: 다른 시드로 만든 가상 개체 3묶음(144대)에서 개체별 첫 경보일과 그때의 RUL 예측을 모은다\n",
        "cases = []\n",
        "for sd in (2, 3, 4):\n",
        "    un, sn, ww = make_runs(seed=sd)\n",
        "    D = add_feat(build_features(sn, ww))\n",
        "    D[\"alarm\"] = latch(D, ae.err(D), ae.thr, k=2)\n",
        "    D[\"rul_hat\"] = gbr.predict(D[F7])\n",
        "    for u in un.itertuples():\n",
        "        g = D[D[\"unit\"] == u.unit]\n",
        "        a = g[g[\"alarm\"]]\n",
        "        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,\n",
        "                      \"rul_hat\": float(a[\"rul_hat\"].iloc[0]) if len(a) else None, \"load\": float(a[\"load\"].iloc[0]) if len(a) else None})\n",
        "C = pd.DataFrame(cases)\n",
        "print(f\"[교육용가상자료] 시험 개체 {len(C)}대(수명 {C['life'].min()}~{C['life'].max()}일), 경보 발생 {int(C['a'].notna().sum())}대\")\n",
        "\n",
        "\n",
        "def cost(mode, life, maint, kind=\"pdm\"):\n",
        "    p = MODES[mode][3] + KIT_COST + INSPECT[kind] * min(maint if maint is not None else life, life) / 30\n",
        "    if maint is None or maint >= life:                       # 정비 전에 고장\n",
        "        return MODES[mode][3] * UNPLANNED[0] + KIT_COST + UNPLANNED[1] + INSPECT[kind] * life / 30, True\n",
        "    return p + PLANNED_FIXED + MODES[mode][3] * (life - maint) / life, False   # 남은 수명 낭비분 포함\n",
        "\n",
        "\n",
        "def commander(row, rng):\n",
        "    \"\"\"[모의 검사관] 잔존 예측이 길면 반려, 10%는 무응답.\"\"\"\n",
        "    if rng.random() < 0.10:\n",
        "        return \"무응답\", \"기한 내 응답 없음 -> 기본 계획(정기 교체) 유지\"\n",
        "    if row[\"rul_hat\"] > 90:\n",
        "        return \"반려\", f\"예측 잔존 {row['rul_hat']:.0f}일로 충분 -> 현 계획 유지\"\n",
        "    if row[\"load\"] > 1.0:\n",
        "        return \"반려\", \"고부하 운용일 신호로 판단 -> 추가 관찰\"\n",
        "    return \"승인\", f\"이상 지속, 예측 잔존 {row['rul_hat']:.0f}일 -> 조기 입고\"\n",
        "\n",
        "\n",
        "def run(mode, delay=2, scale=1.0, seed=0, audit=None):\n",
        "    \"\"\"mode: base(정기 교체) / hitl(제안 -> 검사관 결정) / auto(모든 제안 적용: 비교용 상한, 실제 운용 금지)\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    tot, fails, dec = 0.0, 0, {\"승인\": 0, \"반려\": 0, \"무응답\": 0}\n",
        "    for _, r in C.iterrows():\n",
        "        life = r[\"life\"] / scale\n",
        "        maint, kind = PERIOD, \"base\"\n",
        "        if mode != \"base\" and r[\"a\"] is not None:\n",
        "            a, rh = r[\"a\"] / scale, r[\"rul_hat\"] / scale\n",
        "            inbound = a + max(int(0.4 * rh), 3)\n",
        "            row = {\"rul_hat\": rh, \"load\": r[\"load\"]}\n",
        "            if mode == \"auto\":\n",
        "                d, why = \"승인\", \"비교용 전부 적용\"\n",
        "            else:\n",
        "                d, why = commander(row, rng)\n",
        "                audit.add(event=\"PROPOSE\", id=r[\"id\"], day=round(a), inbound=int(inbound), status=\"PENDING\")\n",
        "                audit.add(event=\"DECIDE\", id=r[\"id\"], decision=d, reason=why, by=\"검사관(모의)\")\n",
        "            dec[d] += 1\n",
        "            if d == \"승인\":\n",
        "                maint, kind = inbound + (delay if mode == \"hitl\" else 0), \"pdm\"\n",
        "        c, f = cost(r[\"mode\"], life, maint, kind)\n",
        "        tot += c; fails += f\n",
        "    return tot, fails, dec\n",
        "\n",
        "\n",
        "b_cost, b_fail, _ = run(\"base\")\n",
        "audit = Audit()\n",
        "h_cost, h_fail, dec = run(\"hitl\", seed=3, audit=audit)\n",
        "a_cost, a_fail, _ = run(\"auto\")\n",
        "print(\"\\n[비용 비교] 방식 | 비용지수 | 사후 고장 건수\")\n",
        "for nm, c, f in [(\"정기 교체(기본 계획)\", b_cost, b_fail), (\"경보 제안 + 검사관 승인\", h_cost, h_fail), (\"제안 전부 적용(비교용)\", a_cost, a_fail)]:\n",
        "    print(f\" {nm:22s} {c:8.0f}  {f:3d}\")\n",
        "cut = 1 - h_cost / b_cost\n",
        "print(f\"HITL 적용 비용 절감률 {cut:.1%} (목표 15% 이상: {'달성' if cut >= 0.15 else '미달'}), 고장 {b_fail} -> {h_fail}건\")\n",
        "print(f\"결정 {sum(dec.values())}건: \" + \", \".join(f\"{k} {v}건\" for k, v in dec.items()))\n",
        "n_prop = sum(e[\"event\"] == \"PROPOSE\" for e in audit.log); n_dec = sum(e[\"event\"] == \"DECIDE\" for e in audit.log)\n",
        "print(f\"제안 {n_prop}건 중 결정 기록 {n_dec}건(로깅률 {n_dec / n_prop:.0%}), 자동 승인 0건, 무응답 시 기본 계획 유지\")\n",
        "rej = [e[\"reason\"] for e in audit.log if e.get(\"decision\") == \"반려\"]\n",
        "print(\"반려 사유 예:\", rej[0] if rej else \"-\")\n",
        "print(f\"감사 로그 무결성: {'정상' if audit.verify() is None else '이상'}\", end=\", \")\n",
        "audit.log[5][\"reason\"] = \"사후 수정\"\n",
        "print(f\"사후 변조 시 seq {audit.verify()}에서 탐지\")\n",
        "\n",
        "print(\"\\n[민감도 1] 승인 후 입고까지 지연(일) -> 비용 절감률\")\n",
        "for d in (0, 2, 7, 14, 21):\n",
        "    c, f, _ = run(\"hitl\", delay=d, seed=3, audit=Audit())\n",
        "    print(f\"  {d:2d}일: {1 - c / b_cost:.1%} (고장 {f}건)\")\n",
        "print(\"[민감도 2] 가동 강도 급증(전시 전환 가정, 열화가 빨라져 수명이 1/배수로 단축) -> 절감률\")\n",
        "for s in (1.0, 1.3, 1.6, 2.0):\n",
        "    bc, bf, _ = run(\"base\", scale=s); hc, hf, _ = run(\"hitl\", scale=s, seed=3, audit=Audit())\n",
        "    print(f\"  x{s}: 정기 교체 비용 {bc:.0f}(고장 {bf}) -> HITL {hc:.0f}(고장 {hf}), 절감 {1 - hc / bc:.1%}\")\n",
        "\n",
        "plt.figure(figsize=(6, 4))\n",
        "plt.bar([\"periodic\", \"HITL\", \"auto(ref)\"], [b_cost, h_cost, a_cost]); plt.ylabel(\"cost index\"); plt.title(\"HITL early-inbound simulator\")\n",
        "plt.tight_layout(); plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 07 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "[교육용가상자료] 시험 개체 144대(수명 101~212일), 경보 발생 144대\n",
        "\n",
        "[비용 비교] 방식 | 비용지수 | 사후 고장 건수\n",
        " 정기 교체(기본 계획)               1640    0\n",
        " 경보 제안 + 검사관 승인             1389    1\n",
        " 제안 전부 적용(비교용)              1364    1\n",
        "HITL 적용 비용 절감률 15.3% (목표 15% 이상: 달성), 고장 0 -> 1건\n",
        "결정 144건: 승인 126건, 반려 7건, 무응답 11건\n",
        "제안 144건 중 결정 기록 144건(로깅률 100%), 자동 승인 0건, 무응답 시 기본 계획 유지\n",
        "반려 사유 예: 고부하 운용일 신호로 판단 -> 추가 관찰\n",
        "감사 로그 무결성: 정상, 사후 변조 시 seq 5에서 탐지\n",
        "\n",
        "[민감도 1] 승인 후 입고까지 지연(일) -> 비용 절감률\n",
        "   0일: 14.8% (고장 1건)\n",
        "   2일: 15.3% (고장 1건)\n",
        "   7일: 14.1% (고장 4건)\n",
        "  14일: 10.5% (고장 10건)\n",
        "  21일: 1.8% (고장 21건)\n",
        "[민감도 2] 가동 강도 급증(전시 전환 가정, 열화가 빨라져 수명이 1/배수로 단축) -> 절감률\n",
        "  x1.0: 정기 교체 비용 1640(고장 0) -> HITL 1389(고장 1), 절감 15.3%\n",
        "  x1.3: 정기 교체 비용 1571(고장 5) -> HITL 1363(고장 1), 절감 13.2%\n",
        "  x1.6: 정기 교체 비용 2005(고장 41) -> HITL 1451(고장 8), 절감 27.7%\n",
        "  x2.0: 정기 교체 비용 3394(고장 128) -> HITL 1572(고장 16), 절감 53.7%\n",
        "```\n",
        "</details>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 08. 장시간 비행/기동 센서 데이터 대용량 스트림 주입 스트레스 테스트\n",
        "**소요시간**: 약 3.0시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "장시간 기동 센서 데이터 대용량 스트림 주입 스트레스 테스트를 수행한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **장시간 기동 센서 데이터를 대용량 스트림으로 주입해 처리 시간, 경보 정확도, 불량 패킷 격리를 확인하는 스트레스 시험을 수행한다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **가상 센서 측정 1,000건에 결측 패킷, 센서 고정, 포화, 범위 이탈을 섞어 한꺼번에 주입하고 정상 처리·격리·유실 건수와 처리 시간을 센다 과정을 진행하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 4개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `src/dashboard/mnt04_app.py (Streamlit/Dash 등)` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 의사결정 지원에 최적화된 시각적 위계 및 사용자 경험(UX)\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: 대시보드 로딩 시간 2초 이내, 경보 발령 정확도 95% 이상\n",
        "\n",
        "\n",
        "- 진행 기준: 웹 기반 분석 대시보드 구동 및 실시간 필터링 기능 검증\n",
        "\n",
        "\n",
        "- 검증 방법: 실제 군수 데이터 연동 화면 테스트 및 UI 동작 검증\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"장시간 기동 센서 데이터를 1,000~100,000건 규모로 한꺼번에 주입해 검증·특징 추출·오토인코더 점수·경보 판정을 배치로 처리하는 스트림 시험을 수행하고, 처리 시간과 상태 화면 구성 시간, 경보 정확도, 불량 패킷(결측·고정값·포화·중복) 격리와 유실 0건 여부를 확인한다.\"\"\"\n",
        "import time\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.stats import kurtosis\n",
        "from sklearn.neural_network import MLPRegressor\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)\n",
        "# ---------------------------------------------------------------------------\n",
        "BASE = [\"rms\", \"kurt\", \"crest\", \"p2p\", \"b50\", \"b100\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "ROLL = [\"rms\", \"kurt\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "\n",
        "\n",
        "def build_features(snap, W):\n",
        "    D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)\n",
        "    D[\"temp_res\"] = D[\"temp\"] - (55 + 22 * D[\"load\"] + 0.5 * D[\"ambient\"])        # 부하·외기 대비 온도 잔차\n",
        "    D[\"press_res\"] = D[\"press\"] - (450 - 12 * D[\"load\"])                           # 부하 대비 압력 잔차\n",
        "    D[\"rms_n\"] = D[\"rms\"] / (0.5 + D[\"load\"])                                      # 부하 정규화 RMS\n",
        "    D[\"hf_ratio\"] = (D[\"b150\"] + D[\"b165\"] + D[\"b200\"]) / (D[\"b50\"] + 1e-6)         # 고주파/기본주파 비\n",
        "    g = D.groupby(\"unit\")\n",
        "    for c in ROLL:\n",
        "        D[c + \"_m3\"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())\n",
        "    for c in (\"rms_n\", \"temp_res\", \"press_res\"):\n",
        "        D[c + \"_d\"] = g[c].diff().fillna(0)\n",
        "    return D\n",
        "\n",
        "\n",
        "DERIVED = [\"temp_res\", \"press_res\", \"rms_n\", \"hf_ratio\"] + [c + \"_m3\" for c in ROLL] + [\"rms_n_d\", \"temp_res_d\", \"press_res_d\"]\n",
        "FEATS = BASE + DERIVED\n",
        "AE_COLS = [\"rms_n\", \"kurt\", \"b150\", \"b165\", \"b200\", \"hf_ratio\", \"temp_res\", \"press_res\"]\n",
        "\n",
        "\n",
        "def split_units(units, n_train=24):\n",
        "    return list(units[\"unit\"][:n_train]), list(units[\"unit\"][n_train:])\n",
        "\n",
        "\n",
        "class AE:\n",
        "    \"\"\"오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다.\"\"\"\n",
        "\n",
        "    def __init__(self, cols=AE_COLS, q=0.99, seed=0):\n",
        "        self.cols, self.q, self.seed = cols, q, seed\n",
        "\n",
        "    def fit(self, healthy):\n",
        "        self.sc = StandardScaler().fit(healthy[self.cols])\n",
        "        Z = self.sc.transform(healthy[self.cols])\n",
        "        self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)\n",
        "        self.thr = float(np.quantile(self.err(healthy), self.q))\n",
        "        return self\n",
        "\n",
        "    def contrib(self, d):\n",
        "        Z = self.sc.transform(d[self.cols])\n",
        "        return (Z - self.net.predict(Z)) ** 2                 # 입력 항목별 재구성 오차(이상 기여도)\n",
        "\n",
        "    def err(self, d):\n",
        "        return self.contrib(d).mean(axis=1)\n",
        "\n",
        "\n",
        "def latch(d, e, thr, k=2):\n",
        "    \"\"\"같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다.\"\"\"\n",
        "    out, e = np.zeros(len(d), bool), np.asarray(e)\n",
        "    for idx in d.reset_index(drop=True).groupby(\"unit\").indices.values():\n",
        "        run, on = 0, False\n",
        "        for i in idx:\n",
        "            run = run + 1 if e[i] > thr else 0\n",
        "            on = on or run >= k\n",
        "            out[i] = on\n",
        "    return out\n",
        "\n",
        "units, snap, W = make_runs()\n",
        "D = build_features(snap, W)\n",
        "tr_u, va_u = split_units(units)\n",
        "tr, va = D[D[\"unit\"].isin(tr_u)].reset_index(drop=True), D[D[\"unit\"].isin(va_u)].reset_index(drop=True)\n",
        "Wva = W[D[\"unit\"].isin(va_u).values]\n",
        "ae = AE().fit(tr[tr[\"deg\"] == 0])\n",
        "KEEP = [\"unit\", \"day\", \"load\", \"ambient\", \"temp\", \"press\"]\n",
        "\n",
        "\n",
        "def score_batch(w, meta):\n",
        "    \"\"\"스트림 한 묶음: 파형 -> 특징 -> 이상 점수(배수). 점수는 개체 연속성(래치)과 무관한 단건 값.\"\"\"\n",
        "    f = extract(w)\n",
        "    d = pd.DataFrame({\"rms_n\": f[\"rms\"] / (0.5 + meta[\"load\"]), \"kurt\": f[\"kurt\"], \"b150\": f[\"b150\"], \"b165\": f[\"b165\"], \"b200\": f[\"b200\"],\n",
        "                      \"hf_ratio\": (f[\"b150\"] + f[\"b165\"] + f[\"b200\"]) / (f[\"b50\"] + 1e-6),\n",
        "                      \"temp_res\": meta[\"temp\"] - (55 + 22 * meta[\"load\"] + 0.5 * meta[\"ambient\"]), \"press_res\": meta[\"press\"] - (450 - 12 * meta[\"load\"])})\n",
        "    return ae.err(d) / ae.thr\n",
        "\n",
        "\n",
        "def valid(w, meta):\n",
        "    \"\"\"불량 패킷 판정: 결측(NaN), 센서 고정(분산 0), 포화(|값|>12), 값 범위 이탈.\"\"\"\n",
        "    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)\n",
        "    bad |= ~meta[\"temp\"].between(20, 160).values | ~meta[\"press\"].between(100, 600).values\n",
        "    return ~bad\n",
        "\n",
        "\n",
        "def make_stream(n, bad_rate, seed=0):\n",
        "    rng = np.random.default_rng(seed)\n",
        "    pick = rng.integers(0, len(va), n)\n",
        "    w = Wva[pick] + rng.normal(0, 0.02, (n, NS))\n",
        "    meta = va.iloc[pick][KEEP + [\"deg\"]].reset_index(drop=True)\n",
        "    kind = np.array([\"\"] * n, dtype=object)\n",
        "    bad = rng.random(n) < bad_rate\n",
        "    k = rng.integers(0, 4, n)\n",
        "    for t, nm in enumerate([\"결측 패킷\", \"센서 고정\", \"포화\", \"범위 이탈\"]):\n",
        "        m = bad & (k == t)\n",
        "        kind[m] = nm\n",
        "        if t == 0: w[m, 10:40] = np.nan\n",
        "        if t == 1: w[m] = 0.5\n",
        "        if t == 2: w[m, ::7] = 50.0\n",
        "        if t == 3: meta.loc[m, \"temp\"] = 400.0\n",
        "    return w, meta, kind\n",
        "\n",
        "\n",
        "# 1) 1,000건 동시 주입\n",
        "t0 = time.perf_counter()\n",
        "w, meta, kind = make_stream(1000, 0.06)\n",
        "ok = valid(w, meta)\n",
        "s = np.full(1000, np.nan)\n",
        "s[ok] = score_batch(w[ok], meta[ok].reset_index(drop=True))\n",
        "sec = time.perf_counter() - t0\n",
        "inj = int((kind != \"\").sum())\n",
        "print(f\"[교육용가상자료] 1,000건 동시 주입: 정상 처리 {int(ok.sum())}건, 격리 {int((~ok).sum())}건, 유실 {1000 - int(ok.sum()) - int((~ok).sum())}건, 처리 시간 {sec:.2f}초 (목표 2초 이내: {'달성' if sec <= 2 else '미달'})\")\n",
        "print(f\"  불량 주입 {inj}건 중 격리 {int(((kind != '') & ~ok).sum())}건 (\" + \", \".join(f\"{k} {int((kind == k).sum())}\" for k in [\"결측 패킷\", \"센서 고정\", \"포화\", \"범위 이탈\"]) + \")\")\n",
        "\n",
        "# 2) 경보 정확도: 정상 구간과 열화가 뚜렷한 구간(deg>=0.3)을 판정 대상으로 하고, 열화 초기 회색 구간은 정확도에서 제외하고 따로 표시\n",
        "al = s > 1.0\n",
        "y = (meta[\"deg\"] >= 0.3).values\n",
        "eval_m = ok & ((meta[\"deg\"] == 0).values | y)\n",
        "acc = (al[eval_m] == y[eval_m]).mean()\n",
        "tp, fp, fn = (al & y & ok).sum(), (al & ~y & eval_m).sum(), (~al & y & ok).sum()\n",
        "print(f\"경보 정확도 {acc:.1%} (목표 95% 이상: {'달성' if acc >= 0.95 else '미달'}), 재현율 {tp / (tp + fn):.1%}, 정밀도 {tp / (tp + fp):.1%}, 경보 {int((al & ok).sum())}건\")\n",
        "gray = ok & (meta[\"deg\"] > 0).values & (meta[\"deg\"] < 0.3).values\n",
        "print(f\"  참고: 열화 초기(0~0.3) 구간 {int(gray.sum())}건 중 경보 {al[gray].mean():.0%}(회색 구간은 정확도에서 제외)\")\n",
        "\n",
        "# 3) 상태 화면 구성: 최근 30일 x 24대 요약(개체별 최댓값·경보 여부) 집계 시간\n",
        "t0 = time.perf_counter()\n",
        "last = va[va[\"day\"] >= va.groupby(\"unit\")[\"day\"].transform(\"max\") - 30]\n",
        "sc = score_batch(Wva[last.index.values], last.reset_index(drop=True))\n",
        "view = last.assign(score=sc).groupby(\"unit\").agg(max_score=(\"score\", \"max\"), last_score=(\"score\", \"last\"))\n",
        "view[\"상태\"] = np.where(view[\"last_score\"] > 1, \"경보\", np.where(view[\"max_score\"] > 1, \"주의\", \"정상\"))\n",
        "tv = time.perf_counter() - t0\n",
        "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()))\n",
        "\n",
        "# 4) 규모 확대: 배치 처리량과 메모리 한계 대응(배치 크기 5,000)\n",
        "print(\"\\n[주입 규모별 처리] 건수 | 시간 | 처리량 | 격리 | 유실\")\n",
        "sizes = (1000, 10000, 50000, 100000)\n",
        "ts = []\n",
        "for n in sizes:\n",
        "    w, meta, kind = make_stream(n, 0.06, seed=n)\n",
        "    t0 = time.perf_counter(); done = quarantine = 0\n",
        "    for i in range(0, n, 5000):\n",
        "        wb, mb = w[i:i + 5000], meta.iloc[i:i + 5000].reset_index(drop=True)\n",
        "        okb = valid(wb, mb)\n",
        "        if okb.any():\n",
        "            score_batch(wb[okb], mb[okb].reset_index(drop=True))\n",
        "        done += int(okb.sum()); quarantine += int((~okb).sum())\n",
        "    t = time.perf_counter() - t0; ts.append(t)\n",
        "    print(f\" {n:7d}  {t:5.2f}초  {n / t:8.0f}건/초  {quarantine:5d}  {n - done - quarantine}\")\n",
        "print(\"운용 기준: 전체 격리율이 20%를 넘으면 센서 계통 점검 알림, 처리 지연 시 FFT 대역 5개만 계산하는 축소 모드로 전환\")\n",
        "\n",
        "plt.figure(figsize=(6, 4))\n",
        "plt.plot(sizes, ts, marker=\"o\"); plt.xscale(\"log\"); plt.xlabel(\"stream size\"); plt.ylabel(\"seconds\"); plt.title(\"Stream stress test\")\n",
        "plt.tight_layout(); plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 08 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "[교육용가상자료] 1,000건 동시 주입: 정상 처리 945건, 격리 55건, 유실 0건, 처리 시간 0.03초 (목표 2초 이내: 달성)\n",
        "  불량 주입 55건 중 격리 55건 (결측 패킷 15, 센서 고정 16, 포화 15, 범위 이탈 9)\n",
        "경보 정확도 98.0% (목표 95% 이상: 달성), 재현율 98.0%, 정밀도 97.3%, 경보 330건\n",
        "  참고: 열화 초기(0~0.3) 구간 231건 중 경보 15%(회색 구간은 정확도에서 제외)\n",
        "상태 화면 구성(24대 x 최근 30일 744건) 0.02초 (목표 2초 이내: 달성), 경보 24대\n",
        "\n",
        "[주입 규모별 처리] 건수 | 시간 | 처리량 | 격리 | 유실\n",
        "    1000   0.02초     48180건/초     68  0\n",
        "   10000   0.21초     48054건/초    618  0\n",
        "   50000   1.23초     40710건/초   2951  0\n",
        "  100000   3.51초     28472건/초   5940  0\n",
        "운용 기준: 전체 격리율이 20%를 넘으면 센서 계통 점검 알림, 처리 지연 시 FFT 대역 5개만 계산하는 축소 모드로 전환\n",
        "```\n",
        "</details>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 09. AAR 기반 예지정비 오경보 분석 및 이상 감지 임계 파라미터 재보정\n",
        "**소요시간**: 약 3.0시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "AAR 기반 예지정비 오경보를 분석하고 임계 파라미터를 재보정한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **AAR 기반으로 예지정비 오경보를 원인별로 분석하고 이상 감지 임계 파라미터를 재보정해 효과를 지표로 산출한다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **AAR 기간 오경보를 고부하 운용일, 외기 급변, 센서 스파이크, 열화 초기 경계 등으로 분류하고, 정상 분위수·연속 일수·고부하일 임계 완화를 AAR 기간 자료만으로 격자 탐색해 재보정 수행하고 검증하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 4개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `reports/mnt04_roi_analysis.md` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 추정 산식의 객관성 및 군 재정/군수 기획 부서 설득력\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: 비용 절감액, 시간 단축률 등 3대 정량 효과 지표 산출\n",
        "\n",
        "\n",
        "- 진행 기준: ROI 분석 보고서 완성 및 주요 지표 민감도 분석 수록\n",
        "\n",
        "\n",
        "- 검증 방법: 회계/기획 기준 모델과의 비교 검증\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"AAR 기간의 예지정비 오경보를 원인(고부하 운용일·외기 급변·센서 스파이크·기타)별로 분석하고, 정상 분위수·연속 일수·고부하일 임계 완화 파라미터를 AAR 기간 자료로만 격자 탐색해 재보정한 뒤 검증 기간에서 오경보·고장 전 확보 시간 변화와 비용 절감·점검 시간 단축·확보 시간 증가의 3대 정량 효과 지표를 가정 단가로 산출한다.\"\"\"\n",
        "import itertools\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "from scipy.stats import kurtosis\n",
        "from sklearn.neural_network import MLPRegressor\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)\n",
        "# ---------------------------------------------------------------------------\n",
        "BASE = [\"rms\", \"kurt\", \"crest\", \"p2p\", \"b50\", \"b100\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "ROLL = [\"rms\", \"kurt\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "\n",
        "\n",
        "def build_features(snap, W):\n",
        "    D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)\n",
        "    D[\"temp_res\"] = D[\"temp\"] - (55 + 22 * D[\"load\"] + 0.5 * D[\"ambient\"])        # 부하·외기 대비 온도 잔차\n",
        "    D[\"press_res\"] = D[\"press\"] - (450 - 12 * D[\"load\"])                           # 부하 대비 압력 잔차\n",
        "    D[\"rms_n\"] = D[\"rms\"] / (0.5 + D[\"load\"])                                      # 부하 정규화 RMS\n",
        "    D[\"hf_ratio\"] = (D[\"b150\"] + D[\"b165\"] + D[\"b200\"]) / (D[\"b50\"] + 1e-6)         # 고주파/기본주파 비\n",
        "    g = D.groupby(\"unit\")\n",
        "    for c in ROLL:\n",
        "        D[c + \"_m3\"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())\n",
        "    for c in (\"rms_n\", \"temp_res\", \"press_res\"):\n",
        "        D[c + \"_d\"] = g[c].diff().fillna(0)\n",
        "    return D\n",
        "\n",
        "\n",
        "DERIVED = [\"temp_res\", \"press_res\", \"rms_n\", \"hf_ratio\"] + [c + \"_m3\" for c in ROLL] + [\"rms_n_d\", \"temp_res_d\", \"press_res_d\"]\n",
        "FEATS = BASE + DERIVED\n",
        "AE_COLS = [\"rms_n\", \"kurt\", \"b150\", \"b165\", \"b200\", \"hf_ratio\", \"temp_res\", \"press_res\"]\n",
        "\n",
        "\n",
        "def split_units(units, n_train=24):\n",
        "    return list(units[\"unit\"][:n_train]), list(units[\"unit\"][n_train:])\n",
        "\n",
        "\n",
        "class AE:\n",
        "    \"\"\"오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다.\"\"\"\n",
        "\n",
        "    def __init__(self, cols=AE_COLS, q=0.99, seed=0):\n",
        "        self.cols, self.q, self.seed = cols, q, seed\n",
        "\n",
        "    def fit(self, healthy):\n",
        "        self.sc = StandardScaler().fit(healthy[self.cols])\n",
        "        Z = self.sc.transform(healthy[self.cols])\n",
        "        self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)\n",
        "        self.thr = float(np.quantile(self.err(healthy), self.q))\n",
        "        return self\n",
        "\n",
        "    def contrib(self, d):\n",
        "        Z = self.sc.transform(d[self.cols])\n",
        "        return (Z - self.net.predict(Z)) ** 2                 # 입력 항목별 재구성 오차(이상 기여도)\n",
        "\n",
        "    def err(self, d):\n",
        "        return self.contrib(d).mean(axis=1)\n",
        "\n",
        "\n",
        "def latch(d, e, thr, k=2):\n",
        "    \"\"\"같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다.\"\"\"\n",
        "    out, e = np.zeros(len(d), bool), np.asarray(e)\n",
        "    for idx in d.reset_index(drop=True).groupby(\"unit\").indices.values():\n",
        "        run, on = 0, False\n",
        "        for i in idx:\n",
        "            run = run + 1 if e[i] > thr else 0\n",
        "            on = on or run >= k\n",
        "            out[i] = on\n",
        "    return out\n",
        "\n",
        "FA_COST, FA_HOURS, LATE_COST = 1.0, 3.0, 0.15      # [교육용가상자료] 오경보 1회 점검 비용지수, 점검 소요시간(h), 확보 시간 1일 손실 비용지수\n",
        "units, snap, W = make_runs(seed=5)\n",
        "D = build_features(snap, W)\n",
        "order = list(units[\"unit\"])\n",
        "tr_u, aar_u, val_u = order[:24], order[24:36], order[36:]\n",
        "tr = D[D[\"unit\"].isin(tr_u)]\n",
        "ae = AE().fit(tr[tr[\"deg\"] == 0])\n",
        "healthy_err = ae.err(tr[tr[\"deg\"] == 0])\n",
        "E = {u: (D[D[\"unit\"] == u].reset_index(drop=True), ae.err(D[D[\"unit\"] == u])) for u in order}\n",
        "\n",
        "\n",
        "def evaluate(us, q, k, relax):\n",
        "    \"\"\"반환: 오경보 에피소드 수, 원인 분류 건수, 고장 전 확보 시간(일) 목록, 미탐 개체 수, 손실 확보 일수\"\"\"\n",
        "    base_thr = float(np.quantile(healthy_err, q))\n",
        "    fa, causes, leads, miss, lost = 0, {\"고부하 운용일\": 0, \"외기 급변·고온\": 0, \"센서 스파이크\": 0, \"열화 초기 경계\": 0, \"기타\": 0}, [], 0, 0.0\n",
        "    for u in us:\n",
        "        g, e = E[u]\n",
        "        thr = base_thr * np.where(g[\"load\"].values > 0.95, 1 + relax, 1.0)\n",
        "        over = e > thr\n",
        "        run = np.zeros(len(g), int)\n",
        "        for i in range(len(g)):\n",
        "            run[i] = run[i - 1] + 1 if over[i] and i else int(over[i])\n",
        "        al = run >= k\n",
        "        an = (g[\"deg\"] >= 0.1).values\n",
        "        start = al & ~np.r_[False, al[:-1]]\n",
        "        for i in np.where(start & ~an)[0]:\n",
        "            fa += 1\n",
        "            r = g.iloc[i]\n",
        "            c = \"고부하 운용일\" if r[\"hi_load\"] else (\"외기 급변·고온\" if r[\"ambient\"] > 24 else (\"센서 스파이크\" if r[\"kurt\"] > 3.2 else (\"열화 초기 경계\" if r[\"deg\"] > 0 else \"기타\")))\n",
        "            causes[c] += 1\n",
        "        first_an = np.argmax(an)\n",
        "        hit = np.where(al & an)[0]\n",
        "        ideal = g[\"rul\"].iloc[first_an]\n",
        "        if len(hit):\n",
        "            leads.append(int(g[\"rul\"].iloc[hit[0]])); lost += ideal - g[\"rul\"].iloc[hit[0]]\n",
        "        else:\n",
        "            miss += 1; lost += ideal\n",
        "    return fa, causes, leads, miss, lost\n",
        "\n",
        "\n",
        "def J(res):\n",
        "    return FA_COST * res[0] + LATE_COST * res[4]\n",
        "\n",
        "\n",
        "cur = dict(q=0.99, k=1, relax=0.0)               # 현행(AAR 이전) 파라미터\n",
        "r_aar = evaluate(aar_u, **cur)\n",
        "print(f\"[교육용가상자료] AAR 기간 {len(aar_u)}대: 현행 파라미터(정상 {cur['q']:.0%} 분위, {cur['k']}일, 고부하 완화 {cur['relax']:.0%}) 오경보 {r_aar[0]}건\")\n",
        "tot = sum(r_aar[1].values())\n",
        "for c, v in sorted(r_aar[1].items(), key=lambda x: -x[1]):\n",
        "    print(f\"  {c:10s} {v:3d}건 ({v / max(tot, 1):.0%})\")\n",
        "\n",
        "# 재보정: AAR 기간 자료만으로 격자 탐색 (미탐 개체가 없는 조합 중 비용 최소)\n",
        "grid = []\n",
        "for q, k, rl in itertools.product((0.98, 0.99, 0.995, 0.999), (1, 2, 3), (0.0, 0.3, 0.6)):\n",
        "    r = evaluate(aar_u, q, k, rl)\n",
        "    if r[3] == 0:\n",
        "        grid.append((J(r), q, k, rl))\n",
        "grid.sort()\n",
        "best = dict(q=grid[0][1], k=grid[0][2], relax=grid[0][3])\n",
        "print(f\"\\n재보정 결과(격자 {4 * 3 * 3}개 중 미탐 0개 조합 {len(grid)}개 탐색): 정상 {best['q']:.1%} 분위, 연속 {best['k']}일, 고부하일 임계 {best['relax']:.0%} 완화\")\n",
        "\n",
        "rows = []\n",
        "for nm, us in ((\"AAR 기간(튜닝)\", aar_u), (\"검증 기간(미사용)\", val_u)):\n",
        "    a, b = evaluate(us, **cur), evaluate(us, **best)\n",
        "    rows.append((nm, a, b))\n",
        "    print(f\" {nm:14s} 오경보 {a[0]:3d} -> {b[0]:3d}건, 고장 전 확보 중앙값 {np.median(a[2]):.0f} -> {np.median(b[2]):.0f}일, 미탐 {a[3]} -> {b[3]}대\")\n",
        "\n",
        "# 3대 정량 효과 지표(검증 3개월 대신 검증 12대 기준, 연환산은 개체 수 비례)\n",
        "a, b = rows[1][1], rows[1][2]\n",
        "jc, jn = J(a), J(b)\n",
        "print(\"\\n[3대 정량 효과 지표] (검증 12대, 가정 단가)\")\n",
        "print(f\" 1) 비용 절감: 비용지수 {jc:.1f} -> {jn:.1f} ({1 - jn / jc:.1%} 절감), 48대 연환산 약 {(jc - jn) * 4:.0f}\")\n",
        "print(f\" 2) 시간 단축: 불필요 점검 {a[0] * FA_HOURS:.0f}h -> {b[0] * FA_HOURS:.0f}h ({1 - b[0] / max(a[0], 1):.1%} 단축)\")\n",
        "print(f\" 3) 확보 시간: 고장 전 첫 경보 확보 중앙값 {np.median(a[2]):.0f}일 -> {np.median(b[2]):.0f}일, 정상 오경보 {a[0]}건 -> {b[0]}건\")\n",
        "print(\"[도입 로드맵] 단기(1~3월): 재보정 임계 적용 + 검사관 확인 병행 / 중기(3~12월): 분기 AAR로 임계 재보정, 부하 보정 항 추가 / 장기(1년~): 군수 정보체계 정비 이력과 연동해 임계 자동 갱신\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 09 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "[교육용가상자료] AAR 기간 12대: 현행 파라미터(정상 99% 분위, 1일, 고부하 완화 0%) 오경보 9건\n",
        "  기타           5건 (56%)\n",
        "  고부하 운용일      2건 (22%)\n",
        "  열화 초기 경계     2건 (22%)\n",
        "  외기 급변·고온     0건 (0%)\n",
        "  센서 스파이크      0건 (0%)\n",
        "\n",
        "재보정 결과(격자 36개 중 미탐 0개 조합 36개 탐색): 정상 99.5% 분위, 연속 1일, 고부하일 임계 30% 완화\n",
        " AAR 기간(튜닝)     오경보   9 ->   4건, 고장 전 확보 중앙값 58 -> 56일, 미탐 0 -> 0대\n",
        " 검증 기간(미사용)     오경보  16 ->  10건, 고장 전 확보 중앙값 60 -> 57일, 미탐 0 -> 0대\n",
        "\n",
        "[3대 정량 효과 지표] (검증 12대, 가정 단가)\n",
        " 1) 비용 절감: 비용지수 24.7 -> 22.0 (10.9% 절감), 48대 연환산 약 11\n",
        " 2) 시간 단축: 불필요 점검 48h -> 30h (37.5% 단축)\n",
        " 3) 확보 시간: 고장 전 첫 경보 확보 중앙값 60일 -> 57일, 정상 오경보 16건 -> 10건\n",
        "[도입 로드맵] 단기(1~3월): 재보정 임계 적용 + 검사관 확인 병행 / 중기(3~12월): 분기 AAR로 임계 재보정, 부하 보정 항 추가 / 장기(1년~): 군수 정보체계 정비 이력과 연동해 임계 자동 갱신\n",
        "```\n",
        "</details>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## 🎯 미션 10. 항공/기갑 장비 예지정비(PdM) 센서 데이터 이상 징후 감지 시스템 통합 완성\n",
        "**소요시간**: 약 3.5시간\n",
        "\n",
        "### 📋 미션 배경 및 수행 지침\n",
        "#### 미션 배경 및 당위성\n",
        "\n",
        "항공/기갑 장비 예지정비(PdM) 센서 데이터 이상 징후 감지 시스템을 완성한다. 실제 전술 현장에서는 단순 탐지에 그치지 않고 신속한 위협 판단과 신뢰도 높은 의사결정이 필수적이므로, **센서 검증부터 이상 감지, RUL 예측, 승인 대기 입고 제안, 수리부속 소요 예고, 재학습까지 한 파일로 통합하고 자체 점검과 브리핑 요약을 출력하는 도구를 완성한다.** 과업을 완수하여 오판을 방지하고 작전 대응력을 극대화해야 합니다.\n",
        "\n",
        "이 문제를 해결하기 위해 본 미션에서는 **지난달·이번 달 실적으로 챔피언 모델과 도전자 모델을 학습해 검증 F1을 비교하고 비슷하거나 나을 때만 교체 수행하고 검증하는** 실습을 체계적으로 수행합니다.\n",
        "\n",
        "학습자는 실습가이드에 제시된 총 4개 세부 활동(기본 기능 구현 ➔ 결과 지표 검증 ➔ 파라미터 튜닝 및 도전과제)을 순차적으로 진행하며 핵심 판정 기준을 달성합니다.\n",
        "\n",
        "**이번 결과물:** `docs/mnt04_operations_guide.md 및 pipeline/mnt04_retrain.py` 형태의 결과물을 완성한다.\n",
        "\n",
        "### 🔍 평가 및 검증 기준 (PASS 조건)\n",
        "#### 미션 판정\n",
        "\n",
        "**평가 구분:** 필수\n",
        "\n",
        "**핵심 평가영역:** 지속 가능한 운영 가능성 및 국방 군수 정보체계와의 연동 용이성\n",
        "\n",
        "#### 미션 PASS 조건\n",
        "\n",
        "- 정량 기준: 자동 재학습 파이프라인 정상 실행, 최종 문서 100% 완비\n",
        "\n",
        "\n",
        "- 진행 기준: 최종 인도물 검수 및 프로젝트 완료 발표자료 완성\n",
        "\n",
        "\n",
        "- 검증 방법: End-to-End 전체 파이프라인 자동화 실행 및 최종 검수\n",
        "\n",
        "\n",
        "#### 채점자 판정 메모\n",
        "\n",
        "> 단순 코드 실행 여부뿐만 아니라 출력된 수치와 근거 자료가 군 실무 요구사항에 부합하는지 면밀히 확인한다.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\"\"\"센서 파형 검증·정제, 오토인코더 이상 감지 모델의 자동 재학습(챔피언·도전자 비교 후 교체), 재보정 임계 적용 경보, 잔존수명(RUL) 예측, 검사관 승인 대기 입고 제안, 수리부속 소요 예고까지 예지정비(PdM) 이상 감지 시스템을 한 파일로 통합 실행하고 자체 점검·운영 문서 완비 점검·정비 담당 장교 브리핑 요약을 출력한다.\"\"\"\n",
        "import time\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.stats import kurtosis\n",
        "from sklearn.ensemble import GradientBoostingRegressor\n",
        "from sklearn.neural_network import MLPRegressor\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# [교육용가상자료] 가상 기동체 48대의 엔진·구동계 진동·온도·압력 센서 자료 생성기 (실제 장비·부대·정비 실적 없음)\n",
        "# 개체마다 정상 -> 열화 -> 고장(수명종료)까지 하루 1회 측정(진동 파형 256표본 + 온도 + 압력 + 부하)을 만든다.\n",
        "# 열화 유형 4가지: bearing(베어링), gear(기어), oil_pump(오일 펌프), cooling(냉각 팬 모터)\n",
        "# ---------------------------------------------------------------------------\n",
        "FS, NS = 1024, 256                              # 가상 표본화 주파수(Hz), 파형 1개당 표본 수\n",
        "TT = np.arange(NS) / FS\n",
        "MODES = {  # 열화 유형: (발생 확률, 평균 수명(일), 교체 수리부속, 가상 단가지수)\n",
        "    \"bearing\": (0.40, 160, \"주 베어링 세트\", 6.0), \"gear\": (0.25, 175, \"구동 기어 세트\", 9.0),\n",
        "    \"oil_pump\": (0.20, 150, \"오일 펌프\", 4.0), \"cooling\": (0.15, 140, \"냉각 팬 모터\", 3.0)}\n",
        "KIT_COST = 1.0                                  # 모든 정비에 함께 소모되는 씰·필터 키트 단가지수\n",
        "BANDS = {\"b50\": (40, 60), \"b100\": (90, 110), \"b150\": (140, 156), \"b165\": (156, 176), \"b200\": (190, 215)}\n",
        "\n",
        "\n",
        "def make_runs(n=48, seed=1):\n",
        "    \"\"\"반환: units(개체 정보), snap(일별 측정 요약), W(진동 파형 [측정 수, 256])\"\"\"\n",
        "    rng = np.random.default_rng(seed)\n",
        "    names = list(MODES)\n",
        "    mode = rng.choice(names, n, p=[MODES[m][0] for m in names])\n",
        "    life = np.clip(np.array([rng.normal(MODES[m][1], 18) for m in mode]), 90, 230).astype(int)\n",
        "    start = np.sort(rng.integers(0, 240, n))\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)})\n",
        "    rows = []\n",
        "    for u in units.itertuples():\n",
        "        d = np.arange(u.life + 1)\n",
        "        onset = 0.45 * u.life\n",
        "        deg = np.clip((d - onset) / (u.life - onset), 0, 1) ** 1.5\n",
        "        load = np.clip(u.base_load + rng.normal(0, 0.06, len(d)), 0.2, 1.0)\n",
        "        hi = rng.random(len(d)) < 0.07                                            # 고부하 운용일(정상이지만 진동·온도 상승)\n",
        "        load = np.where(hi, rng.uniform(0.95, 1.1, len(d)), load)\n",
        "        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}))\n",
        "    snap = pd.concat(rows, ignore_index=True)\n",
        "    S = len(snap)\n",
        "    snap[\"ambient\"] = 15 + 12 * np.sin(2 * np.pi * snap[\"cal\"] / 365) + rng.normal(0, 2, S)\n",
        "    dg, ld, md = snap[\"deg\"].values, snap[\"load\"].values, snap[\"mode\"].values\n",
        "    isb, isg, iso, isc = (md == \"bearing\"), (md == \"gear\"), (md == \"oil_pump\"), (md == \"cooling\")\n",
        "    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)\n",
        "    snap[\"press\"] = 450 - 12 * ld - 85 * dg * iso - 15 * dg * isc + rng.normal(0, 6, S)\n",
        "    ph = rng.uniform(0, 2 * np.pi, (S, 3))\n",
        "    a1, a2 = 1.0 * (0.5 + ld), 0.4 * (0.5 + ld)\n",
        "    W = (a1[:, None] * np.sin(2 * np.pi * 50 * TT + ph[:, :1]) + a2[:, None] * np.sin(2 * np.pi * 100 * TT + ph[:, 1:2])\n",
        "         + rng.normal(0, 0.35, (S, NS)))\n",
        "    W += (0.9 * dg * isb)[:, None] * np.sin(2 * np.pi * 165 * TT + ph[:, 2:3])\n",
        "    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))\n",
        "    W += (0.5 * dg * iso)[:, None] * np.sin(2 * np.pi * 200 * TT + ph[:, 2:3])\n",
        "    k = rng.poisson(10 * dg * isb)                                                 # 베어링 열화: 충격성 펄스(첨도 상승)\n",
        "    pos, amp = rng.integers(0, NS, (S, 30)), rng.normal(0, 3.0, (S, 30)) * dg[:, None]\n",
        "    amp = amp * (np.arange(30)[None, :] < k[:, None])\n",
        "    np.add.at(W, (np.repeat(np.arange(S), 30), pos.ravel()), amp.ravel())\n",
        "    return units, snap, W\n",
        "\n",
        "\n",
        "def extract(W):\n",
        "    \"\"\"진동 파형 -> 특징량: RMS, 첨도, 파고율, 최대-최소, FFT 대역 에너지(50/100/150/165/200Hz 부근)\"\"\"\n",
        "    x = W - W.mean(axis=1, keepdims=True)\n",
        "    rms = np.sqrt((x ** 2).mean(axis=1))\n",
        "    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)}\n",
        "    A = np.abs(np.fft.rfft(x * np.hanning(NS), axis=1)) / NS\n",
        "    f = np.fft.rfftfreq(NS, 1 / FS)\n",
        "    for nm, (lo, hi) in BANDS.items():\n",
        "        out[nm] = np.sqrt((A[:, (f >= lo) & (f < hi)] ** 2).sum(axis=1))\n",
        "    return pd.DataFrame(out)\n",
        "\n",
        "# ---------------------------------------------------------------------------\n",
        "# 공통 분석 조각: 일별 특징량 표(파형 특징 + 온도·압력 + 파생), 개체 단위 학습/검증 분할(앞선 24대 학습, 뒤 24대 검증)\n",
        "# ---------------------------------------------------------------------------\n",
        "BASE = [\"rms\", \"kurt\", \"crest\", \"p2p\", \"b50\", \"b100\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "ROLL = [\"rms\", \"kurt\", \"b150\", \"b165\", \"b200\", \"temp\", \"press\"]\n",
        "\n",
        "\n",
        "def build_features(snap, W):\n",
        "    D = pd.concat([snap.reset_index(drop=True), extract(W)], axis=1)\n",
        "    D[\"temp_res\"] = D[\"temp\"] - (55 + 22 * D[\"load\"] + 0.5 * D[\"ambient\"])        # 부하·외기 대비 온도 잔차\n",
        "    D[\"press_res\"] = D[\"press\"] - (450 - 12 * D[\"load\"])                           # 부하 대비 압력 잔차\n",
        "    D[\"rms_n\"] = D[\"rms\"] / (0.5 + D[\"load\"])                                      # 부하 정규화 RMS\n",
        "    D[\"hf_ratio\"] = (D[\"b150\"] + D[\"b165\"] + D[\"b200\"]) / (D[\"b50\"] + 1e-6)         # 고주파/기본주파 비\n",
        "    g = D.groupby(\"unit\")\n",
        "    for c in ROLL:\n",
        "        D[c + \"_m3\"] = g[c].transform(lambda s: s.rolling(3, min_periods=1).mean())\n",
        "    for c in (\"rms_n\", \"temp_res\", \"press_res\"):\n",
        "        D[c + \"_d\"] = g[c].diff().fillna(0)\n",
        "    return D\n",
        "\n",
        "\n",
        "DERIVED = [\"temp_res\", \"press_res\", \"rms_n\", \"hf_ratio\"] + [c + \"_m3\" for c in ROLL] + [\"rms_n_d\", \"temp_res_d\", \"press_res_d\"]\n",
        "FEATS = BASE + DERIVED\n",
        "AE_COLS = [\"rms_n\", \"kurt\", \"b150\", \"b165\", \"b200\", \"hf_ratio\", \"temp_res\", \"press_res\"]\n",
        "\n",
        "\n",
        "def split_units(units, n_train=24):\n",
        "    return list(units[\"unit\"][:n_train]), list(units[\"unit\"][n_train:])\n",
        "\n",
        "\n",
        "class AE:\n",
        "    \"\"\"오토인코더(입력 8 -> 10 -> 4 -> 10 -> 8): 정상 자료만 재구성하도록 학습하고, 재구성 오차를 이상 점수로 쓴다.\"\"\"\n",
        "\n",
        "    def __init__(self, cols=AE_COLS, q=0.99, seed=0):\n",
        "        self.cols, self.q, self.seed = cols, q, seed\n",
        "\n",
        "    def fit(self, healthy):\n",
        "        self.sc = StandardScaler().fit(healthy[self.cols])\n",
        "        Z = self.sc.transform(healthy[self.cols])\n",
        "        self.net = MLPRegressor(hidden_layer_sizes=(10, 4, 10), max_iter=800, random_state=self.seed).fit(Z, Z)\n",
        "        self.thr = float(np.quantile(self.err(healthy), self.q))\n",
        "        return self\n",
        "\n",
        "    def contrib(self, d):\n",
        "        Z = self.sc.transform(d[self.cols])\n",
        "        return (Z - self.net.predict(Z)) ** 2                 # 입력 항목별 재구성 오차(이상 기여도)\n",
        "\n",
        "    def err(self, d):\n",
        "        return self.contrib(d).mean(axis=1)\n",
        "\n",
        "\n",
        "def latch(d, e, thr, k=2):\n",
        "    \"\"\"같은 개체에서 이상 점수가 k일 연속 기준을 넘으면 경보를 켜고 유지(래치)한다.\"\"\"\n",
        "    out, e = np.zeros(len(d), bool), np.asarray(e)\n",
        "    for idx in d.reset_index(drop=True).groupby(\"unit\").indices.values():\n",
        "        run, on = 0, False\n",
        "        for i in idx:\n",
        "            run = run + 1 if e[i] > thr else 0\n",
        "            on = on or run >= k\n",
        "            out[i] = on\n",
        "    return out\n",
        "\n",
        "CONFIG = {\"q\": 0.995, \"k\": 1, \"relax\": 0.3, \"retrain_tolerance\": 0.02, \"approve_rul\": 90, \"spare_horizon\": 60}   # 미션 9 재보정 결과\n",
        "OPS_GUIDE = {                                                                 # docs/mnt04_operations_guide.md 요약\n",
        "    \"1. 목적과 범위\": \"엔진·구동계 센서로 이상 징후와 잔존수명을 보조 판단. 입고·정비 결정은 검사관.\",\n",
        "    \"2. 입력 데이터\": \"일 1회 진동 파형·온도·압력·부하. 실제 장비 자료 대신 가상 식별코드 사용.\",\n",
        "    \"3. 일일 운용 절차\": \"수집 -> 검증·격리 -> 특징 추출 -> 이상 점수 -> 경보 -> 입고 제안(대기) -> 검사관 결정 -> 결과 기록.\",\n",
        "    \"4. 승인(HITL) 기준\": \"조기 입고 제안은 검사관 승인 후 확정, 무응답 시 기본 계획(정기 교체) 유지, 사유 기록.\",\n",
        "    \"5. 재학습 절차\": \"분기 1회 신규 정상·고장 실적으로 도전자 학습, 검증 F1이 나빠지지 않을 때만 교체.\",\n",
        "    \"6. 장애 대응\": \"모델 오류 시 정기 교체 주기 기준으로 대체, 손상 센서 자료는 수동 확인 목록으로.\",\n",
        "    \"7. 연동 인터페이스\": \"국방 군수 정보체계와 정비 이력·수리부속 소모 표준 필드(장비코드·분류군·일시)를 CSV/API로 교환하는 안을 차기 과제로 둔다.\",\n",
        "}\n",
        "F10 = FEATS + [\"day\", \"hf_max\", \"temp_res_max\"]\n",
        "\n",
        "\n",
        "def prep(units, snap, W):\n",
        "    \"\"\"검증: 파형 결측·고정·포화, 값 범위 이탈 측정은 격리하고 나머지로 특징표를 만든다.\"\"\"\n",
        "    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\n",
        "    ok = ~bad\n",
        "    D = build_features(snap[ok].reset_index(drop=True), W[ok])\n",
        "    D[\"hf_max\"] = D.groupby(\"unit\")[\"hf_ratio\"].cummax(); D[\"temp_res_max\"] = D.groupby(\"unit\")[\"temp_res\"].cummax()\n",
        "    return D, int(bad.sum())\n",
        "\n",
        "\n",
        "def alarms(D, e, thr, cfg=CONFIG):\n",
        "    over = np.asarray(e) > thr * np.where(D[\"load\"].values > 0.95, 1 + cfg[\"relax\"], 1.0)\n",
        "    out = np.zeros(len(D), bool)\n",
        "    for idx in D.reset_index(drop=True).groupby(\"unit\").indices.values():\n",
        "        run = 0\n",
        "        for i in idx:\n",
        "            run = run + 1 if over[i] else 0\n",
        "            out[i] = run >= cfg[\"k\"]\n",
        "    return out\n",
        "\n",
        "\n",
        "def f1(D, al):\n",
        "    y = (D[\"deg\"] >= 0.1).values\n",
        "    tp, fp, fn = (al & y).sum(), (al & ~y).sum(), (~al & y).sum()\n",
        "    return float(2 * tp / (2 * tp + fp + fn))\n",
        "\n",
        "\n",
        "def retrain(champion, D_old, D_new):\n",
        "    \"\"\"챔피언 AE(지난달 자료) vs 도전자 AE(지난달+이번 달 학습 개체) -> 이번 달 검증 개체 F1 비교.\"\"\"\n",
        "    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:])]\n",
        "    ch = AE(seed=1).fit(pd.concat([D_old[D_old[\"deg\"] == 0], new_tr[new_tr[\"deg\"] == 0]]))\n",
        "    thr_c = float(np.quantile(champion.err(D_old[D_old[\"deg\"] == 0]), CONFIG[\"q\"]))\n",
        "    thr_n = float(np.quantile(ch.err(pd.concat([D_old[D_old[\"deg\"] == 0], new_tr[new_tr[\"deg\"] == 0]])), CONFIG[\"q\"]))\n",
        "    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))\n",
        "    swap = f_new >= f_old - CONFIG[\"retrain_tolerance\"]\n",
        "    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)}\n",
        "\n",
        "\n",
        "def predict_safe(model, D):\n",
        "    try:\n",
        "        return np.maximum(model.predict(D[F10]), 0), \"모델\"\n",
        "    except Exception:\n",
        "        return np.maximum(90 - D[\"day\"].values, 0).astype(float), \"정기 주기 기준 대체\"\n",
        "\n",
        "\n",
        "# ---------------- 1. 통합 실행 ----------------\n",
        "t_start = time.perf_counter()\n",
        "Dold, _ = prep(*make_runs(seed=11))                                  # 지난달 실적\n",
        "Dnew, _ = prep(*make_runs(seed=1))                                   # 이번 달 실적\n",
        "u3, s3, w3 = make_runs(seed=21)                                      # 다음 달 운용(신규 측정)\n",
        "w3[[100, 2000], 20:50] = np.nan                                      # 손상 측정 2건 주입(결측 패킷)\n",
        "Dop, n_bad = prep(u3, s3, w3)\n",
        "champion = AE().fit(Dold[Dold[\"deg\"] == 0])\n",
        "champ_thr = float(np.quantile(champion.err(Dold[Dold[\"deg\"] == 0]), CONFIG[\"q\"]))\n",
        "model_ae, thr, rt = retrain(champion, Dold, Dnew)\n",
        "print(f\"[교육용가상자료] 재학습: {rt} -> {'도전자 교체' if rt['교체'] else '챔피언 유지'}\")\n",
        "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))\n",
        "e = model_ae.err(Dop)\n",
        "Dop[\"alarm\"] = alarms(Dop, e, thr)\n",
        "Dop[\"rul_hat\"], src = predict_safe(gbr, Dop)\n",
        "Dop[\"score\"] = e / thr\n",
        "sec = time.perf_counter() - t_start\n",
        "print(f\"운용 측정 {len(s3)}건: 처리 {len(Dop)}건, 격리 {n_bad}건, 경보 측정 {int(Dop['alarm'].sum())}건, RUL 출처 '{src}', 전체 실행 {sec:.1f}초\")\n",
        "\n",
        "# 입고 제안(대기)과 수리부속 소요 예고: 경보 개체의 마지막 측정 기준\n",
        "PART = {\"b165\": \"주 베어링 세트\", \"kurt\": \"주 베어링 세트\", \"b150\": \"구동 기어 세트\", \"b200\": \"오일 펌프\", \"press_res\": \"오일 펌프\", \"temp_res\": \"냉각 팬 모터\"}\n",
        "first = Dop[Dop[\"alarm\"]].groupby(\"unit\").first()\n",
        "cols = model_ae.cols\n",
        "C = model_ae.contrib(Dop)\n",
        "prop = []\n",
        "for u, r in first.iterrows():\n",
        "    i = Dop.index[(Dop[\"unit\"] == u) & (Dop[\"day\"] == r[\"day\"])][0]\n",
        "    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]))], \"미분류\"),\n",
        "                 \"mode\": r[\"mode\"], \"status\": \"PENDING\" if r[\"rul_hat\"] <= CONFIG[\"approve_rul\"] else \"관찰\"})\n",
        "P = pd.DataFrame(prop)\n",
        "P[\"true_part\"] = P[\"mode\"].map({m: MODES[m][2] for m in MODES})\n",
        "# 수리부속 소요 예고: 예측 잔존 30일 이내로 내려온 개체의 의심 부품\n",
        "need = P[P[\"rul_hat\"] <= CONFIG[\"spare_horizon\"]]\n",
        "\n",
        "# ---------------- 2. 자체 점검 ----------------\n",
        "_, bad_src = predict_safe(None, Dop.head(3))\n",
        "y = (Dop[\"deg\"] >= 0.1).values\n",
        "checks = {\n",
        "    \"측정 유실 0건\": len(Dop) + n_bad == len(s3),\n",
        "    \"손상 측정 2건 격리\": n_bad == 2,\n",
        "    \"재학습 파이프라인 실행\": \"교체\" in rt,\n",
        "    \"RUL 예측값 음수 없음\": bool((Dop[\"rul_hat\"] >= 0).all()),\n",
        "    \"모델 장애 시 정기 주기 대체\": bad_src == \"정기 주기 기준 대체\",\n",
        "    \"열화 개체 미탐 0대\": bool(Dop.groupby(\"unit\").apply(lambda g: (g[\"alarm\"] & (g[\"deg\"] >= 0.1)).any() or not (g[\"deg\"] >= 0.1).any()).all()),\n",
        "    \"경보 개체 모두 제안 또는 관찰 분류\": len(P) == int(Dop[Dop[\"alarm\"]][\"unit\"].nunique()),\n",
        "    \"자동 승인 0건(승인 대상은 대기)\": bool(P[\"status\"].isin([\"PENDING\", \"관찰\"]).all()),\n",
        "}\n",
        "print(f\"\\n자체 점검 {sum(checks.values())}/{len(checks)} 통과, 자동 재학습 파이프라인 {'정상 실행' if checks['재학습 파이프라인 실행'] else '실패'}\")\n",
        "for k, v in checks.items():\n",
        "    if not v:\n",
        "        print(\" 실패:\", k)\n",
        "doc = sum(bool(v.strip()) for v in OPS_GUIDE.values()) / len(OPS_GUIDE)\n",
        "print(f\"운영 문서 {len(OPS_GUIDE)}개 장 작성률 {doc:.0%} (목표 100%: {'달성' if doc == 1 else '미달'}) - \" + \", \".join(k.split(\". \")[1] for k in OPS_GUIDE))\n",
        "\n",
        "# ---------------- 3. 화면과 브리핑 ----------------\n",
        "fig, ax = plt.subplots(1, 2, figsize=(11, 4))\n",
        "for u in list(P[\"unit\"])[:3]:\n",
        "    g = Dop[Dop[\"unit\"] == u]\n",
        "    ax[0].plot(g[\"day\"], g[\"score\"], label=u)\n",
        "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()\n",
        "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)\")\n",
        "plt.tight_layout(); plt.show()\n",
        "\n",
        "lead = P[\"rul\"]\n",
        "print(\"\\n[정비 담당 장교 브리핑 요약] [교육용가상자료]\")\n",
        "print(f\"1. 운용 개체 {len(u3)}대 측정 {len(s3)}건 중 손상 {n_bad}건 격리, 열화 개체 {len(P)}대에서 경보 발령, 고장 전 확보 시간 중앙값 {lead.median():.0f}일(최소 {lead.min()}일)\")\n",
        "print(f\"2. 이상 감지 모델: {'신규 실적으로 교체' if rt['교체'] else '기존 유지'}(검증 F1 {max(rt['챔피언 F1'], rt['도전자 F1']):.2f}), 재보정 임계 적용으로 정상 개체 오경보 {int((Dop['alarm'] & (Dop['deg'] == 0)).sum())}건\")\n",
        "print(f\"3. 조기 입고 제안 {int((P['status'] == 'PENDING').sum())}건은 검사관 승인 대기(자동 승인 0건), 예측 잔존 {CONFIG['approve_rul']}일 초과 {int((P['status'] == '관찰').sum())}건은 관찰 유지\")\n",
        "print(\"4. 수리부속 소요 예고(예측 잔존 \" + str(CONFIG[\"spare_horizon\"]) + \"일 이내 개체의 의심 부품): \" + \", \".join(f\"{k} {v}건\" for k, v in need[\"part\"].value_counts().items())\n",
        "      + f\" (첫 경보 시점 의심 부품 정확도 {(P['part'] == P['true_part']).mean():.0%})\")\n",
        "print(f\"5. RUL 예측: 경보 시점 평균 절대오차 {np.mean(np.abs(P['rul_hat'] - np.minimum(P['rul'], 120))):.1f}일\")\n",
        "print(\"6. 건의: 구동 기어·주 베어링 세트 우선 확보(소모 비용 상위), 국방 군수 정보체계와 정비 이력·수리부속 소모 표준 필드 연동 인터페이스 정의\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<details>\n",
        "<summary><b>💡 [참고] 미션 10 예상 실행 결과 (Output)</b></summary>\n",
        "\n",
        "```text\n",
        "[교육용가상자료] 재학습: {'챔피언 F1': 0.868, '도전자 F1': 0.837, '교체': False} -> 챔피언 유지\n",
        "운용 측정 7695건: 처리 7693건, 격리 2건, 경보 측정 2562건, RUL 출처 '모델', 전체 실행 20.6초\n",
        "\n",
        "자체 점검 8/8 통과, 자동 재학습 파이프라인 정상 실행\n",
        "운영 문서 7개 장 작성률 100% (목표 100%: 달성) - 목적과 범위, 입력 데이터, 일일 운용 절차, 승인(HITL) 기준, 재학습 절차, 장애 대응, 연동 인터페이스\n",
        "\n",
        "[정비 담당 장교 브리핑 요약] [교육용가상자료]\n",
        "1. 운용 개체 48대 측정 7695건 중 손상 2건 격리, 열화 개체 48대에서 경보 발령, 고장 전 확보 시간 중앙값 66일(최소 37일)\n",
        "2. 이상 감지 모델: 기존 유지(검증 F1 0.87), 재보정 임계 적용으로 정상 개체 오경보 23건\n",
        "3. 조기 입고 제안 31건은 검사관 승인 대기(자동 승인 0건), 예측 잔존 90일 초과 17건은 관찰 유지\n",
        "4. 수리부속 소요 예고(예측 잔존 60일 이내 개체의 의심 부품): 구동 기어 세트 6건, 주 베어링 세트 3건, 냉각 팬 모터 3건, 오일 펌프 2건 (첫 경보 시점 의심 부품 정확도 52%)\n",
        "5. RUL 예측: 경보 시점 평균 절대오차 10.0일\n",
        "6. 건의: 구동 기어·주 베어링 세트 우선 확보(소모 비용 상위), 국방 군수 정보체계와 정비 이력·수리부속 소모 표준 필드 연동 인터페이스 정의\n",
        "```\n",
        "</details>\n"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}