Files
2026-10-08 16:43:48 +09:00

158 KiB

MNT02 정답 코드 - AI 기반 지능형 부품 단종위험 분석·대응체계

[교육용가상자료] 모든 데이터는 코드 안에서 난수로 만든 가상 자료이며 실제 부대·장비·인원·작전 정보와 무관합니다.

실행 환경

  • Python 3.10 이상, numpy, pandas, scipy, scikit-learn, matplotlib
  • 미션별 코드는 각각 한 파일로 독립 실행됩니다. 외부 파일 없이 코드 안에서 가상 데이터를 생성합니다.
  • 그래프의 글자는 한글 글꼴이 없는 실행 환경에서도 깨지지 않도록 영문으로 표기했습니다. 표·출력 결과는 한글입니다.
  • 난수 시드를 고정해 두어 같은 코드를 실행하면 같은 결과가 나옵니다.

미션 1. 야전 정비창 수리부속 청구/출고 이력 및 장비 가동률 데이터셋 스키마 정의

  • 난이도: 초급 / 예상 소요시간: 1.5시간
  • 코드 설명: 야전 정비창 수리부속 청구/출고 이력과 장비 가동률 데이터의 업무 프로세스 맵, 데이터 스키마(피처 20개 이상)와 데이터 딕셔너리를 정의하고 가상 청구 기록을 검증한다.
"""야전 정비창 수리부속 청구/출고 이력과 장비 가동률 데이터의 업무 프로세스 맵, 데이터 스키마(피처 20개 이상)와 데이터 딕셔너리를 정의하고 가상 청구 기록을 검증한다."""
import json
import re
import numpy as np
import pandas as pd

# 1) 업무 프로세스 맵: 단계별로 어떤 데이터가 생기는지 연결한다
PROCESS = [("1. 정비 소요 발생", "장비 고장·계획정비", "equip_id, equip_type, maint_type"),
           ("2. 수리부속 청구", "정비반 -> 정비창", "req_id, req_date, nsn, qty_req, priority"),
           ("3. 재고 확인·불출", "정비창 재고 담당", "qty_issued, emergency_flag"),
           ("4. 조달 요청", "재고 부족 시 조달부서", "lead_time_m, supplier_count, unit_price"),
           ("5. 단종 위험 검토", "품목 관리 담당", "prod_status, alt_available"),
           ("6. 정비 완료·가동률 집계", "월 단위 집계", "op_hours, availability, downtime_h")]
print("분석 목표: 월별 품목 수요를 예측해 결품·과잉재고를 줄이고, 단종예고·단종 품목의 최종 구매 수량 판단을 지원한다.")
print("[업무 프로세스 맵] 정비 소요 -> 청구 -> 불출 -> 조달 -> 단종 검토 -> 가동률")
for step, who, fields in PROCESS:
    print(f"  {step:18s} | {who:18s} | {fields}")

# 2) 스키마(데이터 딕셔너리): 그룹, 자료형, 규칙, 설명
NSN_RE = r"^\d{4}-\d{2}-\d{3}-\d{4}$"      # NSN = FSC(4) + 국가코드(2) + 일련번호(7), NIIN = 뒤 9자리
SCHEMA = {
    "청구/불출": {"req_id": (str, {"pattern": r"^R\d{6}$"}, "청구 번호"),
              "req_date": (str, {"pattern": r"^\d{4}-\d{2}-\d{2}$"}, "청구 일자"),
              "unit_code": (str, {"pattern": r"^U-\d{2}$"}, "청구 부대(가상 코드)"),
              "equip_id": (str, {"pattern": r"^[ABC]-\d{3}$"}, "장비 번호(가상)"),
              "equip_type": (str, {"enum": ["A", "B", "C"]}, "장비 유형"),
              "maint_type": (str, {"enum": ["계획", "비계획"]}, "정비 구분"),
              "nsn": (str, {"pattern": NSN_RE}, "국가재고번호(가상)"),
              "fsc": (str, {"pattern": r"^\d{4}$"}, "군수품 분류(NSN 앞 4자리)"),
              "niin": (str, {"pattern": r"^\d{9}$"}, "품목 식별번호(NSN 뒤 9자리)"),
              "qty_req": (int, {"min": 1, "max": 500}, "청구 수량"),
              "qty_issued": (int, {"min": 0, "max": 500}, "불출 수량"),
              "unit_of_issue": (str, {"enum": ["EA", "SE", "BX"]}, "불출 단위"),
              "priority": (str, {"enum": ["긴급", "일반"]}, "청구 우선순위"),
              "emergency_flag": (bool, {}, "긴급 불출 여부(이상치로 지우지 않음)")},
    "품목/공급망": {"unit_price": (int, {"min": 0}, "단가(원, 가상)"),
               "lead_time_m": (int, {"min": 0, "max": 36}, "조달 리드타임(월)"),
               "supplier_count": (int, {"min": 0, "max": 20}, "공급 가능 업체 수"),
               "prod_status": (str, {"enum": ["생산중", "단종예고", "단종"]}, "생산(단종) 상태"),
               "alt_available": (bool, {}, "대체품 존재 여부")},
    "가동률(월)": {"ym": (str, {"pattern": r"^\d{4}-\d{2}$"}, "집계 월"),
               "op_hours": (float, {"min": 0, "max": 744}, "장비 운용시간"),
               "availability": (float, {"min": 0, "max": 1}, "장비 가동률"),
               "downtime_h": (float, {"min": 0, "max": 744}, "정비 대기·수리 시간"),
               "exercise_flag": (int, {"enum": [0, 1]}, "대규모 훈련 월 여부")},
}


def check(name, v, typ, rule):
    if typ is float and isinstance(v, int) and not isinstance(v, bool):
        v = float(v)
    if not isinstance(v, typ) or (typ is not bool and isinstance(v, bool)):
        return [f"{name}: 자료형 오류({type(v).__name__})"]
    errs = []
    if "enum" in rule and v not in rule["enum"]:
        errs.append(f"{name}: 허용값 아님 '{v}'")
    if "min" in rule and v < rule["min"]:
        errs.append(f"{name}: 최소값 {rule['min']} 미만")
    if "max" in rule and v > rule["max"]:
        errs.append(f"{name}: 최대값 {rule['max']} 초과")
    if "pattern" in rule and not re.match(rule["pattern"], v):
        errs.append(f"{name}: 형식 불일치 '{v}'")
    return errs


def validate(rec):
    errs = []
    for grp in SCHEMA.values():
        for k, (typ, rule, _) in grp.items():
            errs += [f"{k}: 필수 필드 누락"] if k not in rec else check(k, rec[k], typ, rule)
    if not errs:          # 필드 간 논리 규칙
        if rec["nsn"][:4] != rec["fsc"] or rec["nsn"].replace("-", "")[4:] != rec["niin"]:
            errs.append("nsn/fsc/niin 불일치")
        if rec["qty_issued"] > rec["qty_req"]:
            errs.append("불출 수량이 청구 수량 초과")
    return errs


# 3) [교육용가상자료] 청구 기록 20건 생성
rng = np.random.default_rng(1)
records = []
for i in range(20):
    fsc = str(rng.choice(["2520", "2530", "2940", "6140"]))
    niin = f"00{rng.integers(1000000, 9999999)}"
    q = int(rng.integers(1, 20))
    et = str(rng.choice(["A", "B", "C"]))
    records.append({"req_id": f"R{250001 + i}", "req_date": f"2025-{rng.integers(1, 13):02d}-{rng.integers(1, 29):02d}",
                    "unit_code": f"U-{rng.integers(1, 6):02d}", "equip_id": f"{et}-{rng.integers(1, 200):03d}", "equip_type": et,
                    "maint_type": str(rng.choice(["계획", "비계획"])), "nsn": f"{fsc}-{niin[:2]}-{niin[2:5]}-{niin[5:]}",
                    "fsc": fsc, "niin": niin, "qty_req": q, "qty_issued": int(rng.integers(0, q + 1)),
                    "unit_of_issue": "EA", "priority": str(rng.choice(["긴급", "일반"], p=[0.1, 0.9])),
                    "emergency_flag": bool(rng.random() < 0.05), "unit_price": int(rng.integers(10, 900)) * 1000,
                    "lead_time_m": int(rng.integers(1, 13)), "supplier_count": int(rng.integers(1, 5)),
                    "prod_status": str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1])),
                    "alt_available": bool(rng.random() < 0.5), "ym": "2025-12", "op_hours": float(round(rng.uniform(80, 180), 1)),
                    "availability": float(round(rng.uniform(0.75, 0.97), 3)), "downtime_h": float(round(rng.uniform(5, 60), 1)),
                    "exercise_flag": int(rng.integers(0, 2))})

ok = sum(1 for r in records if not validate(r))
n_feat = sum(len(g) for g in SCHEMA.values())
print(f"\n정의한 피처(필드) {n_feat}개 (목표 20개 이상: {'달성' if n_feat >= 20 else '미달'})")
print(f"정상 샘플 {len(records)}건 중 검증 통과 {ok}건, 오류 {len(records) - ok}건")

# 4) 오류 주입 시험
bad = json.loads(json.dumps(records[0], ensure_ascii=False))
bad["nsn"] = "2520-001234567"            # NSN 형식 오류
bad["qty_req"] = -3                      # 음수 수량
bad["prod_status"] = "미상"              # 허용값 아님
del bad["lead_time_m"]                   # 필수 필드 누락
bad["availability"] = 1.4                # 가동률 범위 초과
print("[오류 주입 결과]")
for e in validate(bad):
    print(" -", e)
bad2 = json.loads(json.dumps(records[1], ensure_ascii=False)); bad2["fsc"] = "9999"
print(" - (논리 규칙)", validate(bad2))

# 5) 데이터 딕셔너리 출력
rows = [(g, k, t.__name__, json.dumps(r, ensure_ascii=False) if r else "-", d) for g, grp in SCHEMA.items() for k, (t, r, d) in grp.items()]
dic = pd.DataFrame(rows, columns=["그룹", "필드", "자료형", "규칙", "설명"])
pd.set_option("display.width", 200); pd.set_option("display.max_colwidth", 34)
print("[데이터 딕셔너리]")
print(dic.to_string(index=False))

실행 결과 예시

분석 목표: 월별 품목 수요를 예측해 결품·과잉재고를 줄이고, 단종예고·단종 품목의 최종 구매 수량 판단을 지원한다.
[업무 프로세스 맵] 정비 소요 -> 청구 -> 불출 -> 조달 -> 단종 검토 -> 가동률
  1. 정비 소요 발생        | 장비 고장·계획정비         | equip_id, equip_type, maint_type
  2. 수리부속 청구         | 정비반 -> 정비창         | req_id, req_date, nsn, qty_req, priority
  3. 재고 확인·불출        | 정비창 재고 담당          | qty_issued, emergency_flag
  4. 조달 요청           | 재고 부족 시 조달부서       | lead_time_m, supplier_count, unit_price
  5. 단종 위험 검토        | 품목 관리 담당           | prod_status, alt_available
  6. 정비 완료·가동률 집계    | 월 단위 집계            | op_hours, availability, downtime_h

정의한 피처(필드) 24개 (목표 20개 이상: 달성)
정상 샘플 20건 중 검증 통과 20건, 오류 0건
[오류 주입 결과]
 - nsn: 형식 불일치 '2520-001234567'
 - qty_req: 최소값 1 미만
 - lead_time_m: 필수 필드 누락
 - prod_status: 허용값 아님 '미상'
 - availability: 최대값 1 초과
 - (논리 규칙) ['nsn/fsc/niin 불일치']
[데이터 딕셔너리]
    그룹             필드   자료형                                           규칙                    설명
 청구/불출         req_id   str                     {"pattern": "^R\\d{6}$"}                 청구 번호
 청구/불출       req_date   str        {"pattern": "^\\d{4}-\\d{2}-\\d{2}$"}                 청구 일자
 청구/불출      unit_code   str                    {"pattern": "^U-\\d{2}$"}          청구 부대(가상 코드)
 청구/불출       equip_id   str                {"pattern": "^[ABC]-\\d{3}$"}             장비 번호(가상)
 청구/불출     equip_type   str                    {"enum": ["A", "B", "C"]}                 장비 유형
 청구/불출     maint_type   str                      {"enum": ["계획", "비계획"]}                 정비 구분
 청구/불출            nsn   str {"pattern": "^\\d{4}-\\d{2}-\\d{3}-\\d{4}$"}            국가재고번호(가상)
 청구/불출            fsc   str                      {"pattern": "^\\d{4}$"}     군수품 분류(NSN 앞 4자리)
 청구/불출           niin   str                      {"pattern": "^\\d{9}$"}    품목 식별번호(NSN 뒤 9자리)
 청구/불출        qty_req   int                       {"min": 1, "max": 500}                 청구 수량
 청구/불출     qty_issued   int                       {"min": 0, "max": 500}                 불출 수량
 청구/불출  unit_of_issue   str                 {"enum": ["EA", "SE", "BX"]}                 불출 단위
 청구/불출       priority   str                       {"enum": ["긴급", "일반"]}               청구 우선순위
 청구/불출 emergency_flag  bool                                            - 긴급 불출 여부(이상치로 지우지 않음)
품목/공급망     unit_price   int                                   {"min": 0}             단가(원, 가상)
품목/공급망    lead_time_m   int                        {"min": 0, "max": 36}            조달 리드타임(월)
품목/공급망 supplier_count   int                        {"min": 0, "max": 20}            공급 가능 업체 수
품목/공급망    prod_status   str              {"enum": ["생산중", "단종예고", "단종"]}             생산(단종) 상태
품목/공급망  alt_available  bool                                            -             대체품 존재 여부
가동률(월)             ym   str               {"pattern": "^\\d{4}-\\d{2}$"}                  집계 월
... (이하 생략)

미션 2. 결측치 보정, 계절성 수요 분리 및 수리부속 품목코드(NSN) 정제 전처리 모듈 구현

  • 난이도: 초급~중급 / 예상 소요시간: 2.0시간
  • 코드 설명: 수리부속 청구/출고 원시 기록에 정제 규칙 10종(NSN 형식 정규화, 중복 제거, 날짜 통일, 단위 환산, 결측 월 계절성 보간 등)을 적용하고, 월별 수요에서 계절성 성분을 분리한 뒤 정합성 검증 통과율을 확인한다.
"""수리부속 청구/출고 원시 기록에 정제 규칙 10종(NSN 형식 정규화, 중복 제거, 날짜 통일, 단위 환산, 결측 월 계절성 보간 등)을 적용하고, 월별 수요에서 계절성 성분을 분리한 뒤 정합성 검증 통과율을 확인한다."""
import re
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 야전 정비창 수리부속 청구/출고·장비 가동률 데이터 생성기
# 실제 부대·품목·재고·단가 정보와 무관하다. 품목 40개 x 48개월(가상 2022-01 ~ 2025-12) 월별 자료를 만든다.
# ---------------------------------------------------------------------------
FSC = {"2520": "power_train", "2530": "brake_steer", "2910": "fuel", "2930": "cooling",
       "2940": "filter", "5330": "gasket", "6140": "battery", "5961": "semiconductor"}
EXERCISE = {3, 8, 10}      # 가상 대규모 훈련 월(훈련 계획으로 미리 알 수 있음)
N_MONTHS = 48


def season(fsc, month):
    """품목군별 계절 계수: 축전지는 겨울, 냉각계통은 여름에 수요가 늘어난다."""
    if fsc == "6140":
        return 1.6 if month in (12, 1, 2) else 0.9
    if fsc == "2930":
        return 1.5 if month in (6, 7, 8) else 0.9
    return 1 + 0.1 * np.sin(2 * np.pi * month / 12)


def make_items(n=40, seed=2):
    """품목 마스터: NSN(가상), 품목군, 장비 유형, 기준 수요, 단가, 조달 리드타임, 공급업체 수, 생산(단종) 상태."""
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n):
        f = str(rng.choice(list(FSC)))
        status = str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1]))
        rows.append({"item_id": f"P{i + 1:03d}", "nsn": f"{f}-00-{rng.integers(100, 1000)}-{rng.integers(1000, 10000)}",
                     "fsc": f, "group": FSC[f], "equip_type": str(rng.choice(["A", "B", "C"])),
                     "base": float(np.round(rng.lognormal(2.7, 0.8), 1)),
                     "price": int(rng.lognormal(12, 1.0)) // 100 * 100,
                     "lead_time": int(rng.integers(1, 7) if status == "생산중" else rng.integers(4, 13)),
                     "suppliers": int(rng.integers(1, 5)) if status == "생산중" else 1,
                     "status": status, "alt_available": bool(rng.random() < (0.6 if status == "생산중" else 0.3))})
    return pd.DataFrame(rows)


def make_history(items, months=N_MONTHS, seed=3):
    """월별 청구 수요(개), 긴급 불출 건수, 장비 유형별 운용시간·가동률."""
    rng = np.random.default_rng(seed)
    t = np.arange(months)
    moy = t % 12 + 1
    ex = np.isin(moy, list(EXERCISE))
    hours = {e: (100 + 20 * rng.random(months)) * np.where(ex, 1.5, 1.0) for e in "ABC"}
    avail = {e: np.clip(0.93 - 0.12 * (hours[e] - 100) / 80 + rng.normal(0, 0.015, months), 0.6, 0.99) for e in "ABC"}
    rows = []
    for it in items.itertuples():
        trend = 1 + rng.uniform(-0.1, 0.3) * t / months
        s = np.array([season(it.fsc, m) for m in moy])
        lam = it.base * s * trend * hours[it.equip_type] / 110
        d = rng.poisson(lam)
        rows.append(pd.DataFrame({"item_id": it.item_id, "t": t, "ym": [f"{2022 + k // 12}-{k % 12 + 1:02d}" for k in t],
                                  "month": moy, "exercise": ex.astype(int), "op_hours": hours[it.equip_type].round(1),
                                  "availability": avail[it.equip_type].round(3), "demand": d,
                                  "emergency": rng.binomial(d, 0.05)}))
    return pd.concat(rows, ignore_index=True)


def make_stock(items, hist, seed=4):
    """현재(가상 2026-01) 재고: 보유량과 미입고 발주량(조달 중인 물량). 단종 품목은 추가 발주가 어려워 미입고량이 적다."""
    rng = np.random.default_rng(seed)
    recent = hist[hist["t"] >= hist["t"].max() - 5].groupby("item_id")["demand"].mean()
    st = items[["item_id", "lead_time", "status"]].copy()
    st["on_hand"] = [int(round(recent[i] * rng.uniform(0.5, 3.0))) for i in st["item_id"]]
    cover = np.where(st["status"] == "생산중", rng.uniform(0.6, 1.3, len(st)), rng.uniform(0.0, 0.6, len(st)))
    st["on_order"] = [int(round(recent[i] * L * c)) for i, L, c in zip(st["item_id"], st["lead_time"], cover)]
    return st[["item_id", "on_hand", "on_order"]]
# ---------------------------------------------------------------------------

items = make_items()
hist = make_history(items)
master = set(items["nsn"])
nsn_of = dict(zip(items["item_id"], items["nsn"]))

# 1) [교육용가상자료] 최근 24개월 월 수요를 개별 청구 기록으로 풀고 일부러 오류를 섞는다
rng = np.random.default_rng(5)
OUTAGE = {"2024-07", "2025-02"}            # 전산 장애로 기록이 올라오지 않은 월(가상)
rows, rid = [], 300000
for r in hist[hist["t"] >= 24].itertuples():
    if r.ym in OUTAGE or r.demand == 0:
        continue
    k = int(np.ceil(r.demand / 5))
    for j, q in enumerate(rng.multinomial(r.demand, np.ones(k) / k)):
        if q == 0:
            continue
        rid += 1
        rows.append({"req_id": f"R{rid}", "date": f"{r.ym}-{rng.integers(1, 29):02d}", "nsn": nsn_of[r.item_id],
                     "qty_req": int(q), "qty_issued": int(q), "uoi": "EA", "emergency": bool(j < r.emergency)})
raw = pd.DataFrame(rows)
n = len(raw)
idx = lambda frac: rng.choice(n, int(n * frac), replace=False)
raw.loc[idx(0.05), "nsn"] = raw["nsn"].str.replace("-", "")                       # 하이픈 누락
raw.loc[idx(0.03), "nsn"] = " " + raw["nsn"] + " "                                  # 앞뒤 공백
raw.loc[idx(0.04), "date"] = raw["date"].str.replace("-", ".")                      # 날짜 형식 혼재
raw.loc[idx(0.02), "qty_req"] = -raw["qty_req"]                                     # 음수(입력 오류)
raw.loc[idx(0.02), "qty_issued"] = raw["qty_req"].abs() + 2                         # 불출 > 청구
bx = idx(0.02); raw.loc[bx, "uoi"] = "BX"; raw.loc[bx, "qty_req"] = 1; raw.loc[bx, "qty_issued"] = 1  # 상자(10EA) 단위
raw.loc[idx(0.01), "qty_req"] = np.nan                                              # 수량 결측
raw.loc[idx(0.005), "date"] = "2027-01-15"                                          # 미래 일자
raw.loc[idx(0.005), "nsn"] = "9999-00-000-0000"                                     # 미등록 NSN
typo = [i for i in idx(0.01) if not raw.at[i, "emergency"]]
raw.loc[typo, "qty_req"] = raw.loc[typo, "qty_req"] * 10                            # 0 하나 더 입력한 오타
raw = pd.concat([raw, raw.sample(60, random_state=1)], ignore_index=True)          # 중복 전송
print(f"원시 청구 기록 {len(raw)}건 (가상), 전산 장애 월 {sorted(OUTAGE)}")


# 2) 정제 규칙 10종
def fmt_nsn(s):
    d = re.sub(r"\D", "", str(s))
    return f"{d[:4]}-{d[4:6]}-{d[6:9]}-{d[9:]}" if len(d) == 13 else str(s).strip()


log = {}
df = raw.copy()
before = len(df); df = df.drop_duplicates("req_id"); log["R1 중복 청구번호 제거"] = before - len(df)
em0 = int(df["emergency"].sum())
fixed = df["nsn"].map(fmt_nsn); log["R2 NSN 형식 정규화(공백·하이픈)"] = int((fixed != df["nsn"]).sum()); df["nsn"] = fixed
unk = ~df["nsn"].isin(master); log["R3 미등록 NSN 격리"] = int(unk.sum()); em_drop = int(df.loc[unk, "emergency"].sum()); df = df[~unk]
d2 = df["date"].str.replace(".", "-", regex=False); log["R4 날짜 형식 통일"] = int((d2 != df["date"]).sum()); df["date"] = pd.to_datetime(d2)
fut = df["date"] > "2025-12-31"; log["R5 미래 일자 제거"] = int(fut.sum()); em_drop += int(df.loc[fut, "emergency"].sum()); df = df[~fut]
neg = df["qty_req"] < 0; log["R6 음수 수량 부호 보정"] = int(neg.sum()); df.loc[neg, "qty_req"] = -df.loc[neg, "qty_req"]
na = df["qty_req"].isna(); log["R7 청구 수량 결측 -> 불출 수량으로 보완"] = int(na.sum()); df.loc[na, "qty_req"] = df.loc[na, "qty_issued"]
bxm = df["uoi"] == "BX"; log["R8 단위 환산(BX=10EA)"] = int(bxm.sum())
df.loc[bxm, ["qty_req", "qty_issued"]] = df.loc[bxm, ["qty_req", "qty_issued"]] * 10; df.loc[bxm, "uoi"] = "EA"
over = df["qty_issued"] > df["qty_req"]; log["R9 불출>청구 -> 청구 수량으로 제한"] = int(over.sum())
df.loc[over, "qty_issued"] = df.loc[over, "qty_req"]
cap = df.groupby("nsn")["qty_req"].transform(lambda s: s.median() * 5)
out = (df["qty_req"] > cap) & ~df["emergency"]
log["R10 오타성 이상치 보정(긴급 불출은 보존)"] = int(out.sum()); df.loc[out, "qty_req"] = (df.loc[out, "qty_req"] / 10).round()
df["qty_req"] = df["qty_req"].astype(int)

