197 KiB
197 KiB
🎯 [AMM02] AI 기반 탄약 저장 최적화 모델
Google Colab 실습 환경 (탄약·특수물자 / Level 3)
📌 1. 프로젝트 개요 및 배경
| 항목 | 내용 |
|---|---|
| 과정 분류 | 탄약·특수물자 (Level 3) |
| 학습 소요 시간 | 10개 미션 (5일 / 총 28.5시간) |
| 환경 사양 | Python 3.10+, numpy, pandas, matplotlib, scipy, scikit-learn |
📖 프로젝트 목표
탄약 저장고는 공간이 한정되어 있고 품목마다 함께 둘 수 없는 제약이 있어, 앞으로 얼마나 소모되고 얼마나 자주 반출될지를 알아야 입고 위치와 재보급 순서를 합리적으로 정할 수 있다. 그러나 작전 단계별 고정 기준 소모율(CSR)만으로 계획하면 교전 강도가 낮은 날은 과잉 재고가 쌓이고 높은 날은 긴급 재보급이 반복되어 저장 공간과 작업 동선이 낭비된다. 본 프로젝트는 교육용 가상 부대·가상 품목·가상 저장고 자료를 코드로 만들어, 탄약·유류 소모율 데이터 정의와 전처리, CSR 기반 Baseline과 탐색적 분석, 적 저항·화력 계획을 반영한 소모량 예측, 보급 수준(DOS) 연계 재보급 산출, 작전 지도 대시보드와 설명가능성, 저장 배치 최적화와 군수참모 결재선(HITL), 7일 고강도 스트레스 시험, AAR 기반 소모 계수 재보정을 단계별로 구현한다. 최종 목표는 예측 소모량을 근거로 재보급 소요와 저장 위치를 함께 제안해 저장 효율과 작업 동선을 개선하고, 불출·대체 탄종·재배치는 군수참모가 근거를 확인해 결정하도록 지원하는 통합 시뮬레이터를 완성하는 것이다.
⚠️ [교육용 가상 데이터셋 원칙]
- 본 노트북의 모든 데이터(제원, 로그, 시나리오 등)는 보안 규정을 준수하여 코드 내 난수 기반 가상 데이터로 생성됩니다.
- 실제 부대 위치, 작전 인가, 실제 전장 정보, 군사 기밀 데이터는 일체 포함되지 않습니다.
In [ ]:
# ============================================================================
# [Google Colab 환경 설정] 필수 라이브러리 설치 및 임포트
# ============================================================================
!pip install -q numpy pandas matplotlib scipy scikit-learn
import os
import sys
import time
import math
import json
import re
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
# 그래프 영문 폰트 기본 설정 (환경 독립성 보장)
plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
print(f"실행 환경 준비 완료 (Python {sys.version.split()[0]})")
In [ ]:
"""작전 단계별 탄약(가상 탄종 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주의: 품목 코드·계수·저장고는 모두 교육용 가상 값이다. 실제 체계 연동 시 국방 표준 군수 코드(재고번호 등) 체계에 맞춰 필드를 대응시킨다.")
In [ ]:
"""부대 규모·교전 강도·이동 거리 변수를 정제·정규화하고, 병력 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()
In [ ]:
"""교전 일수 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에서 변수 반영 모델로 고도화한다.")
In [ ]:
"""적 저항 수준(전일 실적)과 화력 운용 계획(당일 계획 화력 임무 수)을 연계한 파생 피처를 만들고, 일자별 탄약·유류 소모량 시계열 예측 모델을 시계열 교차검증으로 평가해 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()
In [ ]:
"""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()
In [ ]:
"""작전 지도(가상 격자 좌표)에 부대별 탄약·유류 예측 소모량과 보급 수준(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()
In [ ]:
"""가상 저장고의 품목 배치를 할당 문제(정수계획)로 최적화해 작업동선 비용·공간낭비·위치점수를 현행 배치와 비교하고, 재배치·긴급 재보급 우선순위·대체 탄종 불출 제안을 자동 확정하지 않고 군수참모 결재선(승인·수정·반려·무응답)에 올려 사유와 함께 해시 체인 감사 로그에 남기는 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()
In [ ]:
"""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()
In [ ]:
"""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()
In [ ]:
"""데이터 생성·정합성 점검, 소모량 예측 모델 학습과 자동 재학습(기존·신규 모델 비교 후 교체), 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()