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

137 KiB

AMM02 정답 코드 - AI 기반 탄약 저장 최적화 모델

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

실행 환경

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

미션 1. 작전 단계별 탄약(구경별) 및 유류(항공유/경유) 기준 소모율(CSR) 데이터 정의

  • 난이도: 초급 / 예상 소요시간: 1.5시간
  • 코드 설명: 작전 단계별 탄약(가상 탄종 A~D)·유류(항공유/경유) 기준 소모율(CSR) 데이터 규격, 업무 프로세스 맵, 데이터 딕셔너리(필드 20개 이상)를 정의하고 가상 샘플 기록을 검증한다.
"""작전 단계별 탄약(가상 탄종 A~D)·유류(항공유/경유) 기준 소모율(CSR) 데이터 규격, 업무 프로세스 맵, 데이터 딕셔너리(필드 20개 이상)를 정의하고 가상 샘플 기록을 검증한다."""
import json
import re
import numpy as np
import pandas as pd

# 1) 업무 프로세스 맵: 소모 기록 -> 소모 예측 -> 보급 수준(DOS) -> 재보급·저장 배치
PROCESS = ["(1) 일일 소모 보고 수집", "(2) 작전 단계·교전 강도 반영 소모량 산출", "(3) 보유량 대비 보급 수준(DOS) 계산",
           "(4) 재보급 소요·우선순위 결정(군수참모 결재)", "(5) 가상 저장고 입고·배치(위치점수·작업동선)"]