# 3) 월별 집계 + 결측 월 계절성 보간(R11): 작년 같은 달 x (최근 수준 / 작년 수준)
df["ym"] = df["date"].dt.strftime("%Y-%m")
monthly = df.pivot_table(index="ym", columns="nsn", values="qty_req", aggfunc="sum")
all_m = [f"{y}-{m:02d}" for y in (2024, 2025) for m in range(1, 13)]
monthly = monthly.reindex(all_m)
prev = hist[hist["t"].between(12, 23)].assign(nsn=lambda x: x["item_id"].map(nsn_of)).pivot(index="month", columns="nsn", values="demand")
seasonal_fill = monthly.copy()
for ym in OUTAGE:
    i = all_m.index(ym)
    ratio = monthly.iloc[max(0, i - 3):i].mean() / prev.loc[[(i - k) % 12 + 1 for k in (1, 2, 3)]].mean().replace(0, np.nan)
    seasonal_fill.loc[ym] = (prev.loc[int(ym[-2:])] * ratio.fillna(1)).round()
linear_fill = monthly.interpolate()
monthly = monthly.fillna(0)
seasonal_fill = seasonal_fill.fillna(0)
log["R11 결측 월 계절성 보간"] = len(OUTAGE) * monthly.shape[1]
truth = hist[hist["t"] >= 24].assign(nsn=lambda x: x["item_id"].map(nsn_of)).pivot(index="ym", columns="nsn", values="demand")
err_s = np.abs(seasonal_fill.loc[sorted(OUTAGE)] - truth.loc[sorted(OUTAGE)]).values.mean()
err_l = np.abs(linear_fill.loc[sorted(OUTAGE)] - truth.loc[sorted(OUTAGE)]).values.mean()

print("\n[정제 규칙 적용 결과]")
for k, v in log.items():
    print(f"  {k:32s} {v:5d}건")
print(f"결측 월 보간 오차(MAE): 계절성 보간 {err_s:.2f}개 / 단순 선형 보간 {err_l:.2f}개")

# 4) 정합성 검증
checks = {"청구번호 유일": df["req_id"].is_unique, "NSN 형식": df["nsn"].str.match(r"^\d{4}-\d{2}-\d{3}-\d{4}$").all(),
          "마스터 등록 NSN": df["nsn"].isin(master).all(), "수량 양수": (df["qty_req"] > 0).all(),
          "불출<=청구": (df["qty_issued"] <= df["qty_req"]).all(), "단위 EA 통일": (df["uoi"] == "EA").all(),
          "미래 일자 없음": (df["date"] <= "2025-12-31").all(), "결측 월 없음": not seasonal_fill.isna().any().any(),
          "긴급 불출 보존": int(df["emergency"].sum()) == em0 - em_drop,
          "월 합계 = 기록 합계": int(monthly.values.sum()) == int(df["qty_req"].sum())}
rate = sum(checks.values()) / len(checks)
print(f"\n정합성 검증 {sum(checks.values())}/{len(checks)} 통과({rate:.0%}) - 목표 100%: {'달성' if rate == 1 else '미달'}, 적용 규칙 {len(log)}종(목표 10종 이상)")

# 5) 전후 기술통계 비교
stat = pd.DataFrame({"정제 전": pd.to_numeric(raw["qty_req"], errors="coerce").describe(), "정제 후": df["qty_req"].describe()}).round(2)
print("\n[청구 수량 기술통계 전후 비교]")
print(stat.T.to_string())

# 6) 계절성 분리(고전적 분해): 추세 = 12개월 중심 이동평균, 계절지수 = 월별 (실적/추세) 평균
bat = items[items["fsc"] == "6140"].sort_values("base")["nsn"].iloc[-1]      # 수요가 가장 많은 축전지 품목
series = pd.concat([truth.iloc[:0], seasonal_fill])[bat].astype(float)
full = hist[hist["item_id"] == items.set_index("nsn").loc[bat, "item_id"]].set_index("ym")["demand"].astype(float)
full.loc[series.index] = series
trend = full.rolling(12, center=True).mean()
sidx = (full / trend).groupby([int(i[-2:]) for i in full.index]).mean()
sidx = sidx / sidx.mean()
print(f"\n축전지 품목 {bat} 계절지수(1=평균):", " ".join(f"{m}월 {v:.2f}" for m, v in sidx.items()))
print("해석: 겨울(12~2월) 축전지 수요 증가와 훈련 월(3·8·10월) 증가가 함께 나타난다. 훈련 효과는 미션 4에서 별도 피처로 분리한다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(full.values, label="demand"); ax[0].plot(trend.values, label="trend (12M centered MA)")
ax[0].set_title(f"Battery item {bat}: demand vs trend"); ax[0].set_xlabel("month index"); ax[0].legend()
ax[1].bar(sidx.index, sidx.values); ax[1].axhline(1, c="k", ls="--"); ax[1].set_title("Seasonal index by month"); ax[1].set_xlabel("month")
plt.tight_layout(); plt.show()

실행 결과 예시

원시 청구 기록 4918건 (가상), 전산 장애 월 ['2024-07', '2025-02']

[정제 규칙 적용 결과]
  R1 중복 청구번호 제거                       60건
  R2 NSN 형식 정규화(공백·하이픈)              372건
  R3 미등록 NSN 격리                       24건
  R4 날짜 형식 통일                        191건
  R5 미래 일자 제거                         24건
  R6 음수 수량 부호 보정                      94건
  R7 청구 수량 결측 -> 불출 수량으로 보완           47건
  R8 단위 환산(BX=10EA)                   97건
  R9 불출>청구 -> 청구 수량으로 제한              96건
  R10 오타성 이상치 보정(긴급 불출은 보존)           35건
  R11 결측 월 계절성 보간                     80건
결측 월 보간 오차(MAE): 계절성 보간 5.72개 / 단순 선형 보간 7.86개

정합성 검증 10/10 통과(100%) - 목표 100%: 달성, 적용 규칙 11종(목표 10종 이상)

[청구 수량 기술통계 전후 비교]
       count  mean   std   min  25%  50%  75%   max
정제 전  4870.0  4.75  4.76 -12.0  3.0  4.0  6.0  80.0
정제 후  4810.0  4.80  2.17   1.0  3.0  5.0  6.0  20.0

축전지 품목 6140-00-643-1860 계절지수(1=평균): 1월 1.39 2월 1.16 3월 1.13 4월 0.69 5월 0.69 6월 0.66 7월 0.75 8월 1.16 9월 0.85 10월 1.38 11월 0.68 12월 1.46
해석: 겨울(12~2월) 축전지 수요 증가와 훈련 월(3·8·10월) 증가가 함께 나타난다. 훈련 효과는 미션 4에서 별도 피처로 분리한다.

미션 3. 이동평균 및 지수평활법 기반 1차 부품 소요량 예측 Baseline 알고리즘 구현

  • 난이도: 중급 / 예상 소요시간: 2.5시간
  • 코드 설명: 이동평균(3·6개월)과 단순 지수평활법으로 품목별 다음 달 부품 소요량을 예측하는 Baseline을 구현하고, 탐색적 분석(EDA) 차트 8종·상관 검정·ABC(파레토) 분류로 수요 특성과 핵심 관리 품목을 도출한다.
"""이동평균(3·6개월)과 단순 지수평활법으로 품목별 다음 달 부품 소요량을 예측하는 Baseline을 구현하고, 탐색적 분석(EDA) 차트 8종·상관 검정·ABC(파레토) 분류로 수요 특성과 핵심 관리 품목을 도출한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 야전 정비창 수리부속 청구/출고·장비 가동률 데이터 생성기
# 실제 부대·품목·재고·단가 정보와 무관하다. 품목 40개 x 48개월(가상 2022-01 ~ 2025-12) 월별 자료를 만든다.
# ---------------------------------------------------------------------------
FSC = {"2520": "power_train", "2530": "brake_steer", "2910": "fuel", "2930": "cooling",
       "2940": "filter", "5330": "gasket", "6140": "battery", "5961": "semiconductor"}
EXERCISE = {3, 8, 10}      # 가상 대규모 훈련 월(훈련 계획으로 미리 알 수 있음)
N_MONTHS = 48


def season(fsc, month):
    """품목군별 계절 계수: 축전지는 겨울, 냉각계통은 여름에 수요가 늘어난다."""
    if fsc == "6140":
        return 1.6 if month in (12, 1, 2) else 0.9
    if fsc == "2930":
        return 1.5 if month in (6, 7, 8) else 0.9
    return 1 + 0.1 * np.sin(2 * np.pi * month / 12)


def make_items(n=40, seed=2):
    """품목 마스터: NSN(가상), 품목군, 장비 유형, 기준 수요, 단가, 조달 리드타임, 공급업체 수, 생산(단종) 상태."""
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n):
        f = str(rng.choice(list(FSC)))
        status = str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1]))
        rows.append({"item_id": f"P{i + 1:03d}", "nsn": f"{f}-00-{rng.integers(100, 1000)}-{rng.integers(1000, 10000)}",
                     "fsc": f, "group": FSC[f], "equip_type": str(rng.choice(["A", "B", "C"])),
                     "base": float(np.round(rng.lognormal(2.7, 0.8), 1)),
                     "price": int(rng.lognormal(12, 1.0)) // 100 * 100,
                     "lead_time": int(rng.integers(1, 7) if status == "생산중" else rng.integers(4, 13)),
                     "suppliers": int(rng.integers(1, 5)) if status == "생산중" else 1,
                     "status": status, "alt_available": bool(rng.random() < (0.6 if status == "생산중" else 0.3))})
    return pd.DataFrame(rows)


def make_history(items, months=N_MONTHS, seed=3):
    """월별 청구 수요(개), 긴급 불출 건수, 장비 유형별 운용시간·가동률."""
    rng = np.random.default_rng(seed)
    t = np.arange(months)
    moy = t % 12 + 1
    ex = np.isin(moy, list(EXERCISE))
    hours = {e: (100 + 20 * rng.random(months)) * np.where(ex, 1.5, 1.0) for e in "ABC"}
    avail = {e: np.clip(0.93 - 0.12 * (hours[e] - 100) / 80 + rng.normal(0, 0.015, months), 0.6, 0.99) for e in "ABC"}
    rows = []
    for it in items.itertuples():
        trend = 1 + rng.uniform(-0.1, 0.3) * t / months
        s = np.array([season(it.fsc, m) for m in moy])
        lam = it.base * s * trend * hours[it.equip_type] / 110
        d = rng.poisson(lam)
        rows.append(pd.DataFrame({"item_id": it.item_id, "t": t, "ym": [f"{2022 + k // 12}-{k % 12 + 1:02d}" for k in t],
                                  "month": moy, "exercise": ex.astype(int), "op_hours": hours[it.equip_type].round(1),
                                  "availability": avail[it.equip_type].round(3), "demand": d,
                                  "emergency": rng.binomial(d, 0.05)}))
    return pd.concat(rows, ignore_index=True)


def make_stock(items, hist, seed=4):
    """현재(가상 2026-01) 재고: 보유량과 미입고 발주량(조달 중인 물량). 단종 품목은 추가 발주가 어려워 미입고량이 적다."""
    rng = np.random.default_rng(seed)
    recent = hist[hist["t"] >= hist["t"].max() - 5].groupby("item_id")["demand"].mean()
    st = items[["item_id", "lead_time", "status"]].copy()
    st["on_hand"] = [int(round(recent[i] * rng.uniform(0.5, 3.0))) for i in st["item_id"]]
    cover = np.where(st["status"] == "생산중", rng.uniform(0.6, 1.3, len(st)), rng.uniform(0.0, 0.6, len(st)))
    st["on_order"] = [int(round(recent[i] * L * c)) for i, L, c in zip(st["item_id"], st["lead_time"], cover)]
    return st[["item_id", "on_hand", "on_order"]]
# ---------------------------------------------------------------------------

items = make_items()
hist = make_history(items).merge(items[["item_id", "group", "price", "status", "lead_time"]], on="item_id")
TEST_FROM = 42                          # 마지막 6개월(가상 2025-07 ~ 2025-12)을 평가 구간으로 둔다


def ses(y, alpha):
    """단순 지수평활: 다음 달 예측 = alpha x 이번 달 실적 + (1 - alpha) x 이번 달 예측."""
    f = np.empty(len(y)); f[0] = y[:6].mean()
    for i in range(1, len(y)):
        f[i] = alpha * y[i - 1] + (1 - alpha) * f[i - 1]
    return f


# 1) Baseline 예측(한 달 앞): MA3, MA6, SES(품목별 alpha를 학습 구간에서 선택), 작년 같은 달
preds = []
for iid, g in hist.groupby("item_id"):
    y = g.sort_values("t")["demand"].to_numpy(float)
    best = min([0.1, 0.2, 0.3, 0.5], key=lambda a: np.abs(ses(y, a)[12:TEST_FROM] - y[12:TEST_FROM]).mean())
    s = pd.Series(y)
    preds.append(pd.DataFrame({"item_id": iid, "t": np.arange(len(y)), "y": y, "MA3": s.shift(1).rolling(3).mean(),
                               "MA6": s.shift(1).rolling(6).mean(), "SES": ses(y, best), "작년 같은 달": s.shift(12)}))
P = pd.concat(preds)
T = P[P["t"] >= TEST_FROM]
ab = items[items["base"] >= 10]["item_id"]                       # 수요가 충분한 품목(MAPE 해석 가능)
rows = []
for m in ["MA3", "MA6", "SES", "작년 같은 달"]:
    e = T["y"] - T[m]
    big = T[T["item_id"].isin(ab)]
    rows.append({"방법": m, "MAE(개)": round(e.abs().mean(), 2), "WAPE": f"{e.abs().sum() / T['y'].sum():.1%}",
                 "MAPE(수요 많은 품목)": f"{(np.abs(big['y'] - big[m]) / big['y'].clip(lower=1)).mean():.1%}"})
base_tbl = pd.DataFrame(rows)
print("[Baseline 예측 성능 - 평가 구간 6개월]")
print(base_tbl.to_string(index=False))
print("해석: 이동평균은 훈련 월 급증과 계절 변화를 한 달 늦게 따라가 오차가 커진다(미션 4에서 개선).")

# 2) 상관 검정: 품목 평균 대비 수요 비율과 후보 피처
hist = hist.sort_values(["item_id", "t"])
hist["ratio"] = hist["demand"] / hist.groupby("item_id")["demand"].transform("mean")
g = hist.groupby("item_id")["ratio"]
hist["lag1_ratio"] = g.shift(1); hist["lag12_ratio"] = g.shift(12)
hist["winter"] = hist["month"].isin([12, 1, 2]).astype(int)
cand = ["op_hours", "exercise", "availability", "lag12_ratio", "lag1_ratio", "t", "winter", "price"]
res = []
for c in cand:
    d = hist[["ratio", c]].dropna()
    r, p = stats.pearsonr(d[c], d["ratio"])
    res.append({"피처": c, "상관계수": round(r, 3), "p-value": f"{p:.1e}", "유의(p<0.05)": "예" if p < 0.05 else "아니오"})
res = pd.DataFrame(res)
n_sig = (res["유의(p<0.05)"] == "예").sum()
print("\n[수요 비율과 피처의 상관 검정]")
print(res.to_string(index=False))
print(f"유의한 피처 {n_sig}개 (목표 5개 이상: {'달성' if n_sig >= 5 else '미달'})")
t_stat, p_ex = stats.ttest_ind(hist[hist["exercise"] == 1]["ratio"], hist[hist["exercise"] == 0]["ratio"])
print(f"훈련 월과 평시 월 수요 비율 차이 t검정: t={t_stat:.1f}, p={p_ex:.1e}")

# 3) ABC(파레토) 분류: 최근 12개월 수요 금액 기준
val = hist[hist["t"] >= 36].groupby("item_id").apply(lambda x: (x["demand"] * x["price"]).sum(), include_groups=False).sort_values(ascending=False)
cum = val.cumsum() / val.sum()
abc = pd.Series(np.where(cum <= 0.8, "A", np.where(cum <= 0.95, "B", "C")), index=val.index)
top = val.index[: int(np.ceil(len(val) * 0.2))]
print(f"\nABC 분류: A {sum(abc == 'A')}개 / B {sum(abc == 'B')}개 / C {sum(abc == 'C')}개, 상위 20% {len(top)}개 품목이 수요 금액의 {val[top].sum() / val.sum():.0%} 차지")
st = items.set_index("item_id").loc[top]
print("상위 20% 중 단종예고·단종 품목:", ", ".join(f"{i}({s})" for i, s in st["status"].items() if s != "생산중") or "없음")

# 4) EDA 차트 8종
fig, ax = plt.subplots(2, 4, figsize=(18, 8))
tot = hist.groupby("t")["demand"].sum()
ax[0, 0].plot(tot.index, tot.values); ax[0, 0].set_title("1. Total monthly demand")
hist.groupby(["t", "group"])["demand"].sum().unstack().plot(ax=ax[0, 1], legend=False); ax[0, 1].set_title("2. Demand by FSC group")
hist.boxplot(column="ratio", by="month", ax=ax[0, 2], grid=False); ax[0, 2].set_title("3. Seasonality (demand ratio)")
ax[0, 3].scatter(hist["op_hours"], hist["ratio"], s=4, alpha=0.3); ax[0, 3].set_title("4. Operating hours vs demand ratio")
hist.groupby("exercise")["ratio"].mean().plot.bar(ax=ax[1, 0]); ax[1, 0].set_title("5. Exercise month effect")
ax[1, 1].plot(np.arange(1, len(cum) + 1), cum.values, "o-", ms=3); ax[1, 1].axhline(0.8, ls="--", c="gray"); ax[1, 1].set_title("6. Pareto curve (ABC)")
items.assign(st=items["status"].map({"생산중": "active", "단종예고": "EOL notice", "단종": "discontinued"})).boxplot(column="lead_time", by="st", ax=ax[1, 2], grid=False); ax[1, 2].set_title("7. Lead time by production status")
ag = hist.groupby("t").agg(av=("availability", "mean"), d=("demand", "sum"))
ax[1, 3].scatter(ag["av"], ag["d"]); ax[1, 3].set_title("8. Availability vs total demand")
for a in ax.ravel():
    a.set_xlabel(""); a.tick_params(labelsize=7)
plt.suptitle("MNT02 EDA - 8 charts (synthetic data)"); plt.tight_layout(); plt.show()
print("EDA 차트 8종 생성 완료. 핵심 발견: 운용시간·훈련 월이 수요를 끌어올리고, 단종 품목은 리드타임이 길어 결품 위험이 크다.")

실행 결과 예시

[Baseline 예측 성능 - 평가 구간 6개월]
     방법  MAE(개)  WAPE MAPE(수요 많은 품목)
    MA3    7.34 28.7%          32.7%
    MA6    7.09 27.7%          31.6%
    SES    6.53 25.5%          28.9%
작년 같은 달    5.38 21.0%          22.8%
해석: 이동평균은 훈련 월 급증과 계절 변화를 한 달 늦게 따라가 오차가 커진다(미션 4에서 개선).

[수요 비율과 피처의 상관 검정]
          피처   상관계수  p-value 유의(p<0.05)
    op_hours  0.477 8.3e-110          예
    exercise  0.470 4.4e-106          예
availability -0.450  1.6e-96          예
 lag12_ratio  0.339  3.7e-40          예
  lag1_ratio -0.111  1.5e-06          예
           t  0.066  3.8e-03          예
      winter -0.032  1.6e-01        아니오
       price  0.000  1.0e+00        아니오
유의한 피처 6개 (목표 5개 이상: 달성)
훈련 월과 평시 월 수요 비율 차이 t검정: t=23.3, p=4.4e-106

ABC 분류: A 13개 / B 13개 / C 14개, 상위 20% 8개 품목이 수요 금액의 67% 차지
상위 20% 중 단종예고·단종 품목: P016(단종), P028(단종), P032(단종예고), P007(단종예고)
EDA 차트 8종 생성 완료. 핵심 발견: 운용시간·훈련 월이 수요를 끌어올리고, 단종 품목은 리드타임이 길어 결품 위험이 크다.

미션 4. 시계열 머신러닝(Prophet/LSTM) 기반 결전 대비 야전 부품 수요 예측 모델 고도화

  • 난이도: 중급 / 예상 소요시간: 3.0시간
  • 코드 설명: 과거 수요 지연값·이동평균·계절지수·훈련 계획·운용시간 등 파생 피처 15개로 시계열 머신러닝 수요 예측 모델(Prophet/LSTM 대신 포아송 그래디언트 부스팅)을 만들고, 시계열 교차검증으로 MAPE를 확인한 뒤 작전 템포 증가(결전 대비) 시 수요를 예측한다.
"""과거 수요 지연값·이동평균·계절지수·훈련 계획·운용시간 등 파생 피처 15개로 시계열 머신러닝 수요 예측 모델(Prophet/LSTM 대신 포아송 그래디언트 부스팅)을 만들고, 시계열 교차검증으로 MAPE를 확인한 뒤 작전 템포 증가(결전 대비) 시 수요를 예측한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.inspection import permutation_importance

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 야전 정비창 수리부속 청구/출고·장비 가동률 데이터 생성기
# 실제 부대·품목·재고·단가 정보와 무관하다. 품목 40개 x 48개월(가상 2022-01 ~ 2025-12) 월별 자료를 만든다.
# ---------------------------------------------------------------------------
FSC = {"2520": "power_train", "2530": "brake_steer", "2910": "fuel", "2930": "cooling",
       "2940": "filter", "5330": "gasket", "6140": "battery", "5961": "semiconductor"}
EXERCISE = {3, 8, 10}      # 가상 대규모 훈련 월(훈련 계획으로 미리 알 수 있음)
N_MONTHS = 48


def season(fsc, month):
    """품목군별 계절 계수: 축전지는 겨울, 냉각계통은 여름에 수요가 늘어난다."""
    if fsc == "6140":
        return 1.6 if month in (12, 1, 2) else 0.9
    if fsc == "2930":
        return 1.5 if month in (6, 7, 8) else 0.9
    return 1 + 0.1 * np.sin(2 * np.pi * month / 12)


def make_items(n=40, seed=2):
    """품목 마스터: NSN(가상), 품목군, 장비 유형, 기준 수요, 단가, 조달 리드타임, 공급업체 수, 생산(단종) 상태."""
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n):
        f = str(rng.choice(list(FSC)))
        status = str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1]))
        rows.append({"item_id": f"P{i + 1:03d}", "nsn": f"{f}-00-{rng.integers(100, 1000)}-{rng.integers(1000, 10000)}",
                     "fsc": f, "group": FSC[f], "equip_type": str(rng.choice(["A", "B", "C"])),
                     "base": float(np.round(rng.lognormal(2.7, 0.8), 1)),
                     "price": int(rng.lognormal(12, 1.0)) // 100 * 100,
                     "lead_time": int(rng.integers(1, 7) if status == "생산중" else rng.integers(4, 13)),
                     "suppliers": int(rng.integers(1, 5)) if status == "생산중" else 1,
                     "status": status, "alt_available": bool(rng.random() < (0.6 if status == "생산중" else 0.3))})
    return pd.DataFrame(rows)


