310 KiB
310 KiB
🎯 [SUP02] 군수품 시효성 품목 재고통제·잔여일 관리
Google Colab 실습 환경 (보급·재산·시효관리 / Level 2)
📌 1. 프로젝트 개요 및 배경
| 항목 | 내용 |
|---|---|
| 과정 분류 | 보급·재산·시효관리 (Level 2) |
| 학습 소요 시간 | 10개 미션 (5일 / 총 28.5시간) |
| 환경 사양 | Python 3.10+, numpy, pandas, matplotlib, scipy, scikit-learn |
📖 프로젝트 목표
군수품 중에는 정비 주기와 사용 가능 기한이 정해져 있어 시효를 넘기거나 정비 시점을 놓치면 손실이 생기는 품목이 많다. 이 과제는 시효성 품목의 잔여일 관리 문제를 시간 한도가 있는 장비·부품의 잔여 수명(RUL) 관리와 시효성 부품 키트의 선입선출(FEFO) 재고 통제로 연결해 다룬다. 가상 장비 보유 현황의 운용 시간, 환경, 정비 이력으로 생존 분석 기반 수명 모델을 만들고, 정비창 능력·예산·키트 유효기간 제약 아래 12개월 계획정비 일정을 추천한다. 일정 확정과 도태 건의는 군수지휘관 승인 절차로 두고, 대규모 부하 시험, AAR 재보정, 자동 재학습 파이프라인까지 완성해 손실을 줄이는 재고·정비 통제 도구를 만든다.
⚠️ [교육용 가상 데이터셋 원칙]
- 본 노트북의 모든 데이터(제원, 로그, 시나리오 등)는 보안 규정을 준수하여 코드 내 난수 기반 가상 데이터로 생성됩니다.
- 실제 부대 위치, 작전 인가, 실제 전장 정보, 군사 기밀 데이터는 일체 포함되지 않습니다.
In [ ]:
# ============================================================================
# [Google Colab 환경 설정] 필수 라이브러리 설치 및 임포트
# ============================================================================
!pip install -q numpy pandas matplotlib scipy scikit-learn
import os
import sys
import time
import math
import json
import re
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
# 그래프 영문 폰트 기본 설정 (환경 독립성 보장)
plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
print(f"실행 환경 준비 완료 (Python {sys.version.split()[0]})")
In [ ]:
"""장비 운용 시간·기동 거리·정비 이력·창정비 주기와 시효성 부품 키트 로트의 데이터셋 스키마를 설계하고 데이터 사전(피처 20개 이상)을 실제 생성 데이터와 대조하며, 허용값 검증기에 오류 6종을 주입해 검출하고 장비 수명관리 업무 프로세스 맵을 만든다."""
import copy
import numpy as np
import pandas as pd
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
# 1) 데이터 사전 (필드, 구분, 자료형, 허용 범위/값, 설명)
DICT = [
("eq_id", "식별", "str", r"^EQ\d{6}$", "장비 식별코드(가상)"), ("unit", "식별", "str", "U1~U6", "가상 부대 코드"), ("etype", "식별", "str", "truck|generator|trailer|radio", "장비 구분(가상)"),
("env", "환경", "str", "normal|rough|salt|cold", "주 운용 환경(일반/험지/염분/혹한)"), ("rate", "운용", "float", "0.5~20", "일 평균 운용시간(h/일)"),
("hours", "운용", "float", "0~60000", "누적 운용시간(h)"), ("km", "운용", "float", "0~", "누적 기동거리(km)"), ("age_y", "운용", "float", "0.5~25", "사용 연수(년)"),
("n_ovh", "정비", "int", "0~5", "창정비 이력 횟수"), ("h_since", "정비", "float", "0~hours", "마지막 창정비 이후 운용시간(h)"), ("unpl", "정비", "int", "0~", "비계획 정비 건수(긴급 포함)"),
("n_maint", "정비", "int", "0~", "계획 정비 건수"), ("insp", "정비", "float", "0~1.3", "정비 점검 마모 비율(소모 수명/총 수명)"),
("t_fail", "라벨", "float", "0~", "고장까지 운용시간 - 모의 정답(실제 운용에는 없음)"), ("event", "라벨", "int", "0|1", "고장 여부(1=고장 이력, 0=생존)"),
("dur", "라벨", "float", "0~", "관측 운용시간(고장 시점 또는 중도절단)"), ("rul_h", "라벨", "float", "0~", "남은 운용시간 - 모의 정답"), ("rul_d", "라벨", "float", "0~", "남은 일수 - 모의 정답"),
("unpl_per_kh", "파생", "float", "0~", "1000h당 비계획 정비율"), ("km_per_h", "파생", "float", "5~40", "운용시간당 기동거리"), ("k_rough", "파생", "int", "0|1", "험지 더미"),
("k_salt", "파생", "int", "0|1", "염분 더미"), ("k_cold", "파생", "int", "0|1", "혹한 더미"), ("k_truck", "파생", "int", "0|1", "차량 더미"), ("k_gen", "파생", "int", "0|1", "발전기 더미"),
("k_trailer", "파생", "int", "0|1", "트레일러 더미"), ("aging_idx", "파생", "float", "0~2", "환경 가중 노후 지수(1이면 주기 도달)"), ("lin_rul_d", "파생", "float", "-", "선형 감가상각 잔여 일수"),
("hs_frac", "파생", "float", "0~2", "창정비 주기 대비 경과 비율"), ("insp_eta_h", "파생", "float", "0~", "점검 마모로 환산한 총 수명 추정(h)")]
dd = pd.DataFrame(DICT, columns=["필드", "구분", "형", "허용값", "설명"])
print(f"[1] 데이터 사전 {len(dd)}개 필드 " + str(dd["구분"].value_counts().to_dict()))
LOT_SCHEMA = [("lot_id", "str", "시효성 부품 키트 로트 번호"), ("qty", "int", "재고 수량(개)"), ("expiry_d", "int", "유효기한까지 남은 일수"), ("margin", "int", "정비일 기준 최소 유효기간 여유 14일(규칙)")]
print(" 시효성 부품 키트 로트 스키마:", ", ".join(f"{a}({b})" for a, b, _ in LOT_SCHEMA), "- 정비 예정일에 유효기간이 부족한 키트는 쓸 수 없다")
# 2) 실제 생성 데이터와 대조
df = make_fleet(300, seed=1)
feat = add_features(df)
miss = sorted(set(feat.columns) - set(dd["필드"]))
extra = sorted(set(dd["필드"]) - set(feat.columns))
fe = dd[~dd["구분"].isin(["라벨", "식별"])]
print(f"[2] 데이터 사전 대조({len(feat)}행): 사전 누락 컬럼 {miss or '없음'}, 데이터에 없는 필드 {extra or '없음'}")
print(f" 예측 피처 후보(식별·라벨 제외) {len(fe)}개 -> 목표(20개 이상) {'달성' if len(fe) >= 20 else '미달'} / 창정비 규정 주기 {CYCLE:.0f}h(가정)")
# 3) 허용값 검증기와 오류 주입
RULES = {"eq_id": lambda s: s.astype(str).str.match(r"^EQ\d{6}$"), "unit": lambda s: s.isin(UNITS), "env": lambda s: s.isin(ENVS), "rate": lambda s: s.between(0.5, 20),
"hours": lambda s: s.between(0, 60000), "h_since": lambda s: s >= 0, "insp": lambda s: s.between(0, 1.3), "age_y": lambda s: s.between(0.5, 25)}
def validate(d):
bad = {}
for c, f in RULES.items():
n = int((~f(d[c])).sum() + d[c].isna().sum())
if n: bad[c] = n
n = int((d["h_since"] > d["hours"] + 1e-6).sum())
if n: bad["h_since>hours(모순)"] = n
n = int(d["eq_id"].duplicated().sum())
if n: bad["eq_id 중복"] = n
n = int((d["km"] / d["hours"].clip(lower=1) > 60).sum())
if n: bad["km 단위 의심"] = n
return bad
print(f"\n[3] 정상 샘플 {len(df)}건 검증 오류 {sum(validate(df).values())}건")
inj = {"환경 표기 오류": ("env", 0, "험지"), "운용시간 음수": ("hours", 1, -120.0), "미등록 부대": ("unit", 2, "U9"), "점검 마모 비율 범위 초과": ("insp", 3, 2.5),
"창정비 후 시간이 누적보다 큼": ("h_since", 4, 99999.0), "km 단위 오류(m)": ("km", 5, 1.2e7)}
det = 0
for name, (col, i, val) in inj.items():
d2 = copy.deepcopy(df)
d2.loc[i, col] = val
bad = validate(d2)
ok = bool(bad)
det += ok
print(f" - {name}: {'검출' if ok else '미검출'} {bad}")
print(f" 주입 오류 {len(inj)}종 중 {det}종 검출")
# 4) 장비 수명관리 업무 프로세스 맵
proc = [("정비 소요·고장 이력 접수", "정비계획관", "정비 요청, 고장 보고", "고장 이력 표", 30), ("운용시간·기동거리 집계", "운용 담당", "일일 운용 보고", "hours, km, rate", 25),
("정비 이력·점검 입력", "정비 담당", "정비·점검 기록", "n_maint, unpl, insp", 25), ("환경 가중치 적용", "분석 담당", "운용 환경 구분", "aging_idx", 15), ("잔여 수명 예측", "분석 담당", "피처 테이블", "RUL(일), 고장 확률", 40),
("정비창 능력·예산 대조", "정비계획관", "월별 슬롯, 예산, 키트 재고", "정비 일정 후보", 35), ("정비 일정 확정·도태 건의", "군수지휘관", "추천안, 사유", "승인 일정", 30), ("정비 실행·결과 입력", "정비창", "승인 일정", "정비 완료 기록", 20),
("예측 오차 환류(AAR)", "분석 담당", "예측 대 실적", "재보정 보고", 30)]
pm = pd.DataFrame(proc, columns=["단계", "담당", "입력", "산출", "현행(분)"])
pm.index = range(1, len(pm) + 1)
print(f"\n[4] 업무 프로세스 맵 {len(pm)}단계, 현행 합계 {pm['현행(분)'].sum()}분 (가장 긴 단계: {pm.loc[pm['현행(분)'].idxmax(), '단계']})")
print(pm.to_string())
print("\n프로세스 맵 작성 완료(9단계), 피처 목록 20개 이상 -> 목표", "달성" if len(pm) >= 5 and len(fe) >= 20 else "미달")
print("모든 부대·장비·수치는 교육용 가상 자료이며 실제 부대·장비와 무관하다.")
In [ ]:
"""표기·단위·모순·결측·중복 오류를 주입한 가상 장비 운용 기록에 운용 환경 가중치를 반영한 노후화 지표 산출을 포함한 정제 규칙 11종을 적용하고, 긴급 운용·다수 비계획 정비 같은 군수 특수 예외는 삭제하지 않고 보존하며 정합성 검증 통과율과 정제 전후 통계를 확인한다."""
import numpy as np
import pandas as pd
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
truth = make_fleet(220, seed=5)
rng = np.random.default_rng(2)
raw = truth.drop(columns=["t_fail", "event", "dur", "rul_h", "rul_d"]).copy()
raw["urgent_ops"] = False
# 긴급 운용(가동률이 비정상적으로 높지만 실제 사실인 기록): 삭제하면 안 된다
urg = rng.choice(len(raw), 8, replace=False)
raw.loc[urg, "rate"] = rng.uniform(13, 17, 8)
raw.loc[urg, "urgent_ops"] = True
# 오류 주입
errs = {}
i = rng.choice(len(raw), 12, replace=False); raw.loc[i, "env"] = rng.choice(["험지", "ROUGH", "Salt ", "혹한", "일반", "normal "], 12); errs["환경 표기"] = 12
i = rng.choice(len(raw), 8, replace=False); raw.loc[i, "unit"] = rng.choice(["u3", "U 2", "U9", "u5"], 8); errs["부대 코드"] = 8
i = rng.choice(len(raw), 6, replace=False); raw.loc[i, "hours"] = -raw.loc[i, "hours"]; errs["운용시간 음수"] = 6
i = rng.choice(len(raw), 7, replace=False); raw.loc[i, "h_since"] = raw.loc[i, "hours"] * rng.uniform(1.3, 2.0, 7); errs["창정비 후 시간 > 누적"] = 7
i = rng.choice(len(raw), 6, replace=False); raw.loc[i, "km"] = raw.loc[i, "km"] * 1000; errs["km 단위(m)"] = 6
i = rng.choice(len(raw), 9, replace=False); raw.loc[i, "rate"] = np.nan; errs["가동률 결측"] = 9
i = rng.choice(len(raw), 9, replace=False); raw.loc[i, "insp"] = np.nan; errs["점검 마모 결측"] = 9
i = rng.choice(len(raw), 5, replace=False); raw.loc[i, "insp"] = rng.uniform(1.8, 3.0, 5); errs["점검 마모 범위 초과"] = 5
raw = pd.concat([raw, raw.iloc[rng.choice(len(raw), 7, replace=False)]], ignore_index=True); errs["중복 행"] = 7
raw["eq_id"] = raw["eq_id"].where(rng.random(len(raw)) > 0.05, raw["eq_id"].str.lower()); errs["ID 표기"] = int((raw["eq_id"].str.islower()).sum())
print(f"[0] 오염 데이터 {len(raw)}행(원본 {len(truth)}행), 주입 오류: " + ", ".join(f"{k} {v}" for k, v in errs.items()) + f", 긴급 운용 정상 기록 {len(urg)}행")
# ---------------- 정제 규칙 ----------------
log = []
d = raw.copy()
def rule(no, desc, n):
log.append((no, desc, int(n)))
n0 = len(d)
d["eq_id"] = d["eq_id"].str.upper().str.strip()
rule(1, "장비 ID 대문자·공백 통일", (raw["eq_id"] != d["eq_id"]).sum())
d = d.drop_duplicates("eq_id").reset_index(drop=True)
rule(2, "중복 장비 행 제거", n0 - len(d))
d["unit"] = d["unit"].astype(str).str.upper().str.replace(" ", "")
d["env"] = d["env"].astype(str).str.strip()
EMAP = {"험지": "rough", "rough": "rough", "salt": "salt", "혹한": "cold", "cold": "cold", "일반": "normal", "normal": "normal", "염분": "salt"}
env_fix = d["env"].str.lower().map(EMAP)
rule(3, "운용 환경 표기 표준화(한글·대소문자·공백)", (env_fix != d["env"]).sum())
d["env"] = env_fix.fillna("normal")
iso = ~d["unit"].isin(UNITS)
quarantine = d[iso].copy()
d = d[~iso].reset_index(drop=True)
rule(4, "미등록 부대 코드 행 격리(삭제하지 않고 별도 보관)", iso.sum())
neg = d["hours"] < 0
d.loc[neg, "hours"] = -d.loc[neg, "hours"]
rule(5, "누적 운용시간 부호 오류 보정", neg.sum())
bad_hs = d["h_since"] > d["hours"]
d.loc[bad_hs, "h_since"] = d.loc[bad_hs, "hours"]
rule(6, "창정비 후 운용시간 > 누적 운용시간 모순 보정", bad_hs.sum())
unit_km = d["km"] / d["hours"].clip(lower=1) > 100
d.loc[unit_km, "km"] = d.loc[unit_km, "km"] / 1000
rule(7, "기동거리 단위 오류(m -> km) 보정", unit_km.sum())
nan_rate = d["rate"].isna()
d["rate"] = d["rate"].fillna(d.groupby(["etype"])["rate"].transform("median"))
rule(8, "일 평균 운용시간 결측을 장비 구분 중앙값으로 보정", nan_rate.sum())
out_i = d["insp"] > 1.3
d.loc[out_i, "insp"] = np.nan
miss_i = d["insp"].isna()
d["insp_flag"] = np.where(miss_i, "추정", "측정")
d["insp"] = d["insp"].fillna((d["h_since"] / CYCLE * d["env"].map(W0)).clip(0, 1.3))
rule(9, "점검 마모 비율 범위 초과 제거 후 결측을 환경 가중 노후 지수로 대체(추정 표시)", miss_i.sum())
exc = (d["rate"] > 12) | d["urgent_ops"]
d["exception"] = np.where(exc, "긴급 운용 보존", "")
rule(10, "긴급 운용(가동률 12h/일 초과) 예외 표시 후 이상치 삭제 없이 보존", exc.sum())
d["aging_idx"] = d["h_since"] * d["env"].map(W0) / CYCLE
d["eff_hours"] = d["h_since"] * d["env"].map(W0)
rule(11, "환경 가중치 반영 노후화 지표(aging_idx, eff_hours) 산출", len(d))
print("\n[1] 정제 규칙 적용 로그")
print(pd.DataFrame(log, columns=["규칙", "내용", "처리 건수"]).to_string(index=False))
print(f" 정제 규칙 {len(log)}종 적용(목표 10종 이상: {'달성' if len(log) >= 10 else '미달'}), 격리 {len(quarantine)}행은 담당자 확인 대기")
# ---------------- 정합성 검증 ----------------
checks = {"장비 ID 유일": not d["eq_id"].duplicated().any(), "부대 코드 등록": bool(d["unit"].isin(UNITS).all()), "환경 값 4종": bool(d["env"].isin(ENVS).all()),
"운용시간 0 이상": bool((d["hours"] >= 0).all()), "창정비 후 <= 누적": bool((d["h_since"] <= d["hours"] + 1e-6).all()), "기동 속도 5~40km/h": bool((d["km"] / d["hours"].clip(lower=1)).between(5, 40).all()),
"점검 마모 0~1.3": bool(d["insp"].between(0, 1.3).all()), "필수 컬럼 결측 0": int(d[["rate", "insp", "aging_idx"]].isna().sum().sum()) == 0}
print("\n[2] 정합성 검증")
for k, v in checks.items():
print(f" [{'PASS' if v else 'FAIL'}] {k}")
print(f" 통과율 {sum(checks.values()) / len(checks):.0%} (목표 100%)")
# ---------------- 복원 정확도·예외 보존·전후 통계 ----------------
t = truth.set_index("eq_id").loc[d["eq_id"].isin(truth["eq_id"]).pipe(lambda m: d.loc[m, "eq_id"])]
dd = d.set_index("eq_id").loc[t.index]
print(f"\n[3] 복원 정확도(원본과 비교): 누적 운용시간 오차 {np.abs(dd['hours'] - t['hours']).mean():.1f}h, 기동거리 오차 {np.abs(dd['km'] - t['km']).mean():.1f}km, 가동률 오차 {np.abs(dd['rate'] - t['rate'])[dd['exception'] == ''].mean():.3f}h/일(긴급 운용 제외)")
print(f" 긴급 운용 기록 {int((d['exception'] != '').sum())}행 보존(삭제 0행), 점검 마모 추정 대체 {int((d['insp_flag'] == '추정').sum())}행은 추정 표시")
rows = [("운용시간 최솟값(h)", raw["hours"].min(), d["hours"].min(), truth["hours"].min()), ("운용시간 평균(h)", raw["hours"].mean(), d["hours"].mean(), truth["hours"].mean()),
("기동거리 최댓값(km)", raw["km"].max(), d["km"].max(), truth["km"].max()), ("가동률 유효 건수", raw["rate"].notna().sum(), d["rate"].notna().sum(), truth["rate"].notna().sum()),
("가동률 최댓값(h/일)", raw["rate"].max(), d["rate"].max(), truth["rate"].max()), ("점검 마모 최댓값", raw["insp"].max(), d["insp"].max(), truth["insp"].max())]
print("[4] 전후 기술통계(오염 / 정제 후 / 원본)")
print(pd.DataFrame(rows, columns=["항목", "오염", "정제 후", "원본"]).round(2).to_string(index=False))
w = d.groupby("env")["aging_idx"].mean().round(3)
print(" 환경별 평균 노후 지수(가중치 W0 적용):", w.to_dict())
In [ ]:
"""선형 감가상각(규정 창정비 주기 - 경과 운용시간) 기반 잔여 수명(RUL) Baseline을 구현해 오차와 환경별 편향·한계를 확인하고, EDA 차트 8종·상관분석(p-value)·환경별 차이 검정·부대별 위험 파레토로 핵심 노후화 요인을 찾는다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from sklearn.metrics import r2_score
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
df = make_fleet(900, seed=1)
d = add_features(df)
al = d[d.event == 0].copy()
al["err"] = al["lin_rul_d"] - al["rul_d"] # 양수 = Baseline이 수명을 길게 본다(위험한 과대 예측)
big = al["rul_d"] >= 90
print(f"[1] 가상 장비 {len(d)}대(고장 이력 {int(d.event.sum())}대, 생존 {len(al)}대). Baseline RUL(일) = (규정 주기 {CYCLE:.0f}h - 창정비 후 운용시간) / 일 가동률")
print(f" 생존 장비 평가: MAE {mae(al.rul_d, al.lin_rul_d):.0f}일, RMSE {rmse(al.rul_d, al.lin_rul_d):.0f}일, MAPE(참 RUL 90일 이상) {mape(al.rul_d[big], al.lin_rul_d[big]):.1f}%, R2 {r2_score(al.rul_d, al.lin_rul_d):.3f}, 평균 편향 {al.err.mean():+.0f}일")
over = (al["lin_rul_d"] > al["rul_d"] + 30)
print(f" Baseline이 실제보다 30일 넘게 길게 본 장비 {int(over.sum())}대({over.mean():.0%}) - 고장 전에 정비 시점을 놓칠 위험")
fail_before = d[(d.event == 1) & (d.t_fail < CYCLE)]
print(f" 고장 이력 {int(d.event.sum())}대 중 규정 주기({CYCLE:.0f}h) 이전에 고장 난 장비 {len(fail_before)}대({len(fail_before) / d.event.sum():.0%}) -> 고정 주기는 가혹 환경 장비를 지키지 못한다")
# 조건별 오차
conds = {"일반 환경": al.env == "normal", "험지": al.env == "rough", "염분": al.env == "salt", "혹한": al.env == "cold", "비계획 정비율 상위 25%": al.unpl_per_kh >= al.unpl_per_kh.quantile(.75),
"일 가동률 4h 이상": al.rate >= 4, "점검 마모 0.8 이상(노후)": al.insp >= 0.8}
ct = pd.DataFrame([(k, int(m.sum()), mae(al.rul_d[m], al.lin_rul_d[m]), al.err[m].mean()) for k, m in conds.items()], columns=["조건", "대수", "MAE(일)", "평균 편향(일)"])
print("\n[2] 조건별 Baseline 오차(양수 = 과대 예측)\n" + ct.round(0).to_string(index=False))
print(" -> 험지·염분 장비는 수명을 길게(위험), 일반 환경 장비는 짧게(과잉 정비) 보아 환경별 편향 방향이 반대다. 환경 가중치 미반영이 한계다.")
# 상관분석
rows = []
for c in FEATS:
if al[c].nunique() < 2: continue
r, p = stats.spearmanr(al[c], al["err"])
rows.append((c, r, p))
cr = pd.DataFrame(rows, columns=["피처", "스피어만 상관", "p-value"]).assign(abs_r=lambda x: x["스피어만 상관"].abs()).sort_values("abs_r", ascending=False)
sig = cr[(cr["abs_r"] >= 0.2) & (cr["p-value"] < 0.05)]
print(f"\n[3] Baseline 오차(부호 있음: 양수 = 과대 예측)와 유의한 상관(|r|>=0.2, p<0.05) 피처 {len(sig)}개 (목표 5개 이상: {'달성' if len(sig) >= 5 else '미달'})")
print(sig.head(8)[["피처", "스피어만 상관", "p-value"]].round(4).to_string(index=False))
print(" (lin_rul_d·h_since·hs_frac는 Baseline 산식과 직접 겹치므로 해석에서는 환경 더미·점검 마모·기동 강도 같은 독립 요인을 본다)")
hi, lo = al[al.env.isin(["rough", "salt"])]["err"], al[al.env.isin(["normal"])]["err"]
u = stats.mannwhitneyu(hi, lo)
kw = stats.kruskal(*[g["insp_eta_h"].values for _, g in al.groupby("env")])
print(f" 가설검정: 험지·염분 vs 일반 환경 Baseline 편향 평균 {hi.mean():+.0f}일 vs {lo.mean():+.0f}일(Mann-Whitney p={u.pvalue:.1e}), 환경별 점검 환산 총 수명 차이 Kruskal p={kw.pvalue:.1e}")
# 파레토: 부대별 고위험(참 RUL 180일 미만) 장비
risk = al[al.rul_d < 180].groupby("unit").size().sort_values(ascending=False)
cum = risk.cumsum() / risk.sum()
print(f"\n[4] 파레토: 참 RUL 180일 미만 고위험 장비 {int(risk.sum())}대, 상위 {int((cum <= 0.8).sum() + 1)}개 부대가 80% 차지 -> " + ", ".join(f"{u} {n}대" for u, n in risk.items()))
print(" 핵심 발견: 환경 가중치와 정비 이력·점검 마모가 Baseline이 놓치는 요인이며, 이는 M4 피처 설계(환경 더미, 비계획 정비율, 점검 마모)의 근거가 된다.")
fig, ax = plt.subplots(2, 4, figsize=(18, 8))
ax[0, 0].boxplot([d.loc[d.env == e, "insp_eta_h"] for e in ENVS]); ax[0, 0].set_xticklabels(["normal", "rough", "salt", "cold"]); ax[0, 0].set_title("1 Estimated total life by environment (h)")
ax[0, 1].scatter(d.h_since, d.dur, c=d.event, s=6, cmap="coolwarm"); ax[0, 1].set_title("2 Hours since overhaul vs observed (red = failed)")
ax[0, 2].hist(d.rate, 30); ax[0, 2].set_title("3 Daily usage rate (h/day)")
ax[0, 3].scatter(d.h_since, d.unpl, s=6, c=d.env.map({e: i for i, e in enumerate(ENVS)})); ax[0, 3].set_title("4 Unplanned maintenance vs hours")
ax[1, 0].bar(range(len(ct)), ct["평균 편향(일)"]); ax[1, 0].set_xticks(range(len(ct))); ax[1, 0].set_xticklabels([f"C{i + 1}" for i in range(len(ct))]); ax[1, 0].set_title("5 Baseline bias by condition (days)")
ax[1, 1].scatter(al.insp, al.hs_frac, s=6, c=al.env.map({e: i for i, e in enumerate(ENVS)})); ax[1, 1].set_title("6 Inspection wear vs hours fraction")
ax[1, 2].scatter(al.rul_d, al.lin_rul_d, s=6); ax[1, 2].plot([0, 4000], [0, 4000], "k--"); ax[1, 2].set_title("7 True RUL vs Baseline RUL (days)")
ax[1, 3].bar(risk.index, risk.values); ax[1, 3].set_title("8 High-risk equipment by unit (Pareto)")
plt.tight_layout(); plt.show()
print("\nEDA 차트 8종 생성 완료(목표 8종 이상: 달성)")
In [ ]:
"""생존 분석(와이블 AFT 최대우도·Cox 비례위험을 scipy로 직접 구현)으로 환경별 고장 확률 곡선과 기대수명을 추정하고, 환경 가속 계수를 숨은 정답과 비교하며 C-index·기대수명 MAPE·1년 내 고장 확률 보정으로 모델을 검증한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 생존 분석 도구(scipy만 사용): 와이블 AFT(가속 수명) 모형과 Cox 비례위험 모형
# ---------------------------------------------------------------------------
from scipy.optimize import minimize
def aft_x(d):
return np.c_[np.ones(len(d)), d["env"].eq("rough"), d["env"].eq("salt"), d["env"].eq("cold"), d["unpl"] / (d["h_since"] / 1000 + 0.5)].astype(float)
AFT_COLS = ["절편", "험지", "염분", "혹한", "비계획 정비율"]
def fit_aft(d):
"""창정비 후 운용시간 dur(고장이면 event=1, 생존이면 중도절단)에 대한 와이블 AFT 최대우도 추정. 반환: 계수(마지막은 log 형상)."""
x, t, e = aft_x(d), d["dur"].values, d["event"].values
def nll(p):
lam, k = np.exp(x @ p[:-1]), np.exp(p[-1])
z = (t / lam) ** k
return -(e * (np.log(k) - np.log(lam) + (k - 1) * (np.log(t) - np.log(lam))) - z).sum()
r = minimize(nll, np.r_[np.log(8000), np.zeros(x.shape[1] - 1), np.log(5)], method="BFGS")
return r.x
def aft_pred(p, d, q=0.5, h=None):
"""이미 h(운용시간)까지 생존했다는 조건에서 분위수 q 생존 시점(운용시간 h 기준 총 수명)을 돌려준다. q=0.5가 조건부 중앙값."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
h = d["h_since"].values if h is None else h
return lam * ((h / lam) ** k - np.log(q)) ** (1 / k)
def aft_fail_prob(p, d, h0, h1):
"""h0까지 생존한 장비가 h0~h1 사이에 고장 날 조건부 확률."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
s0, s1 = np.exp(-(h0 / lam) ** k), np.exp(-(h1 / lam) ** k)
return 1 - s1 / s0
def fit_cox(d, iters=25):
"""Cox 비례위험(Breslow 근사) 뉴턴-랩슨. 공변량: 험지·염분·혹한·비계획 정비율. 반환: 계수."""
x = aft_x(d)[:, 1:]
o = np.argsort(-d["dur"].values)
x, e = x[o], d["event"].values[o]
b = np.zeros(x.shape[1])
for _ in range(iters):
r = np.exp(x @ b)
S0 = np.cumsum(r)
S1 = np.cumsum(r[:, None] * x, axis=0)
S2 = np.cumsum(r[:, None, None] * x[:, :, None] * x[:, None, :], axis=0)
m1 = S1 / S0[:, None]
g = ((x - m1) * e[:, None]).sum(axis=0)
H = -sum(e[i] * (S2[i] / S0[i] - np.outer(m1[i], m1[i])) for i in np.where(e == 1)[0])
step = np.linalg.solve(H - 1e-6 * np.eye(len(b)), g)
b = b - step
if np.abs(step).max() < 1e-7:
break
return b
def cindex(time, event, risk):
"""일치도 C-index: 위험 점수가 높은 쪽이 먼저 고장 났는가(관측 가능한 쌍만)."""
t, e, r = np.asarray(time), np.asarray(event), np.asarray(risk)
num = den = 0.0
for i in np.where(e == 1)[0]:
m = t > t[i]
den += m.sum()
num += (r[m] < r[i]).sum() + 0.5 * (r[m] == r[i]).sum()
return num / max(den, 1)
# ---------------------------------------------------------------------------
tr = make_fleet(900, seed=1)
te = make_fleet(600, seed=2)
derived = [c for c in FEATS if c in ("unpl_per_kh", "km_per_h", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer", "aging_idx", "lin_rul_d", "hs_frac", "insp_eta_h")]
print(f"[1] 파생 피처 {len(derived)}개(환경·장비 더미, 비계획 정비율, 환경 가중 노후 지수, 점검 환산 수명 등) / 학습 {len(tr)}대(고장 {int(tr.event.sum())}대, 중도절단 {int((1 - tr.event).sum())}대), 평가 {len(te)}대")
# ---------------- 1) 와이블 AFT·Cox 추정 ----------------
p = fit_aft(tr)
b = fit_cox(tr)
k = float(np.exp(p[-1]))
acc = pd.DataFrame({"AFT 수명비 exp(b)": np.exp(p[1:5]), "숨은 정답 수명비": [1 / ENV_W["rough"], 1 / ENV_W["salt"], 1 / ENV_W["cold"], np.nan], "Cox 위험비": np.exp(b)}, index=AFT_COLS[1:])
print(f"[2] 와이블 AFT 형상 k={k:.2f}(1보다 크면 마모형 고장), 기준 수명(일반 환경) {np.exp(p[0]):.0f}h")
print(acc.round(3).to_string())
print(" 환경 가속(수명비 < 1)이 가혹 환경일수록 크며, 숨은 정답(1/가중치)과 비슷하게 복원된다. 비계획 정비율은 수명을 줄이는 방향이다.")
# ---------------- 2) 고장 확률 곡선 ----------------
grid = np.array([5000, 6000, 7000, 8000, 9000, 10000])
typ = pd.DataFrame({"env": ENVS, "unpl": 1.0, "h_since": 1000.0})
typ["unpl"] = 1.0 * (typ["h_since"] / 1000 + 0.5)
lam = np.exp(aft_x(typ) @ p[:-1])
curves = pd.DataFrame({f"{h}h": 1 - np.exp(-(h / lam) ** k) for h in grid}, index=ENVS)
b10 = lam * (-np.log(0.9)) ** (1 / k)
print("\n[3] 창정비 후 운용시간별 누적 고장 확률(비계획 정비율 1.0/1000h 장비)")
print(curves.apply(lambda c: c.map("{:.0%}".format)).to_string())
print(" B10 수명(고장 10% 시점): " + ", ".join(f"{e} {v:.0f}h" for e, v in zip(ENVS, b10)) + f" / 규정 주기 {CYCLE:.0f}h 시점 고장 확률 일반 {curves.loc['normal', '8000h']:.0%}, 험지 {curves.loc['rough', '8000h']:.0%}")
# ---------------- 3) 평가: 기대수명 MAPE, C-index ----------------
zero = np.full(len(te), 1e-6)
life = aft_pred(p, te, 0.5, h=zero)
m1, m0 = mape(te.t_fail, life), mape(te.t_fail, np.full(len(te), CYCLE))
x = aft_x(te)[:, 1:]
ci_cox, ci_aft = cindex(te.dur, te.event, x @ b), cindex(te.dur, te.event, -life)
ci_base = cindex(te.dur, te.event, np.zeros(len(te)))
print(f"\n[4] 평가(숨은 정답 t_fail은 모의 평가에만 사용): 기대수명(중앙값) MAPE {m1:.1f}% (목표 25% 이하: {'달성' if m1 <= 25 else '미달'}), 규정 주기 {CYCLE:.0f}h를 모두에 적용하면 {m0:.1f}%")
print(f" C-index: AFT {ci_aft:.3f}, Cox {ci_cox:.3f}, 환경 무시(무차별) {ci_base:.3f} (0.5 = 무작위)")
by = pd.DataFrame({"환경": te.env, "정답 평균": te.t_fail, "AFT": life, "규정 주기": CYCLE}).groupby("환경").mean().round(0)
print(" 환경별 평균 수명(h):\n" + by.to_string())
# ---------------- 4) 1년 내 고장 확률 보정(생존 장비) ----------------
a = te[te.event == 0].copy()
a["p1y"] = aft_fail_prob(p, a, a["h_since"].values, a["h_since"].values + 365 * a["rate"].values)
a["real"] = (a["rul_d"] < 365).astype(int)
a["bin"] = pd.qcut(a["p1y"], 5, labels=False, duplicates="drop")
cal = a.groupby("bin").agg(예측확률=("p1y", "mean"), 실제비율=("real", "mean"), 대수=("real", "size"))
print("\n[5] 1년 내 고장 확률 보정(생존 장비, 확률 5분위)\n" + cal.round(2).to_string())
brier = float(((a["p1y"] - a["real"]) ** 2).mean())
print(f" Brier 점수 {brier:.3f}(무정보 기준 {float(a['real'].mean() * (1 - a['real'].mean())):.3f}). 점검 마모를 쓰지 않는 생존 모형은 확률의 순위는 맞지만 개별 장비의 시점은 거칠다 -> M5에서 점검 마모 피처로 보강한다.")
# ---------------- 5) 변수 제거 ----------------
def ablate(cols):
keep = [0] + cols
sub_x = aft_x(tr)[:, keep]; t, e = tr["dur"].values, tr["event"].values
from scipy.optimize import minimize
def nll(q):
lm, kk = np.exp(sub_x @ q[:-1]), np.exp(q[-1]); z = (t / lm) ** kk
return -(e * (np.log(kk) - np.log(lm) + (kk - 1) * (np.log(t) - np.log(lm))) - z).sum()
q = minimize(nll, np.r_[np.log(8000), np.zeros(len(keep) - 1), np.log(5)], method="BFGS").x
lm, kk = np.exp(aft_x(te)[:, keep] @ q[:-1]), np.exp(q[-1])
med = lm * np.log(2) ** (1 / kk)
return mape(te.t_fail, med), cindex(te.dur, te.event, -med)
rows = [("전체(환경+비계획 정비율)", *ablate([1, 2, 3, 4])), ("환경 제거", *ablate([4])), ("비계획 정비율 제거", *ablate([1, 2, 3])), ("모두 제거(절편만)", *ablate([]))]
print("\n[6] 변수 제거 실험\n" + pd.DataFrame(rows, columns=["모형", "기대수명 MAPE(%)", "C-index"]).round(3).to_string(index=False))
fig, ax = plt.subplots(1, 3, figsize=(15, 4.3))
hh = np.linspace(2000, 11000, 100)
for e, l in zip(ENVS, lam):
ax[0].plot(hh, 1 - np.exp(-(hh / l) ** k), label=e)
ax[0].axvline(CYCLE, c="k", ls="--"); ax[0].legend(); ax[0].set_title("Failure probability by environment (h since overhaul)")
ax[1].plot(cal["예측확률"], cal["실제비율"], "o-"); ax[1].plot([0, 1], [0, 1], "k--"); ax[1].set_title("1-year failure probability calibration")
ax[2].bar(["rough", "salt", "cold", "unplanned rate"], np.exp(p[1:5])); ax[2].axhline(1, c="k", ls="--"); ax[2].set_title("AFT life ratio vs normal (exp(b))"); ax[2].tick_params(axis="x", rotation=20)
plt.tight_layout(); plt.show()
In [ ]:
"""와이블 AFT와 점검 마모 피처를 쓰는 그래디언트 부스팅을 결합한 잔여 수명(RUL) 예측 엔진을 기준 모델 대비 MAE·RMSE·R2로 검증하고, 정비창 월 능력·예산·시효성 부품 키트 유효기간 제약 아래 정수계획으로 최적 계획정비(창정비) 일정을 추천해 규정 주기 고정 일정과 비용을 비교한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import r2_score
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 생존 분석 도구(scipy만 사용): 와이블 AFT(가속 수명) 모형과 Cox 비례위험 모형
# ---------------------------------------------------------------------------
from scipy.optimize import minimize
def aft_x(d):
return np.c_[np.ones(len(d)), d["env"].eq("rough"), d["env"].eq("salt"), d["env"].eq("cold"), d["unpl"] / (d["h_since"] / 1000 + 0.5)].astype(float)
AFT_COLS = ["절편", "험지", "염분", "혹한", "비계획 정비율"]
def fit_aft(d):
"""창정비 후 운용시간 dur(고장이면 event=1, 생존이면 중도절단)에 대한 와이블 AFT 최대우도 추정. 반환: 계수(마지막은 log 형상)."""
x, t, e = aft_x(d), d["dur"].values, d["event"].values
def nll(p):
lam, k = np.exp(x @ p[:-1]), np.exp(p[-1])
z = (t / lam) ** k
return -(e * (np.log(k) - np.log(lam) + (k - 1) * (np.log(t) - np.log(lam))) - z).sum()
r = minimize(nll, np.r_[np.log(8000), np.zeros(x.shape[1] - 1), np.log(5)], method="BFGS")
return r.x
def aft_pred(p, d, q=0.5, h=None):
"""이미 h(운용시간)까지 생존했다는 조건에서 분위수 q 생존 시점(운용시간 h 기준 총 수명)을 돌려준다. q=0.5가 조건부 중앙값."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
h = d["h_since"].values if h is None else h
return lam * ((h / lam) ** k - np.log(q)) ** (1 / k)
def aft_fail_prob(p, d, h0, h1):
"""h0까지 생존한 장비가 h0~h1 사이에 고장 날 조건부 확률."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
s0, s1 = np.exp(-(h0 / lam) ** k), np.exp(-(h1 / lam) ** k)
return 1 - s1 / s0
def fit_cox(d, iters=25):
"""Cox 비례위험(Breslow 근사) 뉴턴-랩슨. 공변량: 험지·염분·혹한·비계획 정비율. 반환: 계수."""
x = aft_x(d)[:, 1:]
o = np.argsort(-d["dur"].values)
x, e = x[o], d["event"].values[o]
b = np.zeros(x.shape[1])
for _ in range(iters):
r = np.exp(x @ b)
S0 = np.cumsum(r)
S1 = np.cumsum(r[:, None] * x, axis=0)
S2 = np.cumsum(r[:, None, None] * x[:, :, None] * x[:, None, :], axis=0)
m1 = S1 / S0[:, None]
g = ((x - m1) * e[:, None]).sum(axis=0)
H = -sum(e[i] * (S2[i] / S0[i] - np.outer(m1[i], m1[i])) for i in np.where(e == 1)[0])
step = np.linalg.solve(H - 1e-6 * np.eye(len(b)), g)
b = b - step
if np.abs(step).max() < 1e-7:
break
return b
def cindex(time, event, risk):
"""일치도 C-index: 위험 점수가 높은 쪽이 먼저 고장 났는가(관측 가능한 쌍만)."""
t, e, r = np.asarray(time), np.asarray(event), np.asarray(risk)
num = den = 0.0
for i in np.where(e == 1)[0]:
m = t > t[i]
den += m.sum()
num += (r[m] < r[i]).sum() + 0.5 * (r[m] == r[i]).sum()
return num / max(den, 1)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 잔여 수명(RUL, 일) 예측 엔진: 와이블 AFT 조건부 중앙값 + 그래디언트 부스팅(고장 이력 시점 표본 학습)의 평균
# ---------------------------------------------------------------------------
from sklearn.ensemble import GradientBoostingRegressor
def fit_engine(train, seed=0):
p = fit_aft(train)
sn = add_features(snapshots(train, 5, seed))
sn["aft"] = (aft_pred(p, sn) - sn["h_since"]) / sn["rate"]
g = GradientBoostingRegressor(n_estimators=250, max_depth=3, learning_rate=0.06, random_state=0).fit(sn[FEATS + ["aft"]], sn["rul_d"])
return {"p": p, "g": g}
def predict_rul(eng, d):
"""d는 add_features를 거친 표. 반환: (앙상블 RUL, AFT RUL, GBR RUL) - 일 단위."""
aft = ((aft_pred(eng["p"], d) - d["h_since"]) / d["rate"]).values
x = d[FEATS].copy()
x["aft"] = aft
gb = eng["g"].predict(x)
return np.clip(0.5 * aft + 0.5 * gb, 1, None), aft, gb
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 계획정비(창정비) 일정 추천: 월별 정비창 능력·예산·시효성 부품 키트(유효기한) 제약 아래 기대 비용 최소화 (비용 단가는 교육용 가정값, 만원)
# ---------------------------------------------------------------------------
from scipy.optimize import milp, LinearConstraint, Bounds
from scipy.stats import norm
C_PLAN, C_FAIL, C_WASTE_KH, C_KIT = 600.0, 2400.0, 40.0, 120.0 # 계획 창정비, 비계획 고장(긴급 정비 포함), 미사용 수명 1000h당 손실, 시효성 부품 키트
HORIZON, MDAYS = 12, 30 # 계획 12개월, 월 30일
CAP_M, BUDGET_M = 22, 15000.0 # 월별 정비창 능력(대), 월별 예산
LOTS = [(30, 70), (45, 170), (45, 390), (70, 800)] # 시효성 키트 로트(수량, 유효기한까지 남은 일수)
KIT_MARGIN = 14 # 정비일 기준 최소 유효기간 여유(일)
def exp_cost_matrix(pred, rate, month_days):
"""장비 i를 월 m에 정비할 때의 기대 비용. 참 RUL ~ N(pred, sigma), sigma = 35일 + 0.12 x pred (가정). 마지막 열은 계획 기간 밖(미정비)."""
sig = 35 + 0.12 * pred
C = np.zeros((len(pred), len(month_days) + 1))
for j, day in enumerate(month_days):
z = (day - pred) / sig
pf = norm.cdf(z)
unused = np.maximum(pred - day, 0) * rate * 0.6 # 정비를 앞당긴 만큼 버리는 운용시간 h(근사)
C[:, j] = C_PLAN + C_KIT + pf * (C_FAIL - C_PLAN) + C_WASTE_KH * unused / 1000 * (1 - pf)
C[:, -1] = norm.cdf((HORIZON * MDAYS - pred) / sig) * C_FAIL
return C
def plan_milp(pred, rate, cap=CAP_M, budget=BUDGET_M, lots=LOTS, use_kit=True):
"""정수계획으로 장비별 정비 월(0~11) 또는 미정비(-1)를 배정한다. 후보는 예측 RUL이 계획 기간 +90일 이내인 장비."""
md = np.arange(1, HORIZON + 1) * MDAYS
cand = np.where(pred < HORIZON * MDAYS + 90)[0]
n, M = len(cand), HORIZON
C = exp_cost_matrix(pred[cand], rate[cand], md)
c = C.ravel()
A, lo, hi = [], [], []
for i in range(n): # 장비당 한 번만 배정(월 12칸 + 미정비)
row = np.zeros(n * (M + 1)); row[i * (M + 1):(i + 1) * (M + 1)] = 1; A.append(row); lo.append(1); hi.append(1)
for m in range(M): # 월별 능력·예산
row = np.zeros(n * (M + 1)); row[m::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(cap)
A.append(row * (C_PLAN + C_KIT)); lo.append(0); hi.append(budget)
if use_kit: # 시효성 키트: m월 이후 정비 수요는 그 시점에 유효기간이 남은 로트 수량을 넘을 수 없다
for m in range(M):
sup = sum(q for q, e in lots if e >= md[m] + KIT_MARGIN)
row = np.zeros(n * (M + 1))
for mm in range(m, M):
row[mm::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(sup)
r = milp(c, constraints=LinearConstraint(np.array(A), lo, hi), integrality=np.ones(len(c)), bounds=Bounds(0, 1), options={"time_limit": 20})
month = np.full(len(pred), -1)
if r.x is not None:
x = np.round(r.x).reshape(n, M + 1)
assign = np.where(x[:, :M].sum(axis=1) > 0.5, x[:, :M].argmax(axis=1), -1)
month[cand] = assign
return month
def plan_greedy(pred, margin=45, cap=CAP_M, budget=BUDGET_M):
"""간단한 탐욕 배정: 정비 목표일(예측 RUL - margin)이 빠른 순으로 목표일 이전 가장 늦은 빈 월, 없으면 가장 이른 빈 월에 배정."""
cap = min(cap, int(budget // (C_PLAN + C_KIT))) # 월 예산이 허용하는 대수로 능력을 제한
month = np.full(len(pred), -1)
left = np.full(HORIZON, cap)
for i in np.argsort(pred):
if pred[i] > HORIZON * MDAYS + 60: break
tgt = int(np.clip((pred[i] - margin) // MDAYS - 1, 0, HORIZON - 1))
for m in list(range(tgt, -1, -1)) + list(range(tgt + 1, HORIZON)):
if left[m] > 0:
month[i] = m; left[m] -= 1; break
return month
def plan_fixed(lin_rul, cap=CAP_M, budget=BUDGET_M):
"""관행: 규정 주기 도달 순서대로(선형 잔여 일수) 월별 능력 안에서 배정."""
return plan_greedy(np.asarray(lin_rul, float), margin=0, cap=cap, budget=budget)
def evaluate_plan(month, true_rul, rate, cost_kit=True):
"""참 RUL로 계획의 실제 비용을 계산한다: 정비일 이전에 고장 -> 비계획 고장 비용, 아니면 계획 비용 + 미사용 수명 손실."""
day = np.where(month >= 0, (month + 1) * MDAYS, 10 ** 9)
fail = (true_rul < day) & (true_rul < HORIZON * MDAYS)
planned = (month >= 0) & ~fail
waste = np.where(planned, (true_rul - day) * rate * 0.6, 0) * C_WASTE_KH / 1000
cost = fail * C_FAIL + planned * (C_PLAN + (C_KIT if cost_kit else 0)) + waste
return {"total": float(cost.sum()), "fail": int(fail.sum()), "planned": int(planned.sum()), "waste": float(waste.sum()), "unsched_fail": int(((month < 0) & fail).sum())}
def kit_fefo(month, lots=LOTS):
"""유효기한이 빠른 로트부터 쓰는 FEFO로 키트를 배정하고 정비일에 유효기간이 부족한 배정 수와 만료 폐기 수량을 센다."""
q = [[qty, e] for qty, e in sorted(lots, key=lambda t: t[1])]
short = 0
for m in sorted(month[month >= 0]):
day = (m + 1) * MDAYS
for lot in q:
if lot[1] >= day + KIT_MARGIN and lot[0] > 0:
lot[0] -= 1; break
else:
short += 1
expired = sum(l[0] for l in q if l[1] < HORIZON * MDAYS)
return {"short": short, "expired_unused": expired, "left": sum(l[0] for l in q)}
# ---------------------------------------------------------------------------
tr = make_fleet(900, seed=1)
te = add_features(make_fleet(900, seed=3))
eng = fit_engine(tr)
al = te[te.event == 0].reset_index(drop=True)
pred, aft, gb = predict_rul(eng, al)
y = al["rul_d"].values
big = y >= 90
P = {"기준(선형 감가상각)": al["lin_rul_d"].values, "와이블 AFT": aft, "GBR(점검 마모 포함)": gb, "앙상블(제안)": pred}
tab = pd.DataFrame({k: {"MAE(일)": mae(y, v), "RMSE(일)": rmse(y, v), "MAPE(%)": mape(y[big], v[big]), "R2": r2_score(y, v)} for k, v in P.items()}).T
print(f"[1] 생존 장비 {len(al)}대 RUL 예측 성능(학습 {len(tr)}대와 다른 장비, MAPE는 참 RUL 90일 이상)")
print(tab.round(3).to_string())
b, e = tab.loc["기준(선형 감가상각)"], tab.loc["앙상블(제안)"]
ok = 1 - e["MAE(일)"] / b["MAE(일)"] >= .2 and 1 - e["RMSE(일)"] / b["RMSE(일)"] >= .2
print(f" 기준 대비 MAE 개선 {1 - e['MAE(일)'] / b['MAE(일)']:.1%}, RMSE 개선 {1 - e['RMSE(일)'] / b['RMSE(일)']:.1%} (목표 20% 이상: {'달성' if ok else '미달'}), R2 {e['R2']:.3f} (목표 0.85 이상: {'달성' if e['R2'] >= .85 else '미달'})")
# 점검 마모 피처의 효과: 제거 실험
tr_f = add_features(snapshots(tr, 5, 0))
cols_no = [c for c in FEATS if c not in ("insp", "insp_eta_h")]
g_no = GradientBoostingRegressor(n_estimators=250, max_depth=3, learning_rate=0.06, random_state=0).fit(tr_f[cols_no], tr_f["rul_d"])
p_no = g_no.predict(al[cols_no])
print(f" 점검 마모 피처를 빼면 GBR MAE {mae(y, gb):.0f}일 -> {mae(y, p_no):.0f}일, R2 {r2_score(y, gb):.3f} -> {r2_score(y, p_no):.3f} (정비 점검 기록의 가치)")
# 불규칙 상황별 오차
conds = {"일반": al.env == "normal", "험지": al.env == "rough", "염분": al.env == "salt", "혹한": al.env == "cold", "고가동(일 4h 이상)": al.rate >= 4, "비계획 정비율 상위 25%": al.unpl_per_kh >= al.unpl_per_kh.quantile(.75)}
rr = pd.DataFrame([(k, int(m.sum()), mae(y[m.values], P["기준(선형 감가상각)"][m.values]), mae(y[m.values], pred[m.values])) for k, m in conds.items()], columns=["조건", "대수", "기준 MAE", "앙상블 MAE"])
print("\n[2] 조건별 MAE(일)\n" + rr.round(0).to_string(index=False))
# ---------------- 계획정비 일정 추천 ----------------
rate, true_rul = al["rate"].values, al["rul_d"].values
mf, mg = plan_fixed(al["lin_rul_d"].values), plan_greedy(pred)
mm = plan_milp(pred, rate)
pl = {"규정 주기 고정 일정(관행)": mf, "예측 기반 탐욕 배정": mg, "예측 기반 정수계획(제안)": mm}
rows = []
for k, v in pl.items():
r, kf = evaluate_plan(v, true_rul, rate), kit_fefo(v)
rows.append((k, int((v >= 0).sum()), r["total"], r["fail"], r["planned"], kf["short"]))
pt = pd.DataFrame(rows, columns=["일정안", "배정 대수", "실제 총비용(만원)", "비계획 고장", "계획 정비", "키트 유효기간 부족"])
print(f"\n[3] 12개월 계획정비 일정 비교(월 능력 {CAP_M}대, 월 예산 {BUDGET_M:,.0f}만원, 키트 재고 {sum(q for q, _ in LOTS)}개 로트 {len(LOTS)}종, 비용 단가는 교육용 가정값)")
print(pt.round(0).to_string(index=False))
sv = 1 - pt.iloc[2]["실제 총비용(만원)"] / pt.iloc[0]["실제 총비용(만원)"]
print(f" 정수계획 안은 관행 대비 총비용 {sv:.1%} 절감, 비계획 고장 {int(pt.iloc[0]['비계획 고장'])} -> {int(pt.iloc[2]['비계획 고장'])}대")
cnt = np.bincount(mm[mm >= 0], minlength=HORIZON)
print(" 월별 정비 대수(제안):", " ".join(f"{i + 1}월 {c}" for i, c in enumerate(cnt)), f"(능력 {CAP_M}대 이하 확인: {'OK' if cnt.max() <= CAP_M else '초과'}, 월 예산 확인: {'OK' if (cnt * (C_PLAN + C_KIT)).max() <= BUDGET_M else '초과'})")
kf = kit_fefo(mm)
print(f" 시효성 키트 FEFO 배정: 유효기간 부족 {kf['short']}건, 만료 임박 폐기 예상 {kf['expired_unused']}개 -> 유효기한이 빠른 로트부터 쓰도록 정비 월을 맞춤")
print(" 추천 일정은 정비계획관 검토와 군수지휘관 승인(M7) 전까지 확정되지 않는다. 장비·단가는 교육용 가상 자료이다.")
fig, ax = plt.subplots(1, 3, figsize=(15, 4.3))
ax[0].scatter(y, P["기준(선형 감가상각)"], s=6, label="baseline"); ax[0].scatter(y, pred, s=6, label="ensemble"); ax[0].plot([0, 4000], [0, 4000], "k--"); ax[0].legend(); ax[0].set_title("RUL: true vs predicted (days)")
ax[1].bar(range(len(pt)), pt["실제 총비용(만원)"]); ax[1].set_xticks(range(len(pt))); ax[1].set_xticklabels(["fixed", "greedy", "MILP"]); ax[1].set_title("Total cost of 12-month plan (10k KRW)")
ax[2].bar(np.arange(1, HORIZON + 1), cnt); ax[2].axhline(CAP_M, c="r", ls="--"); ax[2].set_title("Overhauls per month (capacity line)")
plt.tight_layout(); plt.show()
In [ ]:
"""부대별 장비 노후도 현황과 권장 정비 캘린더 대시보드를 구현하고, 장비마다 잔여 수명(RUL) 예측의 기여도 상위 3개 요인을 변수군 평균 대체 방식으로 자동 산출해 자연어 설명으로 보여 준다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 생존 분석 도구(scipy만 사용): 와이블 AFT(가속 수명) 모형과 Cox 비례위험 모형
# ---------------------------------------------------------------------------
from scipy.optimize import minimize
def aft_x(d):
return np.c_[np.ones(len(d)), d["env"].eq("rough"), d["env"].eq("salt"), d["env"].eq("cold"), d["unpl"] / (d["h_since"] / 1000 + 0.5)].astype(float)
AFT_COLS = ["절편", "험지", "염분", "혹한", "비계획 정비율"]
def fit_aft(d):
"""창정비 후 운용시간 dur(고장이면 event=1, 생존이면 중도절단)에 대한 와이블 AFT 최대우도 추정. 반환: 계수(마지막은 log 형상)."""
x, t, e = aft_x(d), d["dur"].values, d["event"].values
def nll(p):
lam, k = np.exp(x @ p[:-1]), np.exp(p[-1])
z = (t / lam) ** k
return -(e * (np.log(k) - np.log(lam) + (k - 1) * (np.log(t) - np.log(lam))) - z).sum()
r = minimize(nll, np.r_[np.log(8000), np.zeros(x.shape[1] - 1), np.log(5)], method="BFGS")
return r.x
def aft_pred(p, d, q=0.5, h=None):
"""이미 h(운용시간)까지 생존했다는 조건에서 분위수 q 생존 시점(운용시간 h 기준 총 수명)을 돌려준다. q=0.5가 조건부 중앙값."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
h = d["h_since"].values if h is None else h
return lam * ((h / lam) ** k - np.log(q)) ** (1 / k)
def aft_fail_prob(p, d, h0, h1):
"""h0까지 생존한 장비가 h0~h1 사이에 고장 날 조건부 확률."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
s0, s1 = np.exp(-(h0 / lam) ** k), np.exp(-(h1 / lam) ** k)
return 1 - s1 / s0
def fit_cox(d, iters=25):
"""Cox 비례위험(Breslow 근사) 뉴턴-랩슨. 공변량: 험지·염분·혹한·비계획 정비율. 반환: 계수."""
x = aft_x(d)[:, 1:]
o = np.argsort(-d["dur"].values)
x, e = x[o], d["event"].values[o]
b = np.zeros(x.shape[1])
for _ in range(iters):
r = np.exp(x @ b)
S0 = np.cumsum(r)
S1 = np.cumsum(r[:, None] * x, axis=0)
S2 = np.cumsum(r[:, None, None] * x[:, :, None] * x[:, None, :], axis=0)
m1 = S1 / S0[:, None]
g = ((x - m1) * e[:, None]).sum(axis=0)
H = -sum(e[i] * (S2[i] / S0[i] - np.outer(m1[i], m1[i])) for i in np.where(e == 1)[0])
step = np.linalg.solve(H - 1e-6 * np.eye(len(b)), g)
b = b - step
if np.abs(step).max() < 1e-7:
break
return b
def cindex(time, event, risk):
"""일치도 C-index: 위험 점수가 높은 쪽이 먼저 고장 났는가(관측 가능한 쌍만)."""
t, e, r = np.asarray(time), np.asarray(event), np.asarray(risk)
num = den = 0.0
for i in np.where(e == 1)[0]:
m = t > t[i]
den += m.sum()
num += (r[m] < r[i]).sum() + 0.5 * (r[m] == r[i]).sum()
return num / max(den, 1)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 잔여 수명(RUL, 일) 예측 엔진: 와이블 AFT 조건부 중앙값 + 그래디언트 부스팅(고장 이력 시점 표본 학습)의 평균
# ---------------------------------------------------------------------------
from sklearn.ensemble import GradientBoostingRegressor
def fit_engine(train, seed=0):
p = fit_aft(train)
sn = add_features(snapshots(train, 5, seed))
sn["aft"] = (aft_pred(p, sn) - sn["h_since"]) / sn["rate"]
g = GradientBoostingRegressor(n_estimators=250, max_depth=3, learning_rate=0.06, random_state=0).fit(sn[FEATS + ["aft"]], sn["rul_d"])
return {"p": p, "g": g}
def predict_rul(eng, d):
"""d는 add_features를 거친 표. 반환: (앙상블 RUL, AFT RUL, GBR RUL) - 일 단위."""
aft = ((aft_pred(eng["p"], d) - d["h_since"]) / d["rate"]).values
x = d[FEATS].copy()
x["aft"] = aft
gb = eng["g"].predict(x)
return np.clip(0.5 * aft + 0.5 * gb, 1, None), aft, gb
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 계획정비(창정비) 일정 추천: 월별 정비창 능력·예산·시효성 부품 키트(유효기한) 제약 아래 기대 비용 최소화 (비용 단가는 교육용 가정값, 만원)
# ---------------------------------------------------------------------------
from scipy.optimize import milp, LinearConstraint, Bounds
from scipy.stats import norm
C_PLAN, C_FAIL, C_WASTE_KH, C_KIT = 600.0, 2400.0, 40.0, 120.0 # 계획 창정비, 비계획 고장(긴급 정비 포함), 미사용 수명 1000h당 손실, 시효성 부품 키트
HORIZON, MDAYS = 12, 30 # 계획 12개월, 월 30일
CAP_M, BUDGET_M = 22, 15000.0 # 월별 정비창 능력(대), 월별 예산
LOTS = [(30, 70), (45, 170), (45, 390), (70, 800)] # 시효성 키트 로트(수량, 유효기한까지 남은 일수)
KIT_MARGIN = 14 # 정비일 기준 최소 유효기간 여유(일)
def exp_cost_matrix(pred, rate, month_days):
"""장비 i를 월 m에 정비할 때의 기대 비용. 참 RUL ~ N(pred, sigma), sigma = 35일 + 0.12 x pred (가정). 마지막 열은 계획 기간 밖(미정비)."""
sig = 35 + 0.12 * pred
C = np.zeros((len(pred), len(month_days) + 1))
for j, day in enumerate(month_days):
z = (day - pred) / sig
pf = norm.cdf(z)
unused = np.maximum(pred - day, 0) * rate * 0.6 # 정비를 앞당긴 만큼 버리는 운용시간 h(근사)
C[:, j] = C_PLAN + C_KIT + pf * (C_FAIL - C_PLAN) + C_WASTE_KH * unused / 1000 * (1 - pf)
C[:, -1] = norm.cdf((HORIZON * MDAYS - pred) / sig) * C_FAIL
return C
def plan_milp(pred, rate, cap=CAP_M, budget=BUDGET_M, lots=LOTS, use_kit=True):
"""정수계획으로 장비별 정비 월(0~11) 또는 미정비(-1)를 배정한다. 후보는 예측 RUL이 계획 기간 +90일 이내인 장비."""
md = np.arange(1, HORIZON + 1) * MDAYS
cand = np.where(pred < HORIZON * MDAYS + 90)[0]
n, M = len(cand), HORIZON
C = exp_cost_matrix(pred[cand], rate[cand], md)
c = C.ravel()
A, lo, hi = [], [], []
for i in range(n): # 장비당 한 번만 배정(월 12칸 + 미정비)
row = np.zeros(n * (M + 1)); row[i * (M + 1):(i + 1) * (M + 1)] = 1; A.append(row); lo.append(1); hi.append(1)
for m in range(M): # 월별 능력·예산
row = np.zeros(n * (M + 1)); row[m::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(cap)
A.append(row * (C_PLAN + C_KIT)); lo.append(0); hi.append(budget)
if use_kit: # 시효성 키트: m월 이후 정비 수요는 그 시점에 유효기간이 남은 로트 수량을 넘을 수 없다
for m in range(M):
sup = sum(q for q, e in lots if e >= md[m] + KIT_MARGIN)
row = np.zeros(n * (M + 1))
for mm in range(m, M):
row[mm::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(sup)
r = milp(c, constraints=LinearConstraint(np.array(A), lo, hi), integrality=np.ones(len(c)), bounds=Bounds(0, 1), options={"time_limit": 20})
month = np.full(len(pred), -1)
if r.x is not None:
x = np.round(r.x).reshape(n, M + 1)
assign = np.where(x[:, :M].sum(axis=1) > 0.5, x[:, :M].argmax(axis=1), -1)
month[cand] = assign
return month
def plan_greedy(pred, margin=45, cap=CAP_M, budget=BUDGET_M):
"""간단한 탐욕 배정: 정비 목표일(예측 RUL - margin)이 빠른 순으로 목표일 이전 가장 늦은 빈 월, 없으면 가장 이른 빈 월에 배정."""
cap = min(cap, int(budget // (C_PLAN + C_KIT))) # 월 예산이 허용하는 대수로 능력을 제한
month = np.full(len(pred), -1)
left = np.full(HORIZON, cap)
for i in np.argsort(pred):
if pred[i] > HORIZON * MDAYS + 60: break
tgt = int(np.clip((pred[i] - margin) // MDAYS - 1, 0, HORIZON - 1))
for m in list(range(tgt, -1, -1)) + list(range(tgt + 1, HORIZON)):
if left[m] > 0:
month[i] = m; left[m] -= 1; break
return month
def plan_fixed(lin_rul, cap=CAP_M, budget=BUDGET_M):
"""관행: 규정 주기 도달 순서대로(선형 잔여 일수) 월별 능력 안에서 배정."""
return plan_greedy(np.asarray(lin_rul, float), margin=0, cap=cap, budget=budget)
def evaluate_plan(month, true_rul, rate, cost_kit=True):
"""참 RUL로 계획의 실제 비용을 계산한다: 정비일 이전에 고장 -> 비계획 고장 비용, 아니면 계획 비용 + 미사용 수명 손실."""
day = np.where(month >= 0, (month + 1) * MDAYS, 10 ** 9)
fail = (true_rul < day) & (true_rul < HORIZON * MDAYS)
planned = (month >= 0) & ~fail
waste = np.where(planned, (true_rul - day) * rate * 0.6, 0) * C_WASTE_KH / 1000
cost = fail * C_FAIL + planned * (C_PLAN + (C_KIT if cost_kit else 0)) + waste
return {"total": float(cost.sum()), "fail": int(fail.sum()), "planned": int(planned.sum()), "waste": float(waste.sum()), "unsched_fail": int(((month < 0) & fail).sum())}
def kit_fefo(month, lots=LOTS):
"""유효기한이 빠른 로트부터 쓰는 FEFO로 키트를 배정하고 정비일에 유효기간이 부족한 배정 수와 만료 폐기 수량을 센다."""
q = [[qty, e] for qty, e in sorted(lots, key=lambda t: t[1])]
short = 0
for m in sorted(month[month >= 0]):
day = (m + 1) * MDAYS
for lot in q:
if lot[1] >= day + KIT_MARGIN and lot[0] > 0:
lot[0] -= 1; break
else:
short += 1
expired = sum(l[0] for l in q if l[1] < HORIZON * MDAYS)
return {"short": short, "expired_unused": expired, "left": sum(l[0] for l in q)}
# ---------------------------------------------------------------------------
tr = make_fleet(900, seed=1)
eng = fit_engine(tr)
fl = add_features(make_fleet(900, seed=3))
al = fl[fl.event == 0].reset_index(drop=True)
al["rul_pred"], _, _ = predict_rul(eng, al)
al["month"] = plan_milp(al["rul_pred"].values, al["rate"].values)
al["p90"] = aft_fail_prob(eng["p"], al, al["h_since"].values, al["h_since"].values + 90 * al["rate"].values)
# ---------------- 1) 부대별 노후도 현황 ----------------
g = al.groupby("unit")
board = pd.DataFrame({"장비 수": g.size(), "평균 노후 지수": g["aging_idx"].mean(), "평균 예측 RUL(일)": g["rul_pred"].mean(), "180일 내 고위험": g.apply(lambda x: int((x.rul_pred < 180).sum())),
"험지·염분 비율": g.apply(lambda x: x.env.isin(["rough", "salt"]).mean()), "12개월 내 권장 정비": g.apply(lambda x: int((x.month >= 0).sum()))})
board["고위험 비율"] = board["180일 내 고위험"] / board["장비 수"]
print(f"[1] 부대별 노후도 현황(생존 장비 {len(al)}대, 가상 부대 6개) [교육용가상자료]")
print(board.round(2).to_string())
worst = board["고위험 비율"].idxmax()
print(f" 고위험 비율이 가장 높은 부대 {worst}({board.loc[worst, '고위험 비율']:.0%}), 환경 구성(험지·염분 {board.loc[worst, '험지·염분 비율']:.0%})과 함께 본다.")
# ---------------- 2) 권장 정비 캘린더 ----------------
cal = pd.crosstab(al.loc[al.month >= 0, "unit"], al.loc[al.month >= 0, "month"] + 1).reindex(columns=range(1, HORIZON + 1), fill_value=0)
print("\n[2] 권장 정비 캘린더(부대 x 월, 대수; 월 예산·능력·키트 유효기간 반영)")
print(cal.to_string())
print(f" 월 합계: " + " ".join(f"{m}월 {int(c)}" for m, c in cal.sum().items()) + f" / 정비 일정이 잡히지 않은 고위험(RUL 180일 미만) 장비 {int(((al.month < 0) & (al.rul_pred < 180)).sum())}대")
top = al.sort_values("rul_pred").head(8)
print("\n 위험 상위 8대 [권장 월은 승인 전 추천안]")
print(pd.DataFrame({"장비": top.eq_id, "부대": top.unit, "환경": top.env.map(ENV_KO), "예측 RUL(일)": top.rul_pred, "90일 내 고장 확률": top.p90, "권장 월": top.month.map(lambda m: f"{m + 1}월" if m >= 0 else "-")}).round(2).to_string(index=False))
# ---------------- 3) 기여도 상위 3개(변수군 평균 대체) ----------------
GROUPS = {"운용 이력(경과 시간·거리·연수)": ["h_since", "hours", "km", "age_y"], "운용 환경": ["env"], "정비 이력(비계획·계획 정비)": ["unpl", "n_maint"], "점검 마모 비율": ["insp"], "일 가동률": ["rate"]}
MEAN = tr[["h_since", "hours", "km", "age_y", "unpl", "n_maint", "insp", "rate"]].mean()
raw_cols = ["eq_id", "unit", "env", "etype", "rate", "hours", "km", "age_y", "n_ovh", "h_since", "unpl", "n_maint", "insp"]
base_raw = al[raw_cols].copy()
full = al["rul_pred"].values
contrib = {}
for gname, cols in GROUPS.items():
x = base_raw.copy()
for c in cols:
x[c] = "normal" if c == "env" else MEAN[c]
xf = add_features(x.assign(event=0, t_fail=0.0, dur=x["h_since"]))
contrib[gname] = full - predict_rul(eng, xf)[0] # 음수 = 이 요인 때문에 평균 장비보다 잔여 수명이 짧다
C = pd.DataFrame(contrib)
t3 = C.apply(lambda r: r.reindex(r.abs().sort_values(ascending=False).index)[:3], axis=1)
okn = int(sum(t3.iloc[i].notna().sum() == 3 and np.isfinite(t3.iloc[i].dropna().values).all() for i in range(len(t3))))
print(f"\n[3] 장비 {len(al)}대 중 기여도 상위 3개 요인 자동 산출 {okn}대 ({okn / len(al):.0%}, 목표 100%)")
print(" 요인군별 평균 기여 크기(일):", ", ".join(f"{k} {v:.0f}" for k, v in C.abs().mean().sort_values(ascending=False).items()))
def explain(i):
r, row = al.loc[i], C.loc[i]
ordered = row.reindex(row.abs().sort_values(ascending=False).index)[:3]
main = ordered.index[0]
why = {"운용 이력(경과 시간·거리·연수)": "창정비 이후 운용 시간이 길어", "운용 환경": f"{ENV_KO[r.env]} 환경 가혹도 때문에", "정비 이력(비계획·계획 정비)": "비계획 정비가 잦아", "점검 마모 비율": "점검 마모 비율이 높아", "일 가동률": "일 가동률이 높아"}[main]
return f"{r.eq_id}({r.unit}, {ENV_KO[r.env]}): {why} 잔여 수명 {r.rul_pred:.0f}일 예측. 주요 요인 " + ", ".join(f"{k.split('(')[0]}({v:+.0f}일)" for k, v in ordered.items())
print("[4] 위험 상위 3대 자연어 설명(음수 = 평균 장비보다 잔여 수명이 짧아지는 요인)")
for i in top.index[:3]:
print(" -", explain(i))
print(" 대시보드와 설명은 정비계획관의 판단 보조이며 일정 확정과 장비 도태 건의는 군수지휘관 승인(M7)에 따른다.")
fig, ax = plt.subplots(2, 2, figsize=(13, 8))
ax[0, 0].bar(board.index, board["고위험 비율"] * 100); ax[0, 0].set_title("1 High-risk share by unit (%, RUL < 180 days)")
im = ax[0, 1].imshow(cal.values, aspect="auto", cmap="YlOrRd"); ax[0, 1].set_xticks(range(HORIZON)); ax[0, 1].set_xticklabels(range(1, HORIZON + 1)); ax[0, 1].set_yticks(range(len(cal))); ax[0, 1].set_yticklabels(cal.index)
ax[0, 1].set_title("2 Recommended overhaul calendar (unit x month)"); plt.colorbar(im, ax=ax[0, 1])
ax[1, 0].bar(board.index, board["평균 노후 지수"]); ax[1, 0].axhline(1, c="r", ls="--"); ax[1, 0].set_title("3 Mean aging index by unit")
r0 = C.loc[top.index[0]]
ax[1, 1].barh(["usage history", "environment", "maintenance history", "inspection wear", "daily usage"][::-1], r0.values[::-1]); ax[1, 1].set_title("4 Contribution to RUL, riskiest unit (days)")
plt.tight_layout(); plt.show()
In [ ]:
"""정비창 능력·예산·시효성 키트 제약 아래 예측 기반 계획정비 일정을 규정 주기 고정 일정과 비교해 비용·키트 소진율을 시나리오·백테스트로 검증하고(비용 15% 이상 절감), 정비 일정 확정과 장비 도태 건의를 군수지휘관(도태는 재정 담당 2단계)이 직접 승인·반려하는 HITL 결재 워크플로우(권한·사유·해시 체인 감사 로그)를 구현한다."""
import hashlib
import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 생존 분석 도구(scipy만 사용): 와이블 AFT(가속 수명) 모형과 Cox 비례위험 모형
# ---------------------------------------------------------------------------
from scipy.optimize import minimize
def aft_x(d):
return np.c_[np.ones(len(d)), d["env"].eq("rough"), d["env"].eq("salt"), d["env"].eq("cold"), d["unpl"] / (d["h_since"] / 1000 + 0.5)].astype(float)
AFT_COLS = ["절편", "험지", "염분", "혹한", "비계획 정비율"]
def fit_aft(d):
"""창정비 후 운용시간 dur(고장이면 event=1, 생존이면 중도절단)에 대한 와이블 AFT 최대우도 추정. 반환: 계수(마지막은 log 형상)."""
x, t, e = aft_x(d), d["dur"].values, d["event"].values
def nll(p):
lam, k = np.exp(x @ p[:-1]), np.exp(p[-1])
z = (t / lam) ** k
return -(e * (np.log(k) - np.log(lam) + (k - 1) * (np.log(t) - np.log(lam))) - z).sum()
r = minimize(nll, np.r_[np.log(8000), np.zeros(x.shape[1] - 1), np.log(5)], method="BFGS")
return r.x
def aft_pred(p, d, q=0.5, h=None):
"""이미 h(운용시간)까지 생존했다는 조건에서 분위수 q 생존 시점(운용시간 h 기준 총 수명)을 돌려준다. q=0.5가 조건부 중앙값."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
h = d["h_since"].values if h is None else h
return lam * ((h / lam) ** k - np.log(q)) ** (1 / k)
def aft_fail_prob(p, d, h0, h1):
"""h0까지 생존한 장비가 h0~h1 사이에 고장 날 조건부 확률."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
s0, s1 = np.exp(-(h0 / lam) ** k), np.exp(-(h1 / lam) ** k)
return 1 - s1 / s0
def fit_cox(d, iters=25):
"""Cox 비례위험(Breslow 근사) 뉴턴-랩슨. 공변량: 험지·염분·혹한·비계획 정비율. 반환: 계수."""
x = aft_x(d)[:, 1:]
o = np.argsort(-d["dur"].values)
x, e = x[o], d["event"].values[o]
b = np.zeros(x.shape[1])
for _ in range(iters):
r = np.exp(x @ b)
S0 = np.cumsum(r)
S1 = np.cumsum(r[:, None] * x, axis=0)
S2 = np.cumsum(r[:, None, None] * x[:, :, None] * x[:, None, :], axis=0)
m1 = S1 / S0[:, None]
g = ((x - m1) * e[:, None]).sum(axis=0)
H = -sum(e[i] * (S2[i] / S0[i] - np.outer(m1[i], m1[i])) for i in np.where(e == 1)[0])
step = np.linalg.solve(H - 1e-6 * np.eye(len(b)), g)
b = b - step
if np.abs(step).max() < 1e-7:
break
return b
def cindex(time, event, risk):
"""일치도 C-index: 위험 점수가 높은 쪽이 먼저 고장 났는가(관측 가능한 쌍만)."""
t, e, r = np.asarray(time), np.asarray(event), np.asarray(risk)
num = den = 0.0
for i in np.where(e == 1)[0]:
m = t > t[i]
den += m.sum()
num += (r[m] < r[i]).sum() + 0.5 * (r[m] == r[i]).sum()
return num / max(den, 1)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 잔여 수명(RUL, 일) 예측 엔진: 와이블 AFT 조건부 중앙값 + 그래디언트 부스팅(고장 이력 시점 표본 학습)의 평균
# ---------------------------------------------------------------------------
from sklearn.ensemble import GradientBoostingRegressor
def fit_engine(train, seed=0):
p = fit_aft(train)
sn = add_features(snapshots(train, 5, seed))
sn["aft"] = (aft_pred(p, sn) - sn["h_since"]) / sn["rate"]
g = GradientBoostingRegressor(n_estimators=250, max_depth=3, learning_rate=0.06, random_state=0).fit(sn[FEATS + ["aft"]], sn["rul_d"])
return {"p": p, "g": g}
def predict_rul(eng, d):
"""d는 add_features를 거친 표. 반환: (앙상블 RUL, AFT RUL, GBR RUL) - 일 단위."""
aft = ((aft_pred(eng["p"], d) - d["h_since"]) / d["rate"]).values
x = d[FEATS].copy()
x["aft"] = aft
gb = eng["g"].predict(x)
return np.clip(0.5 * aft + 0.5 * gb, 1, None), aft, gb
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 계획정비(창정비) 일정 추천: 월별 정비창 능력·예산·시효성 부품 키트(유효기한) 제약 아래 기대 비용 최소화 (비용 단가는 교육용 가정값, 만원)
# ---------------------------------------------------------------------------
from scipy.optimize import milp, LinearConstraint, Bounds
from scipy.stats import norm
C_PLAN, C_FAIL, C_WASTE_KH, C_KIT = 600.0, 2400.0, 40.0, 120.0 # 계획 창정비, 비계획 고장(긴급 정비 포함), 미사용 수명 1000h당 손실, 시효성 부품 키트
HORIZON, MDAYS = 12, 30 # 계획 12개월, 월 30일
CAP_M, BUDGET_M = 22, 15000.0 # 월별 정비창 능력(대), 월별 예산
LOTS = [(30, 70), (45, 170), (45, 390), (70, 800)] # 시효성 키트 로트(수량, 유효기한까지 남은 일수)
KIT_MARGIN = 14 # 정비일 기준 최소 유효기간 여유(일)
def exp_cost_matrix(pred, rate, month_days):
"""장비 i를 월 m에 정비할 때의 기대 비용. 참 RUL ~ N(pred, sigma), sigma = 35일 + 0.12 x pred (가정). 마지막 열은 계획 기간 밖(미정비)."""
sig = 35 + 0.12 * pred
C = np.zeros((len(pred), len(month_days) + 1))
for j, day in enumerate(month_days):
z = (day - pred) / sig
pf = norm.cdf(z)
unused = np.maximum(pred - day, 0) * rate * 0.6 # 정비를 앞당긴 만큼 버리는 운용시간 h(근사)
C[:, j] = C_PLAN + C_KIT + pf * (C_FAIL - C_PLAN) + C_WASTE_KH * unused / 1000 * (1 - pf)
C[:, -1] = norm.cdf((HORIZON * MDAYS - pred) / sig) * C_FAIL
return C
def plan_milp(pred, rate, cap=CAP_M, budget=BUDGET_M, lots=LOTS, use_kit=True):
"""정수계획으로 장비별 정비 월(0~11) 또는 미정비(-1)를 배정한다. 후보는 예측 RUL이 계획 기간 +90일 이내인 장비."""
md = np.arange(1, HORIZON + 1) * MDAYS
cand = np.where(pred < HORIZON * MDAYS + 90)[0]
n, M = len(cand), HORIZON
C = exp_cost_matrix(pred[cand], rate[cand], md)
c = C.ravel()
A, lo, hi = [], [], []
for i in range(n): # 장비당 한 번만 배정(월 12칸 + 미정비)
row = np.zeros(n * (M + 1)); row[i * (M + 1):(i + 1) * (M + 1)] = 1; A.append(row); lo.append(1); hi.append(1)
for m in range(M): # 월별 능력·예산
row = np.zeros(n * (M + 1)); row[m::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(cap)
A.append(row * (C_PLAN + C_KIT)); lo.append(0); hi.append(budget)
if use_kit: # 시효성 키트: m월 이후 정비 수요는 그 시점에 유효기간이 남은 로트 수량을 넘을 수 없다
for m in range(M):
sup = sum(q for q, e in lots if e >= md[m] + KIT_MARGIN)
row = np.zeros(n * (M + 1))
for mm in range(m, M):
row[mm::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(sup)
r = milp(c, constraints=LinearConstraint(np.array(A), lo, hi), integrality=np.ones(len(c)), bounds=Bounds(0, 1), options={"time_limit": 20})
month = np.full(len(pred), -1)
if r.x is not None:
x = np.round(r.x).reshape(n, M + 1)
assign = np.where(x[:, :M].sum(axis=1) > 0.5, x[:, :M].argmax(axis=1), -1)
month[cand] = assign
return month
def plan_greedy(pred, margin=45, cap=CAP_M, budget=BUDGET_M):
"""간단한 탐욕 배정: 정비 목표일(예측 RUL - margin)이 빠른 순으로 목표일 이전 가장 늦은 빈 월, 없으면 가장 이른 빈 월에 배정."""
cap = min(cap, int(budget // (C_PLAN + C_KIT))) # 월 예산이 허용하는 대수로 능력을 제한
month = np.full(len(pred), -1)
left = np.full(HORIZON, cap)
for i in np.argsort(pred):
if pred[i] > HORIZON * MDAYS + 60: break
tgt = int(np.clip((pred[i] - margin) // MDAYS - 1, 0, HORIZON - 1))
for m in list(range(tgt, -1, -1)) + list(range(tgt + 1, HORIZON)):
if left[m] > 0:
month[i] = m; left[m] -= 1; break
return month
def plan_fixed(lin_rul, cap=CAP_M, budget=BUDGET_M):
"""관행: 규정 주기 도달 순서대로(선형 잔여 일수) 월별 능력 안에서 배정."""
return plan_greedy(np.asarray(lin_rul, float), margin=0, cap=cap, budget=budget)
def evaluate_plan(month, true_rul, rate, cost_kit=True):
"""참 RUL로 계획의 실제 비용을 계산한다: 정비일 이전에 고장 -> 비계획 고장 비용, 아니면 계획 비용 + 미사용 수명 손실."""
day = np.where(month >= 0, (month + 1) * MDAYS, 10 ** 9)
fail = (true_rul < day) & (true_rul < HORIZON * MDAYS)
planned = (month >= 0) & ~fail
waste = np.where(planned, (true_rul - day) * rate * 0.6, 0) * C_WASTE_KH / 1000
cost = fail * C_FAIL + planned * (C_PLAN + (C_KIT if cost_kit else 0)) + waste
return {"total": float(cost.sum()), "fail": int(fail.sum()), "planned": int(planned.sum()), "waste": float(waste.sum()), "unsched_fail": int(((month < 0) & fail).sum())}
def kit_fefo(month, lots=LOTS):
"""유효기한이 빠른 로트부터 쓰는 FEFO로 키트를 배정하고 정비일에 유효기간이 부족한 배정 수와 만료 폐기 수량을 센다."""
q = [[qty, e] for qty, e in sorted(lots, key=lambda t: t[1])]
short = 0
for m in sorted(month[month >= 0]):
day = (m + 1) * MDAYS
for lot in q:
if lot[1] >= day + KIT_MARGIN and lot[0] > 0:
lot[0] -= 1; break
else:
short += 1
expired = sum(l[0] for l in q if l[1] < HORIZON * MDAYS)
return {"short": short, "expired_unused": expired, "left": sum(l[0] for l in q)}
# ---------------------------------------------------------------------------
tr = make_fleet(900, seed=1)
eng = fit_engine(tr)
def scenario(seed, shift=1.0, cap=CAP_M, budget=BUDGET_M, lots=LOTS, use_kit=True):
fl = add_features(make_fleet(900, seed=seed, shift=shift))
al = fl[fl.event == 0].reset_index(drop=True)
pred, _, _ = predict_rul(eng, al)
mf = plan_fixed(al["lin_rul_d"].values, cap=cap, budget=budget)
mm = plan_milp(pred, al["rate"].values, cap=cap, budget=budget, lots=lots, use_kit=use_kit)
rf, rm = evaluate_plan(mf, al["rul_d"].values, al["rate"].values), evaluate_plan(mm, al["rul_d"].values, al["rate"].values)
kf, km = kit_fefo(mf, lots), kit_fefo(mm, lots)
stock = sum(q for q, _ in lots)
return al, pred, mf, mm, rf, rm, kf, km, stock
# ---------------- 1) 시나리오 6종 ----------------
SC = [("평상", dict(seed=3)), ("정비창 능력 70%(15대/월)", dict(seed=3, cap=15)), ("예산 삭감(월 10,000만원)", dict(seed=3, budget=10000.0)), ("시효성 키트 재고 절반·단기 유효", dict(seed=3, lots=[(15, 70), (20, 170), (20, 240), (30, 330)])),
("가동률 +30%", dict(seed=3, shift=1.3)), ("다른 보유 장비(시드 4)", dict(seed=4))]
rows, cache = [], {}
for name, kw in SC:
al, pred, mf, mm, rf, rm, kf, km, stock = scenario(**kw)
cache[name] = (al, pred, mm)
rows.append((name, len(al), rf["total"], rm["total"], 1 - rm["total"] / rf["total"], rf["fail"], rm["fail"], kf["expired_unused"], km["expired_unused"], km["short"]))
res = pd.DataFrame(rows, columns=["시나리오", "생존 장비", "관행 비용(만원)", "제안 비용(만원)", "절감률", "관행 비계획 고장", "제안 비계획 고장", "관행 키트 만료 폐기", "제안 키트 만료 폐기", "제안 키트 유효기간 부족"])
print("[1] 시나리오별 12개월 비용 비교(관행 = 규정 주기 도달 순 고정 일정, 단가는 교육용 가정값)")
print(res.assign(절감률=res["절감률"].map("{:.1%}".format)).round(0).to_string(index=False))
avg = res["절감률"].mean()
print(f" 시나리오 6종 평균 절감 {avg:.1%} (목표 15% 이상: {'달성' if avg >= .15 else '미달'}), 최소 {res['절감률'].min():.1%}({res.loc[res['절감률'].idxmin(), '시나리오']})")
# ---------------- 2) 백테스트: 서로 다른 보유 장비 5세트 ----------------
bt = []
for sd in (10, 11, 12, 13, 14):
_, _, _, _, rf, rm, _, _, _ = scenario(seed=sd)
bt.append((sd, 1 - rm["total"] / rf["total"], rf["fail"], rm["fail"]))
bt = pd.DataFrame(bt, columns=["장비 세트", "절감률", "관행 고장", "제안 고장"])
print(f"\n[2] 백테스트(서로 다른 보유 장비 5세트): 평균 절감 {bt['절감률'].mean():.1%}, 최소 {bt['절감률'].min():.1%}, 최대 {bt['절감률'].max():.1%}, 비계획 고장 평균 {bt['관행 고장'].mean():.0f} -> {bt['제안 고장'].mean():.0f}대")
# ---------------- 3) 단가 민감도: 비계획 고장 비용 ----------------
sens = []
base_cfail = C_FAIL
for mult in (0.5, 1.0, 1.5, 2.0):
C_FAIL = base_cfail * mult
_, _, _, _, rf, rm, _, _, _ = scenario(seed=3)
sens.append((mult, 1 - rm["total"] / rf["total"]))
C_FAIL = base_cfail
print("[3] 비계획 고장 비용 민감도:", "; ".join(f"x{m} -> 절감 {s:.0%}" for m, s in sens))
print(" 고장 비용이 낮아지면 앞당긴 정비의 낭비가 상대적으로 커져 절감 폭이 줄어든다. 단가는 재정 부서와 실제 값으로 다시 확인해야 한다.")
# ---------------- 4) HITL 결재 워크플로우 ----------------
USERS = {"ai_bot": "AI", "pln_kim": "PLN", "pln_lee": "PLN", "log_park": "LOG", "fin_choi": "FIN", "aud_han": "AUD"}
PERM = {"AI": {"submit"}, "PLN": {"submit", "redraft"}, "LOG": {"approve", "reject"}, "FIN": {"approve", "reject"}, "AUD": {"read"}}
class Flow:
def __init__(self):
self.docs, self.log, self.blocked = {}, [], []
def _audit(self, user, action, doc, detail):
prev = self.log[-1]["hash"] if self.log else "GENESIS"
rec = {"seq": len(self.log) + 1, "user": user, "action": action, "doc": doc, "detail": detail, "prev": prev}
rec["hash"] = hashlib.sha256(json.dumps(rec, ensure_ascii=False, sort_keys=True).encode()).hexdigest()[:12]
self.log.append(rec)
def verify(self):
prev = "GENESIS"
for r in self.log:
body = {k: v for k, v in r.items() if k != "hash"}
if r["prev"] != prev or hashlib.sha256(json.dumps(body, ensure_ascii=False, sort_keys=True).encode()).hexdigest()[:12] != r["hash"]:
return False, r["seq"]
prev = r["hash"]
return True, None
def deny(self, user, act, doc, why):
self.blocked.append((user, act, doc, why)); self._audit(user, "DENIED:" + act, doc, why); return False
def submit(self, user, kind, title, n, parent=None):
need = "redraft" if parent else "submit"
if need not in PERM.get(USERS.get(user), set()):
return self.deny(user, need, "-", "상신 권한 없음")
did = f"MNT{len(self.docs) + 1:02d}"
self.docs[did] = {"kind": kind, "title": title, "line": ["LOG", "FIN"] if kind == "도태 건의" else ["LOG"], "step": 0, "status": "검토대기", "drafter": user, "done": [], "parent": parent}
self._audit(user, "submit", did, f"{kind} {title} {n}대")
return did
def decide(self, user, did, act, reason=""):
role, d = USERS.get(user), self.docs.get(did)
if role is None: return self.deny(user, act, did, "미등록 사용자")
if d is None or d["status"] != "검토대기": return self.deny(user, act, did, "결재 대상 아님(처리 완료)")
if act not in PERM[role]: return self.deny(user, act, did, f"{role} 역할은 {act} 불가")
if d["line"][d["step"]] != role: return self.deny(user, act, did, "결재 순서 위반")
if user == d["drafter"] or user in d["done"]: return self.deny(user, act, did, "기안자·기결재자 재결재 금지")
if len(reason.strip()) < 10: return self.deny(user, act, did, "사유 10자 이상 필수")
if act == "reject":
d["status"] = "반려"
else:
d["done"].append(user); d["step"] += 1
d["status"] = "승인" if d["step"] == len(d["line"]) else "검토대기"
self._audit(user, act, did, reason); return True
def edit(self, user, did):
if self.docs[did]["status"] == "승인": return self.deny(user, "edit", did, "승인 후 일정 변경 금지(재기안)")
return True
al, pred, mm = cache["평상"]
al = al.assign(rul_pred=pred, month=mm)
scrap = al[(al["age_y"] >= 14) & (al["rul_pred"] < 365)] # 도태 건의 후보: 사용 14년 이상이고 1년 내 수명 종료 예측
wf = Flow()
docs = [wf.submit("ai_bot" if q % 2 == 0 else "pln_kim", "정비 일정", f"{3 * q + 1}~{3 * q + 3}월 일정안", int(((al.month >= 3 * q) & (al.month < 3 * q + 3)).sum())) for q in range(2)]
docs += [wf.submit("pln_lee", "도태 건의", f"{u} 부대 노후 장비", int((scrap.unit == u).sum())) for u in scrap["unit"].value_counts().index[:2]]
two = [i for i in docs if len(wf.docs[i]["line"]) == 2]
print(f"\n[4] AI 추천안 {len(docs)}건(정비 일정 2건, 도태 건의 {len(two)}건; 도태 후보 {len(scrap)}대) 전부 '검토대기' 등록, 자동 승인 {sum(wf.docs[i]['status'] == '승인' for i in docs)}건 (목표 0건). 2단계 결재(군수지휘관 + 재정 담당) 대상 {two}")
wf.decide("log_park", docs[0], "reject", "1~3월 정비 순서에 키트 유효기간 재확인 필요, 재기안 요청")
new = wf.submit("pln_kim", "정비 일정", "1~3월 일정안(재기안)", int(((al.month >= 0) & (al.month < 3)).sum()), parent=docs[0])
wf.decide("log_park", new, "approve", "키트 유효기간 확인 완료, 일정 확정")
wf.decide("log_park", docs[1], "approve", "4~6월 일정, 정비창 능력 범위 내 확인함")
if two:
wf.decide("log_park", two[0], "approve", "도태 건의 사유 타당, 재정 담당 검토 상신")
wf.decide("fin_choi", two[0], "approve", "잔존 가치 대비 수리비 초과 확인, 도태 건의 동의")
tgt = two[-1] if two else docs[-1]
attacks = [("pln_kim", tgt, "approve", "상신자가 직접 승인한다고 처리함"), ("ai_bot", tgt, "approve", "AI 계정이 직접 승인하려고 시도함"), ("aud_han", tgt, "approve", "감사 계정으로 승인 시도함"), ("ghost", tgt, "approve", "외부 계정으로 승인 시도함"),
("log_park", tgt, "approve", ""), ("log_park", tgt, "reject", "짧음"), ("fin_choi", tgt, "approve", "재정 담당이 순서를 건너뛰어 승인"), ("log_park", new, "reject", "승인된 일정을 뒤집어 반려함"),
("pln_lee", tgt, "approve", "상신자 계정으로 승인 시도함")]
for u, doc, a, rs in attacks:
wf.decide(u, doc, a, rs)
wf.edit("pln_kim", new)
print(f" 권한 우회 시도 10건 중 차단 {len(wf.blocked)}건 (차단율 {len(wf.blocked) / 10:.0%}, 목표 100%)")
for u, a, doc, why in wf.blocked[:4]:
print(f" 차단 {u:<9} {a:<8} {doc}: {why}")
ok, _ = wf.verify()
snap = json.loads(json.dumps(wf.log))
wf.log[3]["detail"] = "사후에 고쳐 쓴 사유"
ok2, at = wf.verify()
wf.log = snap
print(f" 감사 로그 {len(wf.log)}건 해시 체인 무결성 {'100%' if ok else '훼손'}; 4번 기록 사후 수정 시 {'통과' if ok2 else f'위변조 검출({at}번)'}")
print(f" {docs[0]}(반려) -> {new}(재기안, parent 연결) 이력:", [(r['action'], r['user']) for r in wf.log if r['doc'] in (docs[0], new)])
print(" 문서 상태:", pd.Series({k: v["status"] for k, v in wf.docs.items()}).value_counts().to_dict(), "| 정비 일정 확정과 장비 도태는 군수지휘관 결재(도태는 재정 담당 포함)로만 확정된다.")
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
x = np.arange(len(res)); ax[0].bar(x - .2, res["관행 비용(만원)"], .4, label="fixed-cycle"); ax[0].bar(x + .2, res["제안 비용(만원)"], .4, label="risk-based plan")
ax[0].set_xticks(x); ax[0].set_xticklabels([f"S{i + 1}" for i in x]); ax[0].legend(); ax[0].set_title("12-month cost by scenario (10k KRW)")
ax[1].plot([s[0] for s in sens], [s[1] * 100 for s in sens], "o-"); ax[1].set_title("Saving (%) vs unplanned-failure cost multiplier"); ax[1].set_xlabel("multiplier")
plt.tight_layout(); plt.show()
In [ ]:
"""보유 장비 1천~10만 대의 수명(RUL) 일괄 예측·위험 경보·계획정비 일정 배정·대시보드 집계 시뮬레이션으로 단계별 처리 시간과 대시보드 로딩 시간을 재고, 결측·환경 오기재·단위 오류 주입 조건에서 위험 경보 정확도를 측정하며 정수계획의 규모 한계와 대응(탐욕 배정)을 정리한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 생존 분석 도구(scipy만 사용): 와이블 AFT(가속 수명) 모형과 Cox 비례위험 모형
# ---------------------------------------------------------------------------
from scipy.optimize import minimize
def aft_x(d):
return np.c_[np.ones(len(d)), d["env"].eq("rough"), d["env"].eq("salt"), d["env"].eq("cold"), d["unpl"] / (d["h_since"] / 1000 + 0.5)].astype(float)
AFT_COLS = ["절편", "험지", "염분", "혹한", "비계획 정비율"]
def fit_aft(d):
"""창정비 후 운용시간 dur(고장이면 event=1, 생존이면 중도절단)에 대한 와이블 AFT 최대우도 추정. 반환: 계수(마지막은 log 형상)."""
x, t, e = aft_x(d), d["dur"].values, d["event"].values
def nll(p):
lam, k = np.exp(x @ p[:-1]), np.exp(p[-1])
z = (t / lam) ** k
return -(e * (np.log(k) - np.log(lam) + (k - 1) * (np.log(t) - np.log(lam))) - z).sum()
r = minimize(nll, np.r_[np.log(8000), np.zeros(x.shape[1] - 1), np.log(5)], method="BFGS")
return r.x
def aft_pred(p, d, q=0.5, h=None):
"""이미 h(운용시간)까지 생존했다는 조건에서 분위수 q 생존 시점(운용시간 h 기준 총 수명)을 돌려준다. q=0.5가 조건부 중앙값."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
h = d["h_since"].values if h is None else h
return lam * ((h / lam) ** k - np.log(q)) ** (1 / k)
def aft_fail_prob(p, d, h0, h1):
"""h0까지 생존한 장비가 h0~h1 사이에 고장 날 조건부 확률."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
s0, s1 = np.exp(-(h0 / lam) ** k), np.exp(-(h1 / lam) ** k)
return 1 - s1 / s0
def fit_cox(d, iters=25):
"""Cox 비례위험(Breslow 근사) 뉴턴-랩슨. 공변량: 험지·염분·혹한·비계획 정비율. 반환: 계수."""
x = aft_x(d)[:, 1:]
o = np.argsort(-d["dur"].values)
x, e = x[o], d["event"].values[o]
b = np.zeros(x.shape[1])
for _ in range(iters):
r = np.exp(x @ b)
S0 = np.cumsum(r)
S1 = np.cumsum(r[:, None] * x, axis=0)
S2 = np.cumsum(r[:, None, None] * x[:, :, None] * x[:, None, :], axis=0)
m1 = S1 / S0[:, None]
g = ((x - m1) * e[:, None]).sum(axis=0)
H = -sum(e[i] * (S2[i] / S0[i] - np.outer(m1[i], m1[i])) for i in np.where(e == 1)[0])
step = np.linalg.solve(H - 1e-6 * np.eye(len(b)), g)
b = b - step
if np.abs(step).max() < 1e-7:
break
return b
def cindex(time, event, risk):
"""일치도 C-index: 위험 점수가 높은 쪽이 먼저 고장 났는가(관측 가능한 쌍만)."""
t, e, r = np.asarray(time), np.asarray(event), np.asarray(risk)
num = den = 0.0
for i in np.where(e == 1)[0]:
m = t > t[i]
den += m.sum()
num += (r[m] < r[i]).sum() + 0.5 * (r[m] == r[i]).sum()
return num / max(den, 1)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 잔여 수명(RUL, 일) 예측 엔진: 와이블 AFT 조건부 중앙값 + 그래디언트 부스팅(고장 이력 시점 표본 학습)의 평균
# ---------------------------------------------------------------------------
from sklearn.ensemble import GradientBoostingRegressor
def fit_engine(train, seed=0):
p = fit_aft(train)
sn = add_features(snapshots(train, 5, seed))
sn["aft"] = (aft_pred(p, sn) - sn["h_since"]) / sn["rate"]
g = GradientBoostingRegressor(n_estimators=250, max_depth=3, learning_rate=0.06, random_state=0).fit(sn[FEATS + ["aft"]], sn["rul_d"])
return {"p": p, "g": g}
def predict_rul(eng, d):
"""d는 add_features를 거친 표. 반환: (앙상블 RUL, AFT RUL, GBR RUL) - 일 단위."""
aft = ((aft_pred(eng["p"], d) - d["h_since"]) / d["rate"]).values
x = d[FEATS].copy()
x["aft"] = aft
gb = eng["g"].predict(x)
return np.clip(0.5 * aft + 0.5 * gb, 1, None), aft, gb
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 계획정비(창정비) 일정 추천: 월별 정비창 능력·예산·시효성 부품 키트(유효기한) 제약 아래 기대 비용 최소화 (비용 단가는 교육용 가정값, 만원)
# ---------------------------------------------------------------------------
from scipy.optimize import milp, LinearConstraint, Bounds
from scipy.stats import norm
C_PLAN, C_FAIL, C_WASTE_KH, C_KIT = 600.0, 2400.0, 40.0, 120.0 # 계획 창정비, 비계획 고장(긴급 정비 포함), 미사용 수명 1000h당 손실, 시효성 부품 키트
HORIZON, MDAYS = 12, 30 # 계획 12개월, 월 30일
CAP_M, BUDGET_M = 22, 15000.0 # 월별 정비창 능력(대), 월별 예산
LOTS = [(30, 70), (45, 170), (45, 390), (70, 800)] # 시효성 키트 로트(수량, 유효기한까지 남은 일수)
KIT_MARGIN = 14 # 정비일 기준 최소 유효기간 여유(일)
def exp_cost_matrix(pred, rate, month_days):
"""장비 i를 월 m에 정비할 때의 기대 비용. 참 RUL ~ N(pred, sigma), sigma = 35일 + 0.12 x pred (가정). 마지막 열은 계획 기간 밖(미정비)."""
sig = 35 + 0.12 * pred
C = np.zeros((len(pred), len(month_days) + 1))
for j, day in enumerate(month_days):
z = (day - pred) / sig
pf = norm.cdf(z)
unused = np.maximum(pred - day, 0) * rate * 0.6 # 정비를 앞당긴 만큼 버리는 운용시간 h(근사)
C[:, j] = C_PLAN + C_KIT + pf * (C_FAIL - C_PLAN) + C_WASTE_KH * unused / 1000 * (1 - pf)
C[:, -1] = norm.cdf((HORIZON * MDAYS - pred) / sig) * C_FAIL
return C
def plan_milp(pred, rate, cap=CAP_M, budget=BUDGET_M, lots=LOTS, use_kit=True):
"""정수계획으로 장비별 정비 월(0~11) 또는 미정비(-1)를 배정한다. 후보는 예측 RUL이 계획 기간 +90일 이내인 장비."""
md = np.arange(1, HORIZON + 1) * MDAYS
cand = np.where(pred < HORIZON * MDAYS + 90)[0]
n, M = len(cand), HORIZON
C = exp_cost_matrix(pred[cand], rate[cand], md)
c = C.ravel()
A, lo, hi = [], [], []
for i in range(n): # 장비당 한 번만 배정(월 12칸 + 미정비)
row = np.zeros(n * (M + 1)); row[i * (M + 1):(i + 1) * (M + 1)] = 1; A.append(row); lo.append(1); hi.append(1)
for m in range(M): # 월별 능력·예산
row = np.zeros(n * (M + 1)); row[m::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(cap)
A.append(row * (C_PLAN + C_KIT)); lo.append(0); hi.append(budget)
if use_kit: # 시효성 키트: m월 이후 정비 수요는 그 시점에 유효기간이 남은 로트 수량을 넘을 수 없다
for m in range(M):
sup = sum(q for q, e in lots if e >= md[m] + KIT_MARGIN)
row = np.zeros(n * (M + 1))
for mm in range(m, M):
row[mm::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(sup)
r = milp(c, constraints=LinearConstraint(np.array(A), lo, hi), integrality=np.ones(len(c)), bounds=Bounds(0, 1), options={"time_limit": 20})
month = np.full(len(pred), -1)
if r.x is not None:
x = np.round(r.x).reshape(n, M + 1)
assign = np.where(x[:, :M].sum(axis=1) > 0.5, x[:, :M].argmax(axis=1), -1)
month[cand] = assign
return month
def plan_greedy(pred, margin=45, cap=CAP_M, budget=BUDGET_M):
"""간단한 탐욕 배정: 정비 목표일(예측 RUL - margin)이 빠른 순으로 목표일 이전 가장 늦은 빈 월, 없으면 가장 이른 빈 월에 배정."""
cap = min(cap, int(budget // (C_PLAN + C_KIT))) # 월 예산이 허용하는 대수로 능력을 제한
month = np.full(len(pred), -1)
left = np.full(HORIZON, cap)
for i in np.argsort(pred):
if pred[i] > HORIZON * MDAYS + 60: break
tgt = int(np.clip((pred[i] - margin) // MDAYS - 1, 0, HORIZON - 1))
for m in list(range(tgt, -1, -1)) + list(range(tgt + 1, HORIZON)):
if left[m] > 0:
month[i] = m; left[m] -= 1; break
return month
def plan_fixed(lin_rul, cap=CAP_M, budget=BUDGET_M):
"""관행: 규정 주기 도달 순서대로(선형 잔여 일수) 월별 능력 안에서 배정."""
return plan_greedy(np.asarray(lin_rul, float), margin=0, cap=cap, budget=budget)
def evaluate_plan(month, true_rul, rate, cost_kit=True):
"""참 RUL로 계획의 실제 비용을 계산한다: 정비일 이전에 고장 -> 비계획 고장 비용, 아니면 계획 비용 + 미사용 수명 손실."""
day = np.where(month >= 0, (month + 1) * MDAYS, 10 ** 9)
fail = (true_rul < day) & (true_rul < HORIZON * MDAYS)
planned = (month >= 0) & ~fail
waste = np.where(planned, (true_rul - day) * rate * 0.6, 0) * C_WASTE_KH / 1000
cost = fail * C_FAIL + planned * (C_PLAN + (C_KIT if cost_kit else 0)) + waste
return {"total": float(cost.sum()), "fail": int(fail.sum()), "planned": int(planned.sum()), "waste": float(waste.sum()), "unsched_fail": int(((month < 0) & fail).sum())}
def kit_fefo(month, lots=LOTS):
"""유효기한이 빠른 로트부터 쓰는 FEFO로 키트를 배정하고 정비일에 유효기간이 부족한 배정 수와 만료 폐기 수량을 센다."""
q = [[qty, e] for qty, e in sorted(lots, key=lambda t: t[1])]
short = 0
for m in sorted(month[month >= 0]):
day = (m + 1) * MDAYS
for lot in q:
if lot[1] >= day + KIT_MARGIN and lot[0] > 0:
lot[0] -= 1; break
else:
short += 1
expired = sum(l[0] for l in q if l[1] < HORIZON * MDAYS)
return {"short": short, "expired_unused": expired, "left": sum(l[0] for l in q)}
# ---------------------------------------------------------------------------
tr = make_fleet(900, seed=1)
t0 = time.perf_counter()
eng = fit_engine(tr)
print(f"[0] 예측 엔진 학습 {time.perf_counter() - t0:.2f}초(학습 {len(tr)}대, 시점 표본 {int(tr.event.sum()) * 5}건)")
ALARM_D, MARGIN = 180, 60 # 위험 정의: 참 RUL 180일 미만 / 경보 기준: 예측 RUL이 180+60일 미만(예측 오차를 고려한 여유)
def pipeline(fl):
"""일괄 처리: 피처 -> RUL 예측 -> 경보 -> 일정 배정(탐욕) -> 대시보드 집계. 단계별 시간을 반환한다."""
t = {}
s = time.perf_counter(); d = add_features(fl); t["피처"] = time.perf_counter() - s
al = d[d.event == 0].reset_index(drop=True)
s = time.perf_counter(); pred, _, _ = predict_rul(eng, al); t["예측"] = time.perf_counter() - s
s = time.perf_counter(); month = plan_greedy(pred, cap=max(CAP_M, int(len(al) * 0.03)), budget=1e9); t["일정"] = time.perf_counter() - s
s = time.perf_counter()
al["pred"], al["month"] = pred, month
summ = al.groupby("unit").agg(장비=("pred", "size"), 평균RUL=("pred", "mean"), 위험=("pred", lambda x: int((x < ALARM_D + MARGIN).sum())))
cal = pd.crosstab(al.loc[al.month >= 0, "unit"], al.loc[al.month >= 0, "month"])
low = al.nsmallest(20, "pred")
fg, ax = plt.subplots(1, 2, figsize=(9, 3.5)); ax[0].bar(summ.index, summ["위험"]); ax[1].imshow(cal.values, aspect="auto"); fg.canvas.draw(); plt.close(fg)
t["대시보드"] = time.perf_counter() - s
return al, t
# ---------------- 1) 규모별 시뮬레이션 ----------------
rows = []
for N in (1000, 5000, 10000, 50000, 100000):
s = time.perf_counter(); fl = make_fleet(N, seed=100 + N % 7); tg = time.perf_counter() - s
al, t = pipeline(fl)
tot = t["피처"] + t["예측"] + t["일정"]
rows.append((N, len(al), tg, t["피처"], t["예측"], t["일정"], t["대시보드"], len(al) / tot))
ld = pd.DataFrame(rows, columns=["보유 장비", "생존 장비", "생성(초)", "피처(초)", "예측(초)", "일정(초)", "대시보드 로딩(초)", "처리량(대/초)"])
print("\n[1] 규모별 일괄 처리 시간(생성은 시뮬레이션용이며 실제 운영에는 없음)")
print(ld.round(3).to_string(index=False))
r10 = ld[ld["보유 장비"] == 10000].iloc[0]
print(f" 1만 대 일괄 수명 시뮬레이션: 피처+예측+일정 {r10['피처(초)'] + r10['예측(초)'] + r10['일정(초)']:.2f}초, 대시보드 로딩 {r10['대시보드 로딩(초)']:.3f}초(목표 2초 이내: {'달성' if r10['대시보드 로딩(초)'] <= 2 else '미달'})")
# ---------------- 2) 위험 경보 정확도(1만 대) ----------------
fl10 = make_fleet(10000, seed=5)
d10 = add_features(fl10)
a10 = d10[d10.event == 0].reset_index(drop=True)
pred10, _, _ = predict_rul(eng, a10)
def alarm_metrics(pred, truth_d, margin=MARGIN):
al_, tr_ = pred < ALARM_D + margin, truth_d < ALARM_D
tp, fp, fn, tn = int((al_ & tr_).sum()), int((al_ & ~tr_).sum()), int((~al_ & tr_).sum()), int((~al_ & ~tr_).sum())
return (tp + tn) / len(pred), tp / max(tp + fn, 1), tp / max(tp + fp, 1), fp, fn
acc, rec, prec, fp, fn = alarm_metrics(pred10, a10["rul_d"].values)
sw = pd.DataFrame([(m, *alarm_metrics(pred10, a10["rul_d"].values, m)) for m in (0, 60, 120, 180)], columns=["여유(일)", "정확도", "재현율", "정밀도", "오경보", "미탐"])
print(f"\n[2] 위험 경보 정확도: 생존 {len(a10)}대, 실제 위험(참 RUL {ALARM_D}일 미만) {int((a10.rul_d < ALARM_D).sum())}대({(a10.rul_d < ALARM_D).mean():.1%})")
print(" 경보 기준(예측 RUL < 180 + 여유) 비교\n" + sw.assign(정확도=sw["정확도"].map("{:.1%}".format), 재현율=sw["재현율"].map("{:.1%}".format), 정밀도=sw["정밀도"].map("{:.1%}".format)).to_string(index=False))
print(f" 채택 여유 {MARGIN}일: 정확도 {acc:.1%} (목표 95% 이상: {'달성' if acc >= .95 else '미달'}), 재현율 {rec:.1%}, 정밀도 {prec:.1%}, 오경보 {fp}, 미탐 {fn}")
print(" 위험 장비가 드물어 정확도만으로는 부풀려지므로 재현율·정밀도를 함께 본다. 여유를 늘리면 재현율은 오르지만 오경보가 늘어 정확도가 95% 아래로 내려간다. 놓친 장비는 정기 점검(점검 마모)으로 보완한다.")
# ---------------- 3) 오류 주입 ----------------
rng = np.random.default_rng(7)
W_FIX = {"normal": 1.0, "rough": 1.15, "salt": 1.1, "cold": 1.05}
def corrupt(kind):
x = a10.copy()
if kind == "점검 마모 20% 결측(대체)":
m = rng.random(len(x)) < 0.2
x.loc[m, "insp"] = (x.loc[m, "h_since"] / CYCLE * x.loc[m, "env"].map(W_FIX)).clip(0, 1.3) # 정제 규칙: 노후 지수로 대체
elif kind == "점검 마모 20% 결측(0으로 방치)":
m = rng.random(len(x)) < 0.2
x.loc[m, "insp"] = 0.0
elif kind == "환경 10% 오기재":
m = rng.random(len(x)) < 0.1
x.loc[m, "env"] = rng.choice(ENVS, int(m.sum()))
elif kind == "가동률 5% 단위 오류(x30)":
m = rng.random(len(x)) < 0.05
x.loc[m, "rate"] = x.loc[m, "rate"] * 30
elif kind == "복합(결측 대체+환경 오기재+단위 오류)":
m = rng.random(len(x)) < 0.2; x.loc[m, "insp"] = (x.loc[m, "h_since"] / CYCLE * x.loc[m, "env"].map(W_FIX)).clip(0, 1.3)
m = rng.random(len(x)) < 0.1; x.loc[m, "env"] = rng.choice(ENVS, int(m.sum()))
m = rng.random(len(x)) < 0.05; x.loc[m, "rate"] = x.loc[m, "rate"] * 30
return add_features(x[["eq_id", "unit", "env", "etype", "rate", "hours", "km", "age_y", "n_ovh", "h_since", "unpl", "n_maint", "insp", "t_fail", "event", "dur", "rul_h", "rul_d"]])
rows = [("정상", acc, rec, prec, fp, fn)]
for kind in ("점검 마모 20% 결측(대체)", "점검 마모 20% 결측(0으로 방치)", "환경 10% 오기재", "가동률 5% 단위 오류(x30)", "복합(결측 대체+환경 오기재+단위 오류)"):
x = corrupt(kind)
p, _, _ = predict_rul(eng, x)
rows.append((kind, *alarm_metrics(p, a10["rul_d"].values)))
er = pd.DataFrame(rows, columns=["조건", "정확도", "재현율", "정밀도", "오경보", "미탐"])
print("\n[3] 오류 주입 시 경보 성능(1만 대)\n" + er.assign(정확도=er["정확도"].map("{:.1%}".format), 재현율=er["재현율"].map("{:.1%}".format), 정밀도=er["정밀도"].map("{:.1%}".format)).to_string(index=False))
worst = er.iloc[1:].sort_values("정확도").iloc[0]
z = er[er["조건"].str.contains("0으로")].iloc[0]
print(f" 정확도가 가장 낮은 조건은 '{worst['조건']}'({worst['정확도']:.1%}, 오경보 {int(worst['오경보'])}건). 결측을 0으로 방치하면 재현율이 {z['재현율']:.1%}로 정상({er.iloc[0]['재현율']:.1%})보다 떨어지므로 정제 규칙(M2)의 대체값 처리가 필요하다.")
print(f" 95% 목표는 정상 조건에서만 달성하며 오류가 겹치면 {er.iloc[1:]['정확도'].min():.1%}까지 낮아진다 -> 입력 검증(M1·M2)이 경보 품질의 전제이다.")
# ---------------- 4) 정수계획 규모 한계 ----------------
mt = []
for n in (50, 150, 300, 600):
sub = a10.head(20000).copy()
pr_ = np.sort(pred10)[:n] * 0 + np.linspace(30, 440, n) # 예측 RUL이 계획 기간 안에 분포하는 후보 n대
s = time.perf_counter(); m = plan_milp(pr_, np.full(n, 3.0), cap=max(CAP_M, n // 8), budget=1e9); dt = time.perf_counter() - s
mt.append((n, dt, int((m >= 0).sum())))
print("\n[4] 정수계획 후보 장비 수별 풀이 시간:", ", ".join(f"{n}대 {s:.2f}초" for n, s, _ in mt))
print(" 후보가 수천 대를 넘으면 정수계획이 느려지므로 위험 상위 후보만 정수계획, 나머지는 탐욕 배정으로 나누는 구조를 쓴다. 경보·일정은 판단 보조이며 확정은 군수지휘관 승인(M7)에 따른다.")
print(" 모든 부대·장비·수치는 교육용 가상 자료이다.")
fig, ax = plt.subplots(1, 3, figsize=(15, 4.2))
ax[0].loglog(ld["보유 장비"], ld["피처(초)"] + ld["예측(초)"] + ld["일정(초)"], "o-", label="batch total"); ax[0].loglog(ld["보유 장비"], ld["대시보드 로딩(초)"], "s-", label="dashboard"); ax[0].axhline(2, c="r", ls="--"); ax[0].legend(); ax[0].set_title("Processing time vs fleet size (s)")
ax[1].barh(range(len(er)), er["재현율"] * 100); ax[1].set_yticks(range(len(er))); ax[1].set_yticklabels([f"C{i + 1}" for i in range(len(er))]); ax[1].set_title("Alarm recall by fault case (%)")
ax[2].plot([m[0] for m in mt], [m[1] for m in mt], "o-"); ax[2].set_title("MILP solve time vs candidate units (s)")
plt.tight_layout(); plt.show()
In [ ]:
"""운영 기간 수명 예측 오차를 AAR로 환경별·원인별 분석하고, 환경 가중치를 운영 고장 이력의 AFT 계수로 재보정하고 운영 데이터로 예측 엔진을 재학습해 보류 장비 집단에서 개선을 검증한 뒤 비용 절감액·가동 중단일수 단축률·투자 회수 기간 3대 효과 지표를 산식과 민감도로 산출한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 생존 분석 도구(scipy만 사용): 와이블 AFT(가속 수명) 모형과 Cox 비례위험 모형
# ---------------------------------------------------------------------------
from scipy.optimize import minimize
def aft_x(d):
return np.c_[np.ones(len(d)), d["env"].eq("rough"), d["env"].eq("salt"), d["env"].eq("cold"), d["unpl"] / (d["h_since"] / 1000 + 0.5)].astype(float)
AFT_COLS = ["절편", "험지", "염분", "혹한", "비계획 정비율"]
def fit_aft(d):
"""창정비 후 운용시간 dur(고장이면 event=1, 생존이면 중도절단)에 대한 와이블 AFT 최대우도 추정. 반환: 계수(마지막은 log 형상)."""
x, t, e = aft_x(d), d["dur"].values, d["event"].values
def nll(p):
lam, k = np.exp(x @ p[:-1]), np.exp(p[-1])
z = (t / lam) ** k
return -(e * (np.log(k) - np.log(lam) + (k - 1) * (np.log(t) - np.log(lam))) - z).sum()
r = minimize(nll, np.r_[np.log(8000), np.zeros(x.shape[1] - 1), np.log(5)], method="BFGS")
return r.x
def aft_pred(p, d, q=0.5, h=None):
"""이미 h(운용시간)까지 생존했다는 조건에서 분위수 q 생존 시점(운용시간 h 기준 총 수명)을 돌려준다. q=0.5가 조건부 중앙값."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
h = d["h_since"].values if h is None else h
return lam * ((h / lam) ** k - np.log(q)) ** (1 / k)
def aft_fail_prob(p, d, h0, h1):
"""h0까지 생존한 장비가 h0~h1 사이에 고장 날 조건부 확률."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
s0, s1 = np.exp(-(h0 / lam) ** k), np.exp(-(h1 / lam) ** k)
return 1 - s1 / s0
def fit_cox(d, iters=25):
"""Cox 비례위험(Breslow 근사) 뉴턴-랩슨. 공변량: 험지·염분·혹한·비계획 정비율. 반환: 계수."""
x = aft_x(d)[:, 1:]
o = np.argsort(-d["dur"].values)
x, e = x[o], d["event"].values[o]
b = np.zeros(x.shape[1])
for _ in range(iters):
r = np.exp(x @ b)
S0 = np.cumsum(r)
S1 = np.cumsum(r[:, None] * x, axis=0)
S2 = np.cumsum(r[:, None, None] * x[:, :, None] * x[:, None, :], axis=0)
m1 = S1 / S0[:, None]
g = ((x - m1) * e[:, None]).sum(axis=0)
H = -sum(e[i] * (S2[i] / S0[i] - np.outer(m1[i], m1[i])) for i in np.where(e == 1)[0])
step = np.linalg.solve(H - 1e-6 * np.eye(len(b)), g)
b = b - step
if np.abs(step).max() < 1e-7:
break
return b
def cindex(time, event, risk):
"""일치도 C-index: 위험 점수가 높은 쪽이 먼저 고장 났는가(관측 가능한 쌍만)."""
t, e, r = np.asarray(time), np.asarray(event), np.asarray(risk)
num = den = 0.0
for i in np.where(e == 1)[0]:
m = t > t[i]
den += m.sum()
num += (r[m] < r[i]).sum() + 0.5 * (r[m] == r[i]).sum()
return num / max(den, 1)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 잔여 수명(RUL, 일) 예측 엔진: 와이블 AFT 조건부 중앙값 + 그래디언트 부스팅(고장 이력 시점 표본 학습)의 평균
# ---------------------------------------------------------------------------
from sklearn.ensemble import GradientBoostingRegressor
def fit_engine(train, seed=0):
p = fit_aft(train)
sn = add_features(snapshots(train, 5, seed))
sn["aft"] = (aft_pred(p, sn) - sn["h_since"]) / sn["rate"]
g = GradientBoostingRegressor(n_estimators=250, max_depth=3, learning_rate=0.06, random_state=0).fit(sn[FEATS + ["aft"]], sn["rul_d"])
return {"p": p, "g": g}
def predict_rul(eng, d):
"""d는 add_features를 거친 표. 반환: (앙상블 RUL, AFT RUL, GBR RUL) - 일 단위."""
aft = ((aft_pred(eng["p"], d) - d["h_since"]) / d["rate"]).values
x = d[FEATS].copy()
x["aft"] = aft
gb = eng["g"].predict(x)
return np.clip(0.5 * aft + 0.5 * gb, 1, None), aft, gb
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 계획정비(창정비) 일정 추천: 월별 정비창 능력·예산·시효성 부품 키트(유효기한) 제약 아래 기대 비용 최소화 (비용 단가는 교육용 가정값, 만원)
# ---------------------------------------------------------------------------
from scipy.optimize import milp, LinearConstraint, Bounds
from scipy.stats import norm
C_PLAN, C_FAIL, C_WASTE_KH, C_KIT = 600.0, 2400.0, 40.0, 120.0 # 계획 창정비, 비계획 고장(긴급 정비 포함), 미사용 수명 1000h당 손실, 시효성 부품 키트
HORIZON, MDAYS = 12, 30 # 계획 12개월, 월 30일
CAP_M, BUDGET_M = 22, 15000.0 # 월별 정비창 능력(대), 월별 예산
LOTS = [(30, 70), (45, 170), (45, 390), (70, 800)] # 시효성 키트 로트(수량, 유효기한까지 남은 일수)
KIT_MARGIN = 14 # 정비일 기준 최소 유효기간 여유(일)
def exp_cost_matrix(pred, rate, month_days):
"""장비 i를 월 m에 정비할 때의 기대 비용. 참 RUL ~ N(pred, sigma), sigma = 35일 + 0.12 x pred (가정). 마지막 열은 계획 기간 밖(미정비)."""
sig = 35 + 0.12 * pred
C = np.zeros((len(pred), len(month_days) + 1))
for j, day in enumerate(month_days):
z = (day - pred) / sig
pf = norm.cdf(z)
unused = np.maximum(pred - day, 0) * rate * 0.6 # 정비를 앞당긴 만큼 버리는 운용시간 h(근사)
C[:, j] = C_PLAN + C_KIT + pf * (C_FAIL - C_PLAN) + C_WASTE_KH * unused / 1000 * (1 - pf)
C[:, -1] = norm.cdf((HORIZON * MDAYS - pred) / sig) * C_FAIL
return C
def plan_milp(pred, rate, cap=CAP_M, budget=BUDGET_M, lots=LOTS, use_kit=True):
"""정수계획으로 장비별 정비 월(0~11) 또는 미정비(-1)를 배정한다. 후보는 예측 RUL이 계획 기간 +90일 이내인 장비."""
md = np.arange(1, HORIZON + 1) * MDAYS
cand = np.where(pred < HORIZON * MDAYS + 90)[0]
n, M = len(cand), HORIZON
C = exp_cost_matrix(pred[cand], rate[cand], md)
c = C.ravel()
A, lo, hi = [], [], []
for i in range(n): # 장비당 한 번만 배정(월 12칸 + 미정비)
row = np.zeros(n * (M + 1)); row[i * (M + 1):(i + 1) * (M + 1)] = 1; A.append(row); lo.append(1); hi.append(1)
for m in range(M): # 월별 능력·예산
row = np.zeros(n * (M + 1)); row[m::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(cap)
A.append(row * (C_PLAN + C_KIT)); lo.append(0); hi.append(budget)
if use_kit: # 시효성 키트: m월 이후 정비 수요는 그 시점에 유효기간이 남은 로트 수량을 넘을 수 없다
for m in range(M):
sup = sum(q for q, e in lots if e >= md[m] + KIT_MARGIN)
row = np.zeros(n * (M + 1))
for mm in range(m, M):
row[mm::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(sup)
r = milp(c, constraints=LinearConstraint(np.array(A), lo, hi), integrality=np.ones(len(c)), bounds=Bounds(0, 1), options={"time_limit": 20})
month = np.full(len(pred), -1)
if r.x is not None:
x = np.round(r.x).reshape(n, M + 1)
assign = np.where(x[:, :M].sum(axis=1) > 0.5, x[:, :M].argmax(axis=1), -1)
month[cand] = assign
return month
def plan_greedy(pred, margin=45, cap=CAP_M, budget=BUDGET_M):
"""간단한 탐욕 배정: 정비 목표일(예측 RUL - margin)이 빠른 순으로 목표일 이전 가장 늦은 빈 월, 없으면 가장 이른 빈 월에 배정."""
cap = min(cap, int(budget // (C_PLAN + C_KIT))) # 월 예산이 허용하는 대수로 능력을 제한
month = np.full(len(pred), -1)
left = np.full(HORIZON, cap)
for i in np.argsort(pred):
if pred[i] > HORIZON * MDAYS + 60: break
tgt = int(np.clip((pred[i] - margin) // MDAYS - 1, 0, HORIZON - 1))
for m in list(range(tgt, -1, -1)) + list(range(tgt + 1, HORIZON)):
if left[m] > 0:
month[i] = m; left[m] -= 1; break
return month
def plan_fixed(lin_rul, cap=CAP_M, budget=BUDGET_M):
"""관행: 규정 주기 도달 순서대로(선형 잔여 일수) 월별 능력 안에서 배정."""
return plan_greedy(np.asarray(lin_rul, float), margin=0, cap=cap, budget=budget)
def evaluate_plan(month, true_rul, rate, cost_kit=True):
"""참 RUL로 계획의 실제 비용을 계산한다: 정비일 이전에 고장 -> 비계획 고장 비용, 아니면 계획 비용 + 미사용 수명 손실."""
day = np.where(month >= 0, (month + 1) * MDAYS, 10 ** 9)
fail = (true_rul < day) & (true_rul < HORIZON * MDAYS)
planned = (month >= 0) & ~fail
waste = np.where(planned, (true_rul - day) * rate * 0.6, 0) * C_WASTE_KH / 1000
cost = fail * C_FAIL + planned * (C_PLAN + (C_KIT if cost_kit else 0)) + waste
return {"total": float(cost.sum()), "fail": int(fail.sum()), "planned": int(planned.sum()), "waste": float(waste.sum()), "unsched_fail": int(((month < 0) & fail).sum())}
def kit_fefo(month, lots=LOTS):
"""유효기한이 빠른 로트부터 쓰는 FEFO로 키트를 배정하고 정비일에 유효기간이 부족한 배정 수와 만료 폐기 수량을 센다."""
q = [[qty, e] for qty, e in sorted(lots, key=lambda t: t[1])]
short = 0
for m in sorted(month[month >= 0]):
day = (m + 1) * MDAYS
for lot in q:
if lot[1] >= day + KIT_MARGIN and lot[0] > 0:
lot[0] -= 1; break
else:
short += 1
expired = sum(l[0] for l in q if l[1] < HORIZON * MDAYS)
return {"short": short, "expired_unused": expired, "left": sum(l[0] for l in q)}
# ---------------------------------------------------------------------------
OPS_W = {"normal": 1.00, "rough": 1.55, "salt": 1.15, "cold": 1.30} # 운영 기간의 숨은 환경 가중치(도입 당시와 달라진 가정)
tr = make_fleet(900, seed=1)
eng0 = fit_engine(tr)
A = make_fleet(700, seed=60, env_w=OPS_W, shift=1.2) # AAR 대상 운영 장비
B = make_fleet(700, seed=70, env_w=OPS_W, shift=1.2) # 재보정에 쓰지 않는 보류 장비
# ---------------- 1) AAR: 운영 예측 오차 분석 ----------------
Af = add_features(A)
al = Af[Af.event == 0].reset_index(drop=True)
pred, _, _ = predict_rul(eng0, al)
y = al["rul_d"].values
al["err"] = pred - y # 양수 = 수명을 길게 예측(위험)
big = y >= 90
print(f"[1] AAR 대상 운영 장비 {len(A)}대(생존 {len(al)}대, 고장 {int(A.event.sum())}대): 현행 엔진 MAE {mae(y, pred):.0f}일, RMSE {rmse(y, pred):.0f}일, MAPE {mape(y[big], pred[big]):.1f}%, R2 {r2_score(y, pred):.3f}, 평균 편향 {al.err.mean():+.0f}일")
ct = pd.DataFrame([(k, int(m.sum()), mae(y[m.values], pred[m.values]), al.err[m].mean()) for k, m in {"일반": al.env == "normal", "험지": al.env == "rough", "염분": al.env == "salt", "혹한": al.env == "cold",
"일 가동률 4h 이상": al.rate >= 4, "비계획 정비율 상위 25%": al.unpl_per_kh >= al.unpl_per_kh.quantile(.75)}.items()], columns=["조건", "대수", "MAE(일)", "평균 편향(일)"])
print(ct.round(0).to_string(index=False))
worst = al.reindex(al.err.abs().sort_values(ascending=False).index).head(40)
def cause(r):
if r.env == "rough" and r.err > 0: return "험지 가혹도 과소 반영"
if r.env == "cold" and r.err > 0: return "혹한 가혹도 과소 반영"
if r.env == "salt" and r.err < 0: return "염분 가혹도 과대 반영"
if r.rate >= 4: return "가동률 상승(운용 강도 변화)"
if r.unpl_per_kh >= 1.2: return "비계획 정비 다수"
return "기타(개체 편차)"
cs = worst.apply(cause, axis=1).value_counts()
print(" 오차 큰 상위 40대의 원인 분류:", ", ".join(f"{k} {v}대({v / 40:.0%})" for k, v in cs.items()))
# ---------------- 2) 재보정: 환경 가중치 -> 운영 데이터 재학습 ----------------
pa = fit_aft(A)
life_ratio = np.exp(pa[1:4])
W1 = {"normal": 1.0, "rough": 1 / life_ratio[0], "salt": 1 / life_ratio[1], "cold": 1 / life_ratio[2]}
wt = pd.DataFrame({"도입 시 추정(W0)": W0, "재보정(운영 AFT)": W1, "숨은 정답(운영)": OPS_W}).loc[ENVS]
wt["변화율"] = wt["재보정(운영 AFT)"] / wt["도입 시 추정(W0)"] - 1
print("\n[2] 환경 가중치 매트릭스(클수록 빨리 노후)")
print(wt.round(3).to_string())
flag = wt[wt["변화율"].abs() > 0.2].index.tolist()
print(f" 변화율 20% 초과 환경 {flag or '없음'}는 정비계획관 검토·군수지휘관 승인 뒤 적용한다(자동 반영 금지).")
eng1 = fit_engine(A, seed=1)
Bf0, Bf1 = add_features(B, W0), add_features(B, W1)
bl0, bl1 = Bf0[Bf0.event == 0].reset_index(drop=True), Bf1[Bf1.event == 0].reset_index(drop=True)
p0, _, _ = predict_rul(eng0, bl0)
p1, _, _ = predict_rul(eng1, bl1)
yb = bl0["rul_d"].values
bb = yb >= 90
tab = pd.DataFrame({"재보정 전": [mae(yb, p0), rmse(yb, p0), mape(yb[bb], p0[bb]), r2_score(yb, p0), (p0 - yb).mean()], "재보정 후": [mae(yb, p1), rmse(yb, p1), mape(yb[bb], p1[bb]), r2_score(yb, p1), (p1 - yb).mean()]},
index=["MAE(일)", "RMSE(일)", "MAPE(%)", "R2", "평균 편향(일)"])
print(f"\n[3] 보류 장비 {len(bl0)}대(재보정에 쓰지 않음) 검증\n" + tab.round(3).to_string())
print(f" MAE {1 - tab.loc['MAE(일)', '재보정 후'] / tab.loc['MAE(일)', '재보정 전']:.1%} 감소, 환경별 편향(전 -> 후):", ", ".join(f"{e} {(p0 - yb)[(bl0.env == e).values].mean():+.0f} -> {(p1 - yb)[(bl0.env == e).values].mean():+.0f}" for e in ENVS))
# ---------------- 3) ROI 3대 지표(교육용 가정) ----------------
rate_b, true_b = bl0["rate"].values, yb
m0, m1 = plan_milp(p0, rate_b), plan_milp(p1, rate_b)
r0, r1 = evaluate_plan(m0, true_b, rate_b), evaluate_plan(m1, true_b, rate_b)
D_FAIL, D_PLAN = 45, 12 # 비계획 고장/계획 창정비 가동 중단 일수(가정)
down0, down1 = r0["fail"] * D_FAIL + r0["planned"] * D_PLAN, r1["fail"] * D_FAIL + r1["planned"] * D_PLAN
SCALE, INV, OPEX = 1.0, 12000.0, 1200.0 # 관리 장비 규모 배수(연), 도입 비용, 연 운영비(만원)
def roi(c0, c1, d0, d1, scale=SCALE, inv=INV):
save = (c0 - c1) * scale - OPEX
return save, 1 - d1 / d0, inv / max(save / 12, 1e-9)
sv, tr_, pb = roi(r0["total"], r1["total"], down0, down1)
print(f"\n[4] 12개월 계획 비교(보류 장비): 비용 {r0['total']:,.0f} -> {r1['total']:,.0f}만원, 비계획 고장 {r0['fail']} -> {r1['fail']}대, 계획 정비 {r0['planned']} -> {r1['planned']}대")
print(f" ROI 3대 지표 [교육용 가정: 관리 규모 x{SCALE:.0f}, 도입 {INV:.0f}만원, 운영 {OPEX:.0f}만원/년, 고장 중단 {D_FAIL}일·창정비 {D_PLAN}일]")
print(f" (1) 연간 비용 절감액 {sv:,.0f}만원 (2) 가동 중단일수 단축률 {tr_:.1%}({down0:,}일 -> {down1:,}일) (3) 투자 회수 기간 {pb:.1f}개월")
print(" 산식: 절감액 = (재보정 전 - 후 12개월 비용) x 규모 - 운영비, 단축률 = 1 - 후 중단일수/전 중단일수, 회수 기간 = 도입비 / (월 절감액)")
print(" 관리 규모 민감도:", "; ".join(f"x{s} -> 절감 {roi(r0['total'], r1['total'], down0, down1, scale=s)[0]:,.0f}만원, 회수 {roi(r0['total'], r1['total'], down0, down1, scale=s)[2]:.0f}개월" for s in (0.5, 1, 2, 3)))
half = r0["total"] - (r0["total"] - r1["total"]) / 2
print(f" 비관(재보정 효과 절반만 실현, x1): 절감 {roi(r0['total'], half, down0, down1, scale=1)[0]:,.0f}만원, 회수 {roi(r0['total'], half, down0, down1, scale=1)[2]:.0f}개월")
print(" 가정값이므로 재정·기획 부서와 실제 단가로 다시 산정해야 한다. 재보정 가중치는 군수지휘관 승인 후 운영에 반영한다.")
fig, ax = plt.subplots(1, 3, figsize=(15, 4.3))
ax[0].bar(range(len(ct)), ct["평균 편향(일)"]); ax[0].set_xticks(range(len(ct))); ax[0].set_xticklabels([f"C{i + 1}" for i in range(len(ct))]); ax[0].set_title("AAR: mean bias by condition (days)")
x = np.arange(4); ax[1].bar(x - .27, wt["도입 시 추정(W0)"], .27, label="initial W0"); ax[1].bar(x, wt["재보정(운영 AFT)"], .27, label="recalibrated"); ax[1].bar(x + .27, wt["숨은 정답(운영)"], .27, label="hidden truth")
ax[1].set_xticks(x); ax[1].set_xticklabels(ENVS); ax[1].legend(); ax[1].set_title("Environment weight matrix")
ax[2].scatter(yb, p0, s=6, alpha=.5, label="before"); ax[2].scatter(yb, p1, s=6, alpha=.5, label="after"); ax[2].plot([0, 4000], [0, 4000], "k--"); ax[2].legend(); ax[2].set_title("Hold-out RUL: true vs predicted (days)")
plt.tight_layout(); plt.show()
In [ ]:
"""장비 보유 현황 갱신, 정합성 검증, 분포 변화(PSI) 점검, 자동 재학습(챔피언/챌린저), 잔여 수명(RUL) 예측, 정비창 능력·예산·시효성 키트 제약의 12개월 계획정비 일정 추천과 승인대기·도태 후보 분리까지 한 번에 실행하는 통합 파이프라인을 완성하고 운영 문서 완비율·자체 점검·인터페이스 명세·브리핑 요약을 출력한다."""
import importlib
import json
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score
# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 장비 보유 현황·고장 이력 생성기 (실제 부대·장비·수치 없음)
# 가상 부대 6개(U1~U6)의 장비가 운용 환경(normal/rough/salt/cold)에 따라 다른 속도로 노후화한다. 창정비 후 운용시간(h_since)이
# 장비마다 숨은 수명(t_fail)을 넘으면 고장(event=1) 이력, 넘지 않았으면 현재까지 생존(event=0, 관측 중도절단)이다.
# ---------------------------------------------------------------------------
ENVS = ["normal", "rough", "salt", "cold"]
ENV_W = {"normal": 1.00, "rough": 1.35, "salt": 1.25, "cold": 1.12} # 환경 가중치(숨은 정답, 가정값): 클수록 수명이 짧다
ENV_KO = {"normal": "일반", "rough": "험지", "salt": "염분", "cold": "혹한"}
UNITS = ["U1", "U2", "U3", "U4", "U5", "U6"]
UNIT_ENV = {"U1": [.7, .1, .1, .1], "U2": [.2, .6, .1, .1], "U3": [.2, .1, .6, .1], "U4": [.3, .1, .1, .5], "U5": [.5, .3, .1, .1], "U6": [.4, .1, .4, .1]}
ETYPES = ["truck", "generator", "trailer", "radio"]
CYCLE = 8000.0 # 규정 창정비 주기(운용시간 h, 가정값)
KSHAPE = 12.0 # 수명 분포 와이블 형상(마모형 고장, 클수록 수명이 모여 있다)
INSP_NOISE = 0.025 # 정비 점검 마모 비율 측정 오차(가정값)
def make_fleet(n=900, seed=0, env_w=None, shift=1.0):
"""장비 n대. env_w로 숨은 환경 가중치를 바꿔 운영 분포 변화를 모의할 수 있다. shift>1이면 가동률이 높아진다."""
rng = np.random.default_rng(seed)
w = env_w or ENV_W
unit = rng.choice(UNITS, n)
env = np.array([rng.choice(ENVS, p=UNIT_ENV[u]) for u in unit])
et = rng.choice(ETYPES, n, p=[.4, .2, .2, .2])
rate = rng.lognormal(np.log(3.0 * shift), 0.25, n) # 일 평균 운용시간(h/일)
h_since = rng.uniform(300, 11500, n) # 마지막 창정비 이후 운용시간
n_ovh = rng.integers(0, 3, n)
hours = h_since + n_ovh * rng.uniform(5000, 8000, n) # 누적 운용시간
km = hours * rng.normal(17, 3, n) * np.where(env == "rough", 0.8, 1.0) # 누적 기동거리(km)
age_y = np.clip(hours / (rate * 365) * rng.uniform(0.9, 1.15, n), 0.5, 25)
unpl = rng.poisson(0.35 * h_since / 1000 * np.array([w[e] for e in env]) / 1.3) # 비계획 정비 건수(환경이 가혹할수록 많다)
n_maint = rng.poisson(0.6 * h_since / 1000 + 1) # 계획 정비 건수
scale = 9600 / np.array([w[e] for e in env]) * 0.985 ** unpl * np.exp(rng.normal(0, 0.04, n))
t_fail = scale * rng.weibull(KSHAPE, n) # 창정비 후 고장까지 운용시간(숨은 정답)
event = (t_fail <= h_since).astype(int)
dur = np.where(event == 1, t_fail, h_since)
df = pd.DataFrame({"eq_id": [f"EQ{seed:02d}{i:04d}" for i in range(n)], "unit": unit, "env": env, "etype": et, "rate": rate, "hours": hours, "km": km, "age_y": age_y, "n_ovh": n_ovh,
"h_since": h_since, "unpl": unpl, "n_maint": n_maint, "t_fail": t_fail, "event": event, "dur": dur})
df["insp"] = np.clip(dur / t_fail + rng.normal(0, INSP_NOISE, n), 0, 1.3) # 마지막 정비 점검에서 측정한 마모 비율(소모 수명/총 수명, 오차 포함)
df["rul_h"] = np.where(event == 0, t_fail - h_since, np.nan) # 생존 장비의 남은 운용시간(정답)
df["rul_d"] = df["rul_h"] / df["rate"] # 남은 일수(정답)
return df
def mape(y, p):
y, p = np.asarray(y, float), np.asarray(p, float)
return float(np.mean(np.abs(y - p) / np.abs(y)) * 100)
def mae(y, p):
return float(np.mean(np.abs(np.asarray(y, float) - np.asarray(p, float))))
def rmse(y, p):
return float(np.sqrt(np.mean((np.asarray(y, float) - np.asarray(p, float)) ** 2)))
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 노후화 파생 피처와 학습 표본 만들기. 환경 가중치 W0는 도입 시점의 전문가 추정값이며(실제 가중치와 다름) M9에서 재보정한다.
# ---------------------------------------------------------------------------
W0 = {"normal": 1.00, "rough": 1.15, "salt": 1.10, "cold": 1.05}
FEATS = ["h_since", "hours", "km", "age_y", "rate", "n_maint", "unpl", "unpl_per_kh", "km_per_h", "n_ovh", "k_rough", "k_salt", "k_cold", "k_truck", "k_gen", "k_trailer",
"aging_idx", "lin_rul_d", "hs_frac", "insp", "insp_eta_h"]
def add_features(df, w=None):
d = df.copy()
w = w or W0
d["unpl_per_kh"] = d["unpl"] / (d["h_since"] / 1000 + 0.5)
d["km_per_h"] = d["km"] / d["hours"]
for e in ("rough", "salt", "cold"):
d["k_" + e] = (d["env"] == e).astype(int)
d["k_truck"], d["k_gen"], d["k_trailer"] = (d["etype"] == "truck").astype(int), (d["etype"] == "generator").astype(int), (d["etype"] == "trailer").astype(int)
d["aging_idx"] = d["h_since"] * d["env"].map(w) / CYCLE # 환경 가중 노후 지수(1이면 창정비 주기에 도달)
d["lin_rul_d"] = (CYCLE - d["h_since"]) / d["rate"] # 선형 감가상각 잔여 일수(Baseline)
d["hs_frac"] = d["h_since"] / CYCLE
d["insp_eta_h"] = d["h_since"] / np.maximum(d["insp"], 0.05) # 점검 마모 비율로 환산한 총 수명 추정(h)
return d
def snapshots(df, k=4, seed=0):
"""고장 이력 장비의 과거 시점 표본: 고장 전 임의 시점의 점검값·운용시간과 그때의 남은 일수. (생존 장비의 미래는 알 수 없으므로 고장 이력만 쓴다)"""
rng = np.random.default_rng(seed)
f = df[df["event"] == 1]
out = []
for _ in range(k):
s = f.copy()
h = s["t_fail"] * rng.uniform(0.2, 0.98, len(s))
ratio = h / s["h_since"]
s["hours"] = s["hours"] - (s["h_since"] - h)
s["km"] = s["km"] * ratio.clip(upper=1)
s["unpl"] = np.floor(s["unpl"] * ratio)
s["n_maint"] = np.floor(s["n_maint"] * ratio)
s["insp"] = np.clip(h / s["t_fail"] + rng.normal(0, INSP_NOISE, len(s)), 0, 1.3)
s["h_since"] = h
s["rul_d"] = (s["t_fail"] - h) / s["rate"]
s["event"] = 0
out.append(s)
return pd.concat(out, ignore_index=True)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 생존 분석 도구(scipy만 사용): 와이블 AFT(가속 수명) 모형과 Cox 비례위험 모형
# ---------------------------------------------------------------------------
from scipy.optimize import minimize
def aft_x(d):
return np.c_[np.ones(len(d)), d["env"].eq("rough"), d["env"].eq("salt"), d["env"].eq("cold"), d["unpl"] / (d["h_since"] / 1000 + 0.5)].astype(float)
AFT_COLS = ["절편", "험지", "염분", "혹한", "비계획 정비율"]
def fit_aft(d):
"""창정비 후 운용시간 dur(고장이면 event=1, 생존이면 중도절단)에 대한 와이블 AFT 최대우도 추정. 반환: 계수(마지막은 log 형상)."""
x, t, e = aft_x(d), d["dur"].values, d["event"].values
def nll(p):
lam, k = np.exp(x @ p[:-1]), np.exp(p[-1])
z = (t / lam) ** k
return -(e * (np.log(k) - np.log(lam) + (k - 1) * (np.log(t) - np.log(lam))) - z).sum()
r = minimize(nll, np.r_[np.log(8000), np.zeros(x.shape[1] - 1), np.log(5)], method="BFGS")
return r.x
def aft_pred(p, d, q=0.5, h=None):
"""이미 h(운용시간)까지 생존했다는 조건에서 분위수 q 생존 시점(운용시간 h 기준 총 수명)을 돌려준다. q=0.5가 조건부 중앙값."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
h = d["h_since"].values if h is None else h
return lam * ((h / lam) ** k - np.log(q)) ** (1 / k)
def aft_fail_prob(p, d, h0, h1):
"""h0까지 생존한 장비가 h0~h1 사이에 고장 날 조건부 확률."""
lam, k = np.exp(aft_x(d) @ p[:-1]), np.exp(p[-1])
s0, s1 = np.exp(-(h0 / lam) ** k), np.exp(-(h1 / lam) ** k)
return 1 - s1 / s0
def fit_cox(d, iters=25):
"""Cox 비례위험(Breslow 근사) 뉴턴-랩슨. 공변량: 험지·염분·혹한·비계획 정비율. 반환: 계수."""
x = aft_x(d)[:, 1:]
o = np.argsort(-d["dur"].values)
x, e = x[o], d["event"].values[o]
b = np.zeros(x.shape[1])
for _ in range(iters):
r = np.exp(x @ b)
S0 = np.cumsum(r)
S1 = np.cumsum(r[:, None] * x, axis=0)
S2 = np.cumsum(r[:, None, None] * x[:, :, None] * x[:, None, :], axis=0)
m1 = S1 / S0[:, None]
g = ((x - m1) * e[:, None]).sum(axis=0)
H = -sum(e[i] * (S2[i] / S0[i] - np.outer(m1[i], m1[i])) for i in np.where(e == 1)[0])
step = np.linalg.solve(H - 1e-6 * np.eye(len(b)), g)
b = b - step
if np.abs(step).max() < 1e-7:
break
return b
def cindex(time, event, risk):
"""일치도 C-index: 위험 점수가 높은 쪽이 먼저 고장 났는가(관측 가능한 쌍만)."""
t, e, r = np.asarray(time), np.asarray(event), np.asarray(risk)
num = den = 0.0
for i in np.where(e == 1)[0]:
m = t > t[i]
den += m.sum()
num += (r[m] < r[i]).sum() + 0.5 * (r[m] == r[i]).sum()
return num / max(den, 1)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 잔여 수명(RUL, 일) 예측 엔진: 와이블 AFT 조건부 중앙값 + 그래디언트 부스팅(고장 이력 시점 표본 학습)의 평균
# ---------------------------------------------------------------------------
from sklearn.ensemble import GradientBoostingRegressor
def fit_engine(train, seed=0):
p = fit_aft(train)
sn = add_features(snapshots(train, 5, seed))
sn["aft"] = (aft_pred(p, sn) - sn["h_since"]) / sn["rate"]
g = GradientBoostingRegressor(n_estimators=250, max_depth=3, learning_rate=0.06, random_state=0).fit(sn[FEATS + ["aft"]], sn["rul_d"])
return {"p": p, "g": g}
def predict_rul(eng, d):
"""d는 add_features를 거친 표. 반환: (앙상블 RUL, AFT RUL, GBR RUL) - 일 단위."""
aft = ((aft_pred(eng["p"], d) - d["h_since"]) / d["rate"]).values
x = d[FEATS].copy()
x["aft"] = aft
gb = eng["g"].predict(x)
return np.clip(0.5 * aft + 0.5 * gb, 1, None), aft, gb
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# 계획정비(창정비) 일정 추천: 월별 정비창 능력·예산·시효성 부품 키트(유효기한) 제약 아래 기대 비용 최소화 (비용 단가는 교육용 가정값, 만원)
# ---------------------------------------------------------------------------
from scipy.optimize import milp, LinearConstraint, Bounds
from scipy.stats import norm
C_PLAN, C_FAIL, C_WASTE_KH, C_KIT = 600.0, 2400.0, 40.0, 120.0 # 계획 창정비, 비계획 고장(긴급 정비 포함), 미사용 수명 1000h당 손실, 시효성 부품 키트
HORIZON, MDAYS = 12, 30 # 계획 12개월, 월 30일
CAP_M, BUDGET_M = 22, 15000.0 # 월별 정비창 능력(대), 월별 예산
LOTS = [(30, 70), (45, 170), (45, 390), (70, 800)] # 시효성 키트 로트(수량, 유효기한까지 남은 일수)
KIT_MARGIN = 14 # 정비일 기준 최소 유효기간 여유(일)
def exp_cost_matrix(pred, rate, month_days):
"""장비 i를 월 m에 정비할 때의 기대 비용. 참 RUL ~ N(pred, sigma), sigma = 35일 + 0.12 x pred (가정). 마지막 열은 계획 기간 밖(미정비)."""
sig = 35 + 0.12 * pred
C = np.zeros((len(pred), len(month_days) + 1))
for j, day in enumerate(month_days):
z = (day - pred) / sig
pf = norm.cdf(z)
unused = np.maximum(pred - day, 0) * rate * 0.6 # 정비를 앞당긴 만큼 버리는 운용시간 h(근사)
C[:, j] = C_PLAN + C_KIT + pf * (C_FAIL - C_PLAN) + C_WASTE_KH * unused / 1000 * (1 - pf)
C[:, -1] = norm.cdf((HORIZON * MDAYS - pred) / sig) * C_FAIL
return C
def plan_milp(pred, rate, cap=CAP_M, budget=BUDGET_M, lots=LOTS, use_kit=True):
"""정수계획으로 장비별 정비 월(0~11) 또는 미정비(-1)를 배정한다. 후보는 예측 RUL이 계획 기간 +90일 이내인 장비."""
md = np.arange(1, HORIZON + 1) * MDAYS
cand = np.where(pred < HORIZON * MDAYS + 90)[0]
n, M = len(cand), HORIZON
C = exp_cost_matrix(pred[cand], rate[cand], md)
c = C.ravel()
A, lo, hi = [], [], []
for i in range(n): # 장비당 한 번만 배정(월 12칸 + 미정비)
row = np.zeros(n * (M + 1)); row[i * (M + 1):(i + 1) * (M + 1)] = 1; A.append(row); lo.append(1); hi.append(1)
for m in range(M): # 월별 능력·예산
row = np.zeros(n * (M + 1)); row[m::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(cap)
A.append(row * (C_PLAN + C_KIT)); lo.append(0); hi.append(budget)
if use_kit: # 시효성 키트: m월 이후 정비 수요는 그 시점에 유효기간이 남은 로트 수량을 넘을 수 없다
for m in range(M):
sup = sum(q for q, e in lots if e >= md[m] + KIT_MARGIN)
row = np.zeros(n * (M + 1))
for mm in range(m, M):
row[mm::(M + 1)] = 1
A.append(row); lo.append(0); hi.append(sup)
r = milp(c, constraints=LinearConstraint(np.array(A), lo, hi), integrality=np.ones(len(c)), bounds=Bounds(0, 1), options={"time_limit": 20})
month = np.full(len(pred), -1)
if r.x is not None:
x = np.round(r.x).reshape(n, M + 1)
assign = np.where(x[:, :M].sum(axis=1) > 0.5, x[:, :M].argmax(axis=1), -1)
month[cand] = assign
return month
def plan_greedy(pred, margin=45, cap=CAP_M, budget=BUDGET_M):
"""간단한 탐욕 배정: 정비 목표일(예측 RUL - margin)이 빠른 순으로 목표일 이전 가장 늦은 빈 월, 없으면 가장 이른 빈 월에 배정."""
cap = min(cap, int(budget // (C_PLAN + C_KIT))) # 월 예산이 허용하는 대수로 능력을 제한
month = np.full(len(pred), -1)
left = np.full(HORIZON, cap)
for i in np.argsort(pred):
if pred[i] > HORIZON * MDAYS + 60: break
tgt = int(np.clip((pred[i] - margin) // MDAYS - 1, 0, HORIZON - 1))
for m in list(range(tgt, -1, -1)) + list(range(tgt + 1, HORIZON)):
if left[m] > 0:
month[i] = m; left[m] -= 1; break
return month
def plan_fixed(lin_rul, cap=CAP_M, budget=BUDGET_M):
"""관행: 규정 주기 도달 순서대로(선형 잔여 일수) 월별 능력 안에서 배정."""
return plan_greedy(np.asarray(lin_rul, float), margin=0, cap=cap, budget=budget)
def evaluate_plan(month, true_rul, rate, cost_kit=True):
"""참 RUL로 계획의 실제 비용을 계산한다: 정비일 이전에 고장 -> 비계획 고장 비용, 아니면 계획 비용 + 미사용 수명 손실."""
day = np.where(month >= 0, (month + 1) * MDAYS, 10 ** 9)
fail = (true_rul < day) & (true_rul < HORIZON * MDAYS)
planned = (month >= 0) & ~fail
waste = np.where(planned, (true_rul - day) * rate * 0.6, 0) * C_WASTE_KH / 1000
cost = fail * C_FAIL + planned * (C_PLAN + (C_KIT if cost_kit else 0)) + waste
return {"total": float(cost.sum()), "fail": int(fail.sum()), "planned": int(planned.sum()), "waste": float(waste.sum()), "unsched_fail": int(((month < 0) & fail).sum())}
def kit_fefo(month, lots=LOTS):
"""유효기한이 빠른 로트부터 쓰는 FEFO로 키트를 배정하고 정비일에 유효기간이 부족한 배정 수와 만료 폐기 수량을 센다."""
q = [[qty, e] for qty, e in sorted(lots, key=lambda t: t[1])]
short = 0
for m in sorted(month[month >= 0]):
day = (m + 1) * MDAYS
for lot in q:
if lot[1] >= day + KIT_MARGIN and lot[0] > 0:
lot[0] -= 1; break
else:
short += 1
expired = sum(l[0] for l in q if l[1] < HORIZON * MDAYS)
return {"short": short, "expired_unused": expired, "left": sum(l[0] for l in q)}
# ---------------------------------------------------------------------------
libs = ["numpy", "pandas", "scipy", "sklearn", "matplotlib"]
env_ok = sum(importlib.util.find_spec(m) is not None for m in libs)
print(f"[0] 환경 점검: Python {sys.version_info.major}.{sys.version_info.minor}, 필수 모듈 {env_ok}/{len(libs)}")
# ---------------- 1) 보유 현황 갱신·정합성 검증(미션 1~2 요약) ----------------
OPS_W = {"normal": 1.00, "rough": 1.55, "salt": 1.15, "cold": 1.30}
base = make_fleet(900, seed=1) # 기존 보유 현황
new = make_fleet(500, seed=90, env_w=OPS_W, shift=1.1) # 신규 입력분(운영 환경이 달라짐)
raw = pd.concat([base, new], ignore_index=True)
chk = {"장비 ID 유일": not raw["eq_id"].duplicated().any(), "결측 0": int(raw[["rate", "hours", "h_since", "insp", "unpl"]].isna().sum().sum()) == 0, "환경 값 4종": bool(raw["env"].isin(ENVS).all()),
"부대 코드 등록": bool(raw["unit"].isin(UNITS).all()), "운용시간 0 이상": bool((raw["hours"] >= 0).all()), "창정비 후 <= 누적": bool((raw["h_since"] <= raw["hours"] + 1e-6).all()),
"기동 속도 5~40km/h": bool((raw["km"] / raw["hours"]).between(5, 40).all()), "점검 마모 0~1.3": bool(raw["insp"].between(0, 1.3).all())}
print(f"[1] 신규 {len(new)}대 포함 {len(raw)}대 정합성 검증 {sum(chk.values())}/{len(chk)} 통과")
# ---------------- 2) 분포 변화 점검과 자동 재학습 ----------------
def psi(a, b, bins=8):
q = np.unique(np.quantile(a, np.linspace(0, 1, bins + 1)))
q[0], q[-1] = -np.inf, np.inf
pa, pb = np.histogram(a, q)[0] / len(a) + 1e-4, np.histogram(b, q)[0] / len(b) + 1e-4
return float(np.sum((pa - pb) * np.log(pa / pb)))
fb, fn = add_features(base), add_features(new)
drift = {c: psi(fb[c], fn[c]) for c in ["rate", "unpl_per_kh", "insp_eta_h", "aging_idx", "km_per_h"]}
worst = max(drift, key=drift.get)
print("[2] 분포 변화 PSI(기존 vs 신규):", ", ".join(f"{k} {v:.2f}" for k, v in drift.items()), f"-> 가장 큰 변화 {worst}(0.2 이상이면 재학습 검토)")
new_tr, new_ho = new.iloc[:300].reset_index(drop=True), new.iloc[300:].reset_index(drop=True)
champ = fit_engine(base)
chall = fit_engine(pd.concat([base.sample(300, random_state=0), new_tr], ignore_index=True), seed=1) # 최근 분포를 반영해 기존 일부 + 신규 학습분으로 재학습
ho = add_features(new_ho)
ho = ho[ho.event == 0].reset_index(drop=True)
pc, _, _ = predict_rul(champ, ho)
ph, _, _ = predict_rul(chall, ho)
y = ho["rul_d"].values
big = y >= 90
ac, bc = mae(y, pc), mae(y, ph)
promote = (ac - bc) / ac >= 0.03
print(f" 검증 {len(ho)}대(신규 분포) MAE 챔피언 {ac:.0f}일 vs 챌린저 {bc:.0f}일, MAPE {mape(y[big], pc[big]):.1f}% vs {mape(y[big], ph[big]):.1f}% -> {'챌린저 승격' if promote else '기존 모델 유지'}(MAE 3% 이상 개선일 때만 승격)")
model = chall if promote else champ
# ---------------- 3) RUL 예측과 12개월 정비 일정 ----------------
fl = add_features(new)
al = fl[fl.event == 0].reset_index(drop=True)
al["rul_pred"], _, _ = predict_rul(model, al)
al["month"] = plan_milp(al["rul_pred"].values, al["rate"].values)
mf = plan_fixed(al["lin_rul_d"].values)
r_fix, r_new = evaluate_plan(mf, al["rul_d"].values, al["rate"].values), evaluate_plan(al["month"].values, al["rul_d"].values, al["rate"].values)
kf = kit_fefo(al["month"].values)
print(f"\n[3] 신규 보유 {len(al)}대 12개월 정비 일정 [교육용가상자료]: 정비 배정 {int((al.month >= 0).sum())}대, 월별 " + " ".join(f"{m + 1}월 {c}" for m, c in enumerate(np.bincount(al.loc[al.month >= 0, 'month'], minlength=HORIZON))))
print(f" 규정 주기 고정 일정 대비 실제 비용 {r_fix['total']:,.0f} -> {r_new['total']:,.0f}만원({1 - r_new['total'] / r_fix['total']:.1%} 절감), 비계획 고장 {r_fix['fail']} -> {r_new['fail']}대, 키트 유효기간 부족 {kf['short']}건")
al["처리"] = np.select([al.rul_pred < 90, (al.age_y >= 14) & (al.rul_pred < 365)], ["긴급 점검(현장 확인 후 군수지휘관 승인)", "도태 건의 후보(2단계 결재)"], np.where(al.month >= 0, "정비 일정 승인대기", "정기 점검 유지"))
top = al.sort_values("rul_pred").head(6)
print(pd.DataFrame({"장비": top.eq_id, "부대": top.unit, "환경": top.env.map(ENV_KO), "예측 RUL(일)": top.rul_pred, "권장 월": top.month.map(lambda m: f"{m + 1}월" if m >= 0 else "-"), "처리": top["처리"]}).round(0).to_string(index=False))
vc = al["처리"].value_counts()
print(" 처리 구분:", vc.to_dict(), "/ 자동 확정 0건(모든 일정·도태는 군수지휘관 승인 필요)")
# ---------------- 4) 운영 문서 완비율·인터페이스 명세 ----------------
docs = {"데이터 사전·스키마(M1)": 30, "정제 규칙 11종(M2)": 11, "Baseline·EDA(M3)": 8, "생존 분석 모델(M4)": 12, "RUL 엔진·정비 일정(M5)": 4, "대시보드·설명(M6)": 5, "승인 절차·권한(M7)": 4,
"일괄 시뮬레이션·부하(M8)": 5, "AAR·ROI(M9)": 4, "운영·재학습 절차(M10)": 6, "장애 대응": 5, "인터페이스 가이드": 4}
done = sum(v > 0 for v in docs.values())
iface = {"input": {"fleet": ["eq_id:str", "unit:str", "env:str", "rate:float", "hours:float", "h_since:float", "unpl:int", "insp:float"], "lots": ["lot_id", "qty", "expiry_d"], "format": "json-lines"},
"output": {"plan": ["eq_id", "rul_days", "month", "status"], "approval": "PENDING_ONLY"}, "auth": "role-based(가상)", "note": "상위 군수정보체계는 가상 인터페이스로만 정의"}
print(f"\n[4] 운영 문서 {len(docs)}개 항목 중 작성 완료 {done}개 ({done / len(docs):.0%}) / 인터페이스 명세(JSON):")
print(" ", json.dumps(iface, ensure_ascii=False))
# ---------------- 5) 자체 점검·브리핑 ----------------
pm, _, _ = predict_rul(model, ho)
fin_mae, fin_mape, base_mae = mae(y, pm), mape(y[big], pm[big]), mae(y, ho["lin_rul_d"].values)
checks = [("환경 점검 통과", env_ok == len(libs)), ("데이터 정합성 전 항목 통과", all(chk.values())), ("재학습 파이프라인 정상 실행(챔피언/챌린저 판정 출력)", ac > 0 and bc > 0),
("RUL MAPE 25% 이하(참 RUL 90일 이상)", fin_mape <= 25), ("선형 감가상각 대비 MAE 20% 이상 개선", fin_mae <= 0.8 * base_mae), ("고정 주기 대비 계획 비용 15% 이상 절감", r_new["total"] <= 0.85 * r_fix["total"]),
("자동 확정 0건·키트 유효기간 부족 0건", bool(al["처리"].str.contains("승인|점검|건의").all()) and kf["short"] == 0), ("운영 문서 100% 완비", done == len(docs))]
print("\n[5] 자체 점검")
for n, ok in checks:
print(f" [{'PASS' if ok else 'FAIL'}] {n}")
print("\n[6] 브리핑 요약")
print(f" - 신규 {len(new)}대로 분포 변화를 점검({worst} PSI {drift[worst]:.2f})하고 재학습 판단을 자동 수행했다({'챌린저 승격' if promote else '기존 모델 유지'}).")
print(f" - 신규 분포 RUL MAE {base_mae:.0f}일(선형 감가상각) -> {fin_mae:.0f}일, MAPE {fin_mape:.1f}%. 12개월 계획 비용 {1 - r_new['total'] / r_fix['total']:.1%} 절감, 비계획 고장 {r_fix['fail']} -> {r_new['fail']}대(가정 단가 기준).")
print(" - 정비 일정·긴급 점검·도태 후보는 모두 승인대기로 분리했다. 확정은 군수지휘관(도태는 재정 담당 포함) 결재에 따른다.")
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].bar(list(drift), list(drift.values())); ax[0].axhline(0.2, c="r", ls="--"); ax[0].set_title("Feature drift (PSI), old vs new fleet"); ax[0].tick_params(axis="x", rotation=20)
ax[1].scatter(y, pc, s=6, alpha=.5, label="champion"); ax[1].scatter(y, ph, s=6, alpha=.5, label="challenger"); ax[1].plot([0, 4000], [0, 4000], "k--"); ax[1].legend(); ax[1].set_title("Hold-out RUL: true vs predicted (days)")
plt.suptitle("SUP02 integrated pipeline (synthetic)"); plt.tight_layout(); plt.show()