# 2) 가상 품목·작전 단계·CSR 계수 (실제 탄약 제원·소모율이 아닌 교육용 가상 계수)
PHASES = ["prep", "attack", "defense", "stability"]
ITEMS = {"ammo_A": "탄약", "ammo_B": "탄약", "ammo_C": "탄약", "ammo_D": "탄약", "fuel_jet": "유류", "fuel_diesel": "유류"}
CSR = pd.DataFrame([[1.0, 0.6, 0.3, 0.2, 1.5, 6.0], [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
                    [4.5, 3.0, 1.6, 0.8, 3.0, 7.0], [0.8, 0.5, 0.3, 0.2, 1.0, 5.0]], index=PHASES, columns=list(ITEMS))

# 3) 데이터 딕셔너리: (그룹, 필드, 자료형, 규칙, 설명)
FIELDS = [
    ("식별", "record_id", str, {"pattern": r"^R\d{5}$"}, "기록 번호"),
    ("식별", "unit_id", str, {"pattern": r"^U\d{2}$"}, "가상 부대 코드"),
    ("식별", "day_of_op", int, {"min": 0, "max": 365}, "작전 경과일"),
    ("작전", "phase", str, {"enum": PHASES}, "작전 단계"),
    ("작전", "unit_size", int, {"min": 50, "max": 3000}, "가상 병력 규모(명)"),
    ("작전", "intensity", float, {"min": 0, "max": 1}, "교전 강도(0~1)"),
    ("작전", "resistance", float, {"min": 0, "max": 5}, "적 저항 수준(0~5)"),
    ("작전", "fire_missions", int, {"min": 0, "max": 100}, "계획 화력 임무 수"),
    ("작전", "move_km", float, {"min": 0, "max": 300}, "일일 이동 거리(km)"),
    ("작전", "weather", str, {"enum": ["clear", "rain", "snow"]}, "기상"),
    ("품목", "item_code", str, {"enum": list(ITEMS)}, "가상 품목 코드"),
    ("품목", "item_class", str, {"enum": ["탄약", "유류"]}, "품목 분류"),
    ("품목", "csr_coeff", float, {"min": 0, "max": 50}, "기준 소모율(가상 계수, 100명·1일)"),
    ("품목", "unit_per_pallet", int, {"min": 1, "max": 200}, "팔레트당 단위 수(가상)"),
    ("소모", "consumption_qty", float, {"min": 0, "max": 5000}, "실제 소모량(단위)"),
    ("소모", "emergency_issue", bool, {}, "긴급 불출 여부(예외 보존)"),
    ("재고", "on_hand_qty", float, {"min": 0, "max": 50000}, "보유량(단위)"),
    ("재고", "dos", float, {"min": 0, "max": 60}, "보급 수준(일분)"),
    ("저장", "storage_id", str, {"enum": ["S1", "S2", "S3"]}, "가상 저장고"),
    ("저장", "slot_id", str, {"pattern": r"^[A-C]\d{2}$"}, "저장 위치(열·칸)"),
    ("저장", "slot_distance_m", float, {"min": 0, "max": 500}, "출입구까지 거리(작업동선)"),
    ("저장", "occupancy_rate", float, {"min": 0, "max": 1}, "저장 위치 점유율"),
    ("관리", "source", str, {"enum": ["일일보고", "AAR", "추정"]}, "자료 출처"),
]


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


def validate(rec):
    errs = []
    for _, f, typ, rule, _ in FIELDS:
        errs += [f"{f}: 필수 필드 누락"] if f not in rec else check(f, rec[f], typ, rule)
    if not errs and ITEMS[rec["item_code"]] != rec["item_class"]:
        errs.append("item_class: 품목 코드와 분류 불일치")
    return errs


# 4) [교육용가상자료] 샘플 기록 30건 생성
rng = np.random.default_rng(1)
records = []
for i in range(30):
    ph, it, size = str(rng.choice(PHASES)), str(rng.choice(list(ITEMS))), int(rng.integers(3, 9)) * 100
    qty = float(round(CSR.loc[ph, it] * size / 100 * rng.uniform(0.7, 1.4), 1))
    stock = float(round(qty * rng.uniform(1, 8), 1))
    records.append({"record_id": f"R{i:05d}", "unit_id": f"U{int(rng.integers(1, 9)):02d}", "day_of_op": int(rng.integers(0, 90)),
                    "phase": ph, "unit_size": size, "intensity": float(round(rng.uniform(0, 1), 2)),
                    "resistance": float(round(rng.uniform(0, 5), 1)), "fire_missions": int(rng.integers(0, 12)),
                    "move_km": float(round(rng.uniform(0, 40), 1)), "weather": str(rng.choice(["clear", "rain", "snow"])),
                    "item_code": it, "item_class": ITEMS[it], "csr_coeff": float(CSR.loc[ph, it]),
                    "unit_per_pallet": 20 if ITEMS[it] == "탄약" else 4, "consumption_qty": qty,
                    "emergency_issue": bool(rng.random() < 0.05), "on_hand_qty": stock, "dos": round(stock / qty, 1),
                    "storage_id": str(rng.choice(["S1", "S2", "S3"])), "slot_id": f"{'ABC'[i % 3]}{i % 10:02d}",
                    "slot_distance_m": float(5 + (i % 10) * 6), "occupancy_rate": float(round(rng.uniform(0.2, 1), 2)), "source": "일일보고"})

print("[업무 프로세스 맵] " + " -> ".join(PROCESS))
ok = sum(1 for r in records if not validate(r))
print(f"\n정상 샘플 {len(records)}건 중 검증 통과 {ok}건, 오류 {len(records) - ok}건")

# 5) 오류 주입 시험
bad = json.loads(json.dumps(records[0]))
bad["phase"] = "assault"
bad["intensity"] = 1.7
del bad["dos"]
bad["slot_id"] = "Z-1"
bad["unit_size"] = "600명"
print("\n[오류 주입 결과]")
for e in validate(bad):
    print(" -", e)
bad2 = dict(records[1], item_code="fuel_jet", item_class="탄약")
print(" -", validate(bad2)[0])

# 6) 데이터 딕셔너리와 CSR 기준표
dd = pd.DataFrame([(g, f, t.__name__, ", ".join(f"{k}={v}" for k, v in r.items()) or "-", d) for g, f, t, r, d in FIELDS],
                  columns=["그룹", "필드", "자료형", "규칙", "설명"])
print(f"\n[데이터 딕셔너리] 필드 {len(dd)}개 (목표 20개 이상: {'달성' if len(dd) >= 20 else '미달'})")
for g, s in dd.groupby("그룹", sort=False)["필드"]:
    print(f"  {g}: {', '.join(s)}")
print("\n[작전 단계별 기준 소모율(CSR) - 가상 계수, 100명·1일 기준]")
print(CSR.to_string())
print("\n주의: 품목 코드·계수·저장고는 모두 교육용 가상 값이다. 실제 체계 연동 시 국방 표준 군수 코드(재고번호 등) 체계에 맞춰 필드를 대응시킨다.")

실행 결과 예시

[업무 프로세스 맵] (1) 일일 소모 보고 수집 -> (2) 작전 단계·교전 강도 반영 소모량 산출 -> (3) 보유량 대비 보급 수준(DOS) 계산 -> (4) 재보급 소요·우선순위 결정(군수참모 결재) -> (5) 가상 저장고 입고·배치(위치점수·작업동선)

정상 샘플 30건 중 검증 통과 30건, 오류 0건

[오류 주입 결과]
 - phase: 허용되지 않은 값 'assault'
 - unit_size: 자료형 오류(str)
 - intensity: 최대값 1 초과
 - dos: 필수 필드 누락
 - slot_id: 형식 불일치 'Z-1'
 - item_class: 품목 코드와 분류 불일치

[데이터 딕셔너리] 필드 23개 (목표 20개 이상: 달성)
  식별: record_id, unit_id, day_of_op
  작전: phase, unit_size, intensity, resistance, fire_missions, move_km, weather
  품목: item_code, item_class, csr_coeff, unit_per_pallet
  소모: consumption_qty, emergency_issue
  재고: on_hand_qty, dos
  저장: storage_id, slot_id, slot_distance_m, occupancy_rate
  관리: source

[작전 단계별 기준 소모율(CSR) - 가상 계수, 100명·1일 기준]
           ammo_A  ammo_B  ammo_C  ammo_D  fuel_jet  fuel_diesel
prep          1.0     0.6     0.3     0.2       1.5          6.0
attack        6.0     3.5     2.0     1.2       5.0         12.0
defense       4.5     3.0     1.6     0.8       3.0          7.0
stability     0.8     0.5     0.3     0.2       1.0          5.0

주의: 품목 코드·계수·저장고는 모두 교육용 가상 값이다. 실제 체계 연동 시 국방 표준 군수 코드(재고번호 등) 체계에 맞춰 필드를 대응시킨다.

미션 2. 부대 규모, 교전 강도, 이동 거리 변수 정규화 및 소모량 산출 전처리 모듈 구현

  • 난이도: 초급~중급 / 예상 소요시간: 2.0시간
  • 코드 설명: 부대 규모·교전 강도·이동 거리 변수를 정제·정규화하고, 병력 100명당 환산 소모율을 산출하는 전처리 모듈을 구현한다. 정제 규칙 10종을 적용하고 긴급 불출 같은 군수 예외 기록은 삭제하지 않고 보존하는지 검증한다.
"""부대 규모·교전 강도·이동 거리 변수를 정제·정규화하고, 병력 100명당 환산 소모율을 산출하는 전처리 모듈을 구현한다. 정제 규칙 10종을 적용하고 긴급 불출 같은 군수 예외 기록은 삭제하지 않고 보존하는지 검증한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.impute import KNNImputer

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 부대·가상 품목 일일 소모 기록 생성기 (실제 탄약 제원·소모율과 무관)
# 가상 부대 8개 x 90일. 소모량 = 가상 CSR 계수 x 병력 규모 x (교전 강도·적 저항·화력 계획 효과) x 잡음
# ---------------------------------------------------------------------------
UNITS = [f"U{i:02d}" for i in range(1, 9)]
ITEMS = ["ammo_A", "ammo_B", "ammo_C", "ammo_D", "fuel_jet", "fuel_diesel"]
AMMO, FUEL = ITEMS[:4], ITEMS[4:]
PHASES = ["prep", "attack", "defense", "stability"]
CSR = {  # 가상 계수: 병력 100명·1일 기준 소모 단위(탄약=상자, 유류=드럼)
    "prep":      [1.0, 0.6, 0.3, 0.2, 1.5, 6.0],
    "attack":    [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
    "defense":   [4.5, 3.0, 1.6, 0.8, 3.0, 7.0],
    "stability": [0.8, 0.5, 0.3, 0.2, 1.0, 5.0],
}
PHASE_RESIST = {"prep": 1.0, "attack": 3.2, "defense": 3.6, "stability": 0.8}


def make_ops(n_days=90, seed=7, surge=0):
    """가상 부대별 일일 작전 기록과 품목별 실제 소모량(정답)을 만든다. surge>0이면 마지막 surge일을 고강도 공격 단계로 바꾼다(스트레스 시험용)."""
    rng = np.random.default_rng(seed)
    rows = []
    for uid in UNITS:
        size = int(rng.integers(3, 9)) * 100
        lens = rng.integers(8, 20, 4)
        sched = ["prep"] * lens[0] + ["attack"] * lens[1] + ["defense"] * lens[2] + ["attack"] * lens[3]
        sched = (sched + ["stability"] * n_days)[:n_days]
        if surge:
            sched[-surge:] = ["attack"] * surge
        r = 1.0
        for d, ph in enumerate(sched):
            r = float(np.clip(0.7 * r + 0.3 * (4.8 if surge and d >= n_days - surge else PHASE_RESIST[ph]) + rng.normal(0, 0.35), 0, 5))
            inten = float(np.clip(r / 5 * 0.9 + rng.normal(0, 0.06), 0, 1))
            fire = int(rng.poisson(1 + 8 * inten))
            move = float(rng.gamma(2, {"prep": 4, "attack": 12, "defense": 2, "stability": 5}[ph]))
            weather = str(rng.choice(["clear", "rain", "snow"], p=[0.7, 0.2, 0.1]))
            emerg = bool(rng.random() < 0.02)
            rec = {"unit": uid, "day": d, "phase": ph, "unit_size": size, "intensity": round(inten, 3),
                   "resistance": round(r, 2), "fire_missions": fire, "move_km": round(move, 1),
                   "weather": weather, "emergency": emerg}
            base = size / 100
            for i, it in enumerate(ITEMS):
                c = CSR[ph][i] * base
                if it in AMMO:
                    c *= (0.5 + 0.6 * inten) * (1 + 0.04 * fire * (1.0 if it in ("ammo_C", "ammo_D") else 0.4))
                elif it == "fuel_jet":
                    c *= (0.6 + 0.5 * inten) + 0.03 * fire
                else:
                    c = c * (0.8 + 0.2 * (weather != "clear")) + 0.05 * move * base
                c *= rng.lognormal(0, 0.12) * (2.5 if emerg else 1.0)
                rec[it] = round(float(c), 1)
            rows.append(rec)
    return pd.DataFrame(rows)
# ---------------------------------------------------------------------------


# 1) [교육용가상자료] 원천 기록에 현장 입력 오류를 일부러 섞는다
clean_truth = make_ops(90, seed=7)
raw = clean_truth.copy()
rng = np.random.default_rng(3)
n = len(raw)
idx = lambda k: rng.choice(n, k, replace=False)
raw = raw.astype({"phase": object, "weather": object, "unit_size": float, "fire_missions": float})
raw.loc[idx(15), "phase"] = raw["phase"].str.upper()                          # 표기 불일치
raw.loc[idx(20), "weather"] = "비"                                            # 한글 표기 혼용
raw.loc[idx(12), "intensity"] = raw["intensity"] * 100                         # 백분율로 입력
raw.loc[idx(10), "ammo_A"] = -raw["ammo_A"]                                   # 부호 오류
raw.loc[idx(6), "fire_missions"] = -1
raw.loc[idx(8), "move_km"] = 999.0                                            # 계기 오류
raw.loc[idx(18), "intensity"] = np.nan                                        # 결측
raw.loc[idx(14), "move_km"] = np.nan
raw.loc[idx(5), "unit_size"] = raw["unit_size"] * 10                          # 단위 오기(자릿수)
raw = pd.concat([raw, raw.sample(12, random_state=1)]).sample(frac=1, random_state=2).reset_index(drop=True)  # 중복·순서 섞임
n_emerg = int(clean_truth["emergency"].sum())

# 2) 정제 규칙 10종
log = []
df = raw.copy()
before = len(df); df = df.drop_duplicates(["unit", "day"]); log.append(("R01", "부대·일자 중복 기록 제거", before - len(df)))
m = df["phase"] != df["phase"].str.strip().str.lower(); df["phase"] = df["phase"].str.strip().str.lower(); log.append(("R02", "작전 단계 표기 통일", int(m.sum())))
wmap = {"비": "rain", "눈": "snow", "맑음": "clear"}; m = df["weather"].isin(wmap); df["weather"] = df["weather"].replace(wmap); log.append(("R03", "기상 표기 통일(한글->코드)", int(m.sum())))
m = df["intensity"] > 1; df.loc[m, "intensity"] /= 100; log.append(("R04", "교전 강도 백분율 입력 환산", int(m.sum())))
m = df["unit_size"] > 3000; df.loc[m, "unit_size"] /= 10; log.append(("R05", "병력 규모 자릿수 오기 보정", int(m.sum())))
m = df[ITEMS].lt(0).any(axis=1) | df["fire_missions"].lt(0)
df[ITEMS] = df[ITEMS].mask(df[ITEMS] < 0); df["fire_missions"] = df["fire_missions"].mask(df["fire_missions"] < 0); log.append(("R06", "음수 소모량·화력 임무 결측 처리", int(m.sum())))
m = (df["move_km"] > 150) & ~df["emergency"]; df.loc[m, "move_km"] = np.nan; log.append(("R07", "이동 거리 이상치(150km 초과) 결측 처리", int(m.sum())))
num = ["unit_size", "intensity", "resistance", "fire_missions", "move_km"] + ITEMS
n_missing = int(df[num].isna().any(axis=1).sum())
mu, sd = df[num].mean(), df[num].std()
df[num] = KNNImputer(n_neighbors=5).fit_transform((df[num] - mu) / sd) * sd.values + mu.values
df["fire_missions"] = df["fire_missions"].round(); log.append(("R08", "결측값 KNN 보정(유사 기록 5건)", n_missing))
df = df.sort_values(["unit", "day"]).reset_index(drop=True)
gaps = int(df.groupby("unit")["day"].apply(lambda s: (s.diff().dropna() != 1).sum()).sum()); log.append(("R09", "부대별 일자 정렬·연속성 확인(누락일)", gaps))
kept = int(df["emergency"].sum()); log.append(("R10", "긴급 불출 예외 기록 보존(삭제 금지)", kept))
print(pd.DataFrame(log, columns=["규칙", "내용", "처리 건수"]).to_string(index=False))

# 3) 정규화와 환산 소모율 산출
for c in ["unit_size", "intensity", "move_km"]:
    df[c + "_n"] = (df[c] - df[c].min()) / (df[c].max() - df[c].min())
for it in ITEMS:
    df[it + "_rate"] = df[it] / (df["unit_size"] / 100)          # 병력 100명·1일 환산 소모율

# 4) 정합성 검증
checks = {
    "중복 0건": not df.duplicated(["unit", "day"]).any(),
    "결측 0건": not df[num].isna().any().any(),
    "작전 단계 허용값": df["phase"].isin(PHASES).all(),
    "기상 허용값": df["weather"].isin(["clear", "rain", "snow"]).all(),
    "교전 강도 0~1": df["intensity"].between(0, 1).all(),
    "소모량 0 이상": (df[ITEMS] >= 0).all().all(),
    "병력 규모 복원 일치": (df["unit_size"].values == clean_truth.sort_values(["unit", "day"])["unit_size"].values).all(),
    "정규화 범위 0~1": df[[c + "_n" for c in ["unit_size", "intensity", "move_km"]]].stack().between(0, 1).all(),
    "일자 연속성": gaps == 0,
    "긴급 불출 기록 보존": kept == n_emerg,
}
rate = np.mean(list(checks.values()))
print(f"\n정합성 검증 {sum(checks.values())}/{len(checks)} 통과, 통과율 {rate:.0%} (목표 100%: {'달성' if rate == 1 else '미달'})")
print(f"정제 규칙 {len(log)}종 적용, 원천 {len(raw)}건 -> 정제 {len(df)}건, 긴급 불출 {kept}건 보존")

# 5) 전후 기술통계 비교
cols = ["unit_size", "intensity", "move_km", "ammo_A"]
comp = pd.DataFrame({"정제 전 평균": raw[cols].mean(), "정제 후 평균": df[cols].mean(),
                     "정제 전 최대": raw[cols].max(), "정제 후 최대": df[cols].max()}).round(2)
print("\n[정제 전후 기술통계]")
print(comp.to_string())
print("\n[작전 단계별 환산 소모율 평균(100명·1일, 가상 단위)]")
print(df.groupby("phase")[[it + "_rate" for it in ITEMS]].mean().round(2).rename(columns=lambda c: c[:-5]).to_string())

fig, ax = plt.subplots(1, 3, figsize=(12, 3.5))
for a, c in zip(ax, ["intensity", "move_km", "unit_size"]):
    a.hist(raw[c].dropna(), bins=30, alpha=0.5, label="raw"); a.hist(df[c], bins=30, alpha=0.7, label="cleaned")
    a.set_title(c); a.set_yscale("log"); a.legend()
plt.suptitle("AMM02 preprocessing: raw vs cleaned (synthetic)"); plt.tight_layout(); plt.show()

실행 결과 예시

 규칙                        내용  처리 건수
R01            부대·일자 중복 기록 제거     12
R02               작전 단계 표기 통일     15
R03          기상 표기 통일(한글->코드)     20
R04           교전 강도 백분율 입력 환산     12
R05           병력 규모 자릿수 오기 보정      5
R06        음수 소모량·화력 임무 결측 처리     16
R07 이동 거리 이상치(150km 초과) 결측 처리      8
R08      결측값 KNN 보정(유사 기록 5건)     54
R09     부대별 일자 정렬·연속성 확인(누락일)      0
R10     긴급 불출 예외 기록 보존(삭제 금지)     13

정합성 검증 10/10 통과, 통과율 100% (목표 100%: 달성)
정제 규칙 10종 적용, 원천 732건 -> 정제 720건, 긴급 불출 13건 보존

[정제 전후 기술통계]
           정제 전 평균  정제 후 평균  정제 전 최대  정제 후 최대
unit_size   688.66   650.00   8000.0   800.00
intensity     0.89     0.35     74.4     0.89
move_km      23.52    12.51    999.0    93.10
ammo_A       16.89    17.27    142.1   142.10

[작전 단계별 환산 소모율 평균(100명·1일, 가상 단위)]
           ammo_A  ammo_B  ammo_C  ammo_D  fuel_jet  fuel_diesel
phase                                                           
attack       5.63    3.24    2.07    1.25      5.36        11.93
defense      4.61    3.08    1.82    0.93      3.55         6.58
prep         0.64    0.39    0.20    0.14      1.18         5.66
stability    0.50    0.32    0.20    0.13      0.77         4.86

미션 3. 교전 일수 기반 1차 탄약/유류 기본 소모량 계산 Baseline 알고리즘 구현

  • 난이도: 중급 / 예상 소요시간: 2.5시간
  • 코드 설명: 교전 일수 x 작전 단계별 기준 소모율(CSR) x 병력 규모로 1차 탄약·유류 기본 소모량을 계산하는 Baseline 알고리즘을 구현하고, 탐색적 분석(EDA 차트 8종, 상관·유의성 검정, ABC 파레토 분석)으로 Baseline의 한계와 핵심 변수를 찾는다.
"""교전 일수 x 작전 단계별 기준 소모율(CSR) x 병력 규모로 1차 탄약·유류 기본 소모량을 계산하는 Baseline 알고리즘을 구현하고, 탐색적 분석(EDA 차트 8종, 상관·유의성 검정, ABC 파레토 분석)으로 Baseline의 한계와 핵심 변수를 찾는다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 부대·가상 품목 일일 소모 기록 생성기 (실제 탄약 제원·소모율과 무관)
# 가상 부대 8개 x 90일. 소모량 = 가상 CSR 계수 x 병력 규모 x (교전 강도·적 저항·화력 계획 효과) x 잡음
# ---------------------------------------------------------------------------
UNITS = [f"U{i:02d}" for i in range(1, 9)]
ITEMS = ["ammo_A", "ammo_B", "ammo_C", "ammo_D", "fuel_jet", "fuel_diesel"]
AMMO, FUEL = ITEMS[:4], ITEMS[4:]
PHASES = ["prep", "attack", "defense", "stability"]
CSR = {  # 가상 계수: 병력 100명·1일 기준 소모 단위(탄약=상자, 유류=드럼)
    "prep":      [1.0, 0.6, 0.3, 0.2, 1.5, 6.0],
    "attack":    [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
    "defense":   [4.5, 3.0, 1.6, 0.8, 3.0, 7.0],
    "stability": [0.8, 0.5, 0.3, 0.2, 1.0, 5.0],
}
PHASE_RESIST = {"prep": 1.0, "attack": 3.2, "defense": 3.6, "stability": 0.8}


def make_ops(n_days=90, seed=7, surge=0):
    """가상 부대별 일일 작전 기록과 품목별 실제 소모량(정답)을 만든다. surge>0이면 마지막 surge일을 고강도 공격 단계로 바꾼다(스트레스 시험용)."""
    rng = np.random.default_rng(seed)
    rows = []
    for uid in UNITS:
        size = int(rng.integers(3, 9)) * 100
        lens = rng.integers(8, 20, 4)
        sched = ["prep"] * lens[0] + ["attack"] * lens[1] + ["defense"] * lens[2] + ["attack"] * lens[3]
        sched = (sched + ["stability"] * n_days)[:n_days]
        if surge:
            sched[-surge:] = ["attack"] * surge
        r = 1.0
        for d, ph in enumerate(sched):
            r = float(np.clip(0.7 * r + 0.3 * (4.8 if surge and d >= n_days - surge else PHASE_RESIST[ph]) + rng.normal(0, 0.35), 0, 5))
            inten = float(np.clip(r / 5 * 0.9 + rng.normal(0, 0.06), 0, 1))
            fire = int(rng.poisson(1 + 8 * inten))
            move = float(rng.gamma(2, {"prep": 4, "attack": 12, "defense": 2, "stability": 5}[ph]))
            weather = str(rng.choice(["clear", "rain", "snow"], p=[0.7, 0.2, 0.1]))
            emerg = bool(rng.random() < 0.02)
            rec = {"unit": uid, "day": d, "phase": ph, "unit_size": size, "intensity": round(inten, 3),
                   "resistance": round(r, 2), "fire_missions": fire, "move_km": round(move, 1),
                   "weather": weather, "emergency": emerg}
            base = size / 100
            for i, it in enumerate(ITEMS):
                c = CSR[ph][i] * base
                if it in AMMO:
                    c *= (0.5 + 0.6 * inten) * (1 + 0.04 * fire * (1.0 if it in ("ammo_C", "ammo_D") else 0.4))
                elif it == "fuel_jet":
                    c *= (0.6 + 0.5 * inten) + 0.03 * fire
                else:
                    c = c * (0.8 + 0.2 * (weather != "clear")) + 0.05 * move * base
                c *= rng.lognormal(0, 0.12) * (2.5 if emerg else 1.0)
                rec[it] = round(float(c), 1)
            rows.append(rec)
    return pd.DataFrame(rows)
# ---------------------------------------------------------------------------


def baseline(df):
    """Baseline: 해당 일의 작전 단계 CSR x (병력/100). 교전 강도·적 저항·화력 계획은 반영하지 않는다."""
    return pd.DataFrame({it: [CSR[p][i] * s / 100 for p, s in zip(df["phase"], df["unit_size"])] for i, it in enumerate(ITEMS)}, index=df.index)


df = make_ops(90, seed=7)
base = baseline(df)

# 1) 교전 일수 기반 부대별 작전 기간 총소모량(계획) vs 실제
df["engaged"] = df["phase"].isin(["attack", "defense"])
plan = pd.concat([df[["unit"]], base], axis=1).groupby("unit").sum()
real = df.groupby("unit")[ITEMS].sum()
days = df.groupby("unit")["engaged"].sum()
tot = pd.DataFrame({"교전일수": days, "탄약 계획": plan[AMMO].sum(axis=1).round(), "탄약 실제": real[AMMO].sum(axis=1).round(),
                    "유류 계획": plan[FUEL].sum(axis=1).round(), "유류 실제": real[FUEL].sum(axis=1).round()})
tot["탄약 오차율"] = ((tot["탄약 계획"] - tot["탄약 실제"]) / tot["탄약 실제"]).map("{:+.1%}".format)
print("[부대별 90일 총소모량: Baseline 계획 vs 실제(가상 단위)]")
print(tot.to_string())

# 2) 일일 예측 오차(MAPE)
mape = (np.abs(base - df[ITEMS]) / df[ITEMS]).mean()
print("\nBaseline 일일 MAPE:", ", ".join(f"{k} {v:.1%}" for k, v in mape.items()))

# 3) 상관·유의성 검정: 병력 100명당 탄약 소모율과 각 변수
df["ammo_rate"] = df[AMMO].sum(axis=1) / (df["unit_size"] / 100)
df["bad_weather"] = (df["weather"] != "clear").astype(int)
cand = ["intensity", "resistance", "fire_missions", "move_km", "engaged", "day", "bad_weather", "unit_size"]
rows = []
for c in cand:
    r, p = stats.pearsonr(df[c].astype(float), df["ammo_rate"])
    rows.append({"변수": c, "상관계수": round(r, 3), "p-value": f"{p:.1e}", "유의(p<0.05,|r|>=0.2)": p < 0.05 and abs(r) >= 0.2})
corr = pd.DataFrame(rows).sort_values("상관계수", key=abs, ascending=False)
n_sig = int(corr.iloc[:, 3].sum())
print("\n[탄약 환산 소모율과의 상관·유의성]")
print(corr.to_string(index=False))
print(f"유의미한 상관 변수 {n_sig}개 (목표 5개 이상: {'달성' if n_sig >= 5 else '미달'})")

# 4) ABC 파레토 분석: 부대 x 탄종 조합별 총소모량(상자) 상위 20%
pareto = real[AMMO].stack().sort_values(ascending=False)
top = pareto.iloc[:int(np.ceil(len(pareto) * 0.2))]
print(f"\nABC 분석: 상위 20% 조합({len(top)}개)이 전체 소모량의 {top.sum() / pareto.sum():.1%}를 차지 -> A등급 집중 관리")
print("A등급 예시:", ", ".join(f"{u}-{i}" for u, i in top.index[:5]))

# 5) EDA 차트 8종
fig, ax = plt.subplots(2, 4, figsize=(16, 7))
ax = ax.ravel()
df.groupby("phase")[AMMO].mean().plot.bar(ax=ax[0], title="1. Mean ammo use by phase")
ax[1].scatter(df["intensity"], df["ammo_rate"], s=6, alpha=0.5); ax[1].set_title("2. Intensity vs ammo rate")
ax[2].scatter(df["fire_missions"], df["ammo_rate"], s=6, alpha=0.5, c="tab:red"); ax[2].set_title("3. Fire missions vs ammo rate")
ax[3].scatter(df["move_km"], df["fuel_diesel"] / (df["unit_size"] / 100), s=6, alpha=0.5, c="tab:green"); ax[3].set_title("4. Move km vs diesel rate")
u = df[df["unit"] == "U01"]
ax[4].plot(u["day"], u[AMMO].sum(axis=1), label="actual"); ax[4].plot(u["day"], base.loc[u.index, AMMO].sum(axis=1), label="baseline")
ax[4].legend(); ax[4].set_title("5. U01 daily ammo: actual vs baseline")
im = ax[5].imshow(df[cand + ["ammo_rate"]].astype(float).corr(), cmap="coolwarm", vmin=-1, vmax=1)
ax[5].set_xticks(range(9)); ax[5].set_xticklabels(cand + ["ammo_rate"], rotation=90, fontsize=7); ax[5].set_yticks([])
ax[5].set_title("6. Correlation heatmap")
ax[6].bar(range(len(pareto)), pareto.cumsum() / pareto.sum()); ax[6].axhline(0.8, ls="--", c="k"); ax[6].set_title("7. Pareto (unit x ammo item)")
ax[7].bar(mape.index, mape.values, color="gray"); ax[7].set_title("8. Baseline MAPE by item"); ax[7].tick_params(axis="x", rotation=45)
plt.suptitle("AMM02 EDA - synthetic consumption data"); plt.tight_layout(); plt.show()
print("EDA 차트 8종 생성 (목표 8종 이상: 달성)")
print("핵심 발견: Baseline은 교전 강도·화력 계획을 반영하지 않아 고강도 일에는 과소, 저강도 일에는 과대 산출한다 -> 미션 4에서 변수 반영 모델로 고도화한다.")

실행 결과 예시

[부대별 90일 총소모량: Baseline 계획 vs 실제(가상 단위)]
      교전일수   탄약 계획   탄약 실제   유류 계획   유류 실제  탄약 오차율
unit                                              
U01     48  5114.0  5120.0  7716.0  7959.0   -0.1%
U02     40  3952.0  3885.0  6328.0  6678.0   +1.7%
U03     46  3725.0  3493.0  5715.0  5616.0   +6.6%
U04     41  2893.0  2586.0  4618.0  4504.0  +11.9%
U05     36  2139.0  1971.0  3574.0  3355.0   +8.5%
U06     34  4028.0  3598.0  6820.0  6404.0  +12.0%
U07     26  2945.0  2578.0  5408.0  5282.0  +14.2%
U08     49  4541.0  4124.0  6776.0  6643.0  +10.1%

Baseline 일일 MAPE: ammo_A 43.8%, ammo_B 42.9%, ammo_C 38.1%, ammo_D 37.9%, fuel_jet 25.9%, fuel_diesel 16.8%

[탄약 환산 소모율과의 상관·유의성]
           변수   상관계수  p-value  유의(p<0.05,|r|>=0.2)
      engaged  0.923 1.4e-300                 True
   resistance  0.854 9.2e-206                 True
    intensity  0.838 1.1e-190                 True
fire_missions  0.639  7.9e-84                 True
          day -0.378  6.9e-26                 True
      move_km  0.294  7.5e-16                 True
  bad_weather  0.072  5.2e-02                False
    unit_size  0.017  6.4e-01                False
유의미한 상관 변수 6개 (목표 5개 이상: 달성)

ABC 분석: 상위 20% 조합(7개)이 전체 소모량의 43.0%를 차지 -> A등급 집중 관리
A등급 예시: U01-ammo_A, U08-ammo_A, U02-ammo_A, U06-ammo_A, U03-ammo_A
EDA 차트 8종 생성 (목표 8종 이상: 달성)
핵심 발견: Baseline은 교전 강도·화력 계획을 반영하지 않아 고강도 일에는 과소, 저강도 일에는 과대 산출한다 -> 미션 4에서 변수 반영 모델로 고도화한다.

미션 4. 적 저항 수준 및 화력 운용 계획 연계 일자별 소모량 시계열 예측 모델 고도화

  • 난이도: 중급 / 예상 소요시간: 3.0시간
  • 코드 설명: 적 저항 수준(전일 실적)과 화력 운용 계획(당일 계획 화력 임무 수)을 연계한 파생 피처를 만들고, 일자별 탄약·유류 소모량 시계열 예측 모델을 시계열 교차검증으로 평가해 Baseline(CSR 고정 계수)과 비교한다.
"""적 저항 수준(전일 실적)과 화력 운용 계획(당일 계획 화력 임무 수)을 연계한 파생 피처를 만들고, 일자별 탄약·유류 소모량 시계열 예측 모델을 시계열 교차검증으로 평가해 Baseline(CSR 고정 계수)과 비교한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import TimeSeriesSplit

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 부대·가상 품목 일일 소모 기록 생성기 (실제 탄약 제원·소모율과 무관)
# 가상 부대 8개 x 90일. 소모량 = 가상 CSR 계수 x 병력 규모 x (교전 강도·적 저항·화력 계획 효과) x 잡음
# ---------------------------------------------------------------------------
UNITS = [f"U{i:02d}" for i in range(1, 9)]
ITEMS = ["ammo_A", "ammo_B", "ammo_C", "ammo_D", "fuel_jet", "fuel_diesel"]
AMMO, FUEL = ITEMS[:4], ITEMS[4:]
PHASES = ["prep", "attack", "defense", "stability"]
CSR = {  # 가상 계수: 병력 100명·1일 기준 소모 단위(탄약=상자, 유류=드럼)
    "prep":      [1.0, 0.6, 0.3, 0.2, 1.5, 6.0],
    "attack":    [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
    "defense":   [4.5, 3.0, 1.6, 0.8, 3.0, 7.0],
    "stability": [0.8, 0.5, 0.3, 0.2, 1.0, 5.0],
}
PHASE_RESIST = {"prep": 1.0, "attack": 3.2, "defense": 3.6, "stability": 0.8}


def make_ops(n_days=90, seed=7, surge=0):
    """가상 부대별 일일 작전 기록과 품목별 실제 소모량(정답)을 만든다. surge>0이면 마지막 surge일을 고강도 공격 단계로 바꾼다(스트레스 시험용)."""
    rng = np.random.default_rng(seed)
    rows = []
    for uid in UNITS:
        size = int(rng.integers(3, 9)) * 100
        lens = rng.integers(8, 20, 4)
        sched = ["prep"] * lens[0] + ["attack"] * lens[1] + ["defense"] * lens[2] + ["attack"] * lens[3]
        sched = (sched + ["stability"] * n_days)[:n_days]
        if surge:
            sched[-surge:] = ["attack"] * surge
        r = 1.0
        for d, ph in enumerate(sched):
            r = float(np.clip(0.7 * r + 0.3 * (4.8 if surge and d >= n_days - surge else PHASE_RESIST[ph]) + rng.normal(0, 0.35), 0, 5))
            inten = float(np.clip(r / 5 * 0.9 + rng.normal(0, 0.06), 0, 1))
            fire = int(rng.poisson(1 + 8 * inten))
            move = float(rng.gamma(2, {"prep": 4, "attack": 12, "defense": 2, "stability": 5}[ph]))
            weather = str(rng.choice(["clear", "rain", "snow"], p=[0.7, 0.2, 0.1]))
            emerg = bool(rng.random() < 0.02)
            rec = {"unit": uid, "day": d, "phase": ph, "unit_size": size, "intensity": round(inten, 3),
                   "resistance": round(r, 2), "fire_missions": fire, "move_km": round(move, 1),
                   "weather": weather, "emergency": emerg}
            base = size / 100
            for i, it in enumerate(ITEMS):
                c = CSR[ph][i] * base
                if it in AMMO:
                    c *= (0.5 + 0.6 * inten) * (1 + 0.04 * fire * (1.0 if it in ("ammo_C", "ammo_D") else 0.4))
                elif it == "fuel_jet":
                    c *= (0.6 + 0.5 * inten) + 0.03 * fire
                else:
                    c = c * (0.8 + 0.2 * (weather != "clear")) + 0.05 * move * base
                c *= rng.lognormal(0, 0.12) * (2.5 if emerg else 1.0)
                rec[it] = round(float(c), 1)
            rows.append(rec)
    return pd.DataFrame(rows)
# ---------------------------------------------------------------------------

# ---- 소모량 예측용 파생 피처: 전일까지의 실적 + 당일 작전·화력 계획(사전에 알 수 있는 값만) ----
FEATS = ["lag1", "lag2", "roll3", "roll7", "csr_base", "size100", "fire_plan", "fire_x_engaged", "resist_lag1",
         "resist_trend", "inten_lag1", "move_plan", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep",
         "item_idx", "is_fuel"]
DERIVED = ["lag1", "lag2", "roll3", "roll7", "csr_base", "fire_x_engaged", "resist_lag1", "resist_trend",
           "inten_lag1", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep"]


def build_long(ops):
    """부대x일자 기록을 (부대, 일자, 품목) 행으로 펼치고 파생 피처를 만든다.
    목표값 y = log(실제 소모량 / CSR 기준량): CSR 기준량에 곱할 '보정 계수'를 예측한다. lag·roll 피처도 같은 비율(실적/기준)이다."""
    df = ops.sort_values(["unit", "day"]).reset_index(drop=True)
    g = df.groupby("unit")
    df["resist_lag1"] = g["resistance"].shift(1)
    df["resist_trend"] = df["resist_lag1"] - g["resistance"].shift(3)
    df["inten_lag1"] = g["intensity"].shift(1)
    df["days_in_phase"] = df.groupby((df["phase"] != g["phase"].shift()).cumsum()).cumcount()
    out = []
    for i, it in enumerate(ITEMS):
        csr = np.array([CSR[p][i] for p in df["phase"]])
        ratio = df[it] / (csr * df["unit_size"] / 100)
        prev = ratio.groupby(df["unit"]).shift(1)
        d = pd.DataFrame({"unit": df["unit"], "day": df["day"], "phase": df["phase"], "item": it, "item_idx": i,
                          "is_fuel": int(it in FUEL), "lag1": prev, "lag2": ratio.groupby(df["unit"]).shift(2),
                          "roll3": prev.groupby(df["unit"]).transform(lambda s: s.rolling(3).mean()),
                          "roll7": prev.groupby(df["unit"]).transform(lambda s: s.rolling(7, min_periods=3).mean()),
                          "csr_base": csr, "size100": df["unit_size"] / 100,
                          "fire_plan": df["fire_missions"], "fire_x_engaged": df["fire_missions"] * df["phase"].isin(["attack", "defense"]),
                          "resist_lag1": df["resist_lag1"], "resist_trend": df["resist_trend"], "inten_lag1": df["inten_lag1"],
                          "move_plan": df["move_km"], "days_in_phase": df["days_in_phase"],
                          "engaged": df["phase"].isin(["attack", "defense"]).astype(int),
                          "ph_attack": (df["phase"] == "attack").astype(int), "ph_defense": (df["phase"] == "defense").astype(int),
                          "ph_prep": (df["phase"] == "prep").astype(int), "emergency": df["emergency"],
                          "y": np.log(ratio), "qty": df[it]})
        out.append(d)
    return pd.concat(out, ignore_index=True).dropna(subset=["lag2", "roll3", "resist_trend"]).reset_index(drop=True)


def to_qty(d, y_pred):
    """예측 보정 계수(log)를 소모량으로 되돌린다: 기준량 x exp(예측)."""
    return np.exp(y_pred) * d["csr_base"].values * d["size100"].values


def mape(y, p):
    return float(np.mean(np.abs(np.asarray(p) - np.asarray(y)) / np.asarray(y)))


data = build_long(make_ops(90, seed=7)).sort_values(["day", "unit", "item"]).reset_index(drop=True)
print(f"학습 행 {len(data)}개(부대 8 x 품목 6 x 일자), 입력 피처 {len(FEATS)}개 중 파생 피처 {len(DERIVED)}개 (목표 10개 이상: {'달성' if len(DERIVED) >= 10 else '미달'})")

# 1) 시계열 교차검증: 과거 일자로 학습하고 이후 일자를 예측(미래 정보 누설 방지)
days = np.sort(data["day"].unique())
rows = []
for k, (tr_d, te_d) in enumerate(TimeSeriesSplit(n_splits=5).split(days), start=1):
    tr, te = data[data["day"].isin(days[tr_d])], data[data["day"].isin(days[te_d])]
    m = GradientBoostingRegressor(n_estimators=150, max_depth=3, learning_rate=0.08, subsample=0.8, random_state=0)
    m.fit(tr[FEATS], tr["y"])
    pred = to_qty(te, m.predict(te[FEATS]))
    base = te["csr_base"] * te["size100"]
    rows.append({"폴드": k, "검증 일자": f"{days[te_d][0]}~{days[te_d][-1]}", "Baseline MAPE": mape(te["qty"], base),
                 "모델 MAPE": mape(te["qty"], pred), "전일 비율 유지 MAPE": mape(te["qty"], te["lag1"] * te["csr_base"] * te["size100"])})
cv = pd.DataFrame(rows)
print("\n[시계열 교차검증 결과]")
print(cv.to_string(index=False, float_format=lambda x: f"{x:.1%}"))
avg = cv["모델 MAPE"].mean()
print(f"평균 MAPE: Baseline {cv['Baseline MAPE'].mean():.1%} -> 모델 {avg:.1%} (목표 25% 이하: {'달성' if avg <= 0.25 else '미달'})")

# 2) 최종 모델: 앞 70일 학습, 마지막 18일 검증
tr, te = data[data["day"] < 72], data[data["day"] >= 72]
model = GradientBoostingRegressor(n_estimators=150, max_depth=3, learning_rate=0.08, subsample=0.8, random_state=0).fit(tr[FEATS], tr["y"])
te = te.assign(pred=to_qty(te, model.predict(te[FEATS])), base=te["csr_base"] * te["size100"])
by_item = te.groupby("item")[["qty", "pred", "base"]].apply(lambda g: pd.Series({"Baseline": mape(g["qty"], g["base"]), "모델": mape(g["qty"], g["pred"])}))
print("\n[품목별 검증 MAPE(마지막 18일)]")
print(by_item.T.to_string(float_format=lambda x: f"{x:.1%}"))

# 3) 피처 중요도(차원 축소 후보 확인)
imp = pd.Series(model.feature_importances_, index=FEATS).sort_values(ascending=False)
print("\n피처 중요도 상위 6:", ", ".join(f"{k}({v:.2f})" for k, v in imp.head(6).items()))
print("중요도 0.01 미만(제거 후보):", ", ".join(imp[imp < 0.01].index) or "없음")

# 4) 시각화
full = data.assign(pred=to_qty(data, model.predict(data[FEATS])))
u = full[(full["unit"] == "U01") & (full["item"] == "ammo_A")]
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].plot(u["day"], u["qty"], label="actual"); ax[0].plot(u["day"], u["pred"], label="model")
ax[0].plot(u["day"], u["csr_base"] * u["size100"], ls="--", label="baseline (CSR)"); ax[0].axvline(72, c="gray", ls=":")
ax[0].set_title("U01 ammo_A daily forecast (right of dotted line = hold-out)"); ax[0].set_xlabel("day"); ax[0].legend()
imp.head(10)[::-1].plot.barh(ax=ax[1], title="Feature importance (top 10)")
plt.tight_layout(); plt.show()

실행 결과 예시

학습 행 4176개(부대 8 x 품목 6 x 일자), 입력 피처 19개 중 파생 피처 14개 (목표 10개 이상: 달성)

[시계열 교차검증 결과]
 폴드 검증 일자  Baseline MAPE  모델 MAPE  전일 비율 유지 MAPE
  1 20~33          15.4%    15.4%          23.0%
  2 34~47          17.7%    19.0%          20.4%
  3 48~61          33.4%    13.0%          17.3%
  4 62~75          51.8%    13.6%          21.7%
  5 76~89          48.7%    13.3%          21.8%
평균 MAPE: Baseline 33.4% -> 모델 14.9% (목표 25% 이하: 달성)

[품목별 검증 MAPE(마지막 18일)]
item      ammo_A  ammo_B  ammo_C  ammo_D  fuel_diesel  fuel_jet
Baseline   66.2%   64.7%   60.0%   55.2%        15.4%     36.6%
모델         12.7%   12.3%   13.0%   13.6%        13.5%     12.6%

피처 중요도 상위 6: roll3(0.25), resist_lag1(0.25), csr_base(0.09), fire_x_engaged(0.09), fire_plan(0.07), roll7(0.04)
중요도 0.01 미만(제거 후보): lag2, ph_defense, ph_attack, engaged, ph_prep

미션 5. 부대별 일일 보급 수준(DOS) 연계 긴급 재보급 소요 산출 엔진 개발

  • 난이도: 중급~고급 / 예상 소요시간: 3.5시간
  • 코드 설명: GBR·랜덤포레스트 앙상블로 소모량 예측을 고도화해 Hold-out 기간에서 Baseline 대비 MAE/RMSE 개선율과 R2를 검증하고, 예측 소모량과 보유량으로 부대·품목별 보급 수준(DOS)을 계산해 긴급 재보급 소요를 산출하는 엔진을 구현한다.
"""GBR·랜덤포레스트 앙상블로 소모량 예측을 고도화해 Hold-out 기간에서 Baseline 대비 MAE/RMSE 개선율과 R2를 검증하고, 예측 소모량과 보유량으로 부대·품목별 보급 수준(DOS)을 계산해 긴급 재보급 소요를 산출하는 엔진을 구현한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 부대·가상 품목 일일 소모 기록 생성기 (실제 탄약 제원·소모율과 무관)
# 가상 부대 8개 x 90일. 소모량 = 가상 CSR 계수 x 병력 규모 x (교전 강도·적 저항·화력 계획 효과) x 잡음
# ---------------------------------------------------------------------------
UNITS = [f"U{i:02d}" for i in range(1, 9)]
ITEMS = ["ammo_A", "ammo_B", "ammo_C", "ammo_D", "fuel_jet", "fuel_diesel"]
AMMO, FUEL = ITEMS[:4], ITEMS[4:]
PHASES = ["prep", "attack", "defense", "stability"]
CSR = {  # 가상 계수: 병력 100명·1일 기준 소모 단위(탄약=상자, 유류=드럼)
    "prep":      [1.0, 0.6, 0.3, 0.2, 1.5, 6.0],
    "attack":    [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
    "defense":   [4.5, 3.0, 1.6, 0.8, 3.0, 7.0],
    "stability": [0.8, 0.5, 0.3, 0.2, 1.0, 5.0],
}
PHASE_RESIST = {"prep": 1.0, "attack": 3.2, "defense": 3.6, "stability": 0.8}


def make_ops(n_days=90, seed=7, surge=0):
    """가상 부대별 일일 작전 기록과 품목별 실제 소모량(정답)을 만든다. surge>0이면 마지막 surge일을 고강도 공격 단계로 바꾼다(스트레스 시험용)."""
    rng = np.random.default_rng(seed)
    rows = []
    for uid in UNITS:
        size = int(rng.integers(3, 9)) * 100
        lens = rng.integers(8, 20, 4)
        sched = ["prep"] * lens[0] + ["attack"] * lens[1] + ["defense"] * lens[2] + ["attack"] * lens[3]
        sched = (sched + ["stability"] * n_days)[:n_days]
        if surge:
            sched[-surge:] = ["attack"] * surge
        r = 1.0
        for d, ph in enumerate(sched):
            r = float(np.clip(0.7 * r + 0.3 * (4.8 if surge and d >= n_days - surge else PHASE_RESIST[ph]) + rng.normal(0, 0.35), 0, 5))
            inten = float(np.clip(r / 5 * 0.9 + rng.normal(0, 0.06), 0, 1))
            fire = int(rng.poisson(1 + 8 * inten))
            move = float(rng.gamma(2, {"prep": 4, "attack": 12, "defense": 2, "stability": 5}[ph]))
            weather = str(rng.choice(["clear", "rain", "snow"], p=[0.7, 0.2, 0.1]))
            emerg = bool(rng.random() < 0.02)
            rec = {"unit": uid, "day": d, "phase": ph, "unit_size": size, "intensity": round(inten, 3),
                   "resistance": round(r, 2), "fire_missions": fire, "move_km": round(move, 1),
                   "weather": weather, "emergency": emerg}
            base = size / 100
            for i, it in enumerate(ITEMS):
                c = CSR[ph][i] * base
                if it in AMMO:
                    c *= (0.5 + 0.6 * inten) * (1 + 0.04 * fire * (1.0 if it in ("ammo_C", "ammo_D") else 0.4))
                elif it == "fuel_jet":
                    c *= (0.6 + 0.5 * inten) + 0.03 * fire
                else:
                    c = c * (0.8 + 0.2 * (weather != "clear")) + 0.05 * move * base
                c *= rng.lognormal(0, 0.12) * (2.5 if emerg else 1.0)
                rec[it] = round(float(c), 1)
            rows.append(rec)
    return pd.DataFrame(rows)
# ---------------------------------------------------------------------------

# ---- 소모량 예측용 파생 피처: 전일까지의 실적 + 당일 작전·화력 계획(사전에 알 수 있는 값만) ----
FEATS = ["lag1", "lag2", "roll3", "roll7", "csr_base", "size100", "fire_plan", "fire_x_engaged", "resist_lag1",
         "resist_trend", "inten_lag1", "move_plan", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep",
         "item_idx", "is_fuel"]
DERIVED = ["lag1", "lag2", "roll3", "roll7", "csr_base", "fire_x_engaged", "resist_lag1", "resist_trend",
           "inten_lag1", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep"]


def build_long(ops):
    """부대x일자 기록을 (부대, 일자, 품목) 행으로 펼치고 파생 피처를 만든다.
    목표값 y = log(실제 소모량 / CSR 기준량): CSR 기준량에 곱할 '보정 계수'를 예측한다. lag·roll 피처도 같은 비율(실적/기준)이다."""
    df = ops.sort_values(["unit", "day"]).reset_index(drop=True)
    g = df.groupby("unit")
    df["resist_lag1"] = g["resistance"].shift(1)
    df["resist_trend"] = df["resist_lag1"] - g["resistance"].shift(3)
    df["inten_lag1"] = g["intensity"].shift(1)
    df["days_in_phase"] = df.groupby((df["phase"] != g["phase"].shift()).cumsum()).cumcount()
    out = []
    for i, it in enumerate(ITEMS):
        csr = np.array([CSR[p][i] for p in df["phase"]])
        ratio = df[it] / (csr * df["unit_size"] / 100)
        prev = ratio.groupby(df["unit"]).shift(1)
        d = pd.DataFrame({"unit": df["unit"], "day": df["day"], "phase": df["phase"], "item": it, "item_idx": i,
                          "is_fuel": int(it in FUEL), "lag1": prev, "lag2": ratio.groupby(df["unit"]).shift(2),
                          "roll3": prev.groupby(df["unit"]).transform(lambda s: s.rolling(3).mean()),
                          "roll7": prev.groupby(df["unit"]).transform(lambda s: s.rolling(7, min_periods=3).mean()),
                          "csr_base": csr, "size100": df["unit_size"] / 100,
                          "fire_plan": df["fire_missions"], "fire_x_engaged": df["fire_missions"] * df["phase"].isin(["attack", "defense"]),
                          "resist_lag1": df["resist_lag1"], "resist_trend": df["resist_trend"], "inten_lag1": df["inten_lag1"],
                          "move_plan": df["move_km"], "days_in_phase": df["days_in_phase"],
                          "engaged": df["phase"].isin(["attack", "defense"]).astype(int),
                          "ph_attack": (df["phase"] == "attack").astype(int), "ph_defense": (df["phase"] == "defense").astype(int),
                          "ph_prep": (df["phase"] == "prep").astype(int), "emergency": df["emergency"],
                          "y": np.log(ratio), "qty": df[it]})
        out.append(d)
    return pd.concat(out, ignore_index=True).dropna(subset=["lag2", "roll3", "resist_trend"]).reset_index(drop=True)


def to_qty(d, y_pred):
    """예측 보정 계수(log)를 소모량으로 되돌린다: 기준량 x exp(예측)."""
    return np.exp(y_pred) * d["csr_base"].values * d["size100"].values


def mape(y, p):
    return float(np.mean(np.abs(np.asarray(p) - np.asarray(y)) / np.asarray(y)))

# ---- 소모량 예측 모델(앙상블) 학습·예측 ----
def fit_models(tr, seed=0):
    """GBR + 랜덤포레스트 두 모델을 학습한다(목표값은 CSR 기준량 대비 log 보정 계수)."""
    g = GradientBoostingRegressor(n_estimators=150, max_depth=3, learning_rate=0.08, subsample=0.8, random_state=seed).fit(tr[FEATS], tr["y"])
    r = RandomForestRegressor(n_estimators=120, min_samples_leaf=3, random_state=seed, n_jobs=1).fit(tr[FEATS], tr["y"])
    return g, r


def predict_qty(models, d, w=0.5):
    """두 모델의 log 예측을 가중 평균해 소모량으로 되돌린다."""
    g, r = models
    return to_qty(d, w * g.predict(d[FEATS]) + (1 - w) * r.predict(d[FEATS]))

DOS_ALERT = 2.0      # [교육용가상자료] 긴급 재보급 경보 기준(보급 수준, 일)
DOS_TARGET = 3.0     # 재보급 후 목표 보급 수준(일)

data = build_long(make_ops(90, seed=7))
tr, te = data[data["day"] < 60], data[data["day"] >= 60].copy()
val = tr[tr["day"] >= 48]

# 1) 앙상블 가중치 튜닝: 학습 기간 마지막 12일(검증용)로 고르고, Hold-out(60~89일)은 최종 평가에만 사용
m_inner = fit_models(tr[tr["day"] < 48])
best_w = min(np.arange(0, 1.01, 0.25), key=lambda w: mape(val["qty"], predict_qty(m_inner, val, w)))
models = fit_models(tr)
te["base"] = te["csr_base"] * te["size100"]
te["gbr"] = to_qty(te, models[0].predict(te[FEATS]))
te["ens"] = predict_qty(models, te, best_w)


def scores(col, d):
    """오차율은 CSR 기준량 대비 비율(환산 소모율) 척도로 계산해 품목 간 단위 차이를 없앤다. R2는 소모량 척도."""
    k = d["csr_base"] * d["size100"]
    y, p = d["qty"] / k, d[col] / k
    return {"MAE(비율)": float((y - p).abs().mean()), "RMSE(비율)": float(np.sqrt(((y - p) ** 2).mean())),
            "MAPE": mape(d["qty"], d[col]), "R2(전체)": r2_score(d["qty"], d[col]),
            "R2(정상소모)": r2_score(d.loc[~d["emergency"], "qty"], d.loc[~d["emergency"], col])}


res = pd.DataFrame({"Baseline(CSR)": scores("base", te), "단일 GBR": scores("gbr", te), f"앙상블(w={best_w:.2f})": scores("ens", te)}).T
print("[Hold-out 60~89일 성능, 부대 8 x 품목 6]")
print(res.to_string(float_format=lambda x: f"{x:.3f}"))
b, e = res.iloc[0], res.iloc[2]
imp_mae, imp_rmse = 1 - e["MAE(비율)"] / b["MAE(비율)"], 1 - e["RMSE(비율)"] / b["RMSE(비율)"]
print(f"오차율 개선: MAE {imp_mae:.1%}, RMSE {imp_rmse:.1%} (목표 20% 이상: {'달성' if min(imp_mae, imp_rmse) >= 0.2 else '미달'})")
print(f"R2 전체 {e['R2(전체)']:.3f} ({'달성' if e['R2(전체)'] >= 0.85 else '미달'}), 긴급 불출 제외 정상 소모 {e['R2(정상소모)']:.3f} (목표 0.85 이상: {'달성' if e['R2(정상소모)'] >= 0.85 else '미달'})")
print("-> 전체 R2 미달 원인은 사전 정보가 없는 긴급 불출 급증(2~3배)이다. 예측 모델이 아니라 긴급 경보·예비 재고로 대응한다.")
em = te[te["emergency"]]
print(f"비정형 수요(긴급 불출 {len(em)}건) MAPE: Baseline {mape(em['qty'], em['base']):.1%}, 앙상블 {mape(em['qty'], em['ens']):.1%} -> 급증은 사전 정보가 없어 별도 경보로 보완")

# 2) DOS 연계 긴급 재보급 엔진: 75일 아침 기준, 보유량(가상) / 예측 일일 소모량
rng = np.random.default_rng(11)
day = te[te["day"] == 75].copy()
day["on_hand"] = np.round(day["base"] * rng.uniform(1.0, 5.0, len(day)), 0)
day["dos"] = day["on_hand"] / day["ens"]
day["need"] = np.maximum(0, np.ceil(DOS_TARGET * day["ens"] - day["on_hand"]))
day["flag"] = np.where(day["dos"] < DOS_ALERT, "긴급", np.where(day["dos"] < DOS_TARGET, "주의", "정상"))
pri = {"긴급": 0, "주의": 1, "정상": 2}
out = day.sort_values(["flag", "dos"], key=lambda s: s.map(pri) if s.name == "flag" else s)
print(f"\n[75일 재보급 소요(DOS {DOS_ALERT:.0f}일 미만 긴급, 목표 {DOS_TARGET:.0f}일)] 긴급 {sum(day['flag'] == '긴급')}건, 주의 {sum(day['flag'] == '주의')}건, 정상 {sum(day['flag'] == '정상')}건")
show = out[out["flag"] != "정상"].head(10)[["unit", "item", "phase", "on_hand", "ens", "dos", "need", "flag"]]
show.columns = ["부대", "품목", "작전단계", "보유량", "예측소모/일", "DOS(일)", "재보급량", "판정"]
print(show.to_string(index=False, float_format=lambda x: f"{x:.1f}"))
real_dos = day["on_hand"] / day["qty"]
hit = ((day["dos"] < DOS_ALERT) == (real_dos < DOS_ALERT)).mean()
print(f"긴급 판정 일치율(실제 소모 기준 DOS와 비교): {hit:.1%}")
print("재보급 소요는 군수참모 결재 전 추천값이며 자동 불출하지 않는다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].scatter(te["qty"], te["base"], s=5, alpha=0.4, label="baseline")
ax[0].scatter(te["qty"], te["ens"], s=5, alpha=0.4, label="ensemble")
mx = te["qty"].max(); ax[0].plot([0, mx], [0, mx], "k--", lw=1)
ax[0].set_xlabel("actual"); ax[0].set_ylabel("predicted"); ax[0].set_title("Hold-out actual vs predicted"); ax[0].legend()
piv = day.pivot(index="unit", columns="item", values="dos").clip(upper=6)
im = ax[1].imshow(piv.values, cmap="RdYlGn", vmin=0, vmax=6, aspect="auto")
ax[1].set_xticks(range(len(piv.columns))); ax[1].set_xticklabels(piv.columns, rotation=30)
ax[1].set_yticks(range(len(piv.index))); ax[1].set_yticklabels(piv.index)
ax[1].set_title("Days of supply by unit/item (day 75)"); plt.colorbar(im, ax=ax[1])
plt.tight_layout(); plt.show()

실행 결과 예시

[Hold-out 60~89일 성능, 부대 8 x 품목 6]
               MAE(비율)  RMSE(비율)  MAPE  R2(전체)  R2(정상소모)
Baseline(CSR)    0.313     0.371 0.493   0.828     0.956
단일 GBR           0.110     0.211 0.140   0.827     0.963
앙상블(w=0.25)      0.107     0.210 0.133   0.827     0.966
오차율 개선: MAE 66.0%, RMSE 43.2% (목표 20% 이상: 달성)
R2 전체 0.827 (미달), 긴급 불출 제외 정상 소모 0.966 (목표 0.85 이상: 달성)
-> 전체 R2 미달 원인은 사전 정보가 없는 긴급 불출 급증(2~3배)이다. 예측 모델이 아니라 긴급 경보·예비 재고로 대응한다.
비정형 수요(긴급 불출 24건) MAPE: Baseline 38.7%, 앙상블 58.8% -> 급증은 사전 정보가 없어 별도 경보로 보완

[75일 재보급 소요(DOS 2일 미만 긴급, 목표 3일)] 긴급 4건, 주의 14건, 정상 30건
 부대          품목      작전단계  보유량  예측소모/일  DOS(일)  재보급량 판정
U01 fuel_diesel stability 45.0    35.7     1.3  62.0 긴급
U08 fuel_diesel stability 37.0    29.2     1.3  51.0 긴급
U01    fuel_jet stability 11.0     6.2     1.8   8.0 긴급
U04      ammo_A stability  4.0     2.2     1.8   3.0 긴급
U02 fuel_diesel stability 63.0    31.4     2.0  32.0 주의
U07      ammo_A stability  7.0     3.5     2.0   4.0 주의
U06      ammo_D stability  2.0     1.0     2.1   1.0 주의
U01      ammo_A stability 10.0     4.5     2.2   4.0 주의
U05      ammo_A stability  5.0     2.3     2.2   2.0 주의
U07      ammo_B stability  5.0     2.2     2.3   2.0 주의
긴급 판정 일치율(실제 소모 기준 DOS와 비교): 95.8%
재보급 소요는 군수참모 결재 전 추천값이며 자동 불출하지 않는다.

미션 6. 작전 지도 연계 탄약/유류 소모 예측 및 재보급 소요 가시화 대시보드 구현

  • 난이도: 중급~고급 / 예상 소요시간: 3.0시간
  • 코드 설명: 작전 지도(가상 격자 좌표)에 부대별 탄약·유류 예측 소모량과 보급 수준(DOS)을 표시하는 대시보드를 만들고, 예측 건마다 기여도 상위 3개 인자를 대체값 비교(occlusion) 방식으로 자동 산출해 자연어 설명문을 생성한다(SHAP summary·waterfall 대체 그림 포함).
"""작전 지도(가상 격자 좌표)에 부대별 탄약·유류 예측 소모량과 보급 수준(DOS)을 표시하는 대시보드를 만들고, 예측 건마다 기여도 상위 3개 인자를 대체값 비교(occlusion) 방식으로 자동 산출해 자연어 설명문을 생성한다(SHAP summary·waterfall 대체 그림 포함)."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 부대·가상 품목 일일 소모 기록 생성기 (실제 탄약 제원·소모율과 무관)
# 가상 부대 8개 x 90일. 소모량 = 가상 CSR 계수 x 병력 규모 x (교전 강도·적 저항·화력 계획 효과) x 잡음
# ---------------------------------------------------------------------------
UNITS = [f"U{i:02d}" for i in range(1, 9)]
ITEMS = ["ammo_A", "ammo_B", "ammo_C", "ammo_D", "fuel_jet", "fuel_diesel"]
AMMO, FUEL = ITEMS[:4], ITEMS[4:]
PHASES = ["prep", "attack", "defense", "stability"]
CSR = {  # 가상 계수: 병력 100명·1일 기준 소모 단위(탄약=상자, 유류=드럼)
    "prep":      [1.0, 0.6, 0.3, 0.2, 1.5, 6.0],
    "attack":    [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
    "defense":   [4.5, 3.0, 1.6, 0.8, 3.0, 7.0],
    "stability": [0.8, 0.5, 0.3, 0.2, 1.0, 5.0],
}
PHASE_RESIST = {"prep": 1.0, "attack": 3.2, "defense": 3.6, "stability": 0.8}


def make_ops(n_days=90, seed=7, surge=0):
    """가상 부대별 일일 작전 기록과 품목별 실제 소모량(정답)을 만든다. surge>0이면 마지막 surge일을 고강도 공격 단계로 바꾼다(스트레스 시험용)."""
    rng = np.random.default_rng(seed)
    rows = []
    for uid in UNITS:
        size = int(rng.integers(3, 9)) * 100
        lens = rng.integers(8, 20, 4)
        sched = ["prep"] * lens[0] + ["attack"] * lens[1] + ["defense"] * lens[2] + ["attack"] * lens[3]
        sched = (sched + ["stability"] * n_days)[:n_days]
        if surge:
            sched[-surge:] = ["attack"] * surge
        r = 1.0
        for d, ph in enumerate(sched):
            r = float(np.clip(0.7 * r + 0.3 * (4.8 if surge and d >= n_days - surge else PHASE_RESIST[ph]) + rng.normal(0, 0.35), 0, 5))
            inten = float(np.clip(r / 5 * 0.9 + rng.normal(0, 0.06), 0, 1))
            fire = int(rng.poisson(1 + 8 * inten))
            move = float(rng.gamma(2, {"prep": 4, "attack": 12, "defense": 2, "stability": 5}[ph]))
            weather = str(rng.choice(["clear", "rain", "snow"], p=[0.7, 0.2, 0.1]))
            emerg = bool(rng.random() < 0.02)
            rec = {"unit": uid, "day": d, "phase": ph, "unit_size": size, "intensity": round(inten, 3),
                   "resistance": round(r, 2), "fire_missions": fire, "move_km": round(move, 1),
                   "weather": weather, "emergency": emerg}
            base = size / 100
            for i, it in enumerate(ITEMS):
                c = CSR[ph][i] * base
                if it in AMMO:
                    c *= (0.5 + 0.6 * inten) * (1 + 0.04 * fire * (1.0 if it in ("ammo_C", "ammo_D") else 0.4))
                elif it == "fuel_jet":
                    c *= (0.6 + 0.5 * inten) + 0.03 * fire
                else:
                    c = c * (0.8 + 0.2 * (weather != "clear")) + 0.05 * move * base
                c *= rng.lognormal(0, 0.12) * (2.5 if emerg else 1.0)
                rec[it] = round(float(c), 1)
            rows.append(rec)
    return pd.DataFrame(rows)
# ---------------------------------------------------------------------------

# ---- 소모량 예측용 파생 피처: 전일까지의 실적 + 당일 작전·화력 계획(사전에 알 수 있는 값만) ----
FEATS = ["lag1", "lag2", "roll3", "roll7", "csr_base", "size100", "fire_plan", "fire_x_engaged", "resist_lag1",
         "resist_trend", "inten_lag1", "move_plan", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep",
         "item_idx", "is_fuel"]
DERIVED = ["lag1", "lag2", "roll3", "roll7", "csr_base", "fire_x_engaged", "resist_lag1", "resist_trend",
           "inten_lag1", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep"]


def build_long(ops):
    """부대x일자 기록을 (부대, 일자, 품목) 행으로 펼치고 파생 피처를 만든다.
    목표값 y = log(실제 소모량 / CSR 기준량): CSR 기준량에 곱할 '보정 계수'를 예측한다. lag·roll 피처도 같은 비율(실적/기준)이다."""
    df = ops.sort_values(["unit", "day"]).reset_index(drop=True)
    g = df.groupby("unit")
    df["resist_lag1"] = g["resistance"].shift(1)
    df["resist_trend"] = df["resist_lag1"] - g["resistance"].shift(3)
    df["inten_lag1"] = g["intensity"].shift(1)
    df["days_in_phase"] = df.groupby((df["phase"] != g["phase"].shift()).cumsum()).cumcount()
    out = []
    for i, it in enumerate(ITEMS):
        csr = np.array([CSR[p][i] for p in df["phase"]])
        ratio = df[it] / (csr * df["unit_size"] / 100)
        prev = ratio.groupby(df["unit"]).shift(1)
        d = pd.DataFrame({"unit": df["unit"], "day": df["day"], "phase": df["phase"], "item": it, "item_idx": i,
                          "is_fuel": int(it in FUEL), "lag1": prev, "lag2": ratio.groupby(df["unit"]).shift(2),
                          "roll3": prev.groupby(df["unit"]).transform(lambda s: s.rolling(3).mean()),
                          "roll7": prev.groupby(df["unit"]).transform(lambda s: s.rolling(7, min_periods=3).mean()),
                          "csr_base": csr, "size100": df["unit_size"] / 100,
                          "fire_plan": df["fire_missions"], "fire_x_engaged": df["fire_missions"] * df["phase"].isin(["attack", "defense"]),
                          "resist_lag1": df["resist_lag1"], "resist_trend": df["resist_trend"], "inten_lag1": df["inten_lag1"],
                          "move_plan": df["move_km"], "days_in_phase": df["days_in_phase"],
                          "engaged": df["phase"].isin(["attack", "defense"]).astype(int),
                          "ph_attack": (df["phase"] == "attack").astype(int), "ph_defense": (df["phase"] == "defense").astype(int),
                          "ph_prep": (df["phase"] == "prep").astype(int), "emergency": df["emergency"],
                          "y": np.log(ratio), "qty": df[it]})
        out.append(d)
    return pd.concat(out, ignore_index=True).dropna(subset=["lag2", "roll3", "resist_trend"]).reset_index(drop=True)


def to_qty(d, y_pred):
    """예측 보정 계수(log)를 소모량으로 되돌린다: 기준량 x exp(예측)."""
    return np.exp(y_pred) * d["csr_base"].values * d["size100"].values


def mape(y, p):
    return float(np.mean(np.abs(np.asarray(p) - np.asarray(y)) / np.asarray(y)))

# ---- 소모량 예측 모델(앙상블) 학습·예측 ----
def fit_models(tr, seed=0):
    """GBR + 랜덤포레스트 두 모델을 학습한다(목표값은 CSR 기준량 대비 log 보정 계수)."""
    g = GradientBoostingRegressor(n_estimators=150, max_depth=3, learning_rate=0.08, subsample=0.8, random_state=seed).fit(tr[FEATS], tr["y"])
    r = RandomForestRegressor(n_estimators=120, min_samples_leaf=3, random_state=seed, n_jobs=1).fit(tr[FEATS], tr["y"])
    return g, r


def predict_qty(models, d, w=0.5):
    """두 모델의 log 예측을 가중 평균해 소모량으로 되돌린다."""
    g, r = models
    return to_qty(d, w * g.predict(d[FEATS]) + (1 - w) * r.predict(d[FEATS]))

KOR = {"lag1": "전일 소모 비율", "lag2": "2일 전 소모 비율", "roll3": "최근 3일 평균 소모 비율", "roll7": "최근 7일 평균 소모 비율",
       "csr_base": "작전 단계 CSR 계수", "size100": "병력 규모", "fire_plan": "계획 화력 임무 수", "fire_x_engaged": "교전 중 화력 계획",
       "resist_lag1": "전일 적 저항 수준", "resist_trend": "적 저항 추세", "inten_lag1": "전일 교전 강도", "move_plan": "계획 이동 거리",
       "days_in_phase": "단계 경과일", "engaged": "교전 단계 여부", "ph_attack": "공격 단계", "ph_defense": "방어 단계",
       "ph_prep": "준비 단계", "item_idx": "품목", "is_fuel": "유류 여부"}

data = build_long(make_ops(90, seed=7))
tr = data[data["day"] < 40]
models = fit_models(tr)
D = 41                                     # 대시보드 기준일(교전 단계가 섞인 날)
day = data[data["day"] == D].reset_index(drop=True)
day["pred"] = predict_qty(models, day)
background = tr[FEATS].median()


def logpred(X):
    return 0.5 * models[0].predict(X) + 0.5 * models[1].predict(X)


def explain(X):
    """각 피처를 학습 기간 중앙값으로 바꿨을 때 예측(log 보정 계수)이 얼마나 변하는지 = 그 피처의 기여도."""
    base = logpred(X)
    C = np.zeros((len(X), len(FEATS)))
    for j, f in enumerate(FEATS):
        Z = X.copy(); Z[f] = background[f]
        C[:, j] = base - logpred(Z)
    return pd.DataFrame(C, columns=FEATS)


C = explain(day[FEATS])
top3 = C.abs().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"[기준일 {D}일] 예측 {len(day)}건(부대 8 x 품목 6), 상위 3개 인자 자동 산출률 {cover:.0%} (목표 100%: {'달성' if cover == 1 else '미달'})")


def sentence(i):
    r = day.loc[i]; parts = []
    for f in top3[i]:
        eff = np.exp(C.loc[i, f]) - 1
        parts.append(f"{KOR[f]}({'+' if eff >= 0 else ''}{eff:.0%})")
    return f"{r['unit']} {r['item']} 예측 {r['pred']:.1f}(기준량 대비 {r['pred'] / (r['csr_base'] * r['size100']):.2f}배): 주요 원인 " + ", ".join(parts)


print("\n[부대별 최대 소모 품목 설명(자연어 요약)]")
lead = day.assign(rel=day["pred"] / (day["csr_base"] * day["size100"])).sort_values("rel", ascending=False).groupby("unit").head(1)
for i in lead.sort_values("unit").index:
    print(" -", sentence(i))

imp = C.abs().mean().sort_values(ascending=False)
print("\n[전체 평균 기여도 상위 5(요약 그림)]", ", ".join(f"{KOR[k]} {v:.3f}" for k, v in imp.head(5).items()))
recon = np.corrcoef(C.sum(axis=1), logpred(day[FEATS]) - logpred(pd.DataFrame([background] * len(day))))[0, 1]
print(f"설명 정합성: 기여도 합과 (예측-기준 예측)의 상관 {recon:.2f} (1에 가까울수록 설명이 예측을 잘 분해)")
print("설명은 판단 보조 근거이며 재보급 결정은 군수참모가 근거를 확인해 내린다.")

# 대시보드: 가상 작전 지도 + 기여도 그림
rng = np.random.default_rng(2)
pos = {u: (rng.uniform(0, 50), rng.uniform(0, 30)) for u in UNITS}
ag = day.groupby("unit").apply(lambda g: pd.Series({"ammo": g.loc[g["is_fuel"] == 0, "pred"].sum(), "fuel": g.loc[g["is_fuel"] == 1, "pred"].sum(),
                                                   "phase": g["phase"].iloc[0]}), include_groups=False)
fig, ax = plt.subplots(1, 3, figsize=(16, 4.5))
for u, r in ag.iterrows():
    x, y = pos[u]
    ax[0].scatter(x, y, s=r["ammo"] * 8 + 20, c="tab:red" if r["phase"] in ("attack", "defense") else "tab:blue", alpha=0.6)
    ax[0].annotate(f"{u}\nA{r['ammo']:.0f}/F{r['fuel']:.0f}", (x, y), fontsize=8, ha="center")
ax[0].set_title(f"Virtual op map day {D} (size=ammo forecast, red=engaged)"); ax[0].set_xlabel("grid x"); ax[0].set_ylabel("grid y")
imp.head(10)[::-1].rename(index=lambda k: k).plot.barh(ax=ax[1], title="Mean |contribution| (summary)")
i0 = lead.index[0]; w = C.loc[i0].reindex(C.loc[i0].abs().sort_values(ascending=False).index[:6])
ax[2].barh(w.index[::-1], w.values[::-1], color=["tab:red" if v > 0 else "tab:green" for v in w.values[::-1]])
ax[2].set_title(f"Waterfall-style: {day.loc[i0, 'unit']} {day.loc[i0, 'item']} (log ratio)")
plt.tight_layout(); plt.show()

실행 결과 예시

[기준일 41일] 예측 48건(부대 8 x 품목 6), 상위 3개 인자 자동 산출률 100% (목표 100%: 달성)

[부대별 최대 소모 품목 설명(자연어 요약)]
 - U01 fuel_jet 예측 25.6(기준량 대비 1.07배): 주요 원인 전일 적 저항 수준(-13%), 최근 7일 평균 소모 비율(+5%), 2일 전 소모 비율(-3%)
 - U02 ammo_D 예측 8.9(기준량 대비 1.05배): 주요 원인 전일 적 저항 수준(-14%), 최근 7일 평균 소모 비율(+3%), 최근 3일 평균 소모 비율(-1%)
 - U03 ammo_D 예측 4.9(기준량 대비 1.01배): 주요 원인 전일 적 저항 수준(-12%), 최근 7일 평균 소모 비율(+4%), 단계 경과일(-2%)
 - U04 fuel_diesel 예측 61.6(기준량 대비 1.03배): 주요 원인 계획 이동 거리(+11%), 전일 적 저항 수준(-8%), 품목(+8%)
 - U05 ammo_C 예측 9.1(기준량 대비 1.14배): 주요 원인 전일 적 저항 수준(-11%), 계획 화력 임무 수(+6%), 교전 중 화력 계획(+4%)
 - U06 ammo_D 예측 10.1(기준량 대비 1.05배): 주요 원인 전일 적 저항 수준(-12%), 최근 7일 평균 소모 비율(+3%), 품목(+1%)
 - U07 fuel_diesel 예측 33.8(기준량 대비 0.97배): 주요 원인 품목(+18%), 전일 적 저항 수준(-14%), 작전 단계 CSR 계수(+13%)
 - U08 ammo_D 예측 6.0(기준량 대비 1.08배): 주요 원인 전일 적 저항 수준(-16%), 적 저항 추세(+3%), 최근 7일 평균 소모 비율(+3%)

[전체 평균 기여도 상위 5(요약 그림)] 전일 적 저항 수준 0.158, 작전 단계 CSR 계수 0.038, 품목 0.026, 계획 화력 임무 수 0.023, 최근 7일 평균 소모 비율 0.021
설명 정합성: 기여도 합과 (예측-기준 예측)의 상관 0.91 (1에 가까울수록 설명이 예측을 잘 분해)
설명은 판단 보조 근거이며 재보급 결정은 군수참모가 근거를 확인해 내린다.

미션 7. 보급 우선순위 지정 및 대체 탄종 불출을 위한 군수참모 HITL 결재선 설계

  • 난이도: 고급 / 예상 소요시간: 3.5시간
  • 코드 설명: 가상 저장고의 품목 배치를 할당 문제(정수계획)로 최적화해 작업동선 비용·공간낭비·위치점수를 현행 배치와 비교하고, 재배치·긴급 재보급 우선순위·대체 탄종 불출 제안을 자동 확정하지 않고 군수참모 결재선(승인·수정·반려·무응답)에 올려 사유와 함께 해시 체인 감사 로그에 남기는 HITL 시뮬레이터를 구현한다.
"""가상 저장고의 품목 배치를 할당 문제(정수계획)로 최적화해 작업동선 비용·공간낭비·위치점수를 현행 배치와 비교하고, 재배치·긴급 재보급 우선순위·대체 탄종 불출 제안을 자동 확정하지 않고 군수참모 결재선(승인·수정·반려·무응답)에 올려 사유와 함께 해시 체인 감사 로그에 남기는 HITL 시뮬레이터를 구현한다."""
import hashlib
import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 부대·가상 품목 일일 소모 기록 생성기 (실제 탄약 제원·소모율과 무관)
# 가상 부대 8개 x 90일. 소모량 = 가상 CSR 계수 x 병력 규모 x (교전 강도·적 저항·화력 계획 효과) x 잡음
# ---------------------------------------------------------------------------
UNITS = [f"U{i:02d}" for i in range(1, 9)]
ITEMS = ["ammo_A", "ammo_B", "ammo_C", "ammo_D", "fuel_jet", "fuel_diesel"]
AMMO, FUEL = ITEMS[:4], ITEMS[4:]
PHASES = ["prep", "attack", "defense", "stability"]
CSR = {  # 가상 계수: 병력 100명·1일 기준 소모 단위(탄약=상자, 유류=드럼)
    "prep":      [1.0, 0.6, 0.3, 0.2, 1.5, 6.0],
    "attack":    [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
    "defense":   [4.5, 3.0, 1.6, 0.8, 3.0, 7.0],
    "stability": [0.8, 0.5, 0.3, 0.2, 1.0, 5.0],
}
PHASE_RESIST = {"prep": 1.0, "attack": 3.2, "defense": 3.6, "stability": 0.8}


def make_ops(n_days=90, seed=7, surge=0):
    """가상 부대별 일일 작전 기록과 품목별 실제 소모량(정답)을 만든다. surge>0이면 마지막 surge일을 고강도 공격 단계로 바꾼다(스트레스 시험용)."""
    rng = np.random.default_rng(seed)
    rows = []
    for uid in UNITS:
        size = int(rng.integers(3, 9)) * 100
        lens = rng.integers(8, 20, 4)
        sched = ["prep"] * lens[0] + ["attack"] * lens[1] + ["defense"] * lens[2] + ["attack"] * lens[3]
        sched = (sched + ["stability"] * n_days)[:n_days]
        if surge:
            sched[-surge:] = ["attack"] * surge
        r = 1.0
        for d, ph in enumerate(sched):
            r = float(np.clip(0.7 * r + 0.3 * (4.8 if surge and d >= n_days - surge else PHASE_RESIST[ph]) + rng.normal(0, 0.35), 0, 5))
            inten = float(np.clip(r / 5 * 0.9 + rng.normal(0, 0.06), 0, 1))
            fire = int(rng.poisson(1 + 8 * inten))
            move = float(rng.gamma(2, {"prep": 4, "attack": 12, "defense": 2, "stability": 5}[ph]))
            weather = str(rng.choice(["clear", "rain", "snow"], p=[0.7, 0.2, 0.1]))
            emerg = bool(rng.random() < 0.02)
            rec = {"unit": uid, "day": d, "phase": ph, "unit_size": size, "intensity": round(inten, 3),
                   "resistance": round(r, 2), "fire_missions": fire, "move_km": round(move, 1),
                   "weather": weather, "emergency": emerg}
            base = size / 100
            for i, it in enumerate(ITEMS):
                c = CSR[ph][i] * base
                if it in AMMO:
                    c *= (0.5 + 0.6 * inten) * (1 + 0.04 * fire * (1.0 if it in ("ammo_C", "ammo_D") else 0.4))
                elif it == "fuel_jet":
                    c *= (0.6 + 0.5 * inten) + 0.03 * fire
                else:
                    c = c * (0.8 + 0.2 * (weather != "clear")) + 0.05 * move * base
                c *= rng.lognormal(0, 0.12) * (2.5 if emerg else 1.0)
                rec[it] = round(float(c), 1)
            rows.append(rec)
    return pd.DataFrame(rows)
# ---------------------------------------------------------------------------

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 저장고(H1~H3) 저장칸 배치 모델. 실제 탄약 저장 기준·안전거리와 무관한 가상 계수다.
# 저장칸: 저장고 x 열(row) x 칸(col). 동선 거리 = 반출장(입구)까지의 가상 거리(m).
# 저장군(G1~G3)은 가상 혼적 제한 그룹이며, 저장고마다 받을 수 있는 저장군이 정해져 있다.
# ---------------------------------------------------------------------------
HOUSES = {"H1": {"groups": {"G1", "G2"}, "rows": 2, "cols": 6, "dock": 15.0},
          "H2": {"groups": {"G1"}, "rows": 2, "cols": 5, "dock": 40.0},
          "H3": {"groups": {"G2", "G3"}, "rows": 2, "cols": 6, "dock": 70.0}}
ITEM_GROUP = {"ammo_A": "G1", "ammo_B": "G1", "ammo_C": "G2", "ammo_D": "G3"}
SLOT_CAP = 6                      # 1열 저장칸 가상 팔레트 수용량(2열은 +4)


def make_slots(seed=3):
    """저장칸 목록과 현재 점유(이동 불가 기존 재고) 여부를 만든다."""
    rng = np.random.default_rng(seed)
    rows = []
    for h, c in HOUSES.items():
        for r in range(c["rows"]):
            for k in range(c["cols"]):
                rows.append({"slot": f"{h}-{r + 1}{k + 1:02d}", "house": h, "dist_m": c["dock"] + 6.0 * k + 4.0 * r,
                             "cap": SLOT_CAP if r == 0 else SLOT_CAP + 4, "fixed": bool(rng.random() < 0.15)})
    return pd.DataFrame(rows)


def make_skus(demand, seed=5):
    """품목별 예측 일일 반출량(demand: 품목->상자/일)을 로트 단위 SKU로 나눈다. 반출 빈도 = 일일 반출 횟수(가상)."""
    rng = np.random.default_rng(seed)
    rows = []
    for it, dq in demand.items():
        n_lot = 5
        share = rng.dirichlet(np.ones(n_lot) * 2)
        for j in range(n_lot):
            rows.append({"sku": f"{it}-L{j + 1}", "item": it, "group": ITEM_GROUP[it],
                         "pallets": int(rng.integers(2, SLOT_CAP + 5)), "freq": round(float(dq * share[j] / 4), 2)})
    return pd.DataFrame(rows)


def feasible(sku, slot):
    return (sku["group"] in HOUSES[slot["house"]]["groups"]) and sku["pallets"] <= slot["cap"] and not slot["fixed"]


def travel(assign, skus, slots):
    """작업동선 = 반출 빈도 x 왕복 거리 합(가상 m/일)."""
    d = slots.set_index("slot")["dist_m"]
    return float(sum(2 * skus.loc[i, "freq"] * d[s] for i, s in assign.items()))


def first_fit(skus, slots):
    """현행 방식(대안 비교 기준): 입고 순서대로 비어 있는 첫 번째 적합 칸에 넣는다."""
    used, out = set(), {}
    for i, k in skus.iterrows():
        for _, s in slots.iterrows():
            if s["slot"] not in used and feasible(k, s):
                out[i] = s["slot"]; used.add(s["slot"]); break
    return out


def optimize(skus, slots, lam=8):
    """할당 문제(정수계획 특수형): 비용 = 반출 빈도 x 왕복거리 + lam x 빈 공간(공간낭비). 부적합 칸은 큰 비용."""
    cost = np.full((len(skus), len(slots)), 1e7)
    for i, k in skus.iterrows():
        for j, s in slots.iterrows():
            if feasible(k, s):
                cost[i, j] = 2 * k["freq"] * s["dist_m"] + lam * (s["cap"] - k["pallets"])
    ri, cj = linear_sum_assignment(cost)
    return {int(i): slots.loc[j, "slot"] for i, j in zip(ri, cj) if cost[i, j] < 1e7}, cost


def loc_score(assign, skus, slots):
    """위치점수(0~100): 반출 빈도 상위 SKU일수록 반출장 가까운 칸에 있는 정도(가중 평균)."""
    d = slots.set_index("slot")["dist_m"]; dmin, dmax = d.min(), d.max()
    w = skus["freq"] / skus["freq"].sum()
    return float(100 * sum(w[i] * (1 - (d[s] - dmin) / (dmax - dmin)) for i, s in assign.items()))
# ---------------------------------------------------------------------------

COST_PER_M = 0.02       # [교육용가상자료] 지게차 왕복 1m당 가상 작업비용 지수
SUB_RATE = {"ammo_C": ("ammo_D", 0.8)}   # 가상 대체 탄종과 환산 계수(교육용 가정값)


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


# 1) 반출 수요: 최근 14일 품목별 일일 평균 소모(미션 5 예측값 대신 사용, 가상 상자/일)
ops = make_ops(90, seed=7)
demand = ops[ops["day"] >= 76][AMMO].sum().div(14).round(1).to_dict()
skus, slots = make_skus(demand), make_slots()
print(f"[가상 저장고] 저장칸 {len(slots)}개(현재 점유·이동 불가 {slots['fixed'].sum()}개), SKU {len(skus)}개, 일일 반출 수요 {demand}")

cur = first_fit(skus, slots)
opt, cost = optimize(skus, slots)


def kpi(a):
    used = slots.set_index("slot").loc[list(a.values()), "cap"].sum()
    waste = 1 - skus.loc[list(a.keys()), "pallets"].sum() / used
    t = travel(a, skus, slots)
    return {"작업동선(m/일)": t, "작업비용 지수": t * COST_PER_M, "공간낭비율": waste, "위치점수": loc_score(a, skus, slots), "배치 SKU": len(a)}


kp = pd.DataFrame({"현행(입고순)": kpi(cur), "최적화": kpi(opt)}).T
print(kp.to_string(float_format=lambda x: f"{x:.2f}"))
cut = 1 - kp.loc["최적화", "작업비용 지수"] / kp.loc["현행(입고순)", "작업비용 지수"]
print(f"작업비용 절감률 {cut:.1%} (목표 15% 이상: {'달성' if cut >= 0.15 else '미달'})")

top = skus["freq"].idxmax()
alt = [(slots.loc[j, "slot"], cost[top, j]) for j in np.argsort(cost[top]) if cost[top, j] < 1e7][:3]
print(f"대안배치(반출 최다 {skus.loc[top, 'sku']}): " + ", ".join(f"{s}(비용 {c:.0f})" for s, c in alt) + " -> 1안 점유 시 2·3안 순으로 제시")

# 2) 민감도: 돌발 수요(가상 전시 전환으로 ammo_A 반출 2배)
sk2 = skus.copy(); sk2.loc[sk2["item"] == "ammo_A", "freq"] *= 2
o2, _ = optimize(sk2, slots)
moved = sum(o2[i] != opt[i] for i in opt)
print(f"민감도(ammo_A 반출 2배): 최적 배치 변경 {moved}칸, 작업동선 현행 {travel(cur, sk2, slots):.0f} -> 재최적화 {travel(o2, sk2, slots):.0f} m/일")

# 3) HITL 결재선: 제안은 모두 '결재 대기'로 시작하며 자동 확정하지 않는다
moves = [(skus.loc[i, "sku"], cur.get(i), opt[i], skus.loc[i, "freq"]) for i in opt if cur.get(i) != opt[i]]
moves.sort(key=lambda m: -m[3])
props = [{"id": f"P{n + 1:02d}", "type": "재배치", "what": f"{s} {a}->{b}", "evid": f"반출 {f:.1f}회/일, 동선 감소"} for n, (s, a, b, f) in enumerate(moves[:6])]
rng = np.random.default_rng(4)
dos = {u: round(float(rng.uniform(0.8, 4.0)), 1) for u in UNITS}
urgent = sorted([u for u in UNITS if dos[u] < 2.0], key=lambda u: dos[u])
props += [{"id": f"P{len(props) + k + 1:02d}", "type": "재보급 우선순위", "what": f"{u} 탄약 우선 불출({k + 1}순위)", "evid": f"DOS {dos[u]}일"} for k, u in enumerate(urgent)]
src, (dst, rate) = "ammo_C", SUB_RATE["ammo_C"]
props.append({"id": f"P{len(props) + 1:02d}", "type": "대체 탄종", "what": f"{src} 부족분 40상자 -> {dst} {40 * rate:.0f}상자", "evid": f"가상 환산 계수 {rate}"})

audit = Audit()
for p in props:
    audit.add(t="제안", pid=p["id"], kind=p["type"], what=p["what"], evid=p["evid"], status="결재대기", by="AI")


def staff_decision(p, k):
    """가상 군수참모 판단(시뮬레이션): 근거 확인 후 승인·수정·반려, 일부는 시간 내 무응답."""
    if p["type"] == "대체 탄종":
        return "수정", "대체 수량 절반만 승인, 잔여는 상급부대 청구"
    if k % 5 == 4:
        return "무응답", "결재 시한 초과 -> 현행 유지"
    if p["type"] == "재배치" and k % 3 == 2:
        return "반려", "해당 저장고 점검 일정과 겹침"
    return "승인", "근거 확인"


executed = []
for k, p in enumerate(props):
    d, why = staff_decision(p, k)
    audit.add(t="결재", pid=p["id"], decision=d, reason=why, by="군수참모(가상)")
    if d in ("승인", "수정"):
        executed.append(p["id"]); audit.add(t="집행", pid=p["id"], by="저장관리반(가상)")
tab = pd.DataFrame([{"제안": e["pid"], "결정": e["decision"], "사유": e["reason"]} for e in audit.log if e["t"] == "결재"])
info = pd.DataFrame(props).set_index("id")
tab["유형"], tab["내용"] = tab["제안"].map(info["type"]), tab["제안"].map(info["what"])
print(f"\n[군수참모 결재 결과] 제안 {len(props)}건")
print(tab[["제안", "유형", "내용", "결정", "사유"]].to_string(index=False))
auto = [e for e in audit.log if e["t"] == "집행" and not any(x["t"] == "결재" and x["pid"] == e["pid"] and x["decision"] in ("승인", "수정") for x in audit.log)]
print(f"결재 없이 집행된 건 {len(auto)}건 (자동 확정 금지: {'준수' if not auto else '위반'}), 결재 사유 기록률 {tab['사유'].astype(bool).mean():.0%}")
print(f"감사 로그 {len(audit.log)}건 무결성 검증: {'정상' if audit.verify() is None else '변조 발견'}")
i3 = next(e["seq"] for e in audit.log if e["t"] == "결재" and e["decision"] == "반려")
audit.log[i3]["decision"] = "승인"
print(f"변조 시험({i3}번 반려 기록을 승인으로 수정) -> {audit.verify()}번 기록에서 변조 검출")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
for a, (nm, asg) in zip(ax, [("current (first-fit)", cur), ("optimized", opt)]):
    s = slots.set_index("slot")
    a.scatter(s["dist_m"], [list(HOUSES).index(h) for h in s["house"]], c="lightgray", s=60, marker="s")
    for i, sl in asg.items():
        a.scatter(s.loc[sl, "dist_m"], list(HOUSES).index(s.loc[sl, "house"]), s=40 + 30 * skus.loc[i, "freq"], alpha=0.6)
    a.set_yticks(range(3)); a.set_yticklabels(list(HOUSES)); a.set_xlabel("distance to dock (m)"); a.set_title(f"Placement: {nm} (size=issue freq)")
plt.tight_layout(); plt.show()

실행 결과 예시

[가상 저장고] 저장칸 34개(현재 점유·이동 불가 6개), SKU 20개, 일일 반출 수요 {'ammo_A': 25.8, 'ammo_B': 16.7, 'ammo_C': 10.2, 'ammo_D': 7.1}
         작업동선(m/일)  작업비용 지수  공간낭비율  위치점수  배치 SKU
현행(입고순)    1506.20    30.12   0.29 60.18   20.00
최적화        1195.08    23.90   0.29 71.88   20.00
작업비용 절감률 20.7% (목표 15% 이상: 달성)
대안배치(반출 최다 ammo_A-L3): H1-102(비용 126), H1-201(비용 147), H1-103(비용 157) -> 1안 점유 시 2·3안 순으로 제시
민감도(ammo_A 반출 2배): 최적 배치 변경 7칸, 작업동선 현행 1919 -> 재최적화 1491 m/일

[군수참모 결재 결과] 제안 9건
 제안       유형                             내용  결정                        사유
P01      재배치       ammo_A-L3 H1-104->H1-201  승인                     근거 확인
P02      재배치       ammo_B-L4 H1-206->H1-102  승인                     근거 확인
P03      재배치       ammo_A-L5 H1-106->H1-103  반려          해당 저장고 점검 일정과 겹침
P04      재배치       ammo_A-L2 H1-103->H2-101  승인                     근거 확인
P05      재배치       ammo_B-L3 H1-205->H1-202 무응답         결재 시한 초과 -> 현행 유지
P06      재배치       ammo_B-L5 H2-101->H2-102  반려          해당 저장고 점검 일정과 겹침
P07 재보급 우선순위              U04 탄약 우선 불출(1순위)  승인                     근거 확인
P08 재보급 우선순위              U08 탄약 우선 불출(2순위)  승인                     근거 확인
P09    대체 탄종 ammo_C 부족분 40상자 -> ammo_D 32상자  수정 대체 수량 절반만 승인, 잔여는 상급부대 청구
결재 없이 집행된 건 0건 (자동 확정 금지: 준수), 결재 사유 기록률 100%
감사 로그 24건 무결성 검증: 정상
변조 시험(13번 반려 기록을 승인으로 수정) -> 13번 기록에서 변조 검출

미션 8. 7일간 고강도 전면전 모의 탄약/유류 소모량 극한 스트레스 시뮬레이션

  • 난이도: 고급 / 예상 소요시간: 3.0시간
  • 코드 설명: 7일간 고강도 전면전(가상) 상황에서 부대별 탄약·유류 소모량과 보유량을 일일 모의하고, 보고 누락·보고 지연·긴급 수요 급증·부대 수 10배 부하를 주입해 재보급 경보 대시보드의 갱신 시간(2초 이내)과 경보 발령 정확도(95% 이상)를 검증한다.
"""7일간 고강도 전면전(가상) 상황에서 부대별 탄약·유류 소모량과 보유량을 일일 모의하고, 보고 누락·보고 지연·긴급 수요 급증·부대 수 10배 부하를 주입해 재보급 경보 대시보드의 갱신 시간(2초 이내)과 경보 발령 정확도(95% 이상)를 검증한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 부대·가상 품목 일일 소모 기록 생성기 (실제 탄약 제원·소모율과 무관)
# 가상 부대 8개 x 90일. 소모량 = 가상 CSR 계수 x 병력 규모 x (교전 강도·적 저항·화력 계획 효과) x 잡음
# ---------------------------------------------------------------------------
UNITS = [f"U{i:02d}" for i in range(1, 9)]
ITEMS = ["ammo_A", "ammo_B", "ammo_C", "ammo_D", "fuel_jet", "fuel_diesel"]
AMMO, FUEL = ITEMS[:4], ITEMS[4:]
PHASES = ["prep", "attack", "defense", "stability"]
CSR = {  # 가상 계수: 병력 100명·1일 기준 소모 단위(탄약=상자, 유류=드럼)
    "prep":      [1.0, 0.6, 0.3, 0.2, 1.5, 6.0],
    "attack":    [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
    "defense":   [4.5, 3.0, 1.6, 0.8, 3.0, 7.0],
    "stability": [0.8, 0.5, 0.3, 0.2, 1.0, 5.0],
}
PHASE_RESIST = {"prep": 1.0, "attack": 3.2, "defense": 3.6, "stability": 0.8}


def make_ops(n_days=90, seed=7, surge=0):
    """가상 부대별 일일 작전 기록과 품목별 실제 소모량(정답)을 만든다. surge>0이면 마지막 surge일을 고강도 공격 단계로 바꾼다(스트레스 시험용)."""
    rng = np.random.default_rng(seed)
    rows = []
    for uid in UNITS:
        size = int(rng.integers(3, 9)) * 100
        lens = rng.integers(8, 20, 4)
        sched = ["prep"] * lens[0] + ["attack"] * lens[1] + ["defense"] * lens[2] + ["attack"] * lens[3]
        sched = (sched + ["stability"] * n_days)[:n_days]
        if surge:
            sched[-surge:] = ["attack"] * surge
        r = 1.0
        for d, ph in enumerate(sched):
            r = float(np.clip(0.7 * r + 0.3 * (4.8 if surge and d >= n_days - surge else PHASE_RESIST[ph]) + rng.normal(0, 0.35), 0, 5))
            inten = float(np.clip(r / 5 * 0.9 + rng.normal(0, 0.06), 0, 1))
            fire = int(rng.poisson(1 + 8 * inten))
            move = float(rng.gamma(2, {"prep": 4, "attack": 12, "defense": 2, "stability": 5}[ph]))
            weather = str(rng.choice(["clear", "rain", "snow"], p=[0.7, 0.2, 0.1]))
            emerg = bool(rng.random() < 0.02)
            rec = {"unit": uid, "day": d, "phase": ph, "unit_size": size, "intensity": round(inten, 3),
                   "resistance": round(r, 2), "fire_missions": fire, "move_km": round(move, 1),
                   "weather": weather, "emergency": emerg}
            base = size / 100
            for i, it in enumerate(ITEMS):
                c = CSR[ph][i] * base
                if it in AMMO:
                    c *= (0.5 + 0.6 * inten) * (1 + 0.04 * fire * (1.0 if it in ("ammo_C", "ammo_D") else 0.4))
                elif it == "fuel_jet":
                    c *= (0.6 + 0.5 * inten) + 0.03 * fire
                else:
                    c = c * (0.8 + 0.2 * (weather != "clear")) + 0.05 * move * base
                c *= rng.lognormal(0, 0.12) * (2.5 if emerg else 1.0)
                rec[it] = round(float(c), 1)
            rows.append(rec)
    return pd.DataFrame(rows)
# ---------------------------------------------------------------------------

# ---- 소모량 예측용 파생 피처: 전일까지의 실적 + 당일 작전·화력 계획(사전에 알 수 있는 값만) ----
FEATS = ["lag1", "lag2", "roll3", "roll7", "csr_base", "size100", "fire_plan", "fire_x_engaged", "resist_lag1",
         "resist_trend", "inten_lag1", "move_plan", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep",
         "item_idx", "is_fuel"]
DERIVED = ["lag1", "lag2", "roll3", "roll7", "csr_base", "fire_x_engaged", "resist_lag1", "resist_trend",
           "inten_lag1", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep"]


def build_long(ops):
    """부대x일자 기록을 (부대, 일자, 품목) 행으로 펼치고 파생 피처를 만든다.
    목표값 y = log(실제 소모량 / CSR 기준량): CSR 기준량에 곱할 '보정 계수'를 예측한다. lag·roll 피처도 같은 비율(실적/기준)이다."""
    df = ops.sort_values(["unit", "day"]).reset_index(drop=True)
    g = df.groupby("unit")
    df["resist_lag1"] = g["resistance"].shift(1)
    df["resist_trend"] = df["resist_lag1"] - g["resistance"].shift(3)
    df["inten_lag1"] = g["intensity"].shift(1)
    df["days_in_phase"] = df.groupby((df["phase"] != g["phase"].shift()).cumsum()).cumcount()
    out = []
    for i, it in enumerate(ITEMS):
        csr = np.array([CSR[p][i] for p in df["phase"]])
        ratio = df[it] / (csr * df["unit_size"] / 100)
        prev = ratio.groupby(df["unit"]).shift(1)
        d = pd.DataFrame({"unit": df["unit"], "day": df["day"], "phase": df["phase"], "item": it, "item_idx": i,
                          "is_fuel": int(it in FUEL), "lag1": prev, "lag2": ratio.groupby(df["unit"]).shift(2),
                          "roll3": prev.groupby(df["unit"]).transform(lambda s: s.rolling(3).mean()),
                          "roll7": prev.groupby(df["unit"]).transform(lambda s: s.rolling(7, min_periods=3).mean()),
                          "csr_base": csr, "size100": df["unit_size"] / 100,
                          "fire_plan": df["fire_missions"], "fire_x_engaged": df["fire_missions"] * df["phase"].isin(["attack", "defense"]),
                          "resist_lag1": df["resist_lag1"], "resist_trend": df["resist_trend"], "inten_lag1": df["inten_lag1"],
                          "move_plan": df["move_km"], "days_in_phase": df["days_in_phase"],
                          "engaged": df["phase"].isin(["attack", "defense"]).astype(int),
                          "ph_attack": (df["phase"] == "attack").astype(int), "ph_defense": (df["phase"] == "defense").astype(int),
                          "ph_prep": (df["phase"] == "prep").astype(int), "emergency": df["emergency"],
                          "y": np.log(ratio), "qty": df[it]})
        out.append(d)
    return pd.concat(out, ignore_index=True).dropna(subset=["lag2", "roll3", "resist_trend"]).reset_index(drop=True)


def to_qty(d, y_pred):
    """예측 보정 계수(log)를 소모량으로 되돌린다: 기준량 x exp(예측)."""
    return np.exp(y_pred) * d["csr_base"].values * d["size100"].values


def mape(y, p):
    return float(np.mean(np.abs(np.asarray(p) - np.asarray(y)) / np.asarray(y)))

# ---- 소모량 예측 모델(앙상블) 학습·예측 ----
def fit_models(tr, seed=0):
    """GBR + 랜덤포레스트 두 모델을 학습한다(목표값은 CSR 기준량 대비 log 보정 계수)."""
    g = GradientBoostingRegressor(n_estimators=150, max_depth=3, learning_rate=0.08, subsample=0.8, random_state=seed).fit(tr[FEATS], tr["y"])
    r = RandomForestRegressor(n_estimators=120, min_samples_leaf=3, random_state=seed, n_jobs=1).fit(tr[FEATS], tr["y"])
    return g, r


def predict_qty(models, d, w=0.5):
    """두 모델의 log 예측을 가중 평균해 소모량으로 되돌린다."""
    g, r = models
    return to_qty(d, w * g.predict(d[FEATS]) + (1 - w) * r.predict(d[FEATS]))

SURGE = 7            # 고강도 모의 일수
DOS_ALERT = 2.0      # [교육용가상자료] 경보 기준 보급 수준(일)

data = build_long(make_ops(90, seed=7, surge=SURGE))
start = 90 - SURGE
models = fit_models(data[data["day"] < start])
test = data[data["day"] >= start].copy()


def inject(d, miss=0.0, delay=0.0, spike=0.0, seed=0):
    """오류 주입: 보고 누락(전일 실적 결측 -> CSR 기준 1.0 대입), 보고 지연(전일 대신 2일 전 실적), 긴급 수요 급증(긴급 불출 요청 접수, 실제 소모 x2.5)."""
    rng = np.random.default_rng(seed); d = d.copy()
    m = rng.random(len(d)) < miss
    d.loc[m, ["lag1", "lag2", "roll3", "roll7"]] = 1.0
    l = (rng.random(len(d)) < delay) & ~m
    d.loc[l, "lag1"] = d.loc[l, "lag2"]
    s = rng.random(len(d)) < spike
    d.loc[s, "qty"] *= 2.5
    d.loc[s, "emergency"] = True
    return d


def dashboard(d, stock0, margin=1.0):
    """대시보드 1회 갱신: 예측 -> 보유량 경과 -> DOS·경보 -> 부대별 집계 표. 반환: 경보표, 집계표."""
    # 긴급 불출 요청이 접수된 건은 요청 배수(가상 2.5배)를 예측에 반영한다(요청은 사전에 알 수 있는 정보)
    d = d.assign(pred=predict_qty(models, d) * np.where(d["emergency"], 2.5, 1.0)).sort_values(["unit", "item", "day"])
    resup = d["csr_base"] * d["size100"] * 1.0          # 계획 재보급량(작전 단계 CSR 기준, 가상)
    used = d.groupby(["unit", "item"])["qty"].cumsum() - d["qty"]
    got = resup.groupby([d["unit"], d["item"]]).cumsum() - resup
    d["on_hand"] = d.set_index(["unit", "item"]).index.map(stock0).values + got.values - used.values
    d["dos_pred"] = d["on_hand"].clip(lower=0) / d["pred"]
    d["dos_true"] = d["on_hand"].clip(lower=0) / d["qty"]
    d["alert"], d["truth"] = d["dos_pred"] < DOS_ALERT * margin, d["dos_true"] < DOS_ALERT
    agg = d.groupby(["day", "unit"]).agg(alerts=("alert", "sum"), ammo=("pred", "sum"))
    return d, agg


stock0 = (test[test["day"] == start].set_index(["unit", "item"])["csr_base"] * test[test["day"] == start].set_index(["unit", "item"])["size100"] * 3.0).to_dict()
cases = [("정상 보고", {}), ("보고 누락 10%", {"miss": 0.1}), ("보고 지연 20%", {"delay": 0.2}),
         ("긴급 급증 5%", {"spike": 0.05}), ("복합(누락10+지연20+급증5)", {"miss": 0.1, "delay": 0.2, "spike": 0.05})]
rows = []
for name, kw in cases:
    d = inject(test, seed=1, **kw)
    t0 = time.perf_counter(); out, agg = dashboard(d, stock0); dt = time.perf_counter() - t0
    tp = int((out["alert"] & out["truth"]).sum()); fn = int((~out["alert"] & out["truth"]).sum())
    rows.append({"시나리오": name, "갱신(s)": dt, "경보": int(out["alert"].sum()), "실제 부족": int(out["truth"].sum()),
                 "정확도": (out["alert"] == out["truth"]).mean(), "부족 탐지율": tp / max(tp + fn, 1)})
res = pd.DataFrame(rows)
print(f"[7일 고강도 모의] 부대 8 x 품목 6 x {SURGE}일 = {len(test)}건, 시작 보유량 3일분(가상)")
print(res.to_string(index=False, float_format=lambda x: f"{x:.3f}"))
acc = res["정확도"].min()
print(f"경보 정확도 최저 {acc:.1%} (목표 95% 이상: {'달성' if acc >= 0.95 else '미달'})")

# 부하 시험: 가상 부대 10배(80개) 데이터로 대시보드 갱신 시간 측정
big = pd.concat([test.assign(unit=test["unit"] + f"_{k}") for k in range(10)], ignore_index=True)
st_big = {(u + f"_{k}", it): v for (u, it), v in stock0.items() for k in range(10)}
t0 = time.perf_counter(); out_b, _ = dashboard(big, st_big); dt_b = time.perf_counter() - t0
mx = max(res["갱신(s)"].max(), dt_b)
print(f"부하 시험(부대 80개, {len(big)}건) 갱신 {dt_b:.2f}s, 최대 갱신 {mx:.2f}s (목표 2초 이내: {'달성' if mx <= 2 else '미달'})")

base, _ = dashboard(test, stock0)
first = base[base["truth"]].groupby("item")["day"].min().sub(start - 1).astype(int)
print("품목별 첫 부족 발생(모의 n일차):", ", ".join(f"{k} {v}일차" for k, v in first.items()) or "없음")
ammo_tot = base[base["is_fuel"] == 0].groupby("day")[["qty", "pred"]].sum()
print(f"7일 탄약 총소모(가상 상자): 실제 {ammo_tot['qty'].sum():.0f}, 예측 {ammo_tot['pred'].sum():.0f}, 오차 {ammo_tot['pred'].sum() / ammo_tot['qty'].sum() - 1:+.1%}")
tr_off = []
for mg in (1.0, 1.1, 1.2):
    o, _ = dashboard(test, stock0, mg)
    tr_off.append(f"x{mg}: 정확도 {(o['alert'] == o['truth']).mean():.1%}/탐지율 {(o['alert'] & o['truth']).sum() / o['truth'].sum():.0%}")
print("경보 안전계수 상충관계(정상 보고):", ", ".join(tr_off), "-> 부족 놓침을 줄이려면 안전계수를 올리고 정확도 저하를 감수한다")
print("경보는 재보급 검토 요청이며 불출·우선순위 결정은 군수참모 결재를 거친다.")

fig, ax = plt.subplots(1, 3, figsize=(16, 4))
ammo_tot.plot(ax=ax[0], marker="o", title="7-day surge: total ammo (actual vs forecast)")
piv = base.groupby(["unit", "day"])["alert"].sum().unstack()
im = ax[1].imshow(piv.values, cmap="Reds", aspect="auto"); ax[1].set_yticks(range(len(piv))); ax[1].set_yticklabels(piv.index)
ax[1].set_xticks(range(piv.shape[1])); ax[1].set_xticklabels(range(1, piv.shape[1] + 1)); ax[1].set_xlabel("surge day")
ax[1].set_title("Alerts per unit/day (items with DOS<2)"); plt.colorbar(im, ax=ax[1])
ax[2].bar(range(len(res)), res["정확도"]); ax[2].axhline(0.95, c="r", ls="--"); ax[2].set_ylim(0.8, 1.0)
ax[2].set_xticks(range(len(res))); ax[2].set_xticklabels(["normal", "missing", "delay", "spike", "mixed"]); ax[2].set_title("Alert accuracy by fault scenario")
plt.tight_layout(); plt.show()

실행 결과 예시

[7일 고강도 모의] 부대 8 x 품목 6 x 7일 = 336건, 시작 보유량 3일분(가상)
             시나리오  갱신(s)  경보  실제 부족   정확도  부족 탐지율
            정상 보고  0.041  23     29 0.964   0.690
        보고 누락 10%  0.039  23     29 0.964   0.690
        보고 지연 20%  0.037  22     29 0.961   0.655
         긴급 급증 5%  0.038  67     73 0.964   0.877
복합(누락10+지연20+급증5)  0.041  66     73 0.961   0.863
경보 정확도 최저 96.1% (목표 95% 이상: 달성)
부하 시험(부대 80개, 3360건) 갱신 0.07s, 최대 갱신 0.07s (목표 2초 이내: 달성)
품목별 첫 부족 발생(모의 n일차): ammo_C 4일차, ammo_D 5일차, fuel_jet 4일차
7일 탄약 총소모(가상 상자): 실제 4645, 예측 4633, 오차 -0.2%
경보 안전계수 상충관계(정상 보고): x1.0: 정확도 96.4%/탐지율 69%, x1.1: 정확도 94.3%/탐지율 90%, x1.2: 정확도 87.2%/탐지율 97% -> 부족 놓침을 줄이려면 안전계수를 올리고 정확도 저하를 감수한다
경보는 재보급 검토 요청이며 불출·우선순위 결정은 군수참모 결재를 거친다.

미션 9. AAR 기반 소모량 예측치와 실사격 결과 간 편차 분석 및 소모 계수 재보정

  • 난이도: 고급 / 예상 소요시간: 3.0시간
  • 코드 설명: AAR(사후검토) 자료로 기준 소모율(CSR) 계획치와 실사격·실소모 결과의 편차를 작전 단계·품목·교전 강도별로 분석해 원인을 찾고, 소모 계수를 축소 추정으로 재보정한 뒤 새 훈련 자료에서 검증하며, 비용 절감·작업시간 단축·긴급 재보급 감소 3대 효과 지표와 민감도를 산출한다.
"""AAR(사후검토) 자료로 기준 소모율(CSR) 계획치와 실사격·실소모 결과의 편차를 작전 단계·품목·교전 강도별로 분석해 원인을 찾고, 소모 계수를 축소 추정으로 재보정한 뒤 새 훈련 자료에서 검증하며, 비용 절감·작업시간 단축·긴급 재보급 감소 3대 효과 지표와 민감도를 산출한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 부대·가상 품목 일일 소모 기록 생성기 (실제 탄약 제원·소모율과 무관)
# 가상 부대 8개 x 90일. 소모량 = 가상 CSR 계수 x 병력 규모 x (교전 강도·적 저항·화력 계획 효과) x 잡음
# ---------------------------------------------------------------------------
UNITS = [f"U{i:02d}" for i in range(1, 9)]
ITEMS = ["ammo_A", "ammo_B", "ammo_C", "ammo_D", "fuel_jet", "fuel_diesel"]
AMMO, FUEL = ITEMS[:4], ITEMS[4:]
PHASES = ["prep", "attack", "defense", "stability"]
CSR = {  # 가상 계수: 병력 100명·1일 기준 소모 단위(탄약=상자, 유류=드럼)
    "prep":      [1.0, 0.6, 0.3, 0.2, 1.5, 6.0],
    "attack":    [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
    "defense":   [4.5, 3.0, 1.6, 0.8, 3.0, 7.0],
    "stability": [0.8, 0.5, 0.3, 0.2, 1.0, 5.0],
}
PHASE_RESIST = {"prep": 1.0, "attack": 3.2, "defense": 3.6, "stability": 0.8}


def make_ops(n_days=90, seed=7, surge=0):
    """가상 부대별 일일 작전 기록과 품목별 실제 소모량(정답)을 만든다. surge>0이면 마지막 surge일을 고강도 공격 단계로 바꾼다(스트레스 시험용)."""
    rng = np.random.default_rng(seed)
    rows = []
    for uid in UNITS:
        size = int(rng.integers(3, 9)) * 100
        lens = rng.integers(8, 20, 4)
        sched = ["prep"] * lens[0] + ["attack"] * lens[1] + ["defense"] * lens[2] + ["attack"] * lens[3]
        sched = (sched + ["stability"] * n_days)[:n_days]
        if surge:
            sched[-surge:] = ["attack"] * surge
        r = 1.0
        for d, ph in enumerate(sched):
            r = float(np.clip(0.7 * r + 0.3 * (4.8 if surge and d >= n_days - surge else PHASE_RESIST[ph]) + rng.normal(0, 0.35), 0, 5))
            inten = float(np.clip(r / 5 * 0.9 + rng.normal(0, 0.06), 0, 1))
            fire = int(rng.poisson(1 + 8 * inten))
            move = float(rng.gamma(2, {"prep": 4, "attack": 12, "defense": 2, "stability": 5}[ph]))
            weather = str(rng.choice(["clear", "rain", "snow"], p=[0.7, 0.2, 0.1]))
            emerg = bool(rng.random() < 0.02)
            rec = {"unit": uid, "day": d, "phase": ph, "unit_size": size, "intensity": round(inten, 3),
                   "resistance": round(r, 2), "fire_missions": fire, "move_km": round(move, 1),
                   "weather": weather, "emergency": emerg}
            base = size / 100
            for i, it in enumerate(ITEMS):
                c = CSR[ph][i] * base
                if it in AMMO:
                    c *= (0.5 + 0.6 * inten) * (1 + 0.04 * fire * (1.0 if it in ("ammo_C", "ammo_D") else 0.4))
                elif it == "fuel_jet":
                    c *= (0.6 + 0.5 * inten) + 0.03 * fire
                else:
                    c = c * (0.8 + 0.2 * (weather != "clear")) + 0.05 * move * base
                c *= rng.lognormal(0, 0.12) * (2.5 if emerg else 1.0)
                rec[it] = round(float(c), 1)
            rows.append(rec)
    return pd.DataFrame(rows)
# ---------------------------------------------------------------------------

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 저장고(H1~H3) 저장칸 배치 모델. 실제 탄약 저장 기준·안전거리와 무관한 가상 계수다.
# 저장칸: 저장고 x 열(row) x 칸(col). 동선 거리 = 반출장(입구)까지의 가상 거리(m).
# 저장군(G1~G3)은 가상 혼적 제한 그룹이며, 저장고마다 받을 수 있는 저장군이 정해져 있다.
# ---------------------------------------------------------------------------
HOUSES = {"H1": {"groups": {"G1", "G2"}, "rows": 2, "cols": 6, "dock": 15.0},
          "H2": {"groups": {"G1"}, "rows": 2, "cols": 5, "dock": 40.0},
          "H3": {"groups": {"G2", "G3"}, "rows": 2, "cols": 6, "dock": 70.0}}
ITEM_GROUP = {"ammo_A": "G1", "ammo_B": "G1", "ammo_C": "G2", "ammo_D": "G3"}
SLOT_CAP = 6                      # 1열 저장칸 가상 팔레트 수용량(2열은 +4)


def make_slots(seed=3):
    """저장칸 목록과 현재 점유(이동 불가 기존 재고) 여부를 만든다."""
    rng = np.random.default_rng(seed)
    rows = []
    for h, c in HOUSES.items():
        for r in range(c["rows"]):
            for k in range(c["cols"]):
                rows.append({"slot": f"{h}-{r + 1}{k + 1:02d}", "house": h, "dist_m": c["dock"] + 6.0 * k + 4.0 * r,
                             "cap": SLOT_CAP if r == 0 else SLOT_CAP + 4, "fixed": bool(rng.random() < 0.15)})
    return pd.DataFrame(rows)


def make_skus(demand, seed=5):
    """품목별 예측 일일 반출량(demand: 품목->상자/일)을 로트 단위 SKU로 나눈다. 반출 빈도 = 일일 반출 횟수(가상)."""
    rng = np.random.default_rng(seed)
    rows = []
    for it, dq in demand.items():
        n_lot = 5
        share = rng.dirichlet(np.ones(n_lot) * 2)
        for j in range(n_lot):
            rows.append({"sku": f"{it}-L{j + 1}", "item": it, "group": ITEM_GROUP[it],
                         "pallets": int(rng.integers(2, SLOT_CAP + 5)), "freq": round(float(dq * share[j] / 4), 2)})
    return pd.DataFrame(rows)


def feasible(sku, slot):
    return (sku["group"] in HOUSES[slot["house"]]["groups"]) and sku["pallets"] <= slot["cap"] and not slot["fixed"]


def travel(assign, skus, slots):
    """작업동선 = 반출 빈도 x 왕복 거리 합(가상 m/일)."""
    d = slots.set_index("slot")["dist_m"]
    return float(sum(2 * skus.loc[i, "freq"] * d[s] for i, s in assign.items()))


def first_fit(skus, slots):
    """현행 방식(대안 비교 기준): 입고 순서대로 비어 있는 첫 번째 적합 칸에 넣는다."""
    used, out = set(), {}
    for i, k in skus.iterrows():
        for _, s in slots.iterrows():
            if s["slot"] not in used and feasible(k, s):
                out[i] = s["slot"]; used.add(s["slot"]); break
    return out


def optimize(skus, slots, lam=8):
    """할당 문제(정수계획 특수형): 비용 = 반출 빈도 x 왕복거리 + lam x 빈 공간(공간낭비). 부적합 칸은 큰 비용."""
    cost = np.full((len(skus), len(slots)), 1e7)
    for i, k in skus.iterrows():
        for j, s in slots.iterrows():
            if feasible(k, s):
                cost[i, j] = 2 * k["freq"] * s["dist_m"] + lam * (s["cap"] - k["pallets"])
    ri, cj = linear_sum_assignment(cost)
    return {int(i): slots.loc[j, "slot"] for i, j in zip(ri, cj) if cost[i, j] < 1e7}, cost


def loc_score(assign, skus, slots):
    """위치점수(0~100): 반출 빈도 상위 SKU일수록 반출장 가까운 칸에 있는 정도(가중 평균)."""
    d = slots.set_index("slot")["dist_m"]; dmin, dmax = d.min(), d.max()
    w = skus["freq"] / skus["freq"].sum()
    return float(100 * sum(w[i] * (1 - (d[s] - dmin) / (dmax - dmin)) for i, s in assign.items()))
# ---------------------------------------------------------------------------

PRICE = {"hold": 1.0, "emerg": 25.0}   # [교육용가상자료] 과잉 1단위 보관·수송 비용 지수, 긴급 재보급 1회 비용 지수
SHRINK_K = 20                            # 재보정 축소 강도(표본 n이 작을수록 기존 계수 유지)
SPEED = 80.0                             # 지게차 가상 이동 속도(m/분)


def long_ratio(ops):
    """(부대, 일자, 품목)별 실제 소모량과 CSR 계획량, 비율."""
    out = []
    for i, it in enumerate(ITEMS):
        plan = np.array([CSR[p][i] for p in ops["phase"]]) * ops["unit_size"] / 100
        out.append(pd.DataFrame({"unit": ops["unit"], "day": ops["day"], "phase": ops["phase"], "item": it,
                                 "intensity": ops["intensity"], "plan": plan, "qty": ops[it], "emergency": ops["emergency"]}))
    d = pd.concat(out, ignore_index=True)
    return d.assign(ratio=d["qty"] / d["plan"])


aar = long_ratio(make_ops(90, seed=7))          # 지난 훈련(AAR 대상)
new = long_ratio(make_ops(90, seed=8))          # 다음 훈련(재보정 검증용, 같은 생성 규칙의 다른 가상 부대)

# 1) 편차 분석: 단계 x 품목별 실제/계획 비율 중앙값(1보다 크면 계획 과소)
piv = aar[~aar["emergency"]].groupby(["phase", "item"])["ratio"].median().unstack()[ITEMS]
print("[AAR 편차: 실제/계획 비율 중앙값(긴급 불출 제외)]")
print(piv.to_string(float_format=lambda x: f"{x:.2f}"))
band = pd.cut(aar["intensity"], [0, 0.3, 0.6, 1.0], labels=["low", "mid", "high"], include_lowest=True)
amm = aar[aar["item"].isin(AMMO)].groupby(band, observed=True)["ratio"].median()
print("탄약 비율 vs 교전 강도:", ", ".join(f"{k} {v:.2f}" for k, v in amm.items()),
      "-> 원인: CSR은 단계별 고정값이라 같은 공격 단계라도 강도가 낮은 날은 과대, 높은 날은 과소 계획된다")
em_share = aar.loc[aar["emergency"], "qty"].sum() / aar["qty"].sum()
print(f"긴급 불출 비중 {em_share:.1%} -> 계수 재보정에서 제외하고 예비 재고로 관리")

# 2) 재보정: 단계 x 품목 x 계획 교전 강도 구간별 보정 계수 = 비율 중앙값^(n/(n+K)) (표본이 적은 칸은 기존 계수 유지 쪽으로 축소)
BANDS = [0, 0.3, 0.6, 1.0]
for d in (aar, new):
    d["band"] = pd.cut(d["intensity"], BANDS, labels=["low", "mid", "high"], include_lowest=True).astype(str)
base = aar[~aar["emergency"]]
g = base.groupby(["phase", "item", "band"])["ratio"]
FACT = (g.median() ** (g.size() / (g.size() + SHRINK_K))).to_dict()
print(f"재보정 계수 {len(FACT)}개(단계 4 x 품목 6 x 강도 구간), 예: attack-ammo_A low/mid/high = "
      + "/".join(f"{FACT.get(('attack', 'ammo_A', b), 1):.2f}" for b in ("low", "mid", "high")))


def mape(y, p):
    return float(np.mean(np.abs(p - y) / y))


def apply_csr(d):
    """재보정 계획량 = 기존 CSR 계획량 x 보정 계수(해당 칸이 없으면 1.0)."""
    return d["plan"] * np.array([FACT.get(k, 1.0) for k in zip(d["phase"], d["item"], d["band"])])


aar["plan2"], new["plan2"] = apply_csr(aar), apply_csr(new)


m1, m2 = mape(new["qty"], new["plan"]), mape(new["qty"], new["plan2"])
print(f"\n[재보정 검증(다음 훈련 자료)] 일일 MAPE 기존 CSR {m1:.1%} -> 재보정 {m2:.1%} ({m1 - m2:+.1%}p 개선)")
tot = new.groupby("item")[["qty", "plan", "plan2"]].sum()
tot["기존 오차"], tot["재보정 오차"] = tot["plan"] / tot["qty"] - 1, tot["plan2"] / tot["qty"] - 1
print("품목별 총량 오차:", ", ".join(f"{k} {a:+.1%}->{b:+.1%}" for k, a, b in zip(tot.index, tot["기존 오차"], tot["재보정 오차"])))


# 3) 3대 효과 지표(가상 단가, 90일 x 부대 8 기준). 두 계획 모두 AAR 자료에서 95% 충족 안전계수를 정해 같은 서비스 수준으로 비교한다
SAFE = {c: float(np.quantile(aar.loc[~aar["emergency"], "qty"] / aar.loc[~aar["emergency"], c], 0.95)) for c in ("plan", "plan2")}
print(f"95% 충족 안전계수: 기존 x{SAFE['plan']:.2f}, 재보정 x{SAFE['plan2']:.2f} -> 같은 충족 수준에서 공급 계획량을 비교한다")


def effect(plan_col, p_hold=PRICE["hold"], p_em=PRICE["emerg"]):
    supply = new[plan_col] * SAFE[plan_col]
    excess = (supply - new["qty"]).clip(lower=0).sum()
    short = int((new["qty"] > supply).sum())          # 공급 계획 초과 소모 -> 긴급 재보급 1회
    return excess * p_hold + short * p_em, short


c1, s1 = effect("plan"); c2, s2 = effect("plan2")
ops = make_ops(90, seed=7)
demand = ops[ops["day"] >= 76][AMMO].sum().div(14).to_dict()
skus, slots = make_skus(demand), make_slots()
t_cur, t_opt = travel(first_fit(skus, slots), skus, slots), travel(optimize(skus, slots)[0], skus, slots)
kpis = pd.DataFrame([["비용(과잉 보관·수송+긴급 재보급) 지수", c1, c2, 1 - c2 / c1],
                     ["저장고 반출 작업시간(분/일)", t_cur / SPEED, t_opt / SPEED, 1 - t_opt / t_cur],
                     ["긴급 재보급 횟수(회/90일)", s1, s2, 1 - s2 / max(s1, 1)]], columns=["지표", "현행", "개선", "개선율"])
print("\n[3대 정량 효과 지표(가상 단가 기준 상대 비교값)]")
print(kpis.to_string(index=False, float_format=lambda x: f"{x:.2f}"))
print(f"효과 지표 산출 {len(kpis)}종 (목표 3대 지표: {'달성' if len(kpis) >= 3 else '미달'}) - 긴급 재보급 횟수는 같은 충족 수준에서 비슷하게 유지되는지(차이 {s2 - s1:+d}회) 확인하는 지표다")

sens = []
for ph in (0.8, 1.0, 1.2):
    for pe in (15.0, 25.0, 40.0):
        a, _ = effect("plan", PRICE["hold"] * ph, pe); b, _ = effect("plan2", PRICE["hold"] * ph, pe)
        sens.append((ph, pe, 1 - b / a))
s = pd.DataFrame(sens, columns=["보관단가 배수", "긴급 1회 비용", "비용 절감률"]).pivot(index="보관단가 배수", columns="긴급 1회 비용", values="비용 절감률")
print("\n[민감도: 비용 절감률]")
print(s.to_string(float_format=lambda x: f"{x:.1%}"))
print("도입 로드맵: 단기(1~3개월) 훈련 AAR로 계수 재보정 / 중기(4~9개월) 예측 모델·저장 배치 최적화 시범 운용 / 장기(10개월~) 군수 정보체계 연계와 정기 재학습")
print("금액·시간은 가정 단가로 계산한 상대 비교값이며 예산 근거로 쓰지 않는다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
im = ax[0].imshow(piv.values, cmap="coolwarm", vmin=0.5, vmax=1.5); ax[0].set_xticks(range(6)); ax[0].set_xticklabels(ITEMS, rotation=30)
ax[0].set_yticks(range(4)); ax[0].set_yticklabels(piv.index); ax[0].set_title("AAR actual/plan ratio (median)"); plt.colorbar(im, ax=ax[0])
x = np.arange(len(tot)); ax[1].bar(x - 0.2, tot["기존 오차"], 0.4, label="old CSR"); ax[1].bar(x + 0.2, tot["재보정 오차"], 0.4, label="recalibrated")
ax[1].set_xticks(x); ax[1].set_xticklabels(tot.index, rotation=30); ax[1].axhline(0, c="k", lw=0.8); ax[1].set_title("Total plan error on next exercise"); ax[1].legend()
plt.tight_layout(); plt.show()

실행 결과 예시

[AAR 편차: 실제/계획 비율 중앙값(긴급 불출 제외)]
item       ammo_A  ammo_B  ammo_C  ammo_D  fuel_jet  fuel_diesel
phase                                                           
attack       0.90    0.89    0.99    1.02      1.04         0.94
defense      0.96    0.99    1.05    1.09      1.10         0.85
prep         0.61    0.64    0.62    0.67      0.76         0.90
stability    0.62    0.63    0.67    0.64      0.76         0.94
탄약 비율 vs 교전 강도: low 0.62, mid 0.90, high 1.06 -> 원인: CSR은 단계별 고정값이라 같은 공격 단계라도 강도가 낮은 날은 과대, 높은 날은 과소 계획된다
긴급 불출 비중 6.0% -> 계수 재보정에서 제외하고 예비 재고로 관리
재보정 계수 60개(단계 4 x 품목 6 x 강도 구간), 예: attack-ammo_A low/mid/high = 0.89/0.87/1.02

[재보정 검증(다음 훈련 자료)] 일일 MAPE 기존 CSR 34.0% -> 재보정 14.5% (+19.5%p 개선)
품목별 총량 오차: ammo_A +12.3%->+0.2%, ammo_B +11.6%->+0.2%, ammo_C -0.4%->-4.4%, ammo_D +0.9%->-3.4%, fuel_diesel +2.1%->-4.0%, fuel_jet +0.1%->-2.7%
95% 충족 안전계수: 기존 x1.23, 재보정 x1.26 -> 같은 충족 수준에서 공급 계획량을 비교한다

[3대 정량 효과 지표(가상 단가 기준 상대 비교값)]
                    지표       현행       개선   개선율
비용(과잉 보관·수송+긴급 재보급) 지수 21556.08 19063.11  0.12
      저장고 반출 작업시간(분/일)    18.88    14.97  0.21
      긴급 재보급 횟수(회/90일)   283.00   284.00 -0.00
효과 지표 산출 3종 (목표 3대 지표: 달성) - 긴급 재보급 횟수는 같은 충족 수준에서 비슷하게 유지되는지(차이 +1회) 확인하는 지표다

[민감도: 비용 절감률]
긴급 1회 비용  15.0  25.0  40.0
보관단가 배수                   
0.8      12.6% 10.7%  8.6%
1.0      13.4% 11.6%  9.6%
1.2      13.9% 12.3% 10.4%
도입 로드맵: 단기(1~3개월) 훈련 AAR로 계수 재보정 / 중기(4~9개월) 예측 모델·저장 배치 최적화 시범 운용 / 장기(10개월~) 군수 정보체계 연계와 정기 재학습
금액·시간은 가정 단가로 계산한 상대 비교값이며 예산 근거로 쓰지 않는다.

미션 10. 탄약/유류 소모량 예측 및 전술 재보급 계획 수립 시뮬레이터 통합 완성

  • 난이도: 최고급 / 예상 소요시간: 3.5시간
  • 코드 설명: 데이터 생성·정합성 점검, 소모량 예측 모델 학습과 자동 재학습(기존·신규 모델 비교 후 교체), DOS 기반 전술 재보급 계획, 가상 저장고 배치 최적화, 군수참모 결재 대기열을 한 번에 실행하는 통합 시뮬레이터를 완성하고, 자체 점검표·운영 가이드 목차·체계 연계 메시지 예시·브리핑 요약을 출력한다.
"""데이터 생성·정합성 점검, 소모량 예측 모델 학습과 자동 재학습(기존·신규 모델 비교 후 교체), DOS 기반 전술 재보급 계획, 가상 저장고 배치 최적화, 군수참모 결재 대기열을 한 번에 실행하는 통합 시뮬레이터를 완성하고, 자체 점검표·운영 가이드 목차·체계 연계 메시지 예시·브리핑 요약을 출력한다."""
import json
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
from scipy.optimize import linear_sum_assignment

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 부대·가상 품목 일일 소모 기록 생성기 (실제 탄약 제원·소모율과 무관)
# 가상 부대 8개 x 90일. 소모량 = 가상 CSR 계수 x 병력 규모 x (교전 강도·적 저항·화력 계획 효과) x 잡음
# ---------------------------------------------------------------------------
UNITS = [f"U{i:02d}" for i in range(1, 9)]
ITEMS = ["ammo_A", "ammo_B", "ammo_C", "ammo_D", "fuel_jet", "fuel_diesel"]
AMMO, FUEL = ITEMS[:4], ITEMS[4:]
PHASES = ["prep", "attack", "defense", "stability"]
CSR = {  # 가상 계수: 병력 100명·1일 기준 소모 단위(탄약=상자, 유류=드럼)
    "prep":      [1.0, 0.6, 0.3, 0.2, 1.5, 6.0],
    "attack":    [6.0, 3.5, 2.0, 1.2, 5.0, 12.0],
    "defense":   [4.5, 3.0, 1.6, 0.8, 3.0, 7.0],
    "stability": [0.8, 0.5, 0.3, 0.2, 1.0, 5.0],
}
PHASE_RESIST = {"prep": 1.0, "attack": 3.2, "defense": 3.6, "stability": 0.8}


def make_ops(n_days=90, seed=7, surge=0):
    """가상 부대별 일일 작전 기록과 품목별 실제 소모량(정답)을 만든다. surge>0이면 마지막 surge일을 고강도 공격 단계로 바꾼다(스트레스 시험용)."""
    rng = np.random.default_rng(seed)
    rows = []
    for uid in UNITS:
        size = int(rng.integers(3, 9)) * 100
        lens = rng.integers(8, 20, 4)
        sched = ["prep"] * lens[0] + ["attack"] * lens[1] + ["defense"] * lens[2] + ["attack"] * lens[3]
        sched = (sched + ["stability"] * n_days)[:n_days]
        if surge:
            sched[-surge:] = ["attack"] * surge
        r = 1.0
        for d, ph in enumerate(sched):
            r = float(np.clip(0.7 * r + 0.3 * (4.8 if surge and d >= n_days - surge else PHASE_RESIST[ph]) + rng.normal(0, 0.35), 0, 5))
            inten = float(np.clip(r / 5 * 0.9 + rng.normal(0, 0.06), 0, 1))
            fire = int(rng.poisson(1 + 8 * inten))
            move = float(rng.gamma(2, {"prep": 4, "attack": 12, "defense": 2, "stability": 5}[ph]))
            weather = str(rng.choice(["clear", "rain", "snow"], p=[0.7, 0.2, 0.1]))
            emerg = bool(rng.random() < 0.02)
            rec = {"unit": uid, "day": d, "phase": ph, "unit_size": size, "intensity": round(inten, 3),
                   "resistance": round(r, 2), "fire_missions": fire, "move_km": round(move, 1),
                   "weather": weather, "emergency": emerg}
            base = size / 100
            for i, it in enumerate(ITEMS):
                c = CSR[ph][i] * base
                if it in AMMO:
                    c *= (0.5 + 0.6 * inten) * (1 + 0.04 * fire * (1.0 if it in ("ammo_C", "ammo_D") else 0.4))
                elif it == "fuel_jet":
                    c *= (0.6 + 0.5 * inten) + 0.03 * fire
                else:
                    c = c * (0.8 + 0.2 * (weather != "clear")) + 0.05 * move * base
                c *= rng.lognormal(0, 0.12) * (2.5 if emerg else 1.0)
                rec[it] = round(float(c), 1)
            rows.append(rec)
    return pd.DataFrame(rows)
# ---------------------------------------------------------------------------

# ---- 소모량 예측용 파생 피처: 전일까지의 실적 + 당일 작전·화력 계획(사전에 알 수 있는 값만) ----
FEATS = ["lag1", "lag2", "roll3", "roll7", "csr_base", "size100", "fire_plan", "fire_x_engaged", "resist_lag1",
         "resist_trend", "inten_lag1", "move_plan", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep",
         "item_idx", "is_fuel"]
DERIVED = ["lag1", "lag2", "roll3", "roll7", "csr_base", "fire_x_engaged", "resist_lag1", "resist_trend",
           "inten_lag1", "days_in_phase", "engaged", "ph_attack", "ph_defense", "ph_prep"]


def build_long(ops):
    """부대x일자 기록을 (부대, 일자, 품목) 행으로 펼치고 파생 피처를 만든다.
    목표값 y = log(실제 소모량 / CSR 기준량): CSR 기준량에 곱할 '보정 계수'를 예측한다. lag·roll 피처도 같은 비율(실적/기준)이다."""
    df = ops.sort_values(["unit", "day"]).reset_index(drop=True)
    g = df.groupby("unit")
    df["resist_lag1"] = g["resistance"].shift(1)
    df["resist_trend"] = df["resist_lag1"] - g["resistance"].shift(3)
    df["inten_lag1"] = g["intensity"].shift(1)
    df["days_in_phase"] = df.groupby((df["phase"] != g["phase"].shift()).cumsum()).cumcount()
    out = []
    for i, it in enumerate(ITEMS):
        csr = np.array([CSR[p][i] for p in df["phase"]])
        ratio = df[it] / (csr * df["unit_size"] / 100)
        prev = ratio.groupby(df["unit"]).shift(1)
        d = pd.DataFrame({"unit": df["unit"], "day": df["day"], "phase": df["phase"], "item": it, "item_idx": i,
                          "is_fuel": int(it in FUEL), "lag1": prev, "lag2": ratio.groupby(df["unit"]).shift(2),
                          "roll3": prev.groupby(df["unit"]).transform(lambda s: s.rolling(3).mean()),
                          "roll7": prev.groupby(df["unit"]).transform(lambda s: s.rolling(7, min_periods=3).mean()),
                          "csr_base": csr, "size100": df["unit_size"] / 100,
                          "fire_plan": df["fire_missions"], "fire_x_engaged": df["fire_missions"] * df["phase"].isin(["attack", "defense"]),
                          "resist_lag1": df["resist_lag1"], "resist_trend": df["resist_trend"], "inten_lag1": df["inten_lag1"],
                          "move_plan": df["move_km"], "days_in_phase": df["days_in_phase"],
                          "engaged": df["phase"].isin(["attack", "defense"]).astype(int),
                          "ph_attack": (df["phase"] == "attack").astype(int), "ph_defense": (df["phase"] == "defense").astype(int),
                          "ph_prep": (df["phase"] == "prep").astype(int), "emergency": df["emergency"],
                          "y": np.log(ratio), "qty": df[it]})
        out.append(d)
    return pd.concat(out, ignore_index=True).dropna(subset=["lag2", "roll3", "resist_trend"]).reset_index(drop=True)


def to_qty(d, y_pred):
    """예측 보정 계수(log)를 소모량으로 되돌린다: 기준량 x exp(예측)."""
    return np.exp(y_pred) * d["csr_base"].values * d["size100"].values


def mape(y, p):
    return float(np.mean(np.abs(np.asarray(p) - np.asarray(y)) / np.asarray(y)))

# ---- 소모량 예측 모델(앙상블) 학습·예측 ----
def fit_models(tr, seed=0):
    """GBR + 랜덤포레스트 두 모델을 학습한다(목표값은 CSR 기준량 대비 log 보정 계수)."""
    g = GradientBoostingRegressor(n_estimators=150, max_depth=3, learning_rate=0.08, subsample=0.8, random_state=seed).fit(tr[FEATS], tr["y"])
    r = RandomForestRegressor(n_estimators=120, min_samples_leaf=3, random_state=seed, n_jobs=1).fit(tr[FEATS], tr["y"])
    return g, r


def predict_qty(models, d, w=0.5):
    """두 모델의 log 예측을 가중 평균해 소모량으로 되돌린다."""
    g, r = models
    return to_qty(d, w * g.predict(d[FEATS]) + (1 - w) * r.predict(d[FEATS]))

# ---------------------------------------------------------------------------
# [교육용가상자료] 가상 저장고(H1~H3) 저장칸 배치 모델. 실제 탄약 저장 기준·안전거리와 무관한 가상 계수다.
# 저장칸: 저장고 x 열(row) x 칸(col). 동선 거리 = 반출장(입구)까지의 가상 거리(m).
# 저장군(G1~G3)은 가상 혼적 제한 그룹이며, 저장고마다 받을 수 있는 저장군이 정해져 있다.
# ---------------------------------------------------------------------------
HOUSES = {"H1": {"groups": {"G1", "G2"}, "rows": 2, "cols": 6, "dock": 15.0},
          "H2": {"groups": {"G1"}, "rows": 2, "cols": 5, "dock": 40.0},
          "H3": {"groups": {"G2", "G3"}, "rows": 2, "cols": 6, "dock": 70.0}}
ITEM_GROUP = {"ammo_A": "G1", "ammo_B": "G1", "ammo_C": "G2", "ammo_D": "G3"}
SLOT_CAP = 6                      # 1열 저장칸 가상 팔레트 수용량(2열은 +4)


def make_slots(seed=3):
    """저장칸 목록과 현재 점유(이동 불가 기존 재고) 여부를 만든다."""
    rng = np.random.default_rng(seed)
    rows = []
    for h, c in HOUSES.items():
        for r in range(c["rows"]):
            for k in range(c["cols"]):
                rows.append({"slot": f"{h}-{r + 1}{k + 1:02d}", "house": h, "dist_m": c["dock"] + 6.0 * k + 4.0 * r,
                             "cap": SLOT_CAP if r == 0 else SLOT_CAP + 4, "fixed": bool(rng.random() < 0.15)})
    return pd.DataFrame(rows)


def make_skus(demand, seed=5):
    """품목별 예측 일일 반출량(demand: 품목->상자/일)을 로트 단위 SKU로 나눈다. 반출 빈도 = 일일 반출 횟수(가상)."""
    rng = np.random.default_rng(seed)
    rows = []
    for it, dq in demand.items():
        n_lot = 5
        share = rng.dirichlet(np.ones(n_lot) * 2)
        for j in range(n_lot):
            rows.append({"sku": f"{it}-L{j + 1}", "item": it, "group": ITEM_GROUP[it],
                         "pallets": int(rng.integers(2, SLOT_CAP + 5)), "freq": round(float(dq * share[j] / 4), 2)})
    return pd.DataFrame(rows)


def feasible(sku, slot):
    return (sku["group"] in HOUSES[slot["house"]]["groups"]) and sku["pallets"] <= slot["cap"] and not slot["fixed"]


def travel(assign, skus, slots):
    """작업동선 = 반출 빈도 x 왕복 거리 합(가상 m/일)."""
    d = slots.set_index("slot")["dist_m"]
    return float(sum(2 * skus.loc[i, "freq"] * d[s] for i, s in assign.items()))


def first_fit(skus, slots):
    """현행 방식(대안 비교 기준): 입고 순서대로 비어 있는 첫 번째 적합 칸에 넣는다."""
    used, out = set(), {}
    for i, k in skus.iterrows():
        for _, s in slots.iterrows():
            if s["slot"] not in used and feasible(k, s):
                out[i] = s["slot"]; used.add(s["slot"]); break
    return out


def optimize(skus, slots, lam=8):
    """할당 문제(정수계획 특수형): 비용 = 반출 빈도 x 왕복거리 + lam x 빈 공간(공간낭비). 부적합 칸은 큰 비용."""
    cost = np.full((len(skus), len(slots)), 1e7)
    for i, k in skus.iterrows():
        for j, s in slots.iterrows():
            if feasible(k, s):
                cost[i, j] = 2 * k["freq"] * s["dist_m"] + lam * (s["cap"] - k["pallets"])
    ri, cj = linear_sum_assignment(cost)
    return {int(i): slots.loc[j, "slot"] for i, j in zip(ri, cj) if cost[i, j] < 1e7}, cost


def loc_score(assign, skus, slots):
    """위치점수(0~100): 반출 빈도 상위 SKU일수록 반출장 가까운 칸에 있는 정도(가중 평균)."""
    d = slots.set_index("slot")["dist_m"]; dmin, dmax = d.min(), d.max()
    w = skus["freq"] / skus["freq"].sum()
    return float(100 * sum(w[i] * (1 - (d[s] - dmin) / (dmax - dmin)) for i, s in assign.items()))
# ---------------------------------------------------------------------------

DOS_ALERT, DOS_TARGET, PLAN_DAY = 2.0, 3.0, 85
t_all = time.perf_counter()
checks = []

# 1) 데이터 수집·정합성 점검
ops = make_ops(90, seed=7)
ok_schema = ops["intensity"].between(0, 1).all() and (ops[ITEMS] >= 0).all().all() and not ops.duplicated(["unit", "day"]).any()
checks.append(("원천 자료 정합성(범위·음수·중복)", ok_schema))
data = build_long(ops)
print(f"[1] 가상 작전 기록 {len(ops)}건 -> 학습 행 {len(data)}개, 정합성 {'통과' if ok_schema else '실패'}")

# 2) 운영 모델(v1: 60일 이전 학습) -> 3) 재학습 파이프라인(v2: 80일 이전 학습), 최근 10일로 비교 후 교체 결정
val = data[data["day"] >= 80]
v1 = fit_models(data[data["day"] < 60])
v2 = fit_models(data[data["day"] < 80])
e1, e2 = mape(val["qty"], predict_qty(v1, val)), mape(val["qty"], predict_qty(v2, val))
base = mape(val["qty"], val["csr_base"] * val["size100"])
champ = v2 if e2 <= e1 else v1
print(f"[2] 검증(80~89일) MAPE: CSR Baseline {base:.1%}, 운영 모델 v1 {e1:.1%}, 재학습 v2 {e2:.1%} -> {'v2로 교체' if e2 <= e1 else 'v1 유지'}")
checks += [("자동 재학습 파이프라인 실행·비교", True), ("운영 모델 MAPE 25% 이하", min(e1, e2) <= 0.25)]

# 4) 전술 재보급 계획(85일 기준)
rng = np.random.default_rng(11)
day = data[data["day"] == PLAN_DAY].copy()
day["pred"] = predict_qty(champ, day)
day["on_hand"] = np.round(day["csr_base"] * day["size100"] * rng.uniform(1.0, 5.0, len(day)))
day["dos"] = day["on_hand"] / day["pred"]
day["need"] = np.maximum(0, np.ceil(DOS_TARGET * day["pred"] - day["on_hand"]))
urg = day[day["dos"] < DOS_ALERT].sort_values("dos")
print(f"[3] {PLAN_DAY}일 재보급 계획: 긴급 {len(urg)}건, 재보급 총량 탄약 {day.loc[day['is_fuel'] == 0, 'need'].sum():.0f}상자 / 유류 {day.loc[day['is_fuel'] == 1, 'need'].sum():.0f}드럼(가상)")
for _, r in urg.head(4).iterrows():
    print(f"    {r['unit']} {r['item']}: DOS {r['dos']:.1f}일, 재보급 {r['need']:.0f}")

# 5) 저장 배치 최적화(예측 탄약 반출 수요 기준)
demand = day[day["is_fuel"] == 0].groupby("item")["pred"].sum().round(1).to_dict()
skus, slots = make_skus(demand), make_slots()
cur, (opt, _) = first_fit(skus, slots), optimize(skus, slots)
tc, to = travel(cur, skus, slots), travel(opt, skus, slots)
print(f"[4] 저장 배치: 작업동선 {tc:.0f} -> {to:.0f} m/일 ({1 - to / tc:.1%} 단축), 위치점수 {loc_score(cur, skus, slots):.1f} -> {loc_score(opt, skus, slots):.1f}")
checks.append(("저장 배치 작업동선 15% 이상 단축", 1 - to / tc >= 0.15))

# 6) 결재 대기열(자동 확정 금지)
queue = [{"id": f"Q{k + 1:02d}", "type": "긴급 재보급", "unit": r["unit"], "item": r["item"], "qty": int(r["need"]), "status": "결재대기"} for k, (_, r) in enumerate(urg.iterrows())]
queue += [{"id": f"Q{len(queue) + 1:02d}", "type": "저장 재배치", "moves": int(sum(cur.get(i) != opt[i] for i in opt)), "status": "결재대기"}]
auto = [q for q in queue if q["status"] != "결재대기"]
print(f"[5] 군수참모 결재 대기 {len(queue)}건, 자동 확정 {len(auto)}건")
checks.append(("자동 확정 0건(HITL)", not auto))

# 7) 운영 가이드 목차·연계 인터페이스
GUIDE = {"1. 목적·범위": "가상 부대 탄약·유류 소모 예측과 저장 배치 지원",
         "2. 데이터 규격": f"작전 기록 {ops.shape[1]}개 필드, 품목 {len(ITEMS)}종",
         "3. 전처리 규칙": "범위·음수·중복 점검, 긴급 불출 보존",
         "4. 예측 모델": f"GBR+RF 앙상블, 검증 MAPE {min(e1, e2):.1%}",
         "5. DOS 재보급 엔진": f"경보 {DOS_ALERT:.0f}일 미만, 목표 {DOS_TARGET:.0f}일",
         "6. 저장 배치": f"할당 최적화, 동선 {1 - to / tc:.0%} 단축",
         "7. 결재선·감사": "자동 확정 금지, 사유 기록",
         "8. 재학습": "최근 자료로 재학습 후 검증 MAPE가 낮을 때만 교체",
         "9. 예외 대응": "보고 누락 시 CSR 기준값 대입, 긴급 요청은 즉시 경보",
         "10. 연계 인터페이스": "RESUPPLY_REQ 가상 메시지 규격"}
done = np.mean([bool(v) for v in GUIDE.values()])
print("[6] 운영 가이드 목차:", " / ".join(GUIDE))
checks.append(("운영 가이드 10개 장 완비", done == 1))
msg = {"msg_type": "RESUPPLY_REQ", "version": "1.0", "unit_id": queue[0].get("unit", "U00") if queue else "U00",
       "item_code": queue[0].get("item", "-") if queue else "-", "qty": queue[0].get("qty", 0) if queue else 0,
       "dos_days": round(float(urg["dos"].iloc[0]), 2) if len(urg) else None, "model_ver": "v2" if champ is v2 else "v1",
       "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, ok in checks:
    print(f"  [{'PASS' if ok else 'FAIL'}] {name}")
print(f"점검 {sum(ok for _, ok in checks)}/{len(checks)} 통과, 운영 가이드 완비율 {done:.0%}, 전체 실행 {elapsed:.1f}s")

print("\n[브리핑 요약]")
print(f"- 소모량 예측 MAPE {min(e1, e2):.1%}로 CSR 고정 계수({base:.1%}) 대비 오차를 줄였고, 재학습 모델은 비교 후 교체한다.")
print(f"- {PLAN_DAY}일 기준 긴급 재보급 {len(urg)}건을 DOS 순으로 정렬해 결재 대기열에 올렸다.")
print(f"- 가상 저장고 배치 최적화로 반출 작업동선을 {1 - to / tc:.0%} 줄이는 안을 제시했다(재배치는 결재 후 시행).")
print("- 모든 불출·대체 탄종·재배치는 군수참모 결재 후 집행하며, 수치는 교육용 가상 자료 기준이다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].bar(["CSR baseline", "model v1", "model v2"], [base, e1, e2], color=["gray", "tab:blue", "tab:green"])
ax[0].axhline(0.25, c="r", ls="--"); ax[0].set_title("Validation MAPE (day 80-89)")
ax[1].barh(day["unit"] + "-" + day["item"], day["dos"].clip(upper=6), color=np.where(day["dos"] < DOS_ALERT, "tab:red", "tab:gray"))
ax[1].axvline(DOS_ALERT, c="r", ls="--"); ax[1].set_title(f"Days of supply on day {PLAN_DAY}"); ax[1].tick_params(axis="y", labelsize=5)
plt.tight_layout(); plt.show()

실행 결과 예시

[1] 가상 작전 기록 720건 -> 학습 행 4176개, 정합성 통과
[2] 검증(80~89일) MAPE: CSR Baseline 46.5%, 운영 모델 v1 13.3%, 재학습 v2 14.7% -> v1 유지
[3] 85일 재보급 계획: 긴급 6건, 재보급 총량 탄약 23상자 / 유류 221드럼(가상)
    U08 fuel_diesel: DOS 1.1일, 재보급 61
    U01 fuel_diesel: DOS 1.2일, 재보급 72
    U04 ammo_A: DOS 1.6일, 재보급 4
    U06 ammo_D: DOS 1.6일, 재보급 2
[4] 저장 배치: 작업동선 1588 -> 1260 m/일 (20.6% 단축), 위치점수 60.4 -> 72.0
[5] 군수참모 결재 대기 7건, 자동 확정 0건
[6] 운영 가이드 목차: 1. 목적·범위 / 2. 데이터 규격 / 3. 전처리 규칙 / 4. 예측 모델 / 5. DOS 재보급 엔진 / 6. 저장 배치 / 7. 결재선·감사 / 8. 재학습 / 9. 예외 대응 / 10. 연계 인터페이스
[7] 체계 연계 메시지 예시(가상 규격): {"msg_type": "RESUPPLY_REQ", "version": "1.0", "unit_id": "U08", "item_code": "fuel_diesel", "qty": 61, "dos_days": 1.14, "model_ver": "v1", "approval": "PENDING", "source": "교육용가상자료"}

[자체 점검표]
  [PASS] 원천 자료 정합성(범위·음수·중복)
  [PASS] 자동 재학습 파이프라인 실행·비교
  [PASS] 운영 모델 MAPE 25% 이하
  [PASS] 저장 배치 작업동선 15% 이상 단축
  [PASS] 자동 확정 0건(HITL)
  [PASS] 운영 가이드 10개 장 완비
  [PASS] 통합 실행 60초 이내
점검 7/7 통과, 운영 가이드 완비율 100%, 전체 실행 7.1s

[브리핑 요약]
- 소모량 예측 MAPE 13.3%로 CSR 고정 계수(46.5%) 대비 오차를 줄였고, 재학습 모델은 비교 후 교체한다.
- 85일 기준 긴급 재보급 6건을 DOS 순으로 정렬해 결재 대기열에 올렸다.
- 가상 저장고 배치 최적화로 반출 작업동선을 21% 줄이는 안을 제시했다(재배치는 결재 후 시행).
- 모든 불출·대체 탄종·재배치는 군수참모 결재 후 집행하며, 수치는 교육용 가상 자료 기준이다.