def make_history(items, months=N_MONTHS, seed=3):
    """월별 청구 수요(개), 긴급 불출 건수, 장비 유형별 운용시간·가동률."""
    rng = np.random.default_rng(seed)
    t = np.arange(months)
    moy = t % 12 + 1
    ex = np.isin(moy, list(EXERCISE))
    hours = {e: (100 + 20 * rng.random(months)) * np.where(ex, 1.5, 1.0) for e in "ABC"}
    avail = {e: np.clip(0.93 - 0.12 * (hours[e] - 100) / 80 + rng.normal(0, 0.015, months), 0.6, 0.99) for e in "ABC"}
    rows = []
    for it in items.itertuples():
        trend = 1 + rng.uniform(-0.1, 0.3) * t / months
        s = np.array([season(it.fsc, m) for m in moy])
        lam = it.base * s * trend * hours[it.equip_type] / 110
        d = rng.poisson(lam)
        rows.append(pd.DataFrame({"item_id": it.item_id, "t": t, "ym": [f"{2022 + k // 12}-{k % 12 + 1:02d}" for k in t],
                                  "month": moy, "exercise": ex.astype(int), "op_hours": hours[it.equip_type].round(1),
                                  "availability": avail[it.equip_type].round(3), "demand": d,
                                  "emergency": rng.binomial(d, 0.05)}))
    return pd.concat(rows, ignore_index=True)


def make_stock(items, hist, seed=4):
    """현재(가상 2026-01) 재고: 보유량과 미입고 발주량(조달 중인 물량). 단종 품목은 추가 발주가 어려워 미입고량이 적다."""
    rng = np.random.default_rng(seed)
    recent = hist[hist["t"] >= hist["t"].max() - 5].groupby("item_id")["demand"].mean()
    st = items[["item_id", "lead_time", "status"]].copy()
    st["on_hand"] = [int(round(recent[i] * rng.uniform(0.5, 3.0))) for i in st["item_id"]]
    cover = np.where(st["status"] == "생산중", rng.uniform(0.6, 1.3, len(st)), rng.uniform(0.0, 0.6, len(st)))
    st["on_order"] = [int(round(recent[i] * L * c)) for i, L, c in zip(st["item_id"], st["lead_time"], cover)]
    return st[["item_id", "on_hand", "on_order"]]
# ---------------------------------------------------------------------------

# ---- 수요 예측 공통 함수(파생 피처 + 머신러닝 모델) ----
FEATS = ["lag1", "lag2", "lag3", "lag12", "ma3", "ma6", "std6", "season_idx", "month_sin", "month_cos",
         "exercise", "op_hours", "op_ratio", "item_level", "trend_t"]


def make_features(hist):
    """품목별 과거 수요로 파생 피처를 만든다(모두 예측 시점 이전 정보만 사용)."""
    df = hist.sort_values(["item_id", "t"]).reset_index(drop=True)
    g = df.groupby("item_id")["demand"]
    for k in (1, 2, 3, 12):
        df[f"lag{k}"] = g.shift(k)
    df["ma3"] = g.transform(lambda s: s.shift(1).rolling(3).mean())
    df["ma6"] = g.transform(lambda s: s.shift(1).rolling(6).mean())
    df["std6"] = g.transform(lambda s: s.shift(1).rolling(6).std())
    df["item_level"] = g.transform(lambda s: s.shift(1).rolling(12).mean())
    df["season_idx"] = df["lag12"] / (g.transform(lambda s: s.shift(12).rolling(12).mean()) + 1)
    df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
    df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)
    df["op_ratio"] = df["op_hours"] / df.groupby("item_id")["op_hours"].transform(lambda s: s.shift(1).rolling(12).mean())
    df["trend_t"] = df["t"]
    return df.dropna().reset_index(drop=True)


def new_model():
    """Prophet/LSTM 대신 사용하는 포아송 손실 그래디언트 부스팅(개수 자료에 적합)."""
    return HistGradientBoostingRegressor(loss="poisson", max_iter=300, learning_rate=0.05, max_leaf_nodes=15,
                                         min_samples_leaf=20, random_state=0)


def mape(y, p, min_y=1):
    y, p = np.asarray(y, float), np.asarray(p, float)
    m = y >= min_y
    return float(np.mean(np.abs(y[m] - p[m]) / y[m]))

items = make_items()
hist = make_history(items)
df = make_features(hist)
ab = set(items[items["base"] >= 10]["item_id"])                  # 수요가 충분한 품목(MAPE 평가 대상)
print(f"파생 피처 {len(FEATS)}개 (목표 10개 이상): {', '.join(FEATS)}")

# 1) 시계열 교차검증(확장 창): 과거로 학습 -> 바로 다음 6개월 평가, 3회 반복
rows = []
for cut in (30, 36, 42):
    tr, te = df[df["t"] < cut], df[(df["t"] >= cut) & (df["t"] < cut + 6)]
    p = new_model().fit(tr[FEATS], tr["demand"]).predict(te[FEATS])
    big = te["item_id"].isin(ab).to_numpy()
    rows.append({"학습 구간": f"t<{cut}", "평가 구간": f"t={cut}~{cut + 5}",
                 "ML MAE": round(np.abs(te["demand"] - p).mean(), 2), "MA3 MAE": round(np.abs(te["demand"] - te["ma3"]).mean(), 2),
                 "ML MAPE": mape(te["demand"][big], p[big]), "MA3 MAPE": mape(te["demand"][big], te["ma3"][big])})
cv = pd.DataFrame(rows)
print("\n[시계열 교차검증 결과 - MAPE는 수요 많은 품목 기준]")
print(cv.assign(**{"ML MAPE": cv["ML MAPE"].map("{:.1%}".format), "MA3 MAPE": cv["MA3 MAPE"].map("{:.1%}".format)}).to_string(index=False))
m_mape = cv["ML MAPE"].mean()
print(f"평균 MAPE: ML {m_mape:.1%} / MA3 {cv['MA3 MAPE'].mean():.1%} -> 목표(MAPE 25% 이하): {'달성' if m_mape <= 0.25 else '미달'}")

# 2) 최종 모델(t<42 학습)과 피처 중요도(순열 중요도)
tr, te = df[df["t"] < 42], df[df["t"] >= 42]
model = new_model().fit(tr[FEATS], tr["demand"])
imp = permutation_importance(model, te[FEATS], te["demand"], n_repeats=5, random_state=0, scoring="neg_mean_absolute_error")
imp = pd.Series(imp.importances_mean, index=FEATS).sort_values(ascending=False)
print("\n[피처 중요도 상위 6개(MAE 증가량)]", ", ".join(f"{k} {v:.2f}" for k, v in imp.head(6).items()))
print("중요도가 0에 가까운 피처(차원 축소 후보):", ", ".join(imp[imp < 0.02].index) or "없음")

# 3) 결전 대비 시나리오: 다음 달(가상 2026-01) 운용시간 1.5배 + 훈련 월 지정
future = []
for it in items.itertuples():
    last = hist[(hist["item_id"] == it.item_id) & (hist["t"] == 47)].iloc[0]
    future.append({**last.to_dict(), "t": 48, "ym": "2026-01", "month": 1, "demand": 0, "exercise": 0})
fut = pd.DataFrame(future)
fut["op_hours"] = hist[hist["t"] == 36].set_index("item_id").loc[fut["item_id"], "op_hours"].to_numpy()  # 작년 1월 계획 운용시간
res = {}
for name, mult, ex in [("평시 계획", 1.0, 0), ("작전 템포 1.5배", 1.5, 1)]:
    f2 = fut.assign(op_hours=fut["op_hours"] * mult, exercise=ex)
    X = make_features(pd.concat([hist, f2], ignore_index=True))
    X = X[X["t"] == 48]
    res[name] = pd.Series(model.predict(X[FEATS]), index=X["item_id"])
comp = pd.DataFrame(res).round(1)
comp["증가율"] = (comp["작전 템포 1.5배"] / comp["평시 계획"] - 1).map("{:+.0%}".format)
comp = comp.join(items.set_index("item_id")[["group", "status"]])
print("\n[결전 대비 다음 달 수요 예측 - 증가량 상위 6개 품목]")
print(comp.loc[(comp["작전 템포 1.5배"] - comp["평시 계획"]).sort_values(ascending=False).index[:6]].to_string())
print(f"전체 예측 수요: 평시 {comp['평시 계획'].sum():.0f}개 -> 작전 템포 1.5배 {comp['작전 템포 1.5배'].sum():.0f}개")

# 4) 예측 그림: 수요 많은 품목 1개
iid = items.sort_values("base")["item_id"].iloc[-1]
s = te[te["item_id"] == iid]
full = df[df["item_id"] == iid]
plt.figure(figsize=(9, 4))
plt.plot(full["t"], full["demand"], "k.-", label="actual")
plt.plot(s["t"], model.predict(s[FEATS]), "r-o", label="ML forecast")
plt.plot(s["t"], s["ma3"], "b--", label="MA3 baseline")
plt.axvline(41.5, c="gray", ls=":"); plt.xlabel("month index"); plt.ylabel("demand (EA)")
plt.title(f"{iid} one-month-ahead forecast (synthetic)"); plt.legend(); plt.show()

실행 결과 예시

파생 피처 15개 (목표 10개 이상): lag1, lag2, lag3, lag12, ma3, ma6, std6, season_idx, month_sin, month_cos, exercise, op_hours, op_ratio, item_level, trend_t

[시계열 교차검증 결과 - MAPE는 수요 많은 품목 기준]
학습 구간   평가 구간  ML MAE  MA3 MAE ML MAPE MA3 MAPE
 t<30 t=30~35    5.85     7.26   25.8%    30.9%
 t<36 t=36~41    4.84     7.26   20.5%    30.9%
 t<42 t=42~47    4.57     7.34   20.8%    32.7%
평균 MAPE: ML 22.3% / MA3 31.5% -> 목표(MAPE 25% 이하): 달성

[피처 중요도 상위 6개(MAE 증가량)] lag12 7.83, item_level 7.20, op_hours 0.95, std6 0.27, ma3 0.19, op_ratio 0.13
중요도가 0에 가까운 피처(차원 축소 후보): ma6, trend_t, lag2

[결전 대비 다음 달 수요 예측 - 증가량 상위 6개 품목]
         평시 계획  작전 템포 1.5배   증가율        group status
item_id                                             
P016     138.9       168.2  +21%       gasket     단종
P028      74.1        99.5  +34%  power_train     단종
P033      73.3        94.1  +28%         fuel    생산중
P035      78.5        94.1  +20%  brake_steer    생산중
P030      43.4        56.8  +31%  brake_steer   단종예고
P026      15.6        24.3  +56%      cooling    생산중
전체 예측 수요: 평시 946개 -> 작전 템포 1.5배 1210개

미션 5. 적정 안전재고량 산출 및 결품 위험 지수(Stockout Risk) 엔진 개발

  • 난이도: 중급~고급 / 예상 소요시간: 3.5시간
  • 코드 설명: 포아송 그래디언트 부스팅·랜덤포레스트 앙상블로 수요 예측을 고도화해 이동평균 기준 모델 대비 MAE/RMSE 개선율과 R2를 검증하고, 예측 오차로 품목별 적정 안전재고·재주문점을 산출한 뒤 리드타임 결품 확률과 단종 상태·단일 공급업체·대체품 부재를 결합한 결품 위험 지수 엔진을 만든다.
"""포아송 그래디언트 부스팅·랜덤포레스트 앙상블로 수요 예측을 고도화해 이동평균 기준 모델 대비 MAE/RMSE 개선율과 R2를 검증하고, 예측 오차로 품목별 적정 안전재고·재주문점을 산출한 뒤 리드타임 결품 확률과 단종 상태·단일 공급업체·대체품 부재를 결합한 결품 위험 지수 엔진을 만든다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from sklearn.ensemble import HistGradientBoostingRegressor, RandomForestRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 야전 정비창 수리부속 청구/출고·장비 가동률 데이터 생성기
# 실제 부대·품목·재고·단가 정보와 무관하다. 품목 40개 x 48개월(가상 2022-01 ~ 2025-12) 월별 자료를 만든다.
# ---------------------------------------------------------------------------
FSC = {"2520": "power_train", "2530": "brake_steer", "2910": "fuel", "2930": "cooling",
       "2940": "filter", "5330": "gasket", "6140": "battery", "5961": "semiconductor"}
EXERCISE = {3, 8, 10}      # 가상 대규모 훈련 월(훈련 계획으로 미리 알 수 있음)
N_MONTHS = 48


def season(fsc, month):
    """품목군별 계절 계수: 축전지는 겨울, 냉각계통은 여름에 수요가 늘어난다."""
    if fsc == "6140":
        return 1.6 if month in (12, 1, 2) else 0.9
    if fsc == "2930":
        return 1.5 if month in (6, 7, 8) else 0.9
    return 1 + 0.1 * np.sin(2 * np.pi * month / 12)


def make_items(n=40, seed=2):
    """품목 마스터: NSN(가상), 품목군, 장비 유형, 기준 수요, 단가, 조달 리드타임, 공급업체 수, 생산(단종) 상태."""
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n):
        f = str(rng.choice(list(FSC)))
        status = str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1]))
        rows.append({"item_id": f"P{i + 1:03d}", "nsn": f"{f}-00-{rng.integers(100, 1000)}-{rng.integers(1000, 10000)}",
                     "fsc": f, "group": FSC[f], "equip_type": str(rng.choice(["A", "B", "C"])),
                     "base": float(np.round(rng.lognormal(2.7, 0.8), 1)),
                     "price": int(rng.lognormal(12, 1.0)) // 100 * 100,
                     "lead_time": int(rng.integers(1, 7) if status == "생산중" else rng.integers(4, 13)),
                     "suppliers": int(rng.integers(1, 5)) if status == "생산중" else 1,
                     "status": status, "alt_available": bool(rng.random() < (0.6 if status == "생산중" else 0.3))})
    return pd.DataFrame(rows)


def make_history(items, months=N_MONTHS, seed=3):
    """월별 청구 수요(개), 긴급 불출 건수, 장비 유형별 운용시간·가동률."""
    rng = np.random.default_rng(seed)
    t = np.arange(months)
    moy = t % 12 + 1
    ex = np.isin(moy, list(EXERCISE))
    hours = {e: (100 + 20 * rng.random(months)) * np.where(ex, 1.5, 1.0) for e in "ABC"}
    avail = {e: np.clip(0.93 - 0.12 * (hours[e] - 100) / 80 + rng.normal(0, 0.015, months), 0.6, 0.99) for e in "ABC"}
    rows = []
    for it in items.itertuples():
        trend = 1 + rng.uniform(-0.1, 0.3) * t / months
        s = np.array([season(it.fsc, m) for m in moy])
        lam = it.base * s * trend * hours[it.equip_type] / 110
        d = rng.poisson(lam)
        rows.append(pd.DataFrame({"item_id": it.item_id, "t": t, "ym": [f"{2022 + k // 12}-{k % 12 + 1:02d}" for k in t],
                                  "month": moy, "exercise": ex.astype(int), "op_hours": hours[it.equip_type].round(1),
                                  "availability": avail[it.equip_type].round(3), "demand": d,
                                  "emergency": rng.binomial(d, 0.05)}))
    return pd.concat(rows, ignore_index=True)


def make_stock(items, hist, seed=4):
    """현재(가상 2026-01) 재고: 보유량과 미입고 발주량(조달 중인 물량). 단종 품목은 추가 발주가 어려워 미입고량이 적다."""
    rng = np.random.default_rng(seed)
    recent = hist[hist["t"] >= hist["t"].max() - 5].groupby("item_id")["demand"].mean()
    st = items[["item_id", "lead_time", "status"]].copy()
    st["on_hand"] = [int(round(recent[i] * rng.uniform(0.5, 3.0))) for i in st["item_id"]]
    cover = np.where(st["status"] == "생산중", rng.uniform(0.6, 1.3, len(st)), rng.uniform(0.0, 0.6, len(st)))
    st["on_order"] = [int(round(recent[i] * L * c)) for i, L, c in zip(st["item_id"], st["lead_time"], cover)]
    return st[["item_id", "on_hand", "on_order"]]
# ---------------------------------------------------------------------------

# ---- 수요 예측 공통 함수(파생 피처 + 머신러닝 모델) ----
FEATS = ["lag1", "lag2", "lag3", "lag12", "ma3", "ma6", "std6", "season_idx", "month_sin", "month_cos",
         "exercise", "op_hours", "op_ratio", "item_level", "trend_t"]


def make_features(hist):
    """품목별 과거 수요로 파생 피처를 만든다(모두 예측 시점 이전 정보만 사용)."""
    df = hist.sort_values(["item_id", "t"]).reset_index(drop=True)
    g = df.groupby("item_id")["demand"]
    for k in (1, 2, 3, 12):
        df[f"lag{k}"] = g.shift(k)
    df["ma3"] = g.transform(lambda s: s.shift(1).rolling(3).mean())
    df["ma6"] = g.transform(lambda s: s.shift(1).rolling(6).mean())
    df["std6"] = g.transform(lambda s: s.shift(1).rolling(6).std())
    df["item_level"] = g.transform(lambda s: s.shift(1).rolling(12).mean())
    df["season_idx"] = df["lag12"] / (g.transform(lambda s: s.shift(12).rolling(12).mean()) + 1)
    df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
    df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)
    df["op_ratio"] = df["op_hours"] / df.groupby("item_id")["op_hours"].transform(lambda s: s.shift(1).rolling(12).mean())
    df["trend_t"] = df["t"]
    return df.dropna().reset_index(drop=True)


def new_model():
    """Prophet/LSTM 대신 사용하는 포아송 손실 그래디언트 부스팅(개수 자료에 적합)."""
    return HistGradientBoostingRegressor(loss="poisson", max_iter=300, learning_rate=0.05, max_leaf_nodes=15,
                                         min_samples_leaf=20, random_state=0)


def mape(y, p, min_y=1):
    y, p = np.asarray(y, float), np.asarray(p, float)
    m = y >= min_y
    return float(np.mean(np.abs(y[m] - p[m]) / y[m]))

# ---- 앙상블 수요 예측·안전재고·결품(단종) 위험 지수 공통 함수 ----
STATUS_W = {"생산중": 0.0, "단종예고": 0.5, "단종": 1.0}     # [교육용가상자료] 단종 상태 가중치


def fit_ens(tr):
    """포아송 그래디언트 부스팅(70%) + 랜덤포레스트(30%) 앙상블."""
    h = new_model().fit(tr[FEATS], tr["demand"])
    r = RandomForestRegressor(n_estimators=150, min_samples_leaf=5, random_state=0, n_jobs=1).fit(tr[FEATS], tr["demand"])
    return h, r


def pred_ens(m, X, w=0.7):
    return np.clip(w * m[0].predict(X[FEATS]) + (1 - w) * m[1].predict(X[FEATS]), 0, None)


def next_month(hist, items, t_new=48, op_mult=1.0, ex=0):
    """다음 달(t_new) 예측용 행: 작년 같은 달 운용시간 x 배수, 훈련 여부 지정."""
    last = hist[hist["t"] == t_new - 1].copy()
    prev = hist[hist["t"] == t_new - 12].set_index("item_id")["op_hours"]
    m = (t_new % 12) + 1
    f = last.assign(t=t_new, ym=f"{2022 + t_new // 12}-{m:02d}", month=m, demand=0, exercise=ex,
                    op_hours=last["item_id"].map(prev).to_numpy() * op_mult)
    X = make_features(pd.concat([hist, f], ignore_index=True))
    return X[X["t"] == t_new].set_index("item_id").loc[items["item_id"]].reset_index()


def stock_policy(items, mu, sigma, z=1.65, review=1):
    """안전재고 = z x 월 예측오차 표준편차 x sqrt(리드타임+검토주기), 재주문점 = 월 수요 x (리드타임+검토주기) + 안전재고."""
    L = items["lead_time"].to_numpy() + review
    ss = z * sigma * np.sqrt(L)
    return pd.DataFrame({"item_id": items["item_id"], "mu": mu, "sigma": sigma, "L": L,
                         "ss": np.ceil(ss), "rop": np.ceil(mu * L + ss)})


def risk_index(items, pol, stock):
    """결품 위험 지수(0~100): 리드타임 동안 수요가 가용 재고를 넘을 확률(정규 근사)을 기본으로,
    단종 상태·단일 공급업체·대체품 부재에 가중을 준다(가상 가중치)."""
    st = stock.set_index("item_id").loc[items["item_id"]]
    avail = (st["on_hand"] + st["on_order"]).to_numpy()
    p = 1 - norm.cdf((avail - pol["mu"].to_numpy() * pol["L"].to_numpy()) / (pol["sigma"].to_numpy() * np.sqrt(pol["L"].to_numpy()) + 1e-6))
    sw = items["status"].map(STATUS_W).to_numpy()
    sup = (items["suppliers"].to_numpy() == 1).astype(float)
    alt = (~items["alt_available"].to_numpy()).astype(float)
    score = 100 * np.clip(0.55 * p + 0.25 * sw + 0.10 * sup + 0.10 * alt, 0, 1)
    return pd.DataFrame({"item_id": items["item_id"], "p_short": p, "status_w": sw, "single_sup": sup, "no_alt": alt,
                         "avail": avail, "risk": score.round(1)})

items = make_items()
hist = make_history(items)
df = make_features(hist)
stock = make_stock(items, hist)

# 1) Hold-out(t=42~47) 성능 비교: 기준 모델(MA3) / 단일 모델 / 앙상블. 앙상블 가중치는 t=36~41 검증 구간에서 고른다
inner = fit_ens(df[df["t"] < 36]); vv = df[(df["t"] >= 36) & (df["t"] < 42)]
best_w = min((0.5, 0.6, 0.7, 0.8, 0.9, 1.0), key=lambda w: mean_absolute_error(vv["demand"], pred_ens(inner, vv, w)))
tr, te = df[df["t"] < 42], df[df["t"] >= 42].copy()
ens = fit_ens(tr)
te["ens"], te["hgb"] = pred_ens(ens, te, best_w), np.clip(ens[0].predict(te[FEATS]), 0, None)
rows = {}
for nm, col in [("기준 모델(MA3)", "ma3"), ("단일 포아송 GBM", "hgb"), (f"앙상블(w={best_w})", "ens")]:
    rows[nm] = {"MAE": mean_absolute_error(te["demand"], te[col]), "RMSE": mean_squared_error(te["demand"], te[col]) ** 0.5,
                "R2": r2_score(te["demand"], te[col])}
res = pd.DataFrame(rows).T
print("[Hold-out 6개월(t=42~47), 품목 40개]")
print(res.to_string(float_format=lambda x: f"{x:.3f}"))
b, e = res.iloc[0], res.iloc[2]
im, ir = 1 - e["MAE"] / b["MAE"], 1 - e["RMSE"] / b["RMSE"]
print(f"오차율 개선 MAE {im:.1%}, RMSE {ir:.1%} (목표 20% 이상: {'달성' if min(im, ir) >= 0.2 else '미달'}), R2 {e['R2']:.3f} (목표 0.85 이상: {'달성' if e['R2'] >= 0.85 else '미달'})")
spike = te["exercise"] == 1
print(f"불규칙 수요(훈련 월) MAE: 기준 {mean_absolute_error(te.loc[spike, 'demand'], te.loc[spike, 'ma3']):.2f} -> 앙상블 {mean_absolute_error(te.loc[spike, 'demand'], te.loc[spike, 'ens']):.2f}")

# 2) 안전재고·재주문점: 다음 달 예측과 Hold-out 품목별 오차 표준편차
X = next_month(hist, items)
mu = pred_ens(fit_ens(df), X, best_w)
err = (te["demand"] - te["ens"]).groupby(te["item_id"]).std().reindex(items["item_id"]).to_numpy()
sigma = np.maximum(err, np.sqrt(np.maximum(mu, 1)))          # 오차 표준편차가 포아송 변동보다 작으면 포아송 기준 사용
pol = stock_policy(items, mu, sigma)
risk = risk_index(items, pol, stock)
out = items[["item_id", "group", "status", "lead_time", "suppliers"]].merge(pol, on="item_id").merge(risk, on="item_id")
print(f"\n[다음 달(가상 2026-01) 안전재고 정책] 서비스 수준 95%(z=1.65), 안전재고 합 {out['ss'].sum():.0f}개")
top = out.sort_values("risk", ascending=False).head(8)
show = top[["item_id", "group", "status", "lead_time", "mu", "ss", "rop", "avail", "p_short", "risk"]].copy()
show.columns = ["품목", "품목군", "생산상태", "LT(월)", "월수요", "안전재고", "재주문점", "가용재고", "결품확률", "위험지수"]
print("[결품 위험 지수 상위 8개]")
print(show.to_string(index=False, float_format=lambda x: f"{x:.2f}"))
grade = pd.cut(out["risk"], [-1, 30, 60, 101], labels=["관찰", "주의", "경보"]).value_counts()
print(f"위험 등급: 경보 {grade['경보']}개, 주의 {grade['주의']}개, 관찰 {grade['관찰']}개 / 재주문점 미만 품목 {(out['avail'] < out['rop']).sum()}개")
print("단종·단종예고 품목은 리드타임이 길고 공급업체가 1곳이라 같은 결품 확률에서도 위험 지수가 높게 나온다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].scatter(te["demand"], te["ma3"], s=8, alpha=0.5, label="MA3 baseline"); ax[0].scatter(te["demand"], te["ens"], s=8, alpha=0.5, label="ensemble")
mx = te["demand"].max(); ax[0].plot([0, mx], [0, mx], "k--", lw=1); ax[0].set_xlabel("actual"); ax[0].set_ylabel("forecast"); ax[0].legend()
ax[0].set_title("Hold-out actual vs forecast")
col = out["status"].map({"생산중": "tab:green", "단종예고": "tab:orange", "단종": "tab:red"})
ax[1].scatter(out["p_short"], out["risk"], c=col); ax[1].set_xlabel("shortage probability over lead time"); ax[1].set_ylabel("risk index")
ax[1].set_title("Stockout risk (green=active, orange=EOL notice, red=discontinued)")
plt.tight_layout(); plt.show()

실행 결과 예시

[Hold-out 6개월(t=42~47), 품목 40개]
             MAE   RMSE    R2
기준 모델(MA3) 7.339 11.130 0.839
단일 포아송 GBM 4.572  6.112 0.952
앙상블(w=0.5) 4.541  6.059 0.952
오차율 개선 MAE 38.1%, RMSE 45.6% (목표 20% 이상: 달성), R2 0.952 (목표 0.85 이상: 달성)
불규칙 수요(훈련 월) MAE: 기준 9.88 -> 앙상블 4.94

[다음 달(가상 2026-01) 안전재고 정책] 서비스 수준 95%(z=1.65), 안전재고 합 926개
[결품 위험 지수 상위 8개]
  품목         품목군 생산상태  LT(월)    월수요  안전재고    재주문점  가용재고  결품확률  위험지수
P018     battery   단종      4  17.48 23.00  111.00    71  0.88 93.60
P016      gasket   단종     12 148.02 95.00 2019.00   714  1.00 90.00
P028 power_train   단종     11  68.99 48.00  876.00   515  1.00 90.00
P004      filter 단종예고      6   7.65 16.00   70.00    12  1.00 87.50
P025      filter 단종예고     12  13.94 23.00  204.00    89  1.00 87.50
P030 brake_steer 단종예고     12  45.84 43.00  639.00   162  1.00 87.50
P038        fuel 단종예고      6  11.87 25.00  108.00    40  1.00 87.40
P019      filter 단종예고      8  11.47 30.00  133.00    59  0.99 87.10
위험 등급: 경보 13개, 주의 2개, 관찰 25개 / 재주문점 미만 품목 20개
단종·단종예고 품목은 리드타임이 길고 공급업체가 1곳이라 같은 결품 확률에서도 위험 지수가 높게 나온다.

미션 6. 야전 부품 재고 현황 및 긴급 조달 필요 품목 가시화 대시보드 구현

  • 난이도: 중급~고급 / 예상 소요시간: 3.0시간
  • 코드 설명: 야전 부품 재고 현황(가용재고·재주문점·결품 위험 지수)과 긴급 조달 필요 품목을 보여 주는 대시보드를 만들고, 품목마다 위험 지수를 끌어올린 요인을 기준값 대체(occlusion) 방식으로 분해해 기여도 상위 3개 인자와 자연어 설명문('리드타임 증가 원인은 단종' 등)을 자동 생성한다.
"""야전 부품 재고 현황(가용재고·재주문점·결품 위험 지수)과 긴급 조달 필요 품목을 보여 주는 대시보드를 만들고, 품목마다 위험 지수를 끌어올린 요인을 기준값 대체(occlusion) 방식으로 분해해 기여도 상위 3개 인자와 자연어 설명문('리드타임 증가 원인은 단종' 등)을 자동 생성한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from sklearn.ensemble import HistGradientBoostingRegressor, RandomForestRegressor

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 야전 정비창 수리부속 청구/출고·장비 가동률 데이터 생성기
# 실제 부대·품목·재고·단가 정보와 무관하다. 품목 40개 x 48개월(가상 2022-01 ~ 2025-12) 월별 자료를 만든다.
# ---------------------------------------------------------------------------
FSC = {"2520": "power_train", "2530": "brake_steer", "2910": "fuel", "2930": "cooling",
       "2940": "filter", "5330": "gasket", "6140": "battery", "5961": "semiconductor"}
EXERCISE = {3, 8, 10}      # 가상 대규모 훈련 월(훈련 계획으로 미리 알 수 있음)
N_MONTHS = 48


def season(fsc, month):
    """품목군별 계절 계수: 축전지는 겨울, 냉각계통은 여름에 수요가 늘어난다."""
    if fsc == "6140":
        return 1.6 if month in (12, 1, 2) else 0.9
    if fsc == "2930":
        return 1.5 if month in (6, 7, 8) else 0.9
    return 1 + 0.1 * np.sin(2 * np.pi * month / 12)


def make_items(n=40, seed=2):
    """품목 마스터: NSN(가상), 품목군, 장비 유형, 기준 수요, 단가, 조달 리드타임, 공급업체 수, 생산(단종) 상태."""
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n):
        f = str(rng.choice(list(FSC)))
        status = str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1]))
        rows.append({"item_id": f"P{i + 1:03d}", "nsn": f"{f}-00-{rng.integers(100, 1000)}-{rng.integers(1000, 10000)}",
                     "fsc": f, "group": FSC[f], "equip_type": str(rng.choice(["A", "B", "C"])),
                     "base": float(np.round(rng.lognormal(2.7, 0.8), 1)),
                     "price": int(rng.lognormal(12, 1.0)) // 100 * 100,
                     "lead_time": int(rng.integers(1, 7) if status == "생산중" else rng.integers(4, 13)),
                     "suppliers": int(rng.integers(1, 5)) if status == "생산중" else 1,
                     "status": status, "alt_available": bool(rng.random() < (0.6 if status == "생산중" else 0.3))})
    return pd.DataFrame(rows)


def make_history(items, months=N_MONTHS, seed=3):
    """월별 청구 수요(개), 긴급 불출 건수, 장비 유형별 운용시간·가동률."""
    rng = np.random.default_rng(seed)
    t = np.arange(months)
    moy = t % 12 + 1
    ex = np.isin(moy, list(EXERCISE))
    hours = {e: (100 + 20 * rng.random(months)) * np.where(ex, 1.5, 1.0) for e in "ABC"}
    avail = {e: np.clip(0.93 - 0.12 * (hours[e] - 100) / 80 + rng.normal(0, 0.015, months), 0.6, 0.99) for e in "ABC"}
    rows = []
    for it in items.itertuples():
        trend = 1 + rng.uniform(-0.1, 0.3) * t / months
        s = np.array([season(it.fsc, m) for m in moy])
        lam = it.base * s * trend * hours[it.equip_type] / 110
        d = rng.poisson(lam)
        rows.append(pd.DataFrame({"item_id": it.item_id, "t": t, "ym": [f"{2022 + k // 12}-{k % 12 + 1:02d}" for k in t],
                                  "month": moy, "exercise": ex.astype(int), "op_hours": hours[it.equip_type].round(1),
                                  "availability": avail[it.equip_type].round(3), "demand": d,
                                  "emergency": rng.binomial(d, 0.05)}))
    return pd.concat(rows, ignore_index=True)


def make_stock(items, hist, seed=4):
    """현재(가상 2026-01) 재고: 보유량과 미입고 발주량(조달 중인 물량). 단종 품목은 추가 발주가 어려워 미입고량이 적다."""
    rng = np.random.default_rng(seed)
    recent = hist[hist["t"] >= hist["t"].max() - 5].groupby("item_id")["demand"].mean()
    st = items[["item_id", "lead_time", "status"]].copy()
    st["on_hand"] = [int(round(recent[i] * rng.uniform(0.5, 3.0))) for i in st["item_id"]]
    cover = np.where(st["status"] == "생산중", rng.uniform(0.6, 1.3, len(st)), rng.uniform(0.0, 0.6, len(st)))
    st["on_order"] = [int(round(recent[i] * L * c)) for i, L, c in zip(st["item_id"], st["lead_time"], cover)]
    return st[["item_id", "on_hand", "on_order"]]
# ---------------------------------------------------------------------------

# ---- 수요 예측 공통 함수(파생 피처 + 머신러닝 모델) ----
FEATS = ["lag1", "lag2", "lag3", "lag12", "ma3", "ma6", "std6", "season_idx", "month_sin", "month_cos",
         "exercise", "op_hours", "op_ratio", "item_level", "trend_t"]


def make_features(hist):
    """품목별 과거 수요로 파생 피처를 만든다(모두 예측 시점 이전 정보만 사용)."""
    df = hist.sort_values(["item_id", "t"]).reset_index(drop=True)
    g = df.groupby("item_id")["demand"]
    for k in (1, 2, 3, 12):
        df[f"lag{k}"] = g.shift(k)
    df["ma3"] = g.transform(lambda s: s.shift(1).rolling(3).mean())
    df["ma6"] = g.transform(lambda s: s.shift(1).rolling(6).mean())
    df["std6"] = g.transform(lambda s: s.shift(1).rolling(6).std())
    df["item_level"] = g.transform(lambda s: s.shift(1).rolling(12).mean())
    df["season_idx"] = df["lag12"] / (g.transform(lambda s: s.shift(12).rolling(12).mean()) + 1)
    df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
    df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)
    df["op_ratio"] = df["op_hours"] / df.groupby("item_id")["op_hours"].transform(lambda s: s.shift(1).rolling(12).mean())
    df["trend_t"] = df["t"]
    return df.dropna().reset_index(drop=True)


def new_model():
    """Prophet/LSTM 대신 사용하는 포아송 손실 그래디언트 부스팅(개수 자료에 적합)."""
    return HistGradientBoostingRegressor(loss="poisson", max_iter=300, learning_rate=0.05, max_leaf_nodes=15,
                                         min_samples_leaf=20, random_state=0)


def mape(y, p, min_y=1):
    y, p = np.asarray(y, float), np.asarray(p, float)
    m = y >= min_y
    return float(np.mean(np.abs(y[m] - p[m]) / y[m]))

# ---- 앙상블 수요 예측·안전재고·결품(단종) 위험 지수 공통 함수 ----
STATUS_W = {"생산중": 0.0, "단종예고": 0.5, "단종": 1.0}     # [교육용가상자료] 단종 상태 가중치


def fit_ens(tr):
    """포아송 그래디언트 부스팅(70%) + 랜덤포레스트(30%) 앙상블."""
    h = new_model().fit(tr[FEATS], tr["demand"])
    r = RandomForestRegressor(n_estimators=150, min_samples_leaf=5, random_state=0, n_jobs=1).fit(tr[FEATS], tr["demand"])
    return h, r


def pred_ens(m, X, w=0.7):
    return np.clip(w * m[0].predict(X[FEATS]) + (1 - w) * m[1].predict(X[FEATS]), 0, None)


def next_month(hist, items, t_new=48, op_mult=1.0, ex=0):
    """다음 달(t_new) 예측용 행: 작년 같은 달 운용시간 x 배수, 훈련 여부 지정."""
    last = hist[hist["t"] == t_new - 1].copy()
    prev = hist[hist["t"] == t_new - 12].set_index("item_id")["op_hours"]
    m = (t_new % 12) + 1
    f = last.assign(t=t_new, ym=f"{2022 + t_new // 12}-{m:02d}", month=m, demand=0, exercise=ex,
                    op_hours=last["item_id"].map(prev).to_numpy() * op_mult)
    X = make_features(pd.concat([hist, f], ignore_index=True))
    return X[X["t"] == t_new].set_index("item_id").loc[items["item_id"]].reset_index()


def stock_policy(items, mu, sigma, z=1.65, review=1):
    """안전재고 = z x 월 예측오차 표준편차 x sqrt(리드타임+검토주기), 재주문점 = 월 수요 x (리드타임+검토주기) + 안전재고."""
    L = items["lead_time"].to_numpy() + review
    ss = z * sigma * np.sqrt(L)
    return pd.DataFrame({"item_id": items["item_id"], "mu": mu, "sigma": sigma, "L": L,
                         "ss": np.ceil(ss), "rop": np.ceil(mu * L + ss)})


def risk_index(items, pol, stock):
    """결품 위험 지수(0~100): 리드타임 동안 수요가 가용 재고를 넘을 확률(정규 근사)을 기본으로,
    단종 상태·단일 공급업체·대체품 부재에 가중을 준다(가상 가중치)."""
    st = stock.set_index("item_id").loc[items["item_id"]]
    avail = (st["on_hand"] + st["on_order"]).to_numpy()
    p = 1 - norm.cdf((avail - pol["mu"].to_numpy() * pol["L"].to_numpy()) / (pol["sigma"].to_numpy() * np.sqrt(pol["L"].to_numpy()) + 1e-6))
    sw = items["status"].map(STATUS_W).to_numpy()
    sup = (items["suppliers"].to_numpy() == 1).astype(float)
    alt = (~items["alt_available"].to_numpy()).astype(float)
    score = 100 * np.clip(0.55 * p + 0.25 * sw + 0.10 * sup + 0.10 * alt, 0, 1)
    return pd.DataFrame({"item_id": items["item_id"], "p_short": p, "status_w": sw, "single_sup": sup, "no_alt": alt,
                         "avail": avail, "risk": score.round(1)})

items = make_items()
hist = make_history(items)
df = make_features(hist)
stock = make_stock(items, hist)
X = next_month(hist, items)
mu = pred_ens(fit_ens(df), X)
level = X["item_level"].to_numpy()                       # 최근 12개월 평균 수요(평시 수준)
sigma = np.maximum(0.25 * mu, np.sqrt(np.maximum(mu, 1)))


def score(it, m, s, st):
    return risk_index(it, stock_policy(it, m, s), st)["risk"].to_numpy()


base = score(items, mu, sigma, stock)
pol = stock_policy(items, mu, sigma)
# 각 요인을 '문제 없는 기준값'으로 바꿨을 때 위험 지수가 얼마나 내려가는지 = 그 요인의 기여도
ref_lt = int(items.loc[items["status"] == "생산중", "lead_time"].median())
st_ok = stock.copy(); st_ok["on_hand"] = np.maximum(st_ok["on_hand"], (pol["rop"] - st_ok["on_order"]).clip(lower=0).astype(int).to_numpy())
FACTORS = {
    "조달 리드타임 장기화": lambda: score(items.assign(lead_time=np.minimum(items["lead_time"], ref_lt)), mu, sigma, stock),
    "수요 증가 예측": lambda: score(items, np.minimum(mu, level), sigma, stock),
    "가용재고 부족": lambda: score(items, mu, sigma, st_ok),
    "예측 불확실성": lambda: score(items, mu, np.sqrt(np.maximum(mu, 1)), stock),
    "단종(예고) 상태": lambda: score(items.assign(status="생산중"), mu, sigma, stock),
    "단일 공급업체": lambda: score(items.assign(suppliers=np.maximum(items["suppliers"], 2)), mu, sigma, stock),
    "대체품 부재": lambda: score(items.assign(alt_available=True), mu, sigma, stock),
}
C = pd.DataFrame({k: base - f() for k, f in FACTORS.items()}, index=items["item_id"]).clip(lower=0)
top3 = C.apply(lambda r: list(r.sort_values(ascending=False).index[:3]), axis=1)
cover = np.mean([len(t) == 3 for t in top3])
print(f"[품목 {len(items)}개] 위험 지수 기여도 상위 3개 인자 자동 산출률 {cover:.0%} (목표 100%: {'달성' if cover == 1 else '미달'})")

dash = items[["item_id", "group", "status", "lead_time"]].assign(risk=base, avail=(stock["on_hand"] + stock["on_order"]).to_numpy(),
                                                                 rop=pol["rop"].to_numpy(), mu=mu.round(1))
dash["urgent"] = (dash["risk"] >= 60) & (dash["avail"] < dash["rop"])
CAUSE = {"조달 리드타임 장기화": "조달 리드타임이 길다", "단종(예고) 상태": "생산 중단(예고)으로 추가 조달이 어렵다",
         "가용재고 부족": "보유·미입고 재고가 재주문점보다 적다", "수요 증가 예측": "다음 달 수요가 평시보다 늘 것으로 예측된다",
         "예측 불확실성": "수요 변동이 커 안전재고가 많이 필요하다", "단일 공급업체": "공급업체가 1곳뿐이다", "대체품 부재": "대체품이 없다"}


def explain(iid):
    r = dash.set_index("item_id").loc[iid]
    f = top3[iid]
    parts = ", ".join(f"{k}(+{C.loc[iid, k]:.0f})" for k in f)
    return f"{iid}({r['group']}, {r['status']}) 위험 {r['risk']:.0f}: {parts}. 주된 원인은 {CAUSE[f[0]]}."


urg = dash[dash["urgent"]].sort_values("risk", ascending=False)
print(f"\n[긴급 조달 필요 품목 {len(urg)}개(위험 60 이상 & 가용재고<재주문점)] 자연어 설명")
for iid in urg["item_id"].head(7):
    print(" -", explain(iid))
print("\n[요인별 평균 기여도]", ", ".join(f"{k} {v:.1f}" for k, v in C.mean().sort_values(ascending=False).items()))
g = dash.groupby("status")["risk"].agg(["count", "mean"]).round(1)
print("[생산 상태별 평균 위험 지수]", ", ".join(f"{s} {r['mean']}({int(r['count'])}개)" for s, r in g.iterrows()))
print("설명은 조달·대체품 검토의 근거 자료이며 발주·대체 결정은 정비장교가 확인 후 내린다.")

fig, ax = plt.subplots(1, 3, figsize=(17, 4.8))
d15 = dash.sort_values("risk", ascending=False).head(15)
ax[0].barh(d15["item_id"], d15["avail"], color=np.where(d15["urgent"], "tab:red", "tab:gray"), label="available")
ax[0].scatter(d15["rop"], d15["item_id"], c="k", marker="|", s=200, label="reorder point"); ax[0].invert_yaxis()
ax[0].set_xscale("symlog"); ax[0].set_title("Top-15 risk items: stock vs reorder point (red=urgent)"); ax[0].legend()
hm = dash.pivot_table(index="group", columns="status", values="risk", aggfunc="mean")
im = ax[1].imshow(hm.values, cmap="Reds", vmin=0, vmax=100, aspect="auto"); plt.colorbar(im, ax=ax[1])
ax[1].set_xticks(range(hm.shape[1])); ax[1].set_xticklabels(["active" if c == "생산중" else "EOL notice" if c == "단종예고" else "discontinued" for c in hm.columns])
ax[1].set_yticks(range(len(hm))); ax[1].set_yticklabels(hm.index); ax[1].set_title("Mean risk index by group/status")
i0 = urg["item_id"].iloc[0]; w = C.loc[i0].sort_values()
ax[2].barh(range(len(w)), w.values, color="tab:red"); ax[2].set_yticks(range(len(w)))
EN = {"조달 리드타임 장기화": "lead time", "수요 증가 예측": "demand up", "가용재고 부족": "low stock", "예측 불확실성": "uncertainty",
      "단종(예고) 상태": "EOL status", "단일 공급업체": "single supplier", "대체품 부재": "no alternate"}
ax[2].set_yticklabels([EN[k] for k in w.index])
ax[2].set_title(f"Risk contribution: {i0}")
plt.tight_layout(); plt.show()

실행 결과 예시

[품목 40개] 위험 지수 기여도 상위 3개 인자 자동 산출률 100% (목표 100%: 달성)

[긴급 조달 필요 품목 13개(위험 60 이상 & 가용재고<재주문점)] 자연어 설명
 - P018(battery, 단종) 위험 96: 가용재고 부족(+49), 단종(예고) 상태(+25), 단일 공급업체(+10). 주된 원인은 보유·미입고 재고가 재주문점보다 적다.
 - P028(power_train, 단종) 위험 90: 조달 리드타임 장기화(+55), 가용재고 부족(+52), 단종(예고) 상태(+25). 주된 원인은 조달 리드타임이 길다.
 - P016(gasket, 단종) 위험 90: 가용재고 부족(+52), 단종(예고) 상태(+25), 단일 공급업체(+10). 주된 원인은 보유·미입고 재고가 재주문점보다 적다.
 - P004(filter, 단종예고) 위험 88: 가용재고 부족(+53), 단종(예고) 상태(+12), 단일 공급업체(+10). 주된 원인은 보유·미입고 재고가 재주문점보다 적다.
 - P025(filter, 단종예고) 위험 88: 가용재고 부족(+53), 조달 리드타임 장기화(+39), 단종(예고) 상태(+12). 주된 원인은 보유·미입고 재고가 재주문점보다 적다.
 - P038(fuel, 단종예고) 위험 88: 가용재고 부족(+52), 단종(예고) 상태(+12), 단일 공급업체(+10). 주된 원인은 보유·미입고 재고가 재주문점보다 적다.
 - P030(brake_steer, 단종예고) 위험 88: 가용재고 부족(+52), 단종(예고) 상태(+12), 단일 공급업체(+10). 주된 원인은 보유·미입고 재고가 재주문점보다 적다.

[요인별 평균 기여도] 가용재고 부족 19.5, 단일 공급업체 4.8, 조달 리드타임 장기화 4.5, 단종(예고) 상태 4.4, 대체품 부재 3.8, 수요 증가 예측 0.3, 예측 불확실성 0.3
[생산 상태별 평균 위험 지수] 단종 92.0(3개), 단종예고 84.5(8개), 생산중 13.7(29개)
설명은 조달·대체품 검토의 근거 자료이며 발주·대체 결정은 정비장교가 확인 후 내린다.

미션 7. 부품 발주량 최종 확정을 위한 정비장교 HITL 승인/조정 워크플로우 설계

  • 난이도: 고급 / 예상 소요시간: 3.5시간
  • 코드 설명: 예측 기반 발주 정책(안전재고 포함 목표재고)과 단종 품목 대체품 전환을 12개월 재고 백테스트로 현행(최근 평균+고정 여유) 정책과 비교해 비용·재고 회전율·충족률을 산출하고, 다음 달 발주량을 예산 제약 선형계획(LP)으로 추천한 뒤 정비장교가 승인·조정·반려·보류로 최종 확정하는 HITL 워크플로우를 해시 체인 감사 로그와 함께 구현한다.
"""예측 기반 발주 정책(안전재고 포함 목표재고)과 단종 품목 대체품 전환을 12개월 재고 백테스트로 현행(최근 평균+고정 여유) 정책과 비교해 비용·재고 회전율·충족률을 산출하고, 다음 달 발주량을 예산 제약 선형계획(LP)으로 추천한 뒤 정비장교가 승인·조정·반려·보류로 최종 확정하는 HITL 워크플로우를 해시 체인 감사 로그와 함께 구현한다."""
import hashlib
import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import linprog
from sklearn.ensemble import HistGradientBoostingRegressor, RandomForestRegressor

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 야전 정비창 수리부속 청구/출고·장비 가동률 데이터 생성기
# 실제 부대·품목·재고·단가 정보와 무관하다. 품목 40개 x 48개월(가상 2022-01 ~ 2025-12) 월별 자료를 만든다.
# ---------------------------------------------------------------------------
FSC = {"2520": "power_train", "2530": "brake_steer", "2910": "fuel", "2930": "cooling",
       "2940": "filter", "5330": "gasket", "6140": "battery", "5961": "semiconductor"}
EXERCISE = {3, 8, 10}      # 가상 대규모 훈련 월(훈련 계획으로 미리 알 수 있음)
N_MONTHS = 48


def season(fsc, month):
    """품목군별 계절 계수: 축전지는 겨울, 냉각계통은 여름에 수요가 늘어난다."""
    if fsc == "6140":
        return 1.6 if month in (12, 1, 2) else 0.9
    if fsc == "2930":
        return 1.5 if month in (6, 7, 8) else 0.9
    return 1 + 0.1 * np.sin(2 * np.pi * month / 12)


def make_items(n=40, seed=2):
    """품목 마스터: NSN(가상), 품목군, 장비 유형, 기준 수요, 단가, 조달 리드타임, 공급업체 수, 생산(단종) 상태."""
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n):
        f = str(rng.choice(list(FSC)))
        status = str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1]))
        rows.append({"item_id": f"P{i + 1:03d}", "nsn": f"{f}-00-{rng.integers(100, 1000)}-{rng.integers(1000, 10000)}",
                     "fsc": f, "group": FSC[f], "equip_type": str(rng.choice(["A", "B", "C"])),
                     "base": float(np.round(rng.lognormal(2.7, 0.8), 1)),
                     "price": int(rng.lognormal(12, 1.0)) // 100 * 100,
                     "lead_time": int(rng.integers(1, 7) if status == "생산중" else rng.integers(4, 13)),
                     "suppliers": int(rng.integers(1, 5)) if status == "생산중" else 1,
                     "status": status, "alt_available": bool(rng.random() < (0.6 if status == "생산중" else 0.3))})
    return pd.DataFrame(rows)


def make_history(items, months=N_MONTHS, seed=3):
    """월별 청구 수요(개), 긴급 불출 건수, 장비 유형별 운용시간·가동률."""
    rng = np.random.default_rng(seed)
    t = np.arange(months)
    moy = t % 12 + 1
    ex = np.isin(moy, list(EXERCISE))
    hours = {e: (100 + 20 * rng.random(months)) * np.where(ex, 1.5, 1.0) for e in "ABC"}
    avail = {e: np.clip(0.93 - 0.12 * (hours[e] - 100) / 80 + rng.normal(0, 0.015, months), 0.6, 0.99) for e in "ABC"}
    rows = []
    for it in items.itertuples():
        trend = 1 + rng.uniform(-0.1, 0.3) * t / months
        s = np.array([season(it.fsc, m) for m in moy])
        lam = it.base * s * trend * hours[it.equip_type] / 110
        d = rng.poisson(lam)
        rows.append(pd.DataFrame({"item_id": it.item_id, "t": t, "ym": [f"{2022 + k // 12}-{k % 12 + 1:02d}" for k in t],
                                  "month": moy, "exercise": ex.astype(int), "op_hours": hours[it.equip_type].round(1),
                                  "availability": avail[it.equip_type].round(3), "demand": d,
                                  "emergency": rng.binomial(d, 0.05)}))
    return pd.concat(rows, ignore_index=True)


def make_stock(items, hist, seed=4):
    """현재(가상 2026-01) 재고: 보유량과 미입고 발주량(조달 중인 물량). 단종 품목은 추가 발주가 어려워 미입고량이 적다."""
    rng = np.random.default_rng(seed)
    recent = hist[hist["t"] >= hist["t"].max() - 5].groupby("item_id")["demand"].mean()
    st = items[["item_id", "lead_time", "status"]].copy()
    st["on_hand"] = [int(round(recent[i] * rng.uniform(0.5, 3.0))) for i in st["item_id"]]
    cover = np.where(st["status"] == "생산중", rng.uniform(0.6, 1.3, len(st)), rng.uniform(0.0, 0.6, len(st)))
    st["on_order"] = [int(round(recent[i] * L * c)) for i, L, c in zip(st["item_id"], st["lead_time"], cover)]
    return st[["item_id", "on_hand", "on_order"]]
# ---------------------------------------------------------------------------

# ---- 수요 예측 공통 함수(파생 피처 + 머신러닝 모델) ----
FEATS = ["lag1", "lag2", "lag3", "lag12", "ma3", "ma6", "std6", "season_idx", "month_sin", "month_cos",
         "exercise", "op_hours", "op_ratio", "item_level", "trend_t"]


def make_features(hist):
    """품목별 과거 수요로 파생 피처를 만든다(모두 예측 시점 이전 정보만 사용)."""
    df = hist.sort_values(["item_id", "t"]).reset_index(drop=True)
    g = df.groupby("item_id")["demand"]
    for k in (1, 2, 3, 12):
        df[f"lag{k}"] = g.shift(k)
    df["ma3"] = g.transform(lambda s: s.shift(1).rolling(3).mean())
    df["ma6"] = g.transform(lambda s: s.shift(1).rolling(6).mean())
    df["std6"] = g.transform(lambda s: s.shift(1).rolling(6).std())
    df["item_level"] = g.transform(lambda s: s.shift(1).rolling(12).mean())
    df["season_idx"] = df["lag12"] / (g.transform(lambda s: s.shift(12).rolling(12).mean()) + 1)
    df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
    df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)
    df["op_ratio"] = df["op_hours"] / df.groupby("item_id")["op_hours"].transform(lambda s: s.shift(1).rolling(12).mean())
    df["trend_t"] = df["t"]
    return df.dropna().reset_index(drop=True)


def new_model():
    """Prophet/LSTM 대신 사용하는 포아송 손실 그래디언트 부스팅(개수 자료에 적합)."""
    return HistGradientBoostingRegressor(loss="poisson", max_iter=300, learning_rate=0.05, max_leaf_nodes=15,
                                         min_samples_leaf=20, random_state=0)


def mape(y, p, min_y=1):
    y, p = np.asarray(y, float), np.asarray(p, float)
    m = y >= min_y
    return float(np.mean(np.abs(y[m] - p[m]) / y[m]))

# ---- 앙상블 수요 예측·안전재고·결품(단종) 위험 지수 공통 함수 ----
STATUS_W = {"생산중": 0.0, "단종예고": 0.5, "단종": 1.0}     # [교육용가상자료] 단종 상태 가중치


def fit_ens(tr):
    """포아송 그래디언트 부스팅(70%) + 랜덤포레스트(30%) 앙상블."""
    h = new_model().fit(tr[FEATS], tr["demand"])
    r = RandomForestRegressor(n_estimators=150, min_samples_leaf=5, random_state=0, n_jobs=1).fit(tr[FEATS], tr["demand"])
    return h, r


def pred_ens(m, X, w=0.7):
    return np.clip(w * m[0].predict(X[FEATS]) + (1 - w) * m[1].predict(X[FEATS]), 0, None)


def next_month(hist, items, t_new=48, op_mult=1.0, ex=0):
    """다음 달(t_new) 예측용 행: 작년 같은 달 운용시간 x 배수, 훈련 여부 지정."""
    last = hist[hist["t"] == t_new - 1].copy()
    prev = hist[hist["t"] == t_new - 12].set_index("item_id")["op_hours"]
    m = (t_new % 12) + 1
    f = last.assign(t=t_new, ym=f"{2022 + t_new // 12}-{m:02d}", month=m, demand=0, exercise=ex,
                    op_hours=last["item_id"].map(prev).to_numpy() * op_mult)
    X = make_features(pd.concat([hist, f], ignore_index=True))
    return X[X["t"] == t_new].set_index("item_id").loc[items["item_id"]].reset_index()


def stock_policy(items, mu, sigma, z=1.65, review=1):
    """안전재고 = z x 월 예측오차 표준편차 x sqrt(리드타임+검토주기), 재주문점 = 월 수요 x (리드타임+검토주기) + 안전재고."""
    L = items["lead_time"].to_numpy() + review
    ss = z * sigma * np.sqrt(L)
    return pd.DataFrame({"item_id": items["item_id"], "mu": mu, "sigma": sigma, "L": L,
                         "ss": np.ceil(ss), "rop": np.ceil(mu * L + ss)})


def risk_index(items, pol, stock):
    """결품 위험 지수(0~100): 리드타임 동안 수요가 가용 재고를 넘을 확률(정규 근사)을 기본으로,
    단종 상태·단일 공급업체·대체품 부재에 가중을 준다(가상 가중치)."""
    st = stock.set_index("item_id").loc[items["item_id"]]
    avail = (st["on_hand"] + st["on_order"]).to_numpy()
    p = 1 - norm.cdf((avail - pol["mu"].to_numpy() * pol["L"].to_numpy()) / (pol["sigma"].to_numpy() * np.sqrt(pol["L"].to_numpy()) + 1e-6))
    sw = items["status"].map(STATUS_W).to_numpy()
    sup = (items["suppliers"].to_numpy() == 1).astype(float)
    alt = (~items["alt_available"].to_numpy()).astype(float)
    score = 100 * np.clip(0.55 * p + 0.25 * sw + 0.10 * sup + 0.10 * alt, 0, 1)
    return pd.DataFrame({"item_id": items["item_id"], "p_short": p, "status_w": sw, "single_sup": sup, "no_alt": alt,
                         "avail": avail, "risk": score.round(1)})

# ---- 월 단위 재고 백테스트(가상 비용 지수) ----
HOLD, ORDER_COST, SHORT_MULT = 0.02, 50_000, 0.6   # [교육용가상자료] 월 보관비율(단가 대비), 발주 1회 비용(원), 결품 1개 손실(단가 배수)


ALT_PREMIUM, ALT_LT = 1.2, 3      # [교육용가상자료] 대체품 전환 시 단가 배수, 대체품 조달 리드타임(월)


def simulate(items, demand, fc, policy, z=1.65, months=None, sigma=None, eol=True):
    """품목별 월 재고 모의. demand/fc: (품목 x 월) 실제·예측 수요. policy: 'current'(최근 12개월 평균 + 3개월 여유 고정) 또는 'ai'(예측 + 안전재고).
    단종 품목은 원 품목 발주 불가. 'ai' 정책에서 eol=True이면 단종 품목 중 대체품이 있는 품목은 대체품(단가 x1.2, 리드타임 3개월)으로 전환 발주한다. 반환: 품목별 보관비·발주비·결품손실·불출량·평균재고."""
    n, T = demand.shape
    L = items["lead_time"].to_numpy(); price = items["price"].to_numpy().astype(float)
    can_order = (items["status"] != "단종").to_numpy()
    if policy == "ai" and eol:
        sw = (~can_order) & items["alt_available"].to_numpy()
        can_order = can_order | sw
        L = np.where(sw, ALT_LT, L); price = np.where(sw, price * ALT_PREMIUM, price)
    avg = demand[:, :1].repeat(T, 1) if months is None else months
    on = np.round(avg[:, 0] * 2.0)
    pipe = np.zeros((n, T + int(L.max()) + 1))                 # 도착 예정 물량(품목 x 월)
    for i in range(n):
        pipe[i, :L[i]] = np.round(avg[i, 0] * (1.0 if can_order[i] else 0.3))
    res = np.zeros((n, 6)); rows = np.arange(n)
    for t in range(T):
        on += pipe[:, t]
        issue = np.minimum(on, demand[:, t]); short = demand[:, t] - issue; on -= issue
        ip = on + pipe[:, t + 1:].sum(1)
        if policy == "current":
            target = avg[:, t] * (L + 3)
        else:
            s = sigma if sigma is not None else np.sqrt(np.maximum(fc[:, t], 1))
            target = fc[:, t] * (L + 1) + z * s * np.sqrt(L + 1)
        q = np.where(can_order & (ip < target), np.ceil(target - ip), 0)
        pipe[rows, t + L] += q
        res += np.c_[HOLD * price * on, ORDER_COST * (q > 0), SHORT_MULT * price * short, issue, on, short]
    out = pd.DataFrame(res, columns=["hold", "order", "short", "issued", "inv", "short_qty"], index=items["item_id"])
    out["inv"] /= T
    return out

items = make_items()
hist = make_history(items)
df = make_features(hist)

# 1) 백테스트(t=36~47, 12개월): 한 달 앞 예측은 t<36으로 학습한 앙상블 사용
m = fit_ens(df[df["t"] < 36])
bt = df[df["t"] >= 36].copy(); bt["fc"] = pred_ens(m, bt)
P = lambda c: bt.pivot(index="item_id", columns="t", values=c).loc[items["item_id"]].to_numpy()
D, F, A = P("demand"), P("fc"), P("item_level")
sig = np.maximum(0.25 * F.mean(1), np.sqrt(np.maximum(F.mean(1), 1)))


def summary(r, dem):
    cost = r[["hold", "order", "short"]].sum()
    return {"보관비": cost["hold"] / 1e6, "발주비": cost["order"] / 1e6, "결품손실": cost["short"] / 1e6, "총비용(백만)": cost.sum() / 1e6,
            "충족률": 1 - r["short_qty"].sum() / dem.sum(), "회전율": r["issued"].sum() / r["inv"].sum()}


runs = {"현행(평균+3개월 여유)": simulate(items, D, F, "current", months=A),
        "예측 정책(z=1.65)": simulate(items, D, F, "ai", 1.65, A, sig, eol=False),
        "예측+단종 대체품 전환": simulate(items, D, F, "ai", 1.65, A, sig, eol=True)}
tab = pd.DataFrame({k: summary(v, D) for k, v in runs.items()}).T
print("[12개월 재고 백테스트(가상 비용 지수, 백만 원)]")
print(tab.to_string(formatters={"충족률": "{:.1%}".format, "회전율": "{:.2f}".format}, float_format=lambda x: f"{x:.1f}"))
c0, c2 = tab.iloc[0]["총비용(백만)"], tab.iloc[2]["총비용(백만)"]
t0, t1 = tab.iloc[0]["회전율"], tab.iloc[1]["회전율"]
print(f"총비용 절감 {1 - c2 / c0:.1%}, 예측 정책만의 회전율 향상 {t1 / t0 - 1:.1%} (목표: 비용 15% 절감 또는 회전율 20% 향상 -> {'달성' if (1 - c2 / c0) >= 0.15 or t1 / t0 - 1 >= 0.2 else '미달'})")
act = items["status"] != "단종"
a0, a1 = runs["현행(평균+3개월 여유)"][act.to_numpy()], runs["예측 정책(z=1.65)"][act.to_numpy()]
print(f"효과 분해: 생산중·단종예고 품목 비용 {1 - a1[['hold', 'order', 'short']].sum().sum() / a0[['hold', 'order', 'short']].sum().sum():.1%} 절감(예측 정책), 나머지는 단종 품목 대체품 전환 효과")

# 민감도: 서비스 수준(z)과 돌발 수요(최근 6개월 수요 1.5배, 가상 전시 전환)
D2 = D.copy(); D2[:, 6:] = np.round(D2[:, 6:] * 1.5)
sens = []
for z in (1.28, 1.65, 2.05):
    for nm, dd in (("평시", D), ("돌발 1.5배", D2)):
        r = summary(simulate(items, dd, F, "ai", z, A, sig), dd)
        sens.append({"z": z, "수요": nm, "총비용": r["총비용(백만)"], "충족률": r["충족률"]})
print("[민감도] " + " / ".join(f"z={s['z']} {s['수요']}: {s['총비용']:.0f}백만, 충족 {s['충족률']:.1%}" for s in sens))


# 2) 다음 달 발주 추천: 예산 제약 LP (위험 지수 가중 부족분 충족 최대화)
class Audit:
    def __init__(self):
        self.log = []

    def add(self, **body):
        prev = self.log[-1]["hash"] if self.log else "GENESIS"
        body = {"seq": len(self.log), **body}
        self.log.append({**body, "prev": prev, "hash": hashlib.sha256((prev + json.dumps(body, sort_keys=True, ensure_ascii=False)).encode()).hexdigest()})

    def verify(self):
        prev = "GENESIS"
        for e in self.log:
            b = {k: v for k, v in e.items() if k not in ("prev", "hash")}
            if e["prev"] != prev or hashlib.sha256((prev + json.dumps(b, sort_keys=True, ensure_ascii=False)).encode()).hexdigest() != e["hash"]:
                return e["seq"]
            prev = e["hash"]
        return None


stock = make_stock(items, hist)
X = next_month(hist, items)
mu = pred_ens(fit_ens(df), X)
sg = np.maximum(0.25 * mu, np.sqrt(np.maximum(mu, 1)))
alt_sw = (items["status"] == "단종") & items["alt_available"]
it2 = items.assign(lead_time=np.where(alt_sw, ALT_LT, items["lead_time"]), price=np.where(alt_sw, items["price"] * ALT_PREMIUM, items["price"]))
pol = stock_policy(it2, mu, sg)
risk = risk_index(items, pol, stock)["risk"].to_numpy()
ip = (stock["on_hand"] + stock["on_order"]).to_numpy()
orderable = (items["status"] != "단종").to_numpy() | alt_sw.to_numpy()
gap = np.where(orderable, np.maximum(0, np.ceil(pol["rop"].to_numpy() - ip)), 0)
price = it2["price"].to_numpy().astype(float)
BUDGET = 0.6 * float((gap * price).sum())
idx = np.where(gap > 0)[0]
lp = linprog(c=-(risk[idx] + 10) / gap[idx], A_ub=[price[idx]], b_ub=[BUDGET], bounds=[(0, g) for g in gap[idx]], method="highs")
q = np.zeros(len(items)); q[idx] = np.floor(lp.x)
print(f"\n[다음 달 발주 추천 LP] 부족 품목 {len(idx)}개, 예산 {BUDGET / 1e6:.1f}백만(부족분 전액의 60%, 가상) 사용 {np.dot(q, price) / 1e6:.1f}백만, 추천 {int((q > 0).sum())}건")

audit, rows = Audit(), []
for k, i in enumerate(np.where(q > 0)[0]):
    iid = items.loc[i, "item_id"]; kind = "대체품 전환" if alt_sw[i] else "정규 발주"
    audit.add(t="추천", item=iid, kind=kind, qty=int(q[i]), risk=float(risk[i]), status="결재대기", by="AI")
    # 가상 정비장교 판단(시뮬레이션): 근거 확인 후 승인·조정·반려, 일부는 시한 내 무응답(보류=발주하지 않음)
    if kind == "대체품 전환":
        dec, final, why = "반려", 0, "대체품 호환성 시험 결과 확인 후 재상신"
    elif q[i] * price[i] > 0.15 * BUDGET:
        dec, final, why = "조정", int(q[i] * 0.7), "고액 발주 분할(1차 70%), 다음 달 재검토"
    elif k % 6 == 5:
        dec, final, why = "보류", 0, "결재 시한 초과 -> 발주하지 않음"
    else:
        dec, final, why = "승인", int(q[i]), "위험 지수·재고 근거 확인"
    audit.add(t="결재", item=iid, decision=dec, final_qty=final, reason=why, by="정비장교(가상)")
    rows.append({"품목": iid, "구분": kind, "위험": risk[i], "추천": int(q[i]), "결정": dec, "확정": final, "사유": why})
res = pd.DataFrame(rows).sort_values("위험", ascending=False)
print(res.head(10).to_string(index=False, float_format=lambda x: f"{x:.0f}"))
print("결정 집계:", ", ".join(f"{k} {v}건" for k, v in res["결정"].value_counts().items()))
auto = [e for e in audit.log if e["t"] == "추천" and e["status"] != "결재대기"]
print(f"자동 확정 {len(auto)}건, 결재 사유 기록률 {res['사유'].astype(bool).mean():.0%}, 감사 로그 {len(audit.log)}건 무결성: {'정상' if audit.verify() is None else '변조'}")
j = next(e["seq"] for e in audit.log if e["t"] == "결재" and e["decision"] != "승인")
old = audit.log[j]["decision"]; audit.log[j]["decision"] = "승인"
print(f"변조 시험({j}번 {old} 기록을 승인으로 수정) -> {audit.verify()}번 기록에서 변조 검출")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
tab[["보관비", "발주비", "결품손실"]].plot.bar(stacked=True, ax=ax[0], rot=0)
ax[0].set_xticklabels(["current", "forecast", "forecast+EOL alt"]); ax[0].set_title("12-month cost (million, synthetic)"); ax[0].legend(["holding", "ordering", "shortage"])
sd = pd.DataFrame(sens)
for nm, g in sd.groupby("수요"):
    ax[1].plot(g["z"], g["총비용"], "o-", label="normal" if nm == "평시" else "surge x1.5")
ax[1].set_xlabel("service level z"); ax[1].set_ylabel("total cost (million)"); ax[1].set_title("Sensitivity"); ax[1].legend()
plt.tight_layout(); plt.show()

실행 결과 예시

[12개월 재고 백테스트(가상 비용 지수, 백만 원)]
                보관비  발주비  결품손실  총비용(백만)   충족률  회전율
현행(평균+3개월 여유)  81.6 21.7 234.3    337.6 86.9% 5.42
예측 정책(z=1.65)  77.2 12.0 235.2    324.4 86.8% 6.09
예측+단종 대체품 전환  108.9 12.7   8.3    129.8 98.6% 5.65
총비용 절감 61.5%, 예측 정책만의 회전율 향상 12.2% (목표: 비용 15% 절감 또는 회전율 20% 향상 -> 달성)
효과 분해: 생산중·단종예고 품목 비용 13.0% 절감(예측 정책), 나머지는 단종 품목 대체품 전환 효과
[민감도] z=1.28 평시: 125백만, 충족 98.5% / z=1.28 돌발 1.5배: 230백만, 충족 90.4% / z=1.65 평시: 130백만, 충족 98.6% / z=1.65 돌발 1.5배: 213백만, 충족 91.6% / z=2.05 평시: 136백만, 충족 98.7% / z=2.05 돌발 1.5배: 195백만, 충족 92.8%

[다음 달 발주 추천 LP] 부족 품목 18개, 예산 93.6백만(부족분 전액의 60%, 가상) 사용 93.5백만, 추천 16건
  품목    구분  위험  추천 결정  확정                         사유
P004 정규 발주  88  45 조정  31 고액 발주 분할(1차 70%), 다음 달 재검토
P019 정규 발주  88  62 승인  62             위험 지수·재고 근거 확인
P030 정규 발주  88 491 승인 491             위험 지수·재고 근거 확인
P025 정규 발주  88 116 조정  81 고액 발주 분할(1차 70%), 다음 달 재검토
P038 정규 발주  88  57 승인  57             위험 지수·재고 근거 확인
P027 정규 발주  87  29 조정  20 고액 발주 분할(1차 70%), 다음 달 재검토
P032 정규 발주  73  26 승인  26             위험 지수·재고 근거 확인
P006 정규 발주  64  23 승인  23             위험 지수·재고 근거 확인
P012 정규 발주  61  45 승인  45             위험 지수·재고 근거 확인
P037 정규 발주  41  14 승인  14             위험 지수·재고 근거 확인
결정 집계: 승인 11건, 조정 3건, 보류 2건
자동 확정 0건, 결재 사유 기록률 100%, 감사 로그 32건 무결성: 정상
변조 시험(1번 조정 기록을 승인으로 수정) -> 1번 기록에서 변조 검출

미션 8. 전시 작전 템포 증가 모의 부품 수요 300% 폭증 시뮬레이션 부하 테스트

  • 난이도: 고급 / 예상 소요시간: 3.0시간
  • 코드 설명: 전시 작전 템포 증가로 부품 수요가 평시의 300%로 폭증하는 3개월 상황을 가상으로 모의해, 품목별 결품 경보(다음 달 예측 수요 > 가용재고)의 정확도를 반복 시험하고, 품목 수 50배 부하에서 재고 대시보드 갱신 시간(2초 이내)을 측정한다. 트리 모델이 학습 범위 밖 수요를 외삽하지 못하는 한계를 작전 템포 계수 보정으로 보완한다.
"""전시 작전 템포 증가로 부품 수요가 평시의 300%로 폭증하는 3개월 상황을 가상으로 모의해, 품목별 결품 경보(다음 달 예측 수요 > 가용재고)의 정확도를 반복 시험하고, 품목 수 50배 부하에서 재고 대시보드 갱신 시간(2초 이내)을 측정한다. 트리 모델이 학습 범위 밖 수요를 외삽하지 못하는 한계를 작전 템포 계수 보정으로 보완한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from sklearn.ensemble import HistGradientBoostingRegressor, RandomForestRegressor

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 야전 정비창 수리부속 청구/출고·장비 가동률 데이터 생성기
# 실제 부대·품목·재고·단가 정보와 무관하다. 품목 40개 x 48개월(가상 2022-01 ~ 2025-12) 월별 자료를 만든다.
# ---------------------------------------------------------------------------
FSC = {"2520": "power_train", "2530": "brake_steer", "2910": "fuel", "2930": "cooling",
       "2940": "filter", "5330": "gasket", "6140": "battery", "5961": "semiconductor"}
EXERCISE = {3, 8, 10}      # 가상 대규모 훈련 월(훈련 계획으로 미리 알 수 있음)
N_MONTHS = 48


def season(fsc, month):
    """품목군별 계절 계수: 축전지는 겨울, 냉각계통은 여름에 수요가 늘어난다."""
    if fsc == "6140":
        return 1.6 if month in (12, 1, 2) else 0.9
    if fsc == "2930":
        return 1.5 if month in (6, 7, 8) else 0.9
    return 1 + 0.1 * np.sin(2 * np.pi * month / 12)


def make_items(n=40, seed=2):
    """품목 마스터: NSN(가상), 품목군, 장비 유형, 기준 수요, 단가, 조달 리드타임, 공급업체 수, 생산(단종) 상태."""
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n):
        f = str(rng.choice(list(FSC)))
        status = str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1]))
        rows.append({"item_id": f"P{i + 1:03d}", "nsn": f"{f}-00-{rng.integers(100, 1000)}-{rng.integers(1000, 10000)}",
                     "fsc": f, "group": FSC[f], "equip_type": str(rng.choice(["A", "B", "C"])),
                     "base": float(np.round(rng.lognormal(2.7, 0.8), 1)),
                     "price": int(rng.lognormal(12, 1.0)) // 100 * 100,
                     "lead_time": int(rng.integers(1, 7) if status == "생산중" else rng.integers(4, 13)),
                     "suppliers": int(rng.integers(1, 5)) if status == "생산중" else 1,
                     "status": status, "alt_available": bool(rng.random() < (0.6 if status == "생산중" else 0.3))})
    return pd.DataFrame(rows)


def make_history(items, months=N_MONTHS, seed=3):
    """월별 청구 수요(개), 긴급 불출 건수, 장비 유형별 운용시간·가동률."""
    rng = np.random.default_rng(seed)
    t = np.arange(months)
    moy = t % 12 + 1
    ex = np.isin(moy, list(EXERCISE))
    hours = {e: (100 + 20 * rng.random(months)) * np.where(ex, 1.5, 1.0) for e in "ABC"}
    avail = {e: np.clip(0.93 - 0.12 * (hours[e] - 100) / 80 + rng.normal(0, 0.015, months), 0.6, 0.99) for e in "ABC"}
    rows = []
    for it in items.itertuples():
        trend = 1 + rng.uniform(-0.1, 0.3) * t / months
        s = np.array([season(it.fsc, m) for m in moy])
        lam = it.base * s * trend * hours[it.equip_type] / 110
        d = rng.poisson(lam)
        rows.append(pd.DataFrame({"item_id": it.item_id, "t": t, "ym": [f"{2022 + k // 12}-{k % 12 + 1:02d}" for k in t],
                                  "month": moy, "exercise": ex.astype(int), "op_hours": hours[it.equip_type].round(1),
                                  "availability": avail[it.equip_type].round(3), "demand": d,
                                  "emergency": rng.binomial(d, 0.05)}))
    return pd.concat(rows, ignore_index=True)


def make_stock(items, hist, seed=4):
    """현재(가상 2026-01) 재고: 보유량과 미입고 발주량(조달 중인 물량). 단종 품목은 추가 발주가 어려워 미입고량이 적다."""
    rng = np.random.default_rng(seed)
    recent = hist[hist["t"] >= hist["t"].max() - 5].groupby("item_id")["demand"].mean()
    st = items[["item_id", "lead_time", "status"]].copy()
    st["on_hand"] = [int(round(recent[i] * rng.uniform(0.5, 3.0))) for i in st["item_id"]]
    cover = np.where(st["status"] == "생산중", rng.uniform(0.6, 1.3, len(st)), rng.uniform(0.0, 0.6, len(st)))
    st["on_order"] = [int(round(recent[i] * L * c)) for i, L, c in zip(st["item_id"], st["lead_time"], cover)]
    return st[["item_id", "on_hand", "on_order"]]
# ---------------------------------------------------------------------------

# ---- 수요 예측 공통 함수(파생 피처 + 머신러닝 모델) ----
FEATS = ["lag1", "lag2", "lag3", "lag12", "ma3", "ma6", "std6", "season_idx", "month_sin", "month_cos",
         "exercise", "op_hours", "op_ratio", "item_level", "trend_t"]


def make_features(hist):
    """품목별 과거 수요로 파생 피처를 만든다(모두 예측 시점 이전 정보만 사용)."""
    df = hist.sort_values(["item_id", "t"]).reset_index(drop=True)
    g = df.groupby("item_id")["demand"]
    for k in (1, 2, 3, 12):
        df[f"lag{k}"] = g.shift(k)
    df["ma3"] = g.transform(lambda s: s.shift(1).rolling(3).mean())
    df["ma6"] = g.transform(lambda s: s.shift(1).rolling(6).mean())
    df["std6"] = g.transform(lambda s: s.shift(1).rolling(6).std())
    df["item_level"] = g.transform(lambda s: s.shift(1).rolling(12).mean())
    df["season_idx"] = df["lag12"] / (g.transform(lambda s: s.shift(12).rolling(12).mean()) + 1)
    df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
    df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)
    df["op_ratio"] = df["op_hours"] / df.groupby("item_id")["op_hours"].transform(lambda s: s.shift(1).rolling(12).mean())
    df["trend_t"] = df["t"]
    return df.dropna().reset_index(drop=True)


def new_model():
    """Prophet/LSTM 대신 사용하는 포아송 손실 그래디언트 부스팅(개수 자료에 적합)."""
    return HistGradientBoostingRegressor(loss="poisson", max_iter=300, learning_rate=0.05, max_leaf_nodes=15,
                                         min_samples_leaf=20, random_state=0)


def mape(y, p, min_y=1):
    y, p = np.asarray(y, float), np.asarray(p, float)
    m = y >= min_y
    return float(np.mean(np.abs(y[m] - p[m]) / y[m]))

# ---- 앙상블 수요 예측·안전재고·결품(단종) 위험 지수 공통 함수 ----
STATUS_W = {"생산중": 0.0, "단종예고": 0.5, "단종": 1.0}     # [교육용가상자료] 단종 상태 가중치


def fit_ens(tr):
    """포아송 그래디언트 부스팅(70%) + 랜덤포레스트(30%) 앙상블."""
    h = new_model().fit(tr[FEATS], tr["demand"])
    r = RandomForestRegressor(n_estimators=150, min_samples_leaf=5, random_state=0, n_jobs=1).fit(tr[FEATS], tr["demand"])
    return h, r


def pred_ens(m, X, w=0.7):
    return np.clip(w * m[0].predict(X[FEATS]) + (1 - w) * m[1].predict(X[FEATS]), 0, None)


def next_month(hist, items, t_new=48, op_mult=1.0, ex=0):
    """다음 달(t_new) 예측용 행: 작년 같은 달 운용시간 x 배수, 훈련 여부 지정."""
    last = hist[hist["t"] == t_new - 1].copy()
    prev = hist[hist["t"] == t_new - 12].set_index("item_id")["op_hours"]
    m = (t_new % 12) + 1
    f = last.assign(t=t_new, ym=f"{2022 + t_new // 12}-{m:02d}", month=m, demand=0, exercise=ex,
                    op_hours=last["item_id"].map(prev).to_numpy() * op_mult)
    X = make_features(pd.concat([hist, f], ignore_index=True))
    return X[X["t"] == t_new].set_index("item_id").loc[items["item_id"]].reset_index()


def stock_policy(items, mu, sigma, z=1.65, review=1):
    """안전재고 = z x 월 예측오차 표준편차 x sqrt(리드타임+검토주기), 재주문점 = 월 수요 x (리드타임+검토주기) + 안전재고."""
    L = items["lead_time"].to_numpy() + review
    ss = z * sigma * np.sqrt(L)
    return pd.DataFrame({"item_id": items["item_id"], "mu": mu, "sigma": sigma, "L": L,
                         "ss": np.ceil(ss), "rop": np.ceil(mu * L + ss)})


def risk_index(items, pol, stock):
    """결품 위험 지수(0~100): 리드타임 동안 수요가 가용 재고를 넘을 확률(정규 근사)을 기본으로,
    단종 상태·단일 공급업체·대체품 부재에 가중을 준다(가상 가중치)."""
    st = stock.set_index("item_id").loc[items["item_id"]]
    avail = (st["on_hand"] + st["on_order"]).to_numpy()
    p = 1 - norm.cdf((avail - pol["mu"].to_numpy() * pol["L"].to_numpy()) / (pol["sigma"].to_numpy() * np.sqrt(pol["L"].to_numpy()) + 1e-6))
    sw = items["status"].map(STATUS_W).to_numpy()
    sup = (items["suppliers"].to_numpy() == 1).astype(float)
    alt = (~items["alt_available"].to_numpy()).astype(float)
    score = 100 * np.clip(0.55 * p + 0.25 * sw + 0.10 * sup + 0.10 * alt, 0, 1)
    return pd.DataFrame({"item_id": items["item_id"], "p_short": p, "status_w": sw, "single_sup": sup, "no_alt": alt,
                         "avail": avail, "risk": score.round(1)})

SURGE = 3.0                                                   # 계획상 작전 템포 수요 배수(300%)
GROUP_MULT = {"power_train": 3.6, "brake_steer": 3.3, "fuel": 3.0, "cooling": 2.8,   # [교육용가상자료] 실제 배수(계획에는 미반영)
              "filter": 2.7, "gasket": 3.2, "battery": 2.6, "semiconductor": 3.1}
REPS, MONTHS = 30, 3

items = make_items()
hist = make_history(items)
df = make_features(hist)
stock = make_stock(items, hist)
ens = fit_ens(df)
X = next_month(hist, items, op_mult=SURGE, ex=1)
raw = pred_ens(ens, X)                                         # 운용시간 3배를 넣은 모델 예측(학습 범위 밖)
base = pred_ens(ens, next_month(hist, items))                  # 평시 예측
adj = base * SURGE                                             # 작전 템포 계수 보정 예측
print(f"[폭증 예측] 평시 합계 {base.sum():.0f}개, 모델 직접 예측 {raw.sum():.0f}개({raw.sum() / base.sum():.1f}배), 템포 계수 보정 {adj.sum():.0f}개({SURGE:.0f}배)")
print("-> 트리 모델은 학습한 운용시간 범위를 넘으면 예측이 더 늘지 않는다(외삽 불가). 계획 배수로 보정한다.")


def run(fc, seed, adaptive=False):
    """3개월 모의: 매월 경보(예측 수요 > 월초 재고) vs 실제 결품(실제 수요 > 월초 재고).
    adaptive=True이면 지난달 품목군별 실제/평시 비율로 다음 달 배수를 갱신한다(실적 반영)."""
    rng = np.random.default_rng(seed)
    true_mult = items["group"].map(GROUP_MULT).to_numpy() * rng.lognormal(0, 0.12, len(items))
    on = stock["on_hand"].to_numpy().astype(float)
    arrive = stock["on_order"].to_numpy() / np.maximum(items["lead_time"].to_numpy(), 1)
    al, tr, bd = [], [], []
    grp = items["group"].to_numpy(); f = fc.copy()
    for _ in range(MONTHS):
        actual = rng.poisson(base * true_mult)
        al.append(f > on); tr.append(actual > on); bd.append(np.abs(f - on) <= 1.65 * np.sqrt(np.maximum(f, 1)))
        if adaptive:
            ratio = pd.Series(actual).groupby(grp).sum() / pd.Series(base).groupby(grp).sum()
            f = base * pd.Series(grp).map(ratio).to_numpy()
        on = np.maximum(on - actual, 0) + arrive
    return np.array(al), np.array(tr), np.array(bd)


rows = []
for nm, fc, ad in (("평시 예측(보정 없음)", base, False), ("모델 직접 예측", raw, False), ("템포 계수 보정", adj, False),
                   ("템포 보정+월별 실적 반영", adj, True)):
    acc, rec, n_al = [], [], []
    for s in range(REPS):
        a, t, bd = run(fc, s, ad)
        acc.append((a == t).mean()); rec.append((a & t).sum() / max(t.sum(), 1)); n_al.append(a.sum())
    if ad:
        clear = [(lambda a, t, bd: ((a == t)[~bd].mean(), bd.mean()))(*run(fc, s, ad)) for s in range(REPS)]
    rows.append({"예측 방식": nm, "경보 정확도": np.mean(acc), "최저 정확도": np.min(acc), "결품 탐지율": np.mean(rec), "월평균 경보": np.mean(n_al) / MONTHS})
res = pd.DataFrame(rows)
print(f"\n[경보 시험: 품목 {len(items)}개 x {MONTHS}개월 x 반복 {REPS}회]")
print(res.to_string(index=False, formatters={"경보 정확도": "{:.1%}".format, "최저 정확도": "{:.1%}".format, "결품 탐지율": "{:.1%}".format, "월평균 경보": "{:.1f}".format}))
acc = res.iloc[3]["경보 정확도"]
print(f"템포 보정+실적 반영 경보 정확도 {acc:.1%} (목표 95% 이상: {'달성' if acc >= 0.95 else '미달'}), 템포 보정만 {res.iloc[2]['경보 정확도']:.1%}")
print(f"오경보 원인: 예측과 재고 차이가 수요 변동 폭(1.65 x 표준편차) 이내인 경계 구간이 판정의 {np.mean([c[1] for c in clear]):.0%}를 차지한다. "
      f"경계 구간을 '주의(검토 필요)'로 따로 표시하면 나머지 판정 정확도는 {np.mean([c[0] for c in clear]):.1%}이다.")


# 부하 시험: 품목 50배(2,000개) 대시보드 갱신 = 예측 + 위험 지수 + 집계
def refresh(Xb, it, st):
    mu = pred_ens(ens, Xb) * SURGE
    pol = stock_policy(it, mu, np.sqrt(np.maximum(mu, 1)))
    rk = risk_index(it, pol, st)
    return rk.assign(group=it["group"].to_numpy()).groupby("group")["risk"].agg(["mean", "max"]), rk


K = 50
Xb = pd.concat([X] * K, ignore_index=True)
itb = pd.concat([items.assign(item_id=items["item_id"] + f"_{k}") for k in range(K)], ignore_index=True)
stb = pd.concat([stock.assign(item_id=stock["item_id"] + f"_{k}") for k in range(K)], ignore_index=True)
times = []
for _ in range(3):
    t0 = time.perf_counter(); agg, rk = refresh(Xb, itb, stb); times.append(time.perf_counter() - t0)
t1 = time.perf_counter(); refresh(X, items, stock); t_small = time.perf_counter() - t1
print(f"\n[대시보드 부하 시험] 품목 {len(itb)}개 갱신 평균 {np.mean(times):.2f}s / 최대 {max(times):.2f}s, 품목 40개 {t_small:.3f}s (목표 2초 이내: {'달성' if max(times) <= 2 else '미달'})")
print("폭증 시 품목군별 최대 위험 지수:", ", ".join(f"{g} {r['max']:.0f}" for g, r in agg.iterrows()))
print("경보는 조달 검토 요청이며 긴급 조달·대체품 전환은 정비장교 승인 후 진행한다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].bar(["normal", "raw", "tempo", "tempo+AAR"], res["경보 정확도"], color=["gray", "tab:orange", "tab:green", "tab:blue"])
ax[0].axhline(0.95, c="r", ls="--"); ax[0].set_ylim(0.5, 1); ax[0].set_title("Alert accuracy under 300% demand surge")
ax[1].scatter(base, raw, s=12, label="model raw"); ax[1].scatter(base, adj, s=12, label="tempo adjusted")
ax[1].plot([0, base.max()], [0, base.max() * 3], "k--", lw=1, label="x3 line"); ax[1].set_xlabel("normal forecast"); ax[1].set_ylabel("surge forecast")
ax[1].set_title("Tree model cannot extrapolate"); ax[1].legend()
plt.tight_layout(); plt.show()

실행 결과 예시

[폭증 예측] 평시 합계 948개, 모델 직접 예측 1193개(1.3배), 템포 계수 보정 2843개(3배)
-> 트리 모델은 학습한 운용시간 범위를 넘으면 예측이 더 늘지 않는다(외삽 불가). 계획 배수로 보정한다.

[경보 시험: 품목 40개 x 3개월 x 반복 30회]
         예측 방식 경보 정확도 최저 정확도 결품 탐지율 월평균 경보
  평시 예측(보정 없음)  45.9%  40.8%  40.3%   14.6
      모델 직접 예측  65.0%  61.7%  61.5%   22.3
      템포 계수 보정  93.6%  88.3%  98.5%   37.7
템포 보정+월별 실적 반영  93.8%  88.3%  98.4%   37.5
템포 보정+실적 반영 경보 정확도 93.8% (목표 95% 이상: 미달), 템포 보정만 93.6%
오경보 원인: 예측과 재고 차이가 수요 변동 폭(1.65 x 표준편차) 이내인 경계 구간이 판정의 19%를 차지한다. 경계 구간을 '주의(검토 필요)'로 따로 표시하면 나머지 판정 정확도는 99.6%이다.

[대시보드 부하 시험] 품목 2000개 갱신 평균 0.05s / 최대 0.05s, 품목 40개 0.024s (목표 2초 이내: 달성)
폭증 시 품목군별 최대 위험 지수: battery 100, brake_steer 88, cooling 65, filter 88, fuel 88, gasket 90, power_train 90, semiconductor 78
경보는 조달 검토 요청이며 긴급 조달·대체품 전환은 정비장교 승인 후 진행한다.

미션 9. AAR 기반 부품 과다/과소 재고 예측 오차 분석 및 안전재고 파라미터 재보정

  • 난이도: 고급 / 예상 소요시간: 3.0시간
  • 코드 설명: AAR 관점에서 지난 12개월 재고 운영의 과다·과소 재고와 수요 예측 오차를 품목군·생산 상태·ABC 등급별로 분석해 원인을 찾고, 앞 6개월 자료로 ABC 등급별 안전재고 계수(z)와 품목군별 예측 편향 보정값을 재보정한 뒤 뒤 6개월에서 검증하며, 비용 절감액·정비 대기시간 단축률·재고 회전율 3대 효과 지표와 민감도, 도입 로드맵을 산출한다.
"""AAR 관점에서 지난 12개월 재고 운영의 과다·과소 재고와 수요 예측 오차를 품목군·생산 상태·ABC 등급별로 분석해 원인을 찾고, 앞 6개월 자료로 ABC 등급별 안전재고 계수(z)와 품목군별 예측 편향 보정값을 재보정한 뒤 뒤 6개월에서 검증하며, 비용 절감액·정비 대기시간 단축률·재고 회전율 3대 효과 지표와 민감도, 도입 로드맵을 산출한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from sklearn.ensemble import HistGradientBoostingRegressor, RandomForestRegressor

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 야전 정비창 수리부속 청구/출고·장비 가동률 데이터 생성기
# 실제 부대·품목·재고·단가 정보와 무관하다. 품목 40개 x 48개월(가상 2022-01 ~ 2025-12) 월별 자료를 만든다.
# ---------------------------------------------------------------------------
FSC = {"2520": "power_train", "2530": "brake_steer", "2910": "fuel", "2930": "cooling",
       "2940": "filter", "5330": "gasket", "6140": "battery", "5961": "semiconductor"}
EXERCISE = {3, 8, 10}      # 가상 대규모 훈련 월(훈련 계획으로 미리 알 수 있음)
N_MONTHS = 48


def season(fsc, month):
    """품목군별 계절 계수: 축전지는 겨울, 냉각계통은 여름에 수요가 늘어난다."""
    if fsc == "6140":
        return 1.6 if month in (12, 1, 2) else 0.9
    if fsc == "2930":
        return 1.5 if month in (6, 7, 8) else 0.9
    return 1 + 0.1 * np.sin(2 * np.pi * month / 12)


def make_items(n=40, seed=2):
    """품목 마스터: NSN(가상), 품목군, 장비 유형, 기준 수요, 단가, 조달 리드타임, 공급업체 수, 생산(단종) 상태."""
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n):
        f = str(rng.choice(list(FSC)))
        status = str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1]))
        rows.append({"item_id": f"P{i + 1:03d}", "nsn": f"{f}-00-{rng.integers(100, 1000)}-{rng.integers(1000, 10000)}",
                     "fsc": f, "group": FSC[f], "equip_type": str(rng.choice(["A", "B", "C"])),
                     "base": float(np.round(rng.lognormal(2.7, 0.8), 1)),
                     "price": int(rng.lognormal(12, 1.0)) // 100 * 100,
                     "lead_time": int(rng.integers(1, 7) if status == "생산중" else rng.integers(4, 13)),
                     "suppliers": int(rng.integers(1, 5)) if status == "생산중" else 1,
                     "status": status, "alt_available": bool(rng.random() < (0.6 if status == "생산중" else 0.3))})
    return pd.DataFrame(rows)


def make_history(items, months=N_MONTHS, seed=3):
    """월별 청구 수요(개), 긴급 불출 건수, 장비 유형별 운용시간·가동률."""
    rng = np.random.default_rng(seed)
    t = np.arange(months)
    moy = t % 12 + 1
    ex = np.isin(moy, list(EXERCISE))
    hours = {e: (100 + 20 * rng.random(months)) * np.where(ex, 1.5, 1.0) for e in "ABC"}
    avail = {e: np.clip(0.93 - 0.12 * (hours[e] - 100) / 80 + rng.normal(0, 0.015, months), 0.6, 0.99) for e in "ABC"}
    rows = []
    for it in items.itertuples():
        trend = 1 + rng.uniform(-0.1, 0.3) * t / months
        s = np.array([season(it.fsc, m) for m in moy])
        lam = it.base * s * trend * hours[it.equip_type] / 110
        d = rng.poisson(lam)
        rows.append(pd.DataFrame({"item_id": it.item_id, "t": t, "ym": [f"{2022 + k // 12}-{k % 12 + 1:02d}" for k in t],
                                  "month": moy, "exercise": ex.astype(int), "op_hours": hours[it.equip_type].round(1),
                                  "availability": avail[it.equip_type].round(3), "demand": d,
                                  "emergency": rng.binomial(d, 0.05)}))
    return pd.concat(rows, ignore_index=True)


def make_stock(items, hist, seed=4):
    """현재(가상 2026-01) 재고: 보유량과 미입고 발주량(조달 중인 물량). 단종 품목은 추가 발주가 어려워 미입고량이 적다."""
    rng = np.random.default_rng(seed)
    recent = hist[hist["t"] >= hist["t"].max() - 5].groupby("item_id")["demand"].mean()
    st = items[["item_id", "lead_time", "status"]].copy()
    st["on_hand"] = [int(round(recent[i] * rng.uniform(0.5, 3.0))) for i in st["item_id"]]
    cover = np.where(st["status"] == "생산중", rng.uniform(0.6, 1.3, len(st)), rng.uniform(0.0, 0.6, len(st)))
    st["on_order"] = [int(round(recent[i] * L * c)) for i, L, c in zip(st["item_id"], st["lead_time"], cover)]
    return st[["item_id", "on_hand", "on_order"]]
# ---------------------------------------------------------------------------

# ---- 수요 예측 공통 함수(파생 피처 + 머신러닝 모델) ----
FEATS = ["lag1", "lag2", "lag3", "lag12", "ma3", "ma6", "std6", "season_idx", "month_sin", "month_cos",
         "exercise", "op_hours", "op_ratio", "item_level", "trend_t"]


def make_features(hist):
    """품목별 과거 수요로 파생 피처를 만든다(모두 예측 시점 이전 정보만 사용)."""
    df = hist.sort_values(["item_id", "t"]).reset_index(drop=True)
    g = df.groupby("item_id")["demand"]
    for k in (1, 2, 3, 12):
        df[f"lag{k}"] = g.shift(k)
    df["ma3"] = g.transform(lambda s: s.shift(1).rolling(3).mean())
    df["ma6"] = g.transform(lambda s: s.shift(1).rolling(6).mean())
    df["std6"] = g.transform(lambda s: s.shift(1).rolling(6).std())
    df["item_level"] = g.transform(lambda s: s.shift(1).rolling(12).mean())
    df["season_idx"] = df["lag12"] / (g.transform(lambda s: s.shift(12).rolling(12).mean()) + 1)
    df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
    df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)
    df["op_ratio"] = df["op_hours"] / df.groupby("item_id")["op_hours"].transform(lambda s: s.shift(1).rolling(12).mean())
    df["trend_t"] = df["t"]
    return df.dropna().reset_index(drop=True)


def new_model():
    """Prophet/LSTM 대신 사용하는 포아송 손실 그래디언트 부스팅(개수 자료에 적합)."""
    return HistGradientBoostingRegressor(loss="poisson", max_iter=300, learning_rate=0.05, max_leaf_nodes=15,
                                         min_samples_leaf=20, random_state=0)


def mape(y, p, min_y=1):
    y, p = np.asarray(y, float), np.asarray(p, float)
    m = y >= min_y
    return float(np.mean(np.abs(y[m] - p[m]) / y[m]))

# ---- 앙상블 수요 예측·안전재고·결품(단종) 위험 지수 공통 함수 ----
STATUS_W = {"생산중": 0.0, "단종예고": 0.5, "단종": 1.0}     # [교육용가상자료] 단종 상태 가중치


def fit_ens(tr):
    """포아송 그래디언트 부스팅(70%) + 랜덤포레스트(30%) 앙상블."""
    h = new_model().fit(tr[FEATS], tr["demand"])
    r = RandomForestRegressor(n_estimators=150, min_samples_leaf=5, random_state=0, n_jobs=1).fit(tr[FEATS], tr["demand"])
    return h, r


def pred_ens(m, X, w=0.7):
    return np.clip(w * m[0].predict(X[FEATS]) + (1 - w) * m[1].predict(X[FEATS]), 0, None)


def next_month(hist, items, t_new=48, op_mult=1.0, ex=0):
    """다음 달(t_new) 예측용 행: 작년 같은 달 운용시간 x 배수, 훈련 여부 지정."""
    last = hist[hist["t"] == t_new - 1].copy()
    prev = hist[hist["t"] == t_new - 12].set_index("item_id")["op_hours"]
    m = (t_new % 12) + 1
    f = last.assign(t=t_new, ym=f"{2022 + t_new // 12}-{m:02d}", month=m, demand=0, exercise=ex,
                    op_hours=last["item_id"].map(prev).to_numpy() * op_mult)
    X = make_features(pd.concat([hist, f], ignore_index=True))
    return X[X["t"] == t_new].set_index("item_id").loc[items["item_id"]].reset_index()


def stock_policy(items, mu, sigma, z=1.65, review=1):
    """안전재고 = z x 월 예측오차 표준편차 x sqrt(리드타임+검토주기), 재주문점 = 월 수요 x (리드타임+검토주기) + 안전재고."""
    L = items["lead_time"].to_numpy() + review
    ss = z * sigma * np.sqrt(L)
    return pd.DataFrame({"item_id": items["item_id"], "mu": mu, "sigma": sigma, "L": L,
                         "ss": np.ceil(ss), "rop": np.ceil(mu * L + ss)})


def risk_index(items, pol, stock):
    """결품 위험 지수(0~100): 리드타임 동안 수요가 가용 재고를 넘을 확률(정규 근사)을 기본으로,
    단종 상태·단일 공급업체·대체품 부재에 가중을 준다(가상 가중치)."""
    st = stock.set_index("item_id").loc[items["item_id"]]
    avail = (st["on_hand"] + st["on_order"]).to_numpy()
    p = 1 - norm.cdf((avail - pol["mu"].to_numpy() * pol["L"].to_numpy()) / (pol["sigma"].to_numpy() * np.sqrt(pol["L"].to_numpy()) + 1e-6))
    sw = items["status"].map(STATUS_W).to_numpy()
    sup = (items["suppliers"].to_numpy() == 1).astype(float)
    alt = (~items["alt_available"].to_numpy()).astype(float)
    score = 100 * np.clip(0.55 * p + 0.25 * sw + 0.10 * sup + 0.10 * alt, 0, 1)
    return pd.DataFrame({"item_id": items["item_id"], "p_short": p, "status_w": sw, "single_sup": sup, "no_alt": alt,
                         "avail": avail, "risk": score.round(1)})

# ---- 월 단위 재고 백테스트(가상 비용 지수) ----
HOLD, ORDER_COST, SHORT_MULT = 0.02, 50_000, 0.6   # [교육용가상자료] 월 보관비율(단가 대비), 발주 1회 비용(원), 결품 1개 손실(단가 배수)


ALT_PREMIUM, ALT_LT = 1.2, 3      # [교육용가상자료] 대체품 전환 시 단가 배수, 대체품 조달 리드타임(월)


def simulate(items, demand, fc, policy, z=1.65, months=None, sigma=None, eol=True):
    """품목별 월 재고 모의. demand/fc: (품목 x 월) 실제·예측 수요. policy: 'current'(최근 12개월 평균 + 3개월 여유 고정) 또는 'ai'(예측 + 안전재고).
    단종 품목은 원 품목 발주 불가. 'ai' 정책에서 eol=True이면 단종 품목 중 대체품이 있는 품목은 대체품(단가 x1.2, 리드타임 3개월)으로 전환 발주한다. 반환: 품목별 보관비·발주비·결품손실·불출량·평균재고."""
    n, T = demand.shape
    L = items["lead_time"].to_numpy(); price = items["price"].to_numpy().astype(float)
    can_order = (items["status"] != "단종").to_numpy()
    if policy == "ai" and eol:
        sw = (~can_order) & items["alt_available"].to_numpy()
        can_order = can_order | sw
        L = np.where(sw, ALT_LT, L); price = np.where(sw, price * ALT_PREMIUM, price)
    avg = demand[:, :1].repeat(T, 1) if months is None else months
    on = np.round(avg[:, 0] * 2.0)
    pipe = np.zeros((n, T + int(L.max()) + 1))                 # 도착 예정 물량(품목 x 월)
    for i in range(n):
        pipe[i, :L[i]] = np.round(avg[i, 0] * (1.0 if can_order[i] else 0.3))
    res = np.zeros((n, 6)); rows = np.arange(n)
    for t in range(T):
        on += pipe[:, t]
        issue = np.minimum(on, demand[:, t]); short = demand[:, t] - issue; on -= issue
        ip = on + pipe[:, t + 1:].sum(1)
        if policy == "current":
            target = avg[:, t] * (L + 3)
        else:
            s = sigma if sigma is not None else np.sqrt(np.maximum(fc[:, t], 1))
            target = fc[:, t] * (L + 1) + z * s * np.sqrt(L + 1)
        q = np.where(can_order & (ip < target), np.ceil(target - ip), 0)
        pipe[rows, t + L] += q
        res += np.c_[HOLD * price * on, ORDER_COST * (q > 0), SHORT_MULT * price * short, issue, on, short]
    out = pd.DataFrame(res, columns=["hold", "order", "short", "issued", "inv", "short_qty"], index=items["item_id"])
    out["inv"] /= T
    return out

WAIT_H = 36            # [교육용가상자료] 결품 1개당 장비 정비 대기시간(시간)

items = make_items()
hist = make_history(items)
df = make_features(hist)
m = fit_ens(df[df["t"] < 36])
bt = df[df["t"] >= 36].copy(); bt["fc"] = pred_ens(m, bt)
P = lambda c: bt.pivot(index="item_id", columns="t", values=c).loc[items["item_id"]].to_numpy()
D, F, A = P("demand"), P("fc"), P("item_level")
sig = np.maximum(0.25 * F.mean(1), np.sqrt(np.maximum(F.mean(1), 1)))
val = items["price"].to_numpy() * D.sum(1)
cum = np.cumsum(np.sort(val)[::-1]) / val.sum()
abc = pd.Series(np.select([val >= np.sort(val)[::-1][np.searchsorted(cum, 0.7)], val >= np.sort(val)[::-1][np.searchsorted(cum, 0.9)]], ["A", "B"], "C"), index=items["item_id"])

# 1) AAR 편차 분석: 예측 편향과 과다·과소 재고
an = items[["item_id", "group", "status"]].assign(abc=abc.values, bias=(F[:, :6].sum(1) - D[:, :6].sum(1)) / np.maximum(D[:, :6].sum(1), 1))
r0 = simulate(items, D[:, :6], F[:, :6], "ai", 1.65, A[:, :6], sig)
an["over"] = (r0["inv"].to_numpy() > 3 * np.maximum(D[:, :6].mean(1), 1))       # 평균 재고가 3개월 수요 초과 = 과다
an["under"] = r0["short_qty"].to_numpy() > 0
print("[AAR 편차 분석: 앞 6개월(t=36~41)]")
g = an.groupby("group").agg(예측편향=("bias", "mean"), 과다품목=("over", "sum"), 결품품목=("under", "sum"))
gt = g.T.astype(object)
gt.loc["예측편향"] = [f"{v:+.0%}" for v in g["예측편향"]]
gt.loc[["과다품목", "결품품목"]] = gt.loc[["과다품목", "결품품목"]].map(lambda v: f"{int(v)}")
print(gt.to_string())
print("ABC별 과다/결품 품목:", ", ".join(f"{k} {int(v['over'])}/{int(v['under'])}" for k, v in an.groupby("abc")[["over", "under"]].sum().iterrows()),
      "| 생산 상태별 결품 품목:", ", ".join(f"{k} {int(v)}" for k, v in an.groupby("status")["under"].sum().items()))
ex = bt["exercise"] == 1
hi = g["예측편향"].abs().sort_values(ascending=False).index[:2]
hi_txt = ", ".join(f"{h}({g.loc[h, '예측편향']:+.0%})" for h in hi)
print(f"훈련 월 예측 편향 {(bt.loc[ex, 'fc'].sum() / bt.loc[ex, 'demand'].sum() - 1):+.1%}, 평시 월 {(bt.loc[~ex, 'fc'].sum() / bt.loc[~ex, 'demand'].sum() - 1):+.1%}"
      f" -> 원인: {hi_txt} 품목군 예측 편향, 모든 등급에 같은 z=1.65 적용, 단종 품목 발주 불가")

# 2) 재보정(앞 6개월로 결정): ABC별 z 격자 탐색 + 품목군별 편향 보정
grp_bias = an.groupby("group")["bias"].median().clip(-0.3, 0.3)
corr = 1 / (1 + items["group"].map(grp_bias).to_numpy())
best = {}
for cls in "ABC":
    mask = (abc == cls).to_numpy()
    costs = {}
    for z in (0.84, 1.28, 1.65, 2.05, 2.33):
        zz = np.where(mask, z, 1.65)
        r = simulate(items, D[:, :6], F[:, :6] * corr[:, None], "ai", zz, A[:, :6], sig)
        costs[z] = r[["hold", "order", "short"]].to_numpy()[mask].sum()
    best[cls] = min(costs, key=costs.get)
z_new = abc.map(best).to_numpy()
print(f"\n[재보정 결과] ABC별 z: A {best['A']}, B {best['B']}, C {best['C']} (기존 1.65 일괄), 품목군 편향 보정 범위 {corr.min():.2f}~{corr.max():.2f}배")

# 3) 검증(뒤 6개월, t=42~47) 및 3대 효과 지표
def ev(fc, z):
    r = simulate(items, D[:, 6:], fc, "ai", z, A[:, 6:], sig)
    return r, r[["hold", "order", "short"]].sum().sum()


r_cur = simulate(items, D[:, 6:], F[:, 6:], "current", months=A[:, 6:])
r_old, c_old = ev(F[:, 6:], 1.65)
r_new, c_new = ev(F[:, 6:] * corr[:, None], z_new)
c_cur = r_cur[["hold", "order", "short"]].sum().sum()
turn = lambda r: float(r["issued"].sum() / r["inv"].sum())
over_amt = lambda r: float((np.maximum(r["inv"].to_numpy() - 3 * np.maximum(D[:, 6:].mean(1), 1), 0) * items["price"].to_numpy()).sum())
k = pd.DataFrame([["6개월 총비용(백만 원)", c_cur / 1e6, c_old / 1e6, c_new / 1e6],
                  ["정비 대기시간(시간)", r_cur["short_qty"].sum() * WAIT_H, r_old["short_qty"].sum() * WAIT_H, r_new["short_qty"].sum() * WAIT_H],
                  ["재고 회전율(회/6개월)", turn(r_cur), turn(r_old), turn(r_new)]],
                 columns=["지표", "현행 정책", "예측 정책(기존 파라미터)", "예측 정책(재보정)"]).set_index("지표")
k["현행 대비 개선율"] = 1 - k["예측 정책(재보정)"] / k["현행 정책"]
k.loc["재고 회전율(회/6개월)", "현행 대비 개선율"] *= -1          # 회전율은 클수록 좋음
print("\n[3대 정량 효과 지표: 뒤 6개월 검증, 가상 단가 기준]")
print(k.to_string(float_format=lambda x: f"{x:,.1f}", formatters={"현행 대비 개선율": "{:.1%}".format}))
print(f"재보정 효과(기존 파라미터 대비): 총비용 {1 - c_new / c_old:+.1%}, 연간 환산 절감액(현행 대비) {(c_cur - c_new) * 2 / 1e6:,.0f}백만 원(가상)")
sens = []
for hm in (0.5, 1.0, 1.5):
    for sm in (0.3, 0.6, 1.0):
        HOLD_, SM_ = HOLD * hm, sm
        cc = lambda r, p: (HOLD_ * p * r["inv"].to_numpy() * 6).sum() + r["order"].sum() + (SM_ * p * r["short_qty"].to_numpy()).sum()
        pr = items["price"].to_numpy()
        sens.append(f"보관x{hm}/결품x{sm}: {1 - cc(r_new, pr) / cc(r_cur, pr):.0%}")
print("[민감도: 현행 대비 비용 절감률]", ", ".join(sens))
print("도입 로드맵: 단기(1~3개월) 분기별 AAR로 z·편향 재보정 / 중기(4~9개월) 단종예고 품목 최종 구매·대체품 검증 절차 정착 / 장기(10개월~) 군수 정보체계 연계와 자동 재학습")
print("금액·시간은 가정 단가로 계산한 상대 비교값이며 예산 근거로 쓰지 않는다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
g["예측편향"].plot.bar(ax=ax[0], rot=30, color="tab:purple"); ax[0].axhline(0, c="k", lw=0.8); ax[0].set_title("Forecast bias by group (AAR, first 6 months)")
k.iloc[:, :3].div(k["현행 정책"], axis=0).T.plot.bar(ax=ax[1], rot=0); ax[1].set_xticklabels(["current", "forecast(old)", "forecast(recal)"])
ax[1].set_title("Effect indicators (current = 1.0)"); ax[1].legend(["cost", "wait hours", "turnover"])
plt.tight_layout(); plt.show()

실행 결과 예시

[AAR 편차 분석: 앞 6개월(t=36~41)]
group battery brake_steer cooling filter fuel gasket power_train semiconductor
예측편향     +11%        +18%     +4%    +1%  -6%    -5%         +1%           -3%
과다품목        1           1       1      2    0      0           0             0
결품품목        1           0       0      0    1      0           0             0
ABC별 과다/결품 품목: A 1/0, B 2/0, C 2/2 | 생산 상태별 결품 품목: 단종 1, 단종예고 0, 생산중 1
훈련 월 예측 편향 +1.5%, 평시 월 +0.1% -> 원인: brake_steer(+18%), battery(+11%) 품목군 예측 편향, 모든 등급에 같은 z=1.65 적용, 단종 품목 발주 불가

[재보정 결과] ABC별 z: A 0.84, B 0.84, C 1.28 (기존 1.65 일괄), 품목군 편향 보정 범위 0.89~1.07배

[3대 정량 효과 지표: 뒤 6개월 검증, 가상 단가 기준]
                 현행 정책  예측 정책(기존 파라미터)  예측 정책(재보정) 현행 대비 개선율
지표                                                          
6개월 총비용(백만 원)    117.7            56.7        55.3     53.0%
정비 대기시간(시간)   17,532.0         1,620.0     1,872.0     89.3%
재고 회전율(회/6개월)      3.1             3.0         3.1      1.2%
재보정 효과(기존 파라미터 대비): 총비용 +2.5%, 연간 환산 절감액(현행 대비) 125백만 원(가상)
[민감도: 현행 대비 비용 절감률] 보관x0.5/결품x0.3: 55%, 보관x0.5/결품x0.6: 69%, 보관x0.5/결품x1.0: 77%, 보관x1.0/결품x0.3: 38%, 보관x1.0/결품x0.6: 55%, 보관x1.0/결품x1.0: 67%, 보관x1.5/결품x0.3: 28%, 보관x1.5/결품x0.6: 45%, 보관x1.5/결품x1.0: 58%
도입 로드맵: 단기(1~3개월) 분기별 AAR로 z·편향 재보정 / 중기(4~9개월) 단종예고 품목 최종 구매·대체품 검증 절차 정착 / 장기(10개월~) 군수 정보체계 연계와 자동 재학습
금액·시간은 가정 단가로 계산한 상대 비교값이며 예산 근거로 쓰지 않는다.

미션 10. AI 기반 야전 정비 부품 수요예측 및 재고최적화 시스템 통합 패키지 완성

  • 난이도: 최고급 / 예상 소요시간: 3.5시간
  • 코드 설명: 수리부속 청구 이력 생성·정합성 점검, 수요 예측 모델 자동 재학습(기존·신규 모델 비교 후 교체), 안전재고·결품(단종) 위험 지수, 단종예고 품목 최종 구매량과 단종 품목 대체 검토, 예산 제약 발주 추천과 정비장교 결재 대기열을 한 번에 실행하는 통합 패키지를 완성하고, 자체 점검표·운영 가이드 목차·체계 연계 메시지 예시·브리핑 요약을 출력한다.
"""수리부속 청구 이력 생성·정합성 점검, 수요 예측 모델 자동 재학습(기존·신규 모델 비교 후 교체), 안전재고·결품(단종) 위험 지수, 단종예고 품목 최종 구매량과 단종 품목 대체 검토, 예산 제약 발주 추천과 정비장교 결재 대기열을 한 번에 실행하는 통합 패키지를 완성하고, 자체 점검표·운영 가이드 목차·체계 연계 메시지 예시·브리핑 요약을 출력한다."""
import json
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import linprog
from sklearn.ensemble import HistGradientBoostingRegressor, RandomForestRegressor
from sklearn.metrics import mean_absolute_error

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 야전 정비창 수리부속 청구/출고·장비 가동률 데이터 생성기
# 실제 부대·품목·재고·단가 정보와 무관하다. 품목 40개 x 48개월(가상 2022-01 ~ 2025-12) 월별 자료를 만든다.
# ---------------------------------------------------------------------------
FSC = {"2520": "power_train", "2530": "brake_steer", "2910": "fuel", "2930": "cooling",
       "2940": "filter", "5330": "gasket", "6140": "battery", "5961": "semiconductor"}
EXERCISE = {3, 8, 10}      # 가상 대규모 훈련 월(훈련 계획으로 미리 알 수 있음)
N_MONTHS = 48


def season(fsc, month):
    """품목군별 계절 계수: 축전지는 겨울, 냉각계통은 여름에 수요가 늘어난다."""
    if fsc == "6140":
        return 1.6 if month in (12, 1, 2) else 0.9
    if fsc == "2930":
        return 1.5 if month in (6, 7, 8) else 0.9
    return 1 + 0.1 * np.sin(2 * np.pi * month / 12)


def make_items(n=40, seed=2):
    """품목 마스터: NSN(가상), 품목군, 장비 유형, 기준 수요, 단가, 조달 리드타임, 공급업체 수, 생산(단종) 상태."""
    rng = np.random.default_rng(seed)
    rows = []
    for i in range(n):
        f = str(rng.choice(list(FSC)))
        status = str(rng.choice(["생산중", "단종예고", "단종"], p=[0.7, 0.2, 0.1]))
        rows.append({"item_id": f"P{i + 1:03d}", "nsn": f"{f}-00-{rng.integers(100, 1000)}-{rng.integers(1000, 10000)}",
                     "fsc": f, "group": FSC[f], "equip_type": str(rng.choice(["A", "B", "C"])),
                     "base": float(np.round(rng.lognormal(2.7, 0.8), 1)),
                     "price": int(rng.lognormal(12, 1.0)) // 100 * 100,
                     "lead_time": int(rng.integers(1, 7) if status == "생산중" else rng.integers(4, 13)),
                     "suppliers": int(rng.integers(1, 5)) if status == "생산중" else 1,
                     "status": status, "alt_available": bool(rng.random() < (0.6 if status == "생산중" else 0.3))})
    return pd.DataFrame(rows)


def make_history(items, months=N_MONTHS, seed=3):
    """월별 청구 수요(개), 긴급 불출 건수, 장비 유형별 운용시간·가동률."""
    rng = np.random.default_rng(seed)
    t = np.arange(months)
    moy = t % 12 + 1
    ex = np.isin(moy, list(EXERCISE))
    hours = {e: (100 + 20 * rng.random(months)) * np.where(ex, 1.5, 1.0) for e in "ABC"}
    avail = {e: np.clip(0.93 - 0.12 * (hours[e] - 100) / 80 + rng.normal(0, 0.015, months), 0.6, 0.99) for e in "ABC"}
    rows = []
    for it in items.itertuples():
        trend = 1 + rng.uniform(-0.1, 0.3) * t / months
        s = np.array([season(it.fsc, m) for m in moy])
        lam = it.base * s * trend * hours[it.equip_type] / 110
        d = rng.poisson(lam)
        rows.append(pd.DataFrame({"item_id": it.item_id, "t": t, "ym": [f"{2022 + k // 12}-{k % 12 + 1:02d}" for k in t],
                                  "month": moy, "exercise": ex.astype(int), "op_hours": hours[it.equip_type].round(1),
                                  "availability": avail[it.equip_type].round(3), "demand": d,
                                  "emergency": rng.binomial(d, 0.05)}))
    return pd.concat(rows, ignore_index=True)


def make_stock(items, hist, seed=4):
    """현재(가상 2026-01) 재고: 보유량과 미입고 발주량(조달 중인 물량). 단종 품목은 추가 발주가 어려워 미입고량이 적다."""
    rng = np.random.default_rng(seed)
    recent = hist[hist["t"] >= hist["t"].max() - 5].groupby("item_id")["demand"].mean()
    st = items[["item_id", "lead_time", "status"]].copy()
    st["on_hand"] = [int(round(recent[i] * rng.uniform(0.5, 3.0))) for i in st["item_id"]]
    cover = np.where(st["status"] == "생산중", rng.uniform(0.6, 1.3, len(st)), rng.uniform(0.0, 0.6, len(st)))
    st["on_order"] = [int(round(recent[i] * L * c)) for i, L, c in zip(st["item_id"], st["lead_time"], cover)]
    return st[["item_id", "on_hand", "on_order"]]
# ---------------------------------------------------------------------------

# ---- 수요 예측 공통 함수(파생 피처 + 머신러닝 모델) ----
FEATS = ["lag1", "lag2", "lag3", "lag12", "ma3", "ma6", "std6", "season_idx", "month_sin", "month_cos",
         "exercise", "op_hours", "op_ratio", "item_level", "trend_t"]


def make_features(hist):
    """품목별 과거 수요로 파생 피처를 만든다(모두 예측 시점 이전 정보만 사용)."""
    df = hist.sort_values(["item_id", "t"]).reset_index(drop=True)
    g = df.groupby("item_id")["demand"]
    for k in (1, 2, 3, 12):
        df[f"lag{k}"] = g.shift(k)
    df["ma3"] = g.transform(lambda s: s.shift(1).rolling(3).mean())
    df["ma6"] = g.transform(lambda s: s.shift(1).rolling(6).mean())
    df["std6"] = g.transform(lambda s: s.shift(1).rolling(6).std())
    df["item_level"] = g.transform(lambda s: s.shift(1).rolling(12).mean())
    df["season_idx"] = df["lag12"] / (g.transform(lambda s: s.shift(12).rolling(12).mean()) + 1)
    df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
    df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)
    df["op_ratio"] = df["op_hours"] / df.groupby("item_id")["op_hours"].transform(lambda s: s.shift(1).rolling(12).mean())
    df["trend_t"] = df["t"]
    return df.dropna().reset_index(drop=True)


def new_model():
    """Prophet/LSTM 대신 사용하는 포아송 손실 그래디언트 부스팅(개수 자료에 적합)."""
    return HistGradientBoostingRegressor(loss="poisson", max_iter=300, learning_rate=0.05, max_leaf_nodes=15,
                                         min_samples_leaf=20, random_state=0)


def mape(y, p, min_y=1):
    y, p = np.asarray(y, float), np.asarray(p, float)
    m = y >= min_y
    return float(np.mean(np.abs(y[m] - p[m]) / y[m]))

# ---- 앙상블 수요 예측·안전재고·결품(단종) 위험 지수 공통 함수 ----
STATUS_W = {"생산중": 0.0, "단종예고": 0.5, "단종": 1.0}     # [교육용가상자료] 단종 상태 가중치


def fit_ens(tr):
    """포아송 그래디언트 부스팅(70%) + 랜덤포레스트(30%) 앙상블."""
    h = new_model().fit(tr[FEATS], tr["demand"])
    r = RandomForestRegressor(n_estimators=150, min_samples_leaf=5, random_state=0, n_jobs=1).fit(tr[FEATS], tr["demand"])
    return h, r


def pred_ens(m, X, w=0.7):
    return np.clip(w * m[0].predict(X[FEATS]) + (1 - w) * m[1].predict(X[FEATS]), 0, None)


def next_month(hist, items, t_new=48, op_mult=1.0, ex=0):
    """다음 달(t_new) 예측용 행: 작년 같은 달 운용시간 x 배수, 훈련 여부 지정."""
    last = hist[hist["t"] == t_new - 1].copy()
    prev = hist[hist["t"] == t_new - 12].set_index("item_id")["op_hours"]
    m = (t_new % 12) + 1
    f = last.assign(t=t_new, ym=f"{2022 + t_new // 12}-{m:02d}", month=m, demand=0, exercise=ex,
                    op_hours=last["item_id"].map(prev).to_numpy() * op_mult)
    X = make_features(pd.concat([hist, f], ignore_index=True))
    return X[X["t"] == t_new].set_index("item_id").loc[items["item_id"]].reset_index()


def stock_policy(items, mu, sigma, z=1.65, review=1):
    """안전재고 = z x 월 예측오차 표준편차 x sqrt(리드타임+검토주기), 재주문점 = 월 수요 x (리드타임+검토주기) + 안전재고."""
    L = items["lead_time"].to_numpy() + review
    ss = z * sigma * np.sqrt(L)
    return pd.DataFrame({"item_id": items["item_id"], "mu": mu, "sigma": sigma, "L": L,
                         "ss": np.ceil(ss), "rop": np.ceil(mu * L + ss)})


def risk_index(items, pol, stock):
    """결품 위험 지수(0~100): 리드타임 동안 수요가 가용 재고를 넘을 확률(정규 근사)을 기본으로,
    단종 상태·단일 공급업체·대체품 부재에 가중을 준다(가상 가중치)."""
    st = stock.set_index("item_id").loc[items["item_id"]]
    avail = (st["on_hand"] + st["on_order"]).to_numpy()
    p = 1 - norm.cdf((avail - pol["mu"].to_numpy() * pol["L"].to_numpy()) / (pol["sigma"].to_numpy() * np.sqrt(pol["L"].to_numpy()) + 1e-6))
    sw = items["status"].map(STATUS_W).to_numpy()
    sup = (items["suppliers"].to_numpy() == 1).astype(float)
    alt = (~items["alt_available"].to_numpy()).astype(float)
    score = 100 * np.clip(0.55 * p + 0.25 * sw + 0.10 * sup + 0.10 * alt, 0, 1)
    return pd.DataFrame({"item_id": items["item_id"], "p_short": p, "status_w": sw, "single_sup": sup, "no_alt": alt,
                         "avail": avail, "risk": score.round(1)})

LTB_MONTHS = 24        # [교육용가상자료] 단종예고 품목 최종 구매(LTB) 대상 기간(개월)
t_all = time.perf_counter()
checks = []

# 1) 자료 수집·정합성
items = make_items()
hist = make_history(items)
stock = make_stock(items, hist)
ok = (hist["demand"] >= 0).all() and items["nsn"].str.match(r"^\d{4}-\d{2}-\d{3}-\d{4}$").all() and not hist.duplicated(["item_id", "t"]).any()
checks.append(("원천 자료 정합성(NSN 형식·음수·중복)", ok))
df = make_features(hist)
print(f"[1] 품목 {len(items)}개 x {N_MONTHS}개월 청구 이력, 학습 행 {len(df)}개, 정합성 {'통과' if ok else '실패'}")

# 2) 재학습 파이프라인: 운영 모델 v1(t<36) vs 재학습 v2(t<42), 최근 6개월(t=42~47) 검증 후 교체 판단
va = df[df["t"] >= 42]
big = va["item_id"].isin(items.loc[items["base"] >= 10, "item_id"]).to_numpy()
v1, v2 = fit_ens(df[df["t"] < 36]), fit_ens(df[df["t"] < 42])
p1, p2 = pred_ens(v1, va), pred_ens(v2, va)
e1, e2 = mape(va["demand"][big], p1[big]), mape(va["demand"][big], p2[big])
eb = mape(va["demand"][big], va["ma3"][big])
use = "v2" if e2 <= e1 else "v1"
print(f"[2] 검증 MAPE(수요 많은 품목): MA3 {eb:.1%}, 운영 v1 {e1:.1%}, 재학습 v2 {e2:.1%} -> {use} {'교체' if use == 'v2' else '유지'}")
checks += [("자동 재학습 실행·비교", True), ("운영 모델 MAPE 25% 이하", min(e1, e2) <= 0.25)]

# 3) 다음 달 예측(전체 자료로 최종 학습) -> 안전재고·위험 지수
final = fit_ens(df)
mu = pred_ens(final, next_month(hist, items))
sg = np.maximum(0.25 * mu, np.sqrt(np.maximum(mu, 1)))
pol = stock_policy(items, mu, sg)
rk = risk_index(items, pol, stock)
alert = rk["risk"] >= 60
print(f"[3] 다음 달 예측 수요 {mu.sum():.0f}개, 결품 위험 경보(60 이상) {int(alert.sum())}개 품목, 안전재고 합 {pol['ss'].sum():.0f}개")

# 4) 단종 대응: 단종예고 -> 최종 구매(LTB) 수량, 단종 -> 대체품·재생정비 검토
ip = (stock["on_hand"] + stock["on_order"]).to_numpy()
eol = items.assign(mu=mu, ip=ip, risk=rk["risk"].to_numpy())
ltb = eol[eol["status"] == "단종예고"].assign(ltb=lambda d: np.maximum(0, np.ceil(d["mu"] * LTB_MONTHS - d["ip"])))
dis = eol[eol["status"] == "단종"]
print(f"[4] 단종예고 {len(ltb)}개 품목 최종 구매 추천 합 {ltb['ltb'].sum():.0f}개(향후 {LTB_MONTHS}개월분), 단종 {len(dis)}개 품목 중 대체품 보유 {int(dis['alt_available'].sum())}개")
for _, r in ltb.sort_values("risk", ascending=False).head(3).iterrows():
    print(f"    {r['item_id']}({r['group']}) 월 {r['mu']:.1f}개, 가용 {r['ip']:.0f}개 -> 최종 구매 {r['ltb']:.0f}개")

# 5) 예산 제약 발주 추천(LP, 생산중 품목. 단종예고는 최종 구매로 별도 처리) -> 결재 대기열
gap = np.where(items["status"] == "생산중", np.maximum(0, np.ceil(pol["rop"].to_numpy() - ip)), 0)
price = items["price"].to_numpy().astype(float)
idx = np.where(gap > 0)[0]
budget = 0.6 * float((gap * price).sum())
lp = linprog(c=-(rk["risk"].to_numpy()[idx] + 10) / gap[idx], A_ub=[price[idx]], b_ub=[budget], bounds=[(0, g) for g in gap[idx]], method="highs")
q = np.zeros(len(items)); q[idx] = np.floor(lp.x)
queue = [{"id": f"Q{k + 1:02d}", "item": items.loc[i, "item_id"], "qty": int(q[i]), "risk": float(rk.loc[i, "risk"]), "status": "결재대기"} for k, i in enumerate(np.where(q > 0)[0])]
queue += [{"id": f"Q{len(queue) + k + 1:02d}", "item": r["item_id"], "qty": int(r["ltb"]), "risk": float(r["risk"]), "status": "결재대기", "type": "LTB"} for k, (_, r) in enumerate(ltb[ltb["ltb"] > 0].iterrows())]
auto = [x for x in queue if x["status"] != "결재대기"]
print(f"[5] 발주 추천 {int((q > 0).sum())}건(예산 {budget / 1e6:.1f}백만 중 {np.dot(q, price) / 1e6:.1f}백만) + 최종 구매 {int((ltb['ltb'] > 0).sum())}건 -> 정비장교 결재 대기 {len(queue)}건, 자동 확정 {len(auto)}건")
checks.append(("자동 확정 0건(HITL)", not auto))

# 6) 운영 가이드·연계 메시지
GUIDE = {"1. 목적·범위": "야전 정비 부품 수요 예측·재고 최적화·단종 대응", "2. 데이터 규격": f"품목 마스터 {items.shape[1]}개 필드, NSN 가상 형식",
         "3. 전처리": "NSN 정규화, 결측 월 계절성 보간, 긴급 불출 보존", "4. 예측 모델": f"포아송 GBM+RF 앙상블, MAPE {min(e1, e2):.1%}",
         "5. 안전재고·위험 지수": "z=1.65, 리드타임 결품 확률+단종·공급업체·대체품", "6. 단종 대응": f"LTB {LTB_MONTHS}개월분, 대체품 호환성 검증",
         "7. 결재선·감사": "자동 확정 금지, 승인·조정·반려·보류 사유 기록", "8. 재학습": "분기별 재학습, 검증 MAPE가 낮을 때만 교체",
         "9. 폭증 대응": "작전 템포 계수 보정, 경계 구간 '주의' 표시", "10. 연계 인터페이스": "ORDER_REQ 가상 메시지 규격"}
done = np.mean([bool(v) for v in GUIDE.values()])
checks.append(("운영 가이드 10개 장 완비", done == 1))
print("[6] 운영 가이드 목차:", " / ".join(GUIDE))
top = max(queue, key=lambda x: x["risk"]) if queue else {}
msg = {"msg_type": "ORDER_REQ", "version": "1.0", "item_id": top.get("item"), "nsn": items.set_index("item_id").loc[top.get("item"), "nsn"] if top else None,
       "qty": top.get("qty"), "risk_index": top.get("risk"), "model_ver": use, "approval": "PENDING", "source": "교육용가상자료"}
print("[7] 체계 연계 메시지 예시(가상 규격):", json.dumps(msg, ensure_ascii=False))
elapsed = time.perf_counter() - t_all
checks.append(("통합 실행 60초 이내", elapsed < 60))

print("\n[자체 점검표]")
for name, okk in checks:
    print(f"  [{'PASS' if okk else 'FAIL'}] {name}")
print(f"점검 {sum(o for _, o in checks)}/{len(checks)} 통과, 운영 가이드 완비율 {done:.0%}, 전체 실행 {elapsed:.1f}s")
print("\n[브리핑 요약]")
print(f"- 부품 수요 예측 MAPE {min(e1, e2):.1%}(이동평균 {eb:.1%})이며, 재학습 모델은 검증 후 교체한다({use} 사용).")
print(f"- 결품 위험 경보 {int(alert.sum())}개 품목 중 단종·단종예고가 {int((alert & (items['status'] != '생산중')).sum())}개로, 단종 대응이 재고 위험의 핵심이다.")
print(f"- 단종예고 {len(ltb)}개 품목의 최종 구매 {ltb['ltb'].sum():.0f}개와 발주 {int((q > 0).sum())}건을 결재 대기열에 올렸다.")
print("- 모든 발주·최종 구매·대체품 전환은 정비장교 결재 후 집행하며, 수치는 교육용 가상 자료 기준이다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].bar(["MA3", "model v1", "model v2"], [eb, e1, e2], color=["gray", "tab:blue", "tab:green"]); ax[0].axhline(0.25, c="r", ls="--")
ax[0].set_title("Validation MAPE (t=42-47)")
col = items["status"].map({"생산중": "tab:green", "단종예고": "tab:orange", "단종": "tab:red"})
ax[1].scatter(items["lead_time"], rk["risk"], c=col, s=20 + mu / 2); ax[1].axhline(60, c="r", ls="--")
ax[1].set_xlabel("lead time (months)"); ax[1].set_ylabel("risk index"); ax[1].set_title("Risk vs lead time (size=demand, red=discontinued)")
plt.tight_layout(); plt.show()

실행 결과 예시

[1] 품목 40개 x 48개월 청구 이력, 학습 행 1000개, 정합성 통과
[2] 검증 MAPE(수요 많은 품목): MA3 32.7%, 운영 v1 20.0%, 재학습 v2 20.9% -> v1 유지
[3] 다음 달 예측 수요 948개, 결품 위험 경보(60 이상) 13개 품목, 안전재고 합 1108개
[4] 단종예고 8개 품목 최종 구매 추천 합 2779개(향후 24개월분), 단종 3개 품목 중 대체품 보유 2개
    P004(filter) 월 8.1개, 가용 12개 -> 최종 구매 183개
    P019(filter) 월 11.5개, 가용 59개 -> 최종 구매 219개
    P025(filter) 월 14.0개, 가용 89개 -> 최종 구매 247개
[5] 발주 추천 10건(예산 31.2백만 중 30.9백만) + 최종 구매 8건 -> 정비장교 결재 대기 18건, 자동 확정 0건
[6] 운영 가이드 목차: 1. 목적·범위 / 2. 데이터 규격 / 3. 전처리 / 4. 예측 모델 / 5. 안전재고·위험 지수 / 6. 단종 대응 / 7. 결재선·감사 / 8. 재학습 / 9. 폭증 대응 / 10. 연계 인터페이스
[7] 체계 연계 메시지 예시(가상 규격): {"msg_type": "ORDER_REQ", "version": "1.0", "item_id": "P004", "nsn": "2940-00-559-8893", "qty": 183, "risk_index": 87.5, "model_ver": "v1", "approval": "PENDING", "source": "교육용가상자료"}

[자체 점검표]
  [PASS] 원천 자료 정합성(NSN 형식·음수·중복)
  [PASS] 자동 재학습 실행·비교
  [PASS] 운영 모델 MAPE 25% 이하
  [PASS] 자동 확정 0건(HITL)
  [PASS] 운영 가이드 10개 장 완비
  [PASS] 통합 실행 60초 이내
점검 6/6 통과, 운영 가이드 완비율 100%, 전체 실행 3.5s

[브리핑 요약]
- 부품 수요 예측 MAPE 20.0%(이동평균 32.7%)이며, 재학습 모델은 검증 후 교체한다(v1 사용).
- 결품 위험 경보 13개 품목 중 단종·단종예고가 11개로, 단종 대응이 재고 위험의 핵심이다.
- 단종예고 8개 품목의 최종 구매 2779개와 발주 10건을 결재 대기열에 올렸다.
- 모든 발주·최종 구매·대체품 전환은 정비장교 결재 후 집행하며, 수치는 교육용 가상 자료 기준이다.