Journal Search Engine

ISSN : 2005-0461(Print)
ISSN : 2287-7975(Online)

Journal of Society of Korea Industrial and Systems Engineering Vol.49 No.2 pp.11-22
DOI : https://doi.org/10.11627/jksie.2026.49.2.011

McCormick-Based MILP Solution for Dynamic Weapon-Target Assignment in Multi-Layer Missile Defense Systems

Joonha Jang^*

, Sangmin Lim^**†

, Kwangjin Yang^**

, Kihoon Kwak^***

, Hwajong Jin^****

, Donghyouk Shim^****

^*Department of Computer Science
^**Department of Mechanical Engineering
^***Department of Systems Engineering, Republic of Korea Air Force Academy
^****Missile Systems Air and Missile Defense, System Development and Integration, LIG D&A, South Korea

† Corresponding Author: sangmin.lim@afa.ac.kr

Received 20/04/2026 Finally Revised 22/05/2026 Accepted 25/05/2026

Abstract

This study proposes a McCormick-based mixed-integer linear programming (MILP) model for real-time dynamic weapon-target assignment (DWTA) in multi-layer missile defense systems. The proposed model simultaneously considers engagement time windows (ETWs), time-varying interception effectiveness, and inter-layer dependency based on the shoot-look-shoot (SLS) doctrine. To reflect operational engagement conditions, a K-factor-based effectiveness model is introduced by incorporating trajectory progress, engagement geometry, missile kinematics, and environmental conditions. Since the survival probability formulation results in nonlinear and nonconvex multilinear terms, the problem is reformulated into an MILP framework using binary-tree decomposition and McCormick linearization. In addition, a rolling horizon strategy and warm-start mechanism are applied to support real-time decision making under dynamically changing battlefield conditions. Numerical experiments on various threat scenarios demonstrate that the proposed approach significantly reduces computation time compared to nonlinear optimization methods while maintaining high solution quality. In particular, the proposed method achieves near optimal solutions within sub-second computation time even for large-scale instances. The results indicate that the proposed McCormick-based MILP framework provides an effective and practical solution for real-time DWTA problems in multi-layer missile defense systems.

Key Words : Ballistic Missile Defense , Dynamic Weapon-Target Assignment , McCormick Linearization , Mixed Integer Linear Programming , Multi-layer Defense

McCormick 기반 혼합정수선형계획법을 이용한 다층 미사일 방어체계의 동적 무기-표적 할당 해법

장준하^*, 임상민^**†, 양광진^**, 곽기훈^***, 진화종^****, 심동혁^****

^*공군사관학교 컴퓨터과학과
^**공군사관학교 기계공학과
^***공군사관학교 시스템공학과
^****LIG 디펜스&에어로스페이스 미사일 시스템 대공체계연구소

초록

키워드 :

This article has been cited by 0 article in crossref

Cited-By

Funding:

1. 서 론

최근 탄도미사일, 순항미사일 및 무인기 기반 공격이 급격히 증가함에 따라, 이를 효과적으로 대응하기 위한 미사일 방어체계의 중요성이 전 세계적으로 부각되고 있다. 특히 현대 전쟁에서는 단거리부터 대륙간에 이르는 다양한 유형의 미사일 위협이 동시에 존재하며, 극초음속 활공체와 같은 고속·기동형 무기체계의 등장으로 기존 방어체계의 대응이 점점 어려워지고 있다. 극초음속 무기는 높은 속도와 기동성을 동시에 가지기 때문에 탐지, 추적 및 요격이 매우 어려운 것으로 알려져 있다[24]. 또한 최근 전장 환경에서는 다수의 미사일 및 무인기를 동시에 투입하여 방어체계의 대응 능력을 초과시키는 포화 공격 전술이 확산되고 있으며, 이는 방어 자원의 급격한 소모와 의사결정 부담 증가를 유발한다[3, 5]. 이러한 변화는 단일 요격체계 중심의 대응 방식에서 벗어나, 다층 미사일 방어체계를 기반으로 한 통합적이고 동적인 자원 배분 전략이 필수적임을 보여준다.

다층 미사일 방어체계는 서로 다른 고도 및 사거리에서 작동하는 방어체계를 계층적으로 구성하여, 다양한 위협에 대해 중첩된 방어를 수행하는 구조이다. 이러한 체계는 상층에서의 요격 실패가 하층 교전에 직접적인 영향을 미치는 구조적 특성을 가지며, 교전 가능 시간 또한 표적의 궤적에 따라 지속적으로 변화하는 특징이 있다. 따라서 다층 방어 환경에서는 단순한 정적 자원 배분이 아닌, 시간에 따라 변화하는 교전 조건과 계층 간 상호의존성을 동시에 고려하는 동적 의사결정이 요구된다.

이러한 문제는 일반적으로 무기-표적 할당(Weapon-Target Assignment, WTA) 문제로 모델링되며, 제한된 무기 자원을 다수의 표적에 최적으로 배분하는 조합최적화 문제로 정의된다. 특히 표적 수가 증가함에 따라 탐색 공간이 지수적으로 증가하는 특성을 가지므로, 대규모 WTA 문제를 실시간으로 해결하는 것은 매우 어렵다.

기존 연구는 주로 단일 시점에서의 자원 배분을 가정하는 정적 무기-표적 할당(Static Weapon-Target Assignment, SWTA) 문제에 집중했다. 그러나 실제 전장 환경에서는 표적의 이동, 교전 가능 시간의 변화, 그리고 교전 결과에 따른 상태 변화가 지속적으로 발생하므로, 이를 반영한 동적 무기-표적 할당(Dynamic Weapon-Target Assignment, DWTA) 문제로의 확장이 필수적이다. 이에 따라 메타휴리스틱 및 비선형 최적화 기반 접근이 제안되어 왔으나, 메타휴리스틱 기법은 전역 최적성을 보장하지 못하고, 비선형 최적화 기법은 계산 시간이 급격히 증가하여 실시간 적용이 어렵다는 한계를 가진다. 또한 기존 연구는 교전 가능 시간 범위, 시간 의존적 요격 성능 변화, 그리고 다층 구조에서의 계층 간 상호의존성을 동시에 반영하는 통합 모형을 충분히 제시하지 못하였다.

이러한 한계를 극복하기 위해 본 연구에서는 다층 미사일 방어체계 환경에서의 DWTA 문제를 대상으로, 시간 의존적 교전 특성과 계층 간 상호작용을 동시에 반영하는 통합 수리모형을 제안한다. 특히 표적의 비행 궤적에 따라 변화하는 교전 가능 시간 범위(Engagement Time Window, ETW)를 반영하고, 교전 시점에 따라 요격 성능이 변화하는 교전 품질 계수(K-factor) 기반 성능 모델을 도입하였다. 또한 상층 교전 결과가 하층 자원 배분에 영향을 미치는 발사-확인-재발사 구조(Shoot-Look-Shoot, SLS)를 포함하여 다층 방어체계의 의사결정 구조를 반영하였다.

본 논문의 구성은 다음과 같다. 제2장에서는 WTA 분야의 연구를 조사한 후 기존 WTA 해법의 한계를 고찰하고, 제3장에서는 탄도미사일의 비행 특성과 DWTA 수리모형 및 McCormick 기반 MILP 해법을 기술한다. 제4장에서 실험 결과를 분석하고, 제5장에서 결론 및 향후 연구 과제를 제시한다.

2. 이론적 고찰

2.1 WTA 분야 관련 연구

WTA 문제는 제한된 자원을 다수의 표적에 배분하여 기대 피해를 최소화하거나 파괴 효과를 극대화하는 대표적인 조합최적화 문제로, Manne에 의해 최초로 정식화되었으며[20] 이후 Lloyd와 Witsenhausen는 해당 문제가 NP-complete임을 증명하였다[18]. 최근 연구에서는 WTA 문제를 정적 WTA와 동적 WTA로 구분하고, 동적 문제를 다시 다단계 구조와 발사-확인-재발사 구조로 세분화하고 있다[1]. 또한 수리모형 측면에서는 기본형 모형, 확률 기반 모형, 시간 제약 모형, 적대적 환경 모형으로 분류하였다. 이러한 분류는 문제에 포함되는 확률 구조, 시간적 변화, 그리고 상호작용 구조의 복잡성을 보여준다.

초기 WTA 연구는 비교적 단순한 정적 환경에서의 기본형 또는 확률기반 모형에 집중되어 발전하였다[20, 23]. 이후 최적해 알고리즘(exact algorithm) 계열에서는 분지한계법(Branch-and-Bound), 정수계획법(integer programming), 열생성(column generation) 등을 통해 정적 WTA의 전역 최적해를 구하려는 시도가 이루어졌으며, 대표적으로 Ahuja et al.[1]은 최적해 알고리즘 및 휴리스틱 알고리즘을 통해 중간 규모 WTA를 효과적으로 해결할 수 있음을 보였다. Lu and Chen[19]은 열 생성 기반의 새로운 최적해 알고리즘을 제안하여 SWTA 문제의 계산 효율을 크게 향상시켰고, Andersen et al.[2]은 기존 WTA 선형화 기법과 최적해/근사 알고리즘을 체계적으로 비교·정리하였다. 최근에는 Bertsimas et al.[4]은 분지가격절단법(Branch-and-Price-and-Cut) 기반 접근을 통해 대규모 WTA를 초단위로 해결하는 프레임워크를 제안함으로써, 최적해 방법의 확장 가능성을 한층 높였다.

한편, 계산 시간 문제를 극복하기 위해 휴리스틱 및 메타휴리스틱 기반 접근도 활발히 연구되었다. Lee et al.[16]은 유전 알고리즘을 이용하여 일반 WTA를 효율적으로 해결하는 방안을 제시하였고, Zhou et al.[25]은 제약이 있는 SWTA 문제에서 이산 입자 군집 최적화(Particle Swarm Optimization, PSO)를 적용하였다. DWTA 및 다단계 WTA에 대해서도 다양한 메타휴리스틱이 제안되었는데, Kong et al.[13]은 향상된 다목적 PSO를 통해 DWTA를 해결하였고, Chang et al.[17]은 다단계 WTA에 대해 적응형 대규모 이웃 탐색(adaptive large neighborhood search) 알고리즘을 제안하여 공격 단계별 유연성을 반영하였다.

최근에는 기계학습 기반 방법론이 등장하고 있다. Na et al.[22]은 포인터 네트워크(pointer network) 기반 강화학습을 이용하여 WTA를 해결하는 방법을 제안하였으며, 발사간격과 교전 가능 시간 범위를 포함한 문제를 학습기반 정책으로 처리하였다. 또한 일부 최근 연구들은 심층강화학습 또는 Q-러닝을 결합한 혼합형 프레임워크를 통해 다단계 무기 할당 문제를 해결하는 방법을 제시하였다[17].

또한 최근 WTA 연구는 보다 현실적인 운용 제약을 반영하는 방향으로 확장되고 있다. Kim et al.[10]은 발사기의 방향오차가 표적 요격 성능에 미치는 영향을 고려한 WTA 모델을 제시하였으며, 이 결과는 실제 운용 환경에서는 단순한 거리·확률뿐 아니라 자세 오차, 접근 방향, 체계 정렬 상태와 같은 기하학적 제약이 중요함을 보여준다. 같은 맥락에서 시간제약 WTA, 대립형 WTA 관련 연구들은 위협 평가, 시간 범위, 적대적 상호작용, 방호자산 가치 등을 포함하는 방향으로 발전해 오고 있다[1, 6].

2.2 기존 연구의 한계 및 해법

기존 WTA 해법은 크게 최적해 알고리즘 계열과 휴리스틱·메타휴리스틱 계열로 구분될 수 있다. 최적해 알고리즘 계열은 전역 최적해 탐색이 가능하다는 장점이 있으나, 문제 규모 증가에 따라 계산 시간이 급격히 증가한다[9]. 반면 휴리스틱 및 메타휴리스틱 계열은 비교적 빠른 계산이 가능하지만, 일반적으로 전역 최적성을 보장하지 못한다[7, 16]. 특히 DWTA 문제는 시간 의존 제약과 비선형 생존확률 구조를 포함하므로, 계산 효율성과 해의 품질을 동시에 확보하기가 더욱 어렵기 때문에 다음과 같은 제약이 존재한다. 첫째, 교전 가능 시간 범위의 동적 변화이다. 탄도미사일 비행 단계에 따라 교전 가능 구간이 실시간으로 변하지만, 많은 기존 연구는 이를 이산 단계 또는 고정 이진값 수준으로 단순화한다[13, 17]. 둘째, 교전 품질 계수의 연속적 열화이다. 교전창 내에서도 할당 시점이 늦어질수록 요격 확률이 감소할 수 있으나, 이를 명시적으로 반영한 DWTA 연구는 제한적이다. 셋째, 다층 체계의 계층 간 상호의존성이다. 상층 요격 결과가 하층 자원 배분에 직접 영향을 미치는 SLS 구조를 명시적으로 포함한 MILP 수준의 통합 모형은 드물고 연구가 제한적이다[23].

이에 본 연구는 다층 미사일방어체계 환경에서의 DWTA 문제를 대상으로, 교전 가능 시간 범위, K-factor, 그리고 계층 간 SLS 상호의존성을 통합 수리모형으로 구성하고, McCormick 선형화를 통해 비선형·non-convex(비볼록) 구조를 MILP 형태로 변환함으로써 전역 최적화 기반의 실시간 해법을 제시한다.

3. 다층 미사일방어체계의 동적 WTA 모형 및 해법

3.1 탄도미사일의 비행특성

탄도미사일의 비행단계는 <Figure 1>과 같이 부스트단계, 중간단계, 종말단계의 3단계로 구분된다. 추력 단계는 수직상승·프로그램 선회·등자세 선회의 3단계로 구분할 수 있다. 탄도미사일은 추력, 질량, 직경 등 기본 제원과 피치 프로그램에 따라 비행거리, 정점고도, 비행시간, 속도가 결정된다[11]. 따라서, 동일한 발사-목표 이격거리에서도 비행궤적에 따라 요격체계의 교전가능시간이 달라진다[9].

다층 DWTA 수리모형은 탄도미사일 비행궤적에 기초한 동적 모형 수립이 필수적이다. 본 연구의 대응 대상 미사일 유형 및 방어체계 제원은 <Table 1>, <Table 2>에 제시되어 있다. 각 표의 수치는 시뮬레이션 가정값이며, 동시 교전 능력은 화력통제 운용 제약을 반영한 수치이다.

3.2 DWTA 모형

본 모형의 가정은 다음과 같다. 첫째, 탄도미사일의 위험도는 예상 낙하지역의 중요도에 따라 설정한다. 둘째, 적 탄도미사일의 시간과 규모 등 공격계획은 알 수 없다. 셋째, 탄도미사일마다 상·하층 방어체계별로 요격미사일을 2기씩 할당하여 Salvo 교전을 수행하며, Shoot-Look-Shoot 방식에 따라 운용한다. 넷째, 방어체계의 격추확률, 동시교전 능력 및 유도탄 보유량은 사전에 알려져 있다. 집합 및 인덱스, 파라미터, 결정변수는 아래와 같이 정의하였다.

◦ 집합 및 인덱스
$u \in U$ : 상층방어체계
$l \in L$ : 하층방어체계
$j \in J (t)$ : 시점 t에서 아직 제거되지 않았고, 교전 대상이 되는 탄도미사일의 집합
◦ 파라미터
$a_{u j} (t), a_{l j} (t)$ : 시점 t 기준, 상/하층 방어체계 u, l이 탄도미사일 j를 교전하기 위해 추적을 시작할 수 있는 가장 이른 시각
$b_{u j} (t), b_{l j} (t)$ : 시점 t 기준, 상/하층 방어체계 u, l이 탄도미사일 j를 교전(명중)할 수 있는 최종 시각
$W_{u j} (t), W_{l j} (t)$ : 시점 t 기준, 상/하층 방어체계 u, l의 탄도미사일 j에 대한 교전 가능 시간 범위; [[math:W_uj(t)=[a_uj(t),b_uj(t)], W_lj(t)=[a_lj(t),b_lj(t)]]]
$M_{u} (t), M_{l} (t)$ : 시점 t에서 상/하층 방어체계 u, l이 보유하고 있는 유도탄 수량
$C_{u}, C_{l}$ : 상/하층 방어체계 u, l의 동시교전 능력
$T_{j}$ : 탄도미사일 j의 위험도
$P_{u}, P_{l}$ : 탄도미사일에 대한 상/하층 방어체계 u/l의 일반 격추확률
$K_{u j}, K_{l j}$ : 탄도미사일 j에 대한 상/하층 방어체계 u/l의 교전 품질 계수; $K_{u j}, K_{l j} \in [0.6, 1]$
◦ 결정변수
$x_{u j} (t), x_{l j} (t)$ : 상/하층 방어체계 u, l이 교전가능 시간 범위 $W_{u j} (t) / W_{l j} (t)$ 를 현재 시점 t의 교전계획에 포함하면 1, 아니면 0

수리모형

min \sum_{j \in J (t)} r T_{j} {\prod_{u \in U} (1 - x_{u j} (t) K_{u j} P_{u}) \prod_{l \in L} (1 - x_{l j} (t) K_{l j} P_{l})}

(1)

subject to

\sum_{u \in U} x_{u j} (t) \leq 1, \forall j \in J (t)

(2)

\sum_{l \in L} x_{l j} (t) \leq 1, \forall j \in J (t)

(3)

\sum_{j \in J (t)} x_{u j} (t) \leq M_{u} (t), \forall u \in U

(4)

\sum_{j \in J (t)} x_{l j} (t) \leq M_{l} (t), \forall l \in L

(5)

\sum_{j : θ \in W_{u j} (t)} x_{u j} (t) \leq C_{u}, \forall u \in U, \forall θ \in ℝ

(6)

\sum_{j : θ \in W_{l j} (t)} x_{l j} (t) \leq C_{l}, \forall l \in L, \forall θ \in ℝ

(7)

목적함수는 식(1)로 표현되며, 각 탄도미사일의 위험도에 방어 실패 확률을 곱한 기대값 손실의 합을 최소화한 것이다. 방어 실패 확률은 병렬신뢰도 개념을 적용하여 상·하층 살보 격추확률( $P_{u}, P_{l}$ )에 교전 품질 계수( $K_{u j}, K_{l j}$ )를 곱한 유효 격추확률을 반영하도록 설계하였다. Kij는 매 최적화 시점마다 적 탄도미사일의 현재 위치·속도로부터 계산하는 고정 파라미터로, 상세 산출방법은 3.3절에서 서술한다.

목적함수는 이진 변수 $x_{u j}, x_{l j} \in 0, 1$ 의 다항곱 $\prod_{u \in U} (1 - x_{u j} (t) K_{{u j P}_{u}}) \prod_{l \in L} (1 - x_{l j} (t) K_{{l j P}_{l}})$ 구조로 인해 비선형·비볼록 형태를 가진다. 그렇기에 최적화를 위한 선형화는 필수적이다. 식(2), (3)은 탄도미사일당 상·하층 각 최대 1개 방어체계를 할당하고, 식(4), (5)는 유도탄 재고 제약을 명시한다. 식(6), (7)은 레이다의 물리적 한계를 반영한 동시 교전능력 제약으로, 특정 시각 θ에서 유효한 교전창의 합이 포대 동시 교전 능력( $C_{u} / C_{l} = 6$ )을 초과할 수 없다[6]. 식(1)~(7)은 특정 시점 t에서의 단일 할당 문제를 정의하며, 동적 할당을 위해 본 모형은 롤링호라이즌(Rolling Horizon) 방식으로 매 타임스텝마다 재실행된다.

3.3 모형의 구현

위 수리모형에서 교전창 $W_{i j} (t)$ 와 교전 품질 계수 $K_{i j}$ 는 탄도미사일의 실제 비행궤적으로부터 산출되는 파라미터이다. 이를 위해 적 탄도미사일의 비행단계별(부스트·중간·종말) 동역학적 특성을 반영한 탄도미사일 궤적 시뮬레이터를 구현하였다. <Figure 2>는 근거리 탄도미사일과 중거리 탄도미사일의 제원을 적용한 대표 궤적을 나타낸다.

상/하층 방어체계의 교전 가능 시간 범위는 탄도미사일 궤적과 방어영역의 교차 구간을 계산하여 산출하며, 탄도미사일 탐지 및 요격미사일 비행시간을 반영한다[8]. 이때, 유효 교전창은 교전가능사거리, 교전가능고도, 궤적 진행률(Progress≥60%)의 세 조건을 동시에 충족하는 구간으로 결정된다. 유효 격추확률은 기본 살보 격추확률에 교전 품질 계수 $K_{i j}$ 를 곱하여 산출한다. $K_{i j}$ 는 시간적·기하학적·운동학적·환경적 4개 요소의 가중 기하 평균으로 산출하며 [0.6, 1.0]으로 경계를 가정한다.

K_{i j} = T^{c_{1}} \cdot G^{c_{2}} \cdot K^{c_{3}} \cdot E^{c_{4}}

$K_{i j}$ 를 구성하는 4개 요소는 다음과 같다. 시간적 요소( $T$ )는 교전창 내 잔여 시간을 반영하며, 할당 시점이 늦어질수록 감소한다.

T = k_{m i n} + (1 - k_{m i n}) \cdot \sqrt{min (t_{d u r a t i o n} / t_{o p t i m a l}, 1)}

기하학적 요소( $G$ )는 탄도미사일 비행 경로와 방어영역의 기하학적 접근각 및 교전 거리를 반영하며, 정면 요격·최적 거리(50 km)에서 최대값을 가지고 이를 벗어날수록 감소한다. 운동학적 요소( $K$ )는 표적의 속도·고도에 따른 요격 난이도를 나타내며, 고도 30 km·속도 2 km/s의 조건에서 최대값을 가지고 초고속 또는 저고도인 표적일수록 감소한다. 환경적 요소( $E$ )는 고도에 따른 대기 밀도 변화를 포함하며, 기상·전자전 등 외부 환경 조건을 반영한다. 본 연구에서는 이상 조건을 가정하여, $E = 1.0$ 으로 고정한다. 요소들의 세부 수치 기준 및 가중치( $c_{1} = 0.25, c_{2} = 0.3, c_{3} = 0.25, c_{4} = 0.20$ )는 전문가 판단에 기반한 휴리스틱 값이다. 해당 가중치 설정이 특정 환경에 과도하게 의존하는지 확인하기 위한 실험 결과는 4.6절에서 서술한다. 또한, 최적 가중치 탐색은 독립적인 연구 과제로 본 연구에서는 다루지 않는다.

<Figure 3>은 교전 가능 시간 범위에 따른 격추확률의 변화를 보여준다. 교전 가능 시간 범위가 긴 경우가 짧은 경우보다 높은 유효 격추확률을 가진다. 유효 격추확률은 기본 살보 격추확률 $P_{u}, P_{l}$ 에 교전창 품질 계수 $K_{i j}$ 를 곱하여 산출한다. 각 탄도미사일에 대한 상·하층 방어체계의 교전 가능 시간 범위 및 유효 격추확률 계산 결과는 <Table 3>과 같은 행렬 자료구조에 저장하여 수리모형에 적용한다.

3.4 해법 개발

3.4.1 선형화의 필요성

식(1)의 목적함수는 $K_{i j}$ 가 파라미터임에도 불구하고 비선형·비볼록 구조를 포함하여 표준 LP 기반 Branch-and-Bound로 직접 풀 수 없다. 이 비선형성의 발생 원천과 선형화 필요성은 아래와 같다.

$K_{i j}$ 는 매 롤링호라이즌 호출 시 계산되어 고정되는 파라미터이므로, 유효 격추확률 $p_{i j} = K_{i j} \cdot P_{i}^{(m)}$ 은 상수이다. 따라서 단일 교전의 생존 기여 인자 $q_{i j} (t) = 1 - p_{i j} (t) x_{i j} (t)$ 는 이진 변수 $x_{i j}$ 의 선형함수이다.

그러나 다층 방어체계의 생존확률은 이 $q_{i j}$ 들의 다항곱으로 표현된다.

S_{j} = \prod_{u \in U} q_{u j} (t) \cdot \prod_{l \in L} q_{l j} (t)

이 곱구조를 순차 보조 변수 $z_{k}$ 로 전개하면 상층 생존확률 $S_{j}^{U}$ (|U|=3 예시)는 <Table 4>와 같다.

$z_{1} = q_{1} j$ 는 $x_{1} j$ 의 선형 함수이므로 비선형 항이 없다. 그러나, z2 이후 각 단계에서 연속 변수 $z_{k}$ 와 연속 변수 $q_{(k + 1) j}$ 의 곱, 즉 Bilinear Term(쌍선형 형식)이 발생한다. 이 구조는 $K_{i j}$ 가 파라미터로 고정되어 있음에도 다층 방어 곱구조에 의해 필연적으로 발생한다. LP 완화에서 $x_{i j} \in [0, 1$ ]로 완화하면 $q_{i j} \in [1 - p_{i j}, 1$ ]이 연속 구간이 되어 각 단계의 곱이 비볼록 구조를 형성하여 하한이 느슨해져 탐색 트리가 급격히 증가한다. 따라서, IPOPT·BARON 등 비선형 솔버를 직접 적용하면 15개 위협 기준 30초 이상 소요되어 실시간 운용이 불가능하다[12].

3.4.2 McCormick 선형화

McCormick 선형화는 3.4.1절에서 제시한 이진 트리의 각 단계 쌍선형 형식에 적용된다. $p_{i j} = K_{i j} \cdot P_{i}^{(m)}$ 은 상수이며, 생존 기여 인자 $q_{i j} (t) = 1 - p_{i j} (t) x_{i j} (t)$ 는 $x_{i j} (t)$ 의 선형 표현식이다. 위협 j의 상층 생존확률은 보조 변수로 분해된다:

z_{1} = q_{1 j}, ⋮ z_{k + 1} = z_{k} \cdot q_{(k + 1) j}, k = 1, \dots, | U | - 1

(8)

각 단계 $z_{(k + 1)} = z_{k} \cdot q_{(k + 1)} j$ 는 두 연속변수의 곱이다. 이진변수를 $x_{i j} \in [0, 1$ ]로 완화하면 $q_{i j} \in [1 - p_{i j}, 1$ ]이 연속 구간이 되므로, McCormick[21] 선형화를 경계값 $z_{k} \in [z_{k}^{L}, z_{k}^{U}$ ]와 $q_{(k + 1) j} \in [1 - p_{(k + 1) j}, 1]$ ]에 적용하여 4개 선형 부등식으로 등가 변환할 수 있다:

ζ_{k + 1} \geq z_{k}^{L} q_{(k + 1) j} + (1 - p_{(k + 1) j}) z_{k} - z_{k}^{L} (1 - p_{(k + 1) j})

(9)

ζ_{k + 1} \geq z_{k}^{U} q_{(k + 1) j} + z_{k} - z_{k}^{U}

(10)

ζ_{k + 1} \leq z_{k}^{U} q_{(k + 1) j} + (1 - p_{(k + 1) j}) z_{k} - z_{k}^{U} (1 - p_{(k + 1) j})

(11)

ζ_{k + 1} \leq z_{k}^{L} q_{(k + 1) j} + z_{k} - z_{k}^{L}

(12)

경계값은 귀납적으로 $z_{1}^{L} = 1 - p_{1} j, z_{1}^{U} = 1; z_{k}^{L} = z_{(k - 1)}^{L} \cdot (1 - p_{k j}), z_{k}^{U} = z_{(k - 1)}^{U}$ 로 계산된다. 하층 L도 동일하게 처리한다. 등가성은 다음과 같이 확인된다. $x_{i j} \in 0, 1$ 이면 $q_{i j}$ 의 값이 정확히 결정되고, 각 단계의 McCormick 부등식(9)~(12)가 등호로 활성화(tight)되어 $ζ_{(k + 1)} = z_{k} \cdot q_{(k + 1) j}$ ]의 유계조건에서 4개 부등식이 이중선형 항의 볼록 다각형(Convex Hull)을 정확히 표현하므로[21], LP 완화 하한이 강화되어 분지한계법 탐색 트리 크기가 대폭 감소하고 전역최적해를 보장하면서도 비선형 솔버 대비 처리 속도를 크게 향상시킨다.

3.4.3 MILP 모형 구성 및 분지한계법 기반 최적화

다항곱 $\prod_{u \in U} q_{u j} \cdot \prod_{l \in L} q_{l j}$ 는 이진 트리구조로 분해하여 선형화할 수 있다. 위협 j에 대한 상층 생존 확률 $S_{j}^{U} = \prod_{u} (q_{u j})$ 를 식(8)과 같이 순차적으로 보조 변수 zk를 도입하여 처리하며, 각 단계의 쌍선형 형식 $z_{k} \cdot q_{(k + 1) j}$ 에 선형화를 반복 적용한다. 하층도 동일하게 처리한다. 최종 MILP는 이진 변수 $x_{i j}$ , 연속변수 $q_{i j}, z_{k}$ 와 선형제약(9)~(12)만으로 표준 형태가 된다. 각 타임스텝의 MILP 변수·제약 규모는 교전 가능 쌍수 E(배터리-위협 쌍 중 실제 교전 가능한 것), 배터리수 B, 위협수 J, 계층수 ℓ에 의해 결정된다. 희소 변수 생성으로 이진 변수 $x_{i j} O (E)$ , 이진 트리 보조 변수 $q_{i j} \cdot z_{k}$ 및 McCormick 부등식(9)~(12)는 각각 $O (E \cdot ℓ)$ , 용량·배정 제약은 $O (B + J)$ 로, 전체 규모는 $O (E \cdot ℓ + B + J)$ 의 준선형 구조를 유지한다. 최종 MILP 모형은 CBC 솔버(PuLP 2.x)를 이용한 분지한계법으로 풀이하며, 솔버 설정은 최적성 갭 허용치 1%, 시간 제한 10초, 병렬 스레드 4개이다. 최적성 갭 허용치 1%는 도출된 해의 목적함수 값이 전역 최적값의 1% 이내임을 보장하며, 이를 통해 실시간 운용 제약 하에서 고품질 근사 최적해를 안정적으로 획득한다. <Figure 4>는 전체 해법 흐름을 나타낸다. 비선형 DWTA 목적함수에서 출발하여 McCormick 선형화로 이중선형 항을 처리한다. 이후 이진 트리 분해로 다중선형곱을 MILP로 변환한 뒤, 롤링호라이즌 기법에 따라 1초 단위의 타임스텝(t)마다 현재 관측된 위협 상태를 기반으로 MILP 모형이 호출되어 최적해를 도출한 뒤 다음 타임스텝에서 이를 반복하는 방식으로 전체 해법흐름을 설명할 수 있다.

3.4.4 실시간 처리를 위한 구현 최적화

실시간 처리 최적화를 위하여 네 가지 방법을 활용하였다.

첫째, 양방향 인덱스를 활용하였다. 알고리즘은 시뮬레이션 초기에 포대→위협목록, 위협→포대목록의 양방향 인덱스를 구축하여, 매 타임스텝의 제약 생성 시 전체 쌍을 순회하지 않고 해당 쌍만 직접 참조하면서 불필요한 반복 연산을 제거하고 메모리 접근 효율을 높였다. 둘째, 희소 변수를 생성하여 변수를 최소화하는 방식을 택했다. 솔버는 교전 가능성 행렬을 기반으로 실제 교전이 가능한(포대, 위협) 쌍에 대해서만 결정변수를 생성한다. 교전 불가 쌍의 변수와 제약을 원천적으로 제외함으로써 문제 규모를 축소시켰다. 셋째, 교전 품질 계수 사전계산표를 활용하였다. 시뮬레이션 시작 전 모든(포대-위협) 쌍의 교전 품질 계수를 궤적 기반으로 사전 계산하여 2차원 테이블에 저장한다. $K_{i j}$ 는 탄도탄의 현재 위치·속도 상태에 따라 매 타임스텝마다 갱신된다. 솔버는 위협 수나 포대 수와 관계없이 테이블에서 단순 인덱스 참조 $O (1)$ 로 해당 값을 즉시 읽어오며, 테이블 구축 비용은 $O (B \times T)$ 이다. 마지막으로, 초기해 재활용(Warm-start)을 적용하였다. 이전 타임스텝 최적해를 초기해로 재활용하여 풀이 시간을 콜드 스타트(Cold-start) 대비 효과적으로 단축하였다. 이상 네 가지 구현 최적화 중 양방향 인덱스·희소 변수 생성에 따른 문제 규모 축소 효과는 4.2절에서, 사전계산표·초기해 재활용에 따른 풀이 시간 단축 효과는 4.4절에서 각각 분석한다.

4. 실험결과 및 분석

4.1 실험 환경 및 시나리오 설계

제안 알고리즘은 위협 규모 3개에서 100개까지 8개 시나리오에서 모의실험을 수행하였다. 실험은 Python 3.x/PuLP 2.x(CBC 솔버) 환경에서 한반도를 가상 전장으로 설정하여 수행하였으며, 실험 장비 사양은 Intel Core i7-12700(3.6 GHz, 12코어), 32 GB RAM이다.

방어체계는 <Table 5>와 같이 L-SAM 3개 포대(상층)·M-SAM 3개 포대(하층), 총 6개 포대로 모든 시나리오에서 동일하게 구성하였다. L-SAM 포대는 발사대 4기·포대당 탑재 미사일 20발, M-SAM 포대는 발사대 6기·포대당 탑재 미사일 30발이며, 포대 배치는 한반도 전구 방어 구역(서울·중부·남부)을 기준으로 하였다.

위협은 고도 55 km 이상의 탄도 궤적을 비행하는 NODONG(L-SAM 담당)과 고도 25 km 수준의 저고도 표적인 SCUD-B(M-SAM 담당)로 구성되며, 전 시나리오에서 NODONG 40%·SCUD-B 60% 비율을 유지하여 방어체계의 교전 부하 균형을 맞추었다. 시나리오별 위협 구성은 <Table 6>과 같으며, 위협 규모 3, 5, 8, 10, 15, 20, 40, 100개의 8개 그룹을 각 1,000회 반복 시뮬레이션한 평균값이다. 1,000회 반복 실행은 무작위로 생성되는 탄도미사일의 발사 원점과 낙탄점의 기하학적 분포에 따른 몬테카를로(Monte Carlo) 시뮬레이션의 통계적 유의성을 유지하기 위함이다.

4.2 문제 규모 분석

전체 가능한 조합을 모두 생성하면 변수 수는 O(포대수×위협수×방어계층수)로 증가하지만, 교전가능 시간범위를 만족하는 쌍에 한정하여 변수를 생성하면 McCormick 선형화 후 MILP의 변수·제약 수는 교전 가능 쌍의 수에 비례한다. 이에 따라 <Table 6>과 같이 위협 규모에 준선형으로 성장한다. 교전가능한(포대-위협) 쌍에 대해서만 변수를 생성하는 희소 인덱싱을 적용한 결과, STRESS100(위협 100개·포대 6개)에서 변수 수 930개, 제약 수 2,254개로 문제 규모를 축소하여 전체 조합 대비 계산 복잡도를 효과적으로 감소시킨다. 희소 변수 인덱싱이 없는 경우 동일 시나리오에서 변수 수는 포대×위협×층=6×100×2=1,200개, 제약 수는 이를 초과하여 실시간 풀이가 어려워진다. 역방향 인덱스와 결합하면 제약조건 생성 비용이 $O (n \times k)$ 에서 단일 조회로 감소하여 문제 규모 증가에도 솔버 호출 준비 시간이 일정하게 유지된다.

4.3 시나리오 분석 및 확장성 평가

4.3.1 BASELINE15 시나리오 분석

BASELINE15 시나리오는 Nodong 6발·Scud-B 9발과 L-SAM·M-SAM 각 3개 포대로 구성된다. 제안 모형은 매 초 활성 위협 $J (t)$ 에 대해 <Algorithm 1> 단계 ②~⑤로 할당을 갱신하고, 단계 ⑥에서 시간창이 도래한 위협을 상층 우선으로 교전한다. <Figure 5>는 매 타임스텝의 위협별 실행 결과(발사 ▼, 할당 ◇, 교전창, 요격 ●)를 정리한 Gantt 차트이며, 시간축 상단에 6개 국면을 함께 표기하였다. ◇ 내부의 포대 라벨(U1-U3, L1-L3)은 단계 ⑤에서 결정된 해당 시점의 할당을 가리킨다.

<Algorithm 1>

Rolling-Horizon DWTA with Warm-Start

Input:Threats $J$ , batteries $B = U \cup L$ , timestep $Δ t = 1$ s

Output:Per-threat assignment $x^{*}$ and intercept times

1: $t \leftarrow 0, x_{prev} \leftarrow \emptyset$
2:
3: $J (t) \leftarrow {j \in J : j active at t}$ ▷ step ②
4: $[a_{i j} (t), b_{i j} (t)]$ for $(i, j) \in B \times J (t)$ ▷ step ③
5: Construct MILP with $w_{i j} = x_{i j} \cdot K_{i j} (t) \cdot P_{i j}$ ▷ step ④
6: $min \sum_{j \in J (t)} T_{j} (1 - \sum_{u \in U} w_{u j}) (1 - \sum_{l \in L} w_{l j})$ …eq. (1)
7: $x^{*} \leftarrow MILP.solve (warmstart = x_{prev})$ ▷ step ⑤
8:
9: $fire (u, j)$ ; on miss→ $fire (l, j)$ [SLS]
10: $x_{prev} \leftarrow x^{*}$ ; $t \leftarrow t + Δ t$ ▷ step ①
11:
12:

실행은 6개 국면([1]~[6])으로 전개된다. [1]( $t = 5 - 52 s$ )에서는 15발의 위협이 순차 발사되며, 솔버는 각 위협이 $J (t)$ 에 진입한 직후 단계 ②~⑤로 초기 할당을 산출한다. [2]( $t = 52 - 114 s$ )는 모든 위협이 활성이나 어느 위협도 교전 영역에 진입하지 않아 시간창이 미개방된 상태이며, 이 구간에서 솔버는 단계 ②~⑤를 매 1초 재실행하면서 향후 교전을 사전 계획한다. [3]( $t = 114 - 219 s$ )에서 단거리 Scud-B의 시간창이 먼저 개방되어 교전 대기로 전환되고, [4]( $t = 219 - 233 s$ )에서 하층 M-SAM이 단계 ⑥에 따라 Scud-B 8발(T07-T12, T14, T15)을 시간창 도래 순서로 14초간 요격하여 활성 위협이 15→7발로 감소한다. [5]( $t = 233 - 289 s$ )에서는 Nodong의 시간창 개방을 대기하며 재계획이 빈번히 발생한다. 또한 불확실성 모델링에 따른 궤적 이탈로 T13은 교전 대상에서 제외된다. [6]( $t = 289 - 327 s$ )에서 상층 L-SAM이 Nodong 6발(T01-T06)을 38초간 요격하여 실 위협 14발 전량을 처리하고 시나리오가 종료된다.

상층(Nodong)과 하층(Scud-B)은 <Figure 5>의 ◇ 시퀀스에서 대조적 양상을 보인다. 상층에서는 L-SAM 3개 포대로 6발을 처리하는 자원 경합이 MILP의 매 타임스텝 재최적화를 유발하여 위협당 평균 4.0회의 재할당(T01: 6회, T03·T05: 5회)을 산출한다. 반면 하층에서는 M-SAM의 풍부한 자원과 넓은 교전창이 warm-start의 빠른 수렴을 가능케하여, 솔버는 발사 시점의 할당을 종단까지 유지한다. <Algorithm 1>은 별도의 분담 규칙 없이 각 위협의 시간창 차이만으로 단거리 Scud를 하층 M-SAM에, 장거리 Nodong을 상층 L-SAM에 자동 분담한다. 단일 계층이 1회 교전으로 모든 위협 요격에 성공함에 따라 SLS 폴백(단계 ⑥의 상층 요격 실패 시 하층의 재교전)은 작동하지 않았다. 솔버는 평균 0.38초(최대 0.85초)에 매 타임스텝 풀이를 완료하여 실시간 운용 요건을 충족한다.

4.3.2 위협 규모별 확장성 평가

확장성 검증을 위해 3개(SMALL3)에서 100개(STRESS100)에 이르는 8개 위협 규모 시나리오에 대해 성능을 평가한 결과, 모든 시나리오에서 Timeout 없이 안정적으로 동작하며, SMALL3에서 LARGE40까지 7개 시나리오에서는 요격률 100%를 달성하였다(<Table 7>). 요격률은 전체 실제 위협(궤적 이탈 제외) 중 요격에 성공한 위협의 비율이다. 특정 타임스텝에서 솔버가 해를 반환하지 못하면 직전 해를 유지하고 다음 타임스텝에서 재시도한다.

평균 풀이 시간은 위협 수 증가에 비례하여 증가하지 않으며, MEDIUM20(0.25초)에서는 오히려 감소한다. 위협 수가 증가할수록 초기해 재활용 적용률이 98% 이상으로 안정되어, 솔버는 더 나은 초기해에서 탐색을 시작한다. SMALL8(0.62초)에서 평균 풀이 시간이 일시적으로 증가하는 이유는 초기해 재활용 적용률이 98.4%로 안정 직전 단계에 있으면서도 8개 위협이 생성하는 교전창 조합의 복잡도가 상대적으로 높기 때문이다.

STRESS100은 총 100발의 위협에 대해 6개 포대(L-SAM 3·M-SAM 3)로 대응하는 포화 공격 시나리오이다. 살보 교전(2발) 기준 가용 교전 횟수는 L-SAM 30회(60발÷2)와 M-SAM 45회(90발÷2)를 합산한 75회가 물리적 최대 교전 한계이며, 알고리즘의 최적성과 무관하게 탄약 제약이 최대 75발 요격을 이론적 상한으로 결정한다. 1,000회 실행 모두 정확히 75발 요격에 수렴함으로써 알고리즘이 이 상한에 도달함을 확인한다.

요격률이 75.0%가 아닌 78.7%로 보고되는 이유는 분모 산정 방식에 기인한다. 불확실성 모델링에 따른 궤도 이탈로 자산 위협 구역을 벗어난 일부 탄도미사일을 DWTA 개입이 불필요한 비위협으로 분류하며, 평균 이탈 수 4.75발을 제외한 실제 자산 위협 수는 약 95.25발이다. 포화 공격 시나리오에서 이와 같이 물리적 상한에 일관되게 수렴하는 결과는 제안 알고리즘이 가용 자원을 최적으로 배분함을 실험적으로 뒷받침한다.

4.4 McCormick 선형화 및 초기해 재활용 효과 검증

McCormick 선형화와 초기해 재활용의 효과를 정량적으로 검증하기 위해, BASELINE15(15개 위협) 시나리오에서 동일 하드웨어 조건으로 3가지 방법을 비교하였다(<Table 8>).

McCormick 선형화는 선형계획법 완화 하한 강화를 통해 분지한계법 탐색 트리를 대폭 감소시키며, 3.4.2절의 등가 변환 증명에 의해 전역 최적해를 보전한다.

계산 비용 절감은 두 단계에서 이루어진다. 첫째, 교전 품질 계수 사전계산표 적용이다. $K_{i j}$ 를 STRESS100 시나리오에서 매 타임스텝마다 계산할 경우, 궤적·고도·운동학적 요소 산출에 약 400 ms가 소요된다. 반면 포대-위협 쌍별 사전계산표(Look-Up Table, LUT)를 구축하면 타임스텝당 조회 비용이 0.8 ms로 단축되어, 본 실험환경 기준 약 500배의 계산 비용 절감을 달성한다. 둘째, 초기해를 재활용한다. 이전 타임스텝의 최적해를 현재 타임스텝의 초기해로 재사용해서, <Table 8>과 같이 Cold-start 1.37초에서 0.38초로 약 3.6배 단축한다. 특히 위협 20개에서 40개 시나리오에서 초기해 재활용 적용률 98% 이상과 함께 최대 풀이 시간이 0.95초 이하로 유지되어 가장 큰 성능 이득을 달성한다.

4.5 McCormick 완화 강도 및 Root Node Gap 분석

McCormick 선형화의 완화 강도를 정량적으로 평가하기 위해 각 시나리오에서 Root Node Gap을 측정하였다. Root Node Gap은 분지한계법의 루트 노드에서 얻은 선형 계획(LP) 완화 목적값 $Z_{L P}^{r} o o t$ 와 최종 정수해 $Z *$ 목적값 사이의 상대 차이로 정의한다. Gap이 솔버 허용 오차(본 실험: $δ_{r e l} = 1$ ) 이하이면 CBC는 루트 노드에서 분기 없이 즉시 종료하며 전역 최적해를 반환한다. 3.4.2절에서 서술한 바와 같이 정수점 $x_{i j} \in 0, 1$ 에서 MILP는 식(1)의 전역 최적해와 동일하다. 각 시나리오별로 롤링호라이즌 전 구간에 걸쳐 Root Node Gap을 측정하였다(gapRel=1%, McCormick 완화 강도만을 분리 측정하기 위해 Warm-start 비활성화). gapRel=0%에서도 동일 실험을 진행하여, gapRel=1% 설정과 관계없이 최적해를 도출함을 검증하였다. 실험 결과는 <Table 9>와 같다.

SMALL3-MEDIUM20에서 Root Node Gap이 0.0000%로 선형완화 단계에서 이미 최적 정수해를 찾아 분기 탐색 없이 즉시 종료된다. 3.4.2절에서 논증한 바와 같이, 정수점에서 McCormick 선형화 오차가 사라지므로 MILP 최적해가 원래 비선형 문제의 전역 최적해와 동일하다. 위협 수가 증가하는 LARGE40과 STRESS100에서도 gap은 각각 0.0001%, 0.0053%로 허용 상한(1%)을 크게 하회한다.

4.6 K-factor 가중치 민감도 분석

교전 품질 계수의 가중치 설정이 특정 환경에 과도하게 의존하는지 확인하기 위해 기준 가중치와 네 가지 대안 가중치 조합을 비교하였다. 모든 경우 가중치 합은 1이 되도록 정규화하였으며, 동일한 위협 시나리오와 동일한 난수 시드를 사용하였다. 실험은 <Table 10>과 같이 설계하였다. 이때, 각 조건에서 시나리오를 500회 반복 수행하여 평균 목적함수 값, 평균 요격률, 평균 풀이 시간, 할당 패턴 변화율(연속 최적화 스텝 간 할당 집합이 바뀐 비율)을 측정하였다. 결과는 <Table 11>과 같다.

전체 평균 요격률은 94.7%로 가중치 구성 간 차이가 없다. 요격률이 100%에 미치지 못한 것은 STRESS100이 용량 포화 시나리오(<Table 7> 참조)이기 때문이며, STRESS100을 제외하면 모든 가중치에서의 요격률 100%를 달성하였다. 따라서 K-factor 가중치 변화는 요격률에 영향을 미치지 않으며, 민감도는 목적함수 값과 할당 패턴 변화율에서 나타난다. 목적함수 값은 기하학적 요소 강조( $c_{1}$ ) 구성에서 가장 높고, 운동학적 요소 강조 구성( $c_{2}$ )에서 가장 낮다. 이는 교전 거리 기반 K-factor를 과도하게 강조하면 원거리 위협에 대한 교전 효율이 저하되어 기대 손실이 증가함을 의미한다. 할당 패턴 변화율 역시 기하학적 요소 강조( $c_{1}$ ) 구성에서 73.0%로 가장 높아, 기하학적 요소가 타임 스텝 간 할당 결정의 변동성을 높임을 확인하였다. 풀이 시간은 모든 구성에서 1초 이내의 실시간 운용 요건을 충족한다. 이를 통해 제안 모형이 특정 K-factor 가중치 조합에 과도하게 의존하지 않으며, 가중치 선택과 무관하게 안정적인 성능을 유지함을 확인할 수 있다.

5. 결 론

본 연구에서는 다층 미사일방어체계 환경에서의 동적 무기-표적 할당 문제를 해결하기 위한 통합 DWTA 프레임워크를 제안하였다. 제안된 방법은 교전 가능 시간 범위, 교전 품질 계수 기반 요격 성능 변화, 그리고 계층 간 발사-확인-재발사 상호의존성을 동시에 고려함으로써 기존 WTA 및 DWTA 연구에서 충분히 반영되지 못했던 현실적 요소들을 통합적으로 모델링하였다. 비선형·비볼록 구조를 가지는 생존확률 모델을 처리하기 위해 McCormick 완화 기법을 적용하여 문제를 MILP 형태로 변환하였으며, 이진 트리 기반 분해 기법을 통해 다항 곱 구조를 체계적으로 선형화하였다. 또한 롤링호라이즌 기반 재계획 구조와 초기해 재활용 초기화, 희소 변수 생성, 양방향 인덱스, 교전 품질 계수 사전계산표와 같은 구현 최적화를 결합하여 동적으로 변화하는 전장 환경에서도 실시간 의사결정이 가능한 계산 효율성을 확보하였다.

계산 실험 결과, 제안한 방법은 3~100개 위협 시나리오에서 평균 0.22초~0.68초의 풀이 시간을 보였으며, 직접 비선형 해법 대비 최대 79배의 속도 향상을 달성하였다. 또한, 높은 해의 품질과 안정적인 계산 성능을 유지함으로써, 실시간 동적 무기-표적 할당 문제에 대한 실용적 적용 가능성을 확인하였다.

그러나 본 연구에는 몇 가지 한계가 존재한다. 첫째, 요격 확률과 교전 품질 계수가 결정론적으로 주어진 것으로 가정되었으며, 센서 오차나 불확실성을 명시적으로 반영하지 않았다. 둘째, 롤링호라이즌 구조로 인해 장기적 최적성이 일부 제한될 가능성이 있다. 셋째, 실험은 시뮬레이션 기반 시나리오를 사용하였으므로 실제 방어체계 운용 데이터에 기반한 추가 검증이 필요하다. 향후 연구에서는 강건 최적화, 확률적 위협 평가, 적대적 환경 모델링, 분산형 방어체계 확장 등을 고려할 필요가 있다.

결론적으로, 본 연구에서 제안한 McCormick 기반 MILP 해법은 다층 미사일방어 환경에서의 동적 WTA 문제를 효율적으로 해결할 수 있는 실용적인 접근법으로, 계산 효율성과 해의 품질을 동시에 확보하면서 실시간 운용체계에 대한 적용 가능성을 가진다.

Acknowledgement

This paper is the result of research supported by LIG D&A Co., Ltd.

Figure

<Figure 1>.

Ballistic Missile Flight Phases

<Figure 2>.

Ballistic Missile Trajectory Simulation

<Figure 3>.

Conceptual Diagram for Interception Probability Calculation

<Figure 4>.

McCormick-based MILP Solution Framework for DWTA

<Figure 5>.

Per-threat execution timeline of rolling-horizon DWTA for Baseline 15

Table

<Table 1>.

Classification of Ballistic Missiles by Range [11]

Type	Classification	Representative	Range (km)	Flight Time (s)	Max. Speed (Mach)
CRBM	Short range	Scud-B	300	240	4.5~5.0
MRBM	Medium range	Nodong	1,300	480~600	9.0+
IRBM	Intermediate range	—	3,000~5,500	900~1,200	12.0~15.0+

<Table 2>.

Interceptor System Specifications

System	Altitude (km)	Range (km)	Simultaneous Engagement	Pk (Single)	Salvo Pk
CRBM	40~70	150~300	6 targets/battery	0.85	0.977
MRBM	15~40	5~50	6 targets/battery	0.80	0.96

<Table 3>.

An Example of Time Window and Interception Probability Matrix for Weapon system vs. Ballistic Missiles

Target Weapon	Target Ballistic Missiles
WeaponSystem	BatteryU1	t: 5-45 sec Kij: 0.85 Pk: 0.83	t: 52-78 sec Kij: 0.88	No Window
BatteryU2	t: 8-35 sec Kij: 0.65	t: 48-82 sec Kij: 0.88 Pk: 0.86	No Window
BatteryL1	No Window	No Window	t: 55-95 sec Kij: 0.95 Pk: 0.91
BatteryL2	No Window	t: 50-76 sec Kij: 0.78	t: 60-90 sec Kij: 0.82
BatteryL3	No Window	No Window	t: 63-93 sec Kij: 0.68

<Table 4>.

Stage-wise Nonlinear Structure in Binary Tree Decomposition

Stage	Expression	Nonlinear Structure
z1=q1j=1−p1jx1j	Linear in x1j (p1j constant)	None
z2=z1⋅q2j	Continuous z1 × Continuous q2j	Bilinear Term
z3=z2⋅q3j	Continuous z2 × Continuous q3j	Bilinear Term

<Table 5>.

The Configuration and Scale of the Defense System for Experimental Scenarios

Component	Batteries	Parameter	Value
L-SAM Battery (Upper Layer)	3	# of Launchers	4
Missiles per Battery	20
M-SAM Battery (Lower Layer)	3	# of Launchers	6
Missiles per Battery	30

<Table 6>.

Problem Instances for each Scenario

Scenario	Threats	Batteries	Avg. Variables	Avg. Constraints
SMALL3	3	6	29	84
SMALL5	5	6	48	130
SMALL8	8	6	81	211
MEDIUM10	10	6	103	266
BASELINE15	15	6	147	375
MEDIUM20	20	6	185	473
LARGE40	40	6	630	1,515
STRESS100	100	6	930	2,254

<Table 7>.

Scalability Results for Each Scenario

Scenario	Threats	Intercept Rate (%)	Avg. Time (s)	Max. Time (s)	Warm start (%)
SMALL3	3	100	0.22	0.3	95.5
SMALL5	5	100	0.38	0.4	98.3
SMALL8	8	100	0.62	0.42	98.4
MEDIUM10	10	100	0.44	0.66	98.4
BASELINE15	15	100	0.38	0.85	98.5
MEDIUM20	20	100	0.25	0.95	98.5
LARGE40	40	100	0.45	0.82	98.6
STRESS100	100	78.7	0.68	1.35	98.6

<Table 8>.

Solver Performance Comparison

Method	Avg. Time (s)	Speedup vs. IPOPT
Nonlinear solver (IPOPT, direct)	~30.00	—
McCormick MILP, Cold-start	1.37	×21.9
McCormick MILP, Warm-start (ours)	0.38	×79

<Table 9>.

Root Node Gap for each Scenario

Scenario	Z_LP^root	Z* (Best integer)	Root Node Gap (%)	Final Mip Gap (%)
SMALL3	Z*	0.97	0.00	0.00
SMALL5	Z*	1.467	0.00	0.00
SMALL8	Z*	2.333	0.00	0.00
MEDIUM10	Z*	2.890	0.00	0.00
BASELINE15	501.312	501.312	0.00	0.00
MEDIUM20	262.592	262.592	0.00	0.00
LARGE40	836.599	836.599	0.0001	0.01
STRESS100	3869.46	3869.46	0.0053	0.05

<Table 10>.

K-factor Weight Parametric Analysis Scenario Instances

Case	c1	c2	c3	c4	Description
Base	0.3	0.25	0.25	0.2	Baseline (this study)
Time-heavy	0.45	0.2	0.2	0.15	cT emphasized
Geo-heavy	0.2	0.4	0.25	0.15	cK emphasized
Kinetic-heavy	0.2	0.25	0.4	0.15	cE emphasized
Equal	0.25	0.25	0.25	0.25	Equal Weights

<Table 11>.

K-factor Weight Parametric Analysis Results

Case	Intercept Rate (%)	Avg. Objective Value	Avg. Time (s)	Assignment Change Rate (%)
Base	94.7	620.6	0.25	66.9
Time-heavy	94.7	626.28	0.2	70.6
Geo-heavy	94.7	651.26	0.25	73.0
Motion-heavy	94.7	616.98	0.4	68.5
Equal	94.7	634.45	0.25	72.2

Reference

Ahuja, R.K., Kumar, A., Jha, K.C., and Orlin, J.B., Exact and Heuristic Algorithms for the Weapon-Target Assignment Problem, Operations Research, 2007, Vol. 55, No. 6, pp. 1136-1146.
Andersen, A.C., Pavlikov, K., and Toffolo, T.A.M., Weapon-Target Assignment Problem: Exact and Approximate Solution Algorithms, Annals of Operations Research, 2022, Vol. 312, No. 2, pp. 581-606.
Arquilla, J., The End of War as We Knew It: Insurgency, Counterinsurgency and Lessons from Saturation Attack, Naval Postgraduate School Review, 2010.
Bertsimas, D., Farias, V.F., and Trichakis, N., Solving Large-Scale Weapon Target Assignment Problems in Seconds Using Branch-Price-and-Cut, Naval Research Logistics, 2025.
Boban, I., Requirements of Modern Air Defence Systems in Contemporary Conflicts, Journal of Defense & Security Analysis, 2025, Vol. 12, No. 1, pp. 45-62.
Cha, S.B., Park, S.H., Jung, J.H., Kim, W.C., and Choi, I.O., Efficient Dwell Time Allocation Method for Multi-Function Radar Resource Management, Journal of KIIT, 2023, Vol. 21, No. 5, pp. 91-100.
Chang, X., Shi, J., Luo, Z., and Liu, Y., Adaptive Large Neighborhood Search Algorithm for Multi-Stage Weapon Target Assignment Problem, Computers & Industrial Engineering, 2023, Vol. 181, 109303.
Hong, S.W., Song, J.Y., and Chang, Y.K., Analysis on Time Performance of Intercept System for Engagement Plan of Missile Defense System, Journal of the KIMST, 2019, Vol. 22, No. 1, pp. 93-105.
Im, J.S., Yoo, B.C., Kim, J.H., and Choi, B.W., A Study of Multi-to-Majority Response on Threat Assessment and Weapon Assignment Algorithm: by Adjusting Ballistic Missiles and Long-Range Artillery Treat, Journal of Korean Society of Industrial and Systems Engineering, 2021, Vol. 44, No. 4, pp. 43-52.
Kim, J.E., Lee, C.H., and Yi, M.Y., A Study on the Weapon-Target Assignment Problem Considering Heading Error, International Journal of Aeronautical and Space Sciences, 2024, Vol. 25, No. 3, pp. 1105-1120.
Kim, J.W. and Kwon, Y.S., Analysis of Flight Trajectory Characteristics of the MRBM by Adjusting the Angle of a Flight Path, Journal of the KIMST, 2015, Vol. 18, No. 2, pp. 173-180.
Kline, A., Ahner, D., and Hill, R., The Weapon-Target Assignment Problem, Computers & Operations Research, 2019, Vol. 105, pp. 226-236.
Kong, L., Chang, J., Yang, J., and Zhang, X., Solving the Dynamic Weapon Target Assignment Problem by an Improved Multiobjective Particle Swarm Optimization Algorithm, Applied Sciences, 2021, Vol. 11, No. 19, 9254.
Kwak, K.H., Lee, J.Y., and Jung, C.Y., The Optimal Allocation Model for SAM Using Multi-Heuristic Algorithm: Focused on Aircraft Defense, Journal of the Korean Operations Research and Management Science Society, 2009, Vol. 34, No. 4, pp. 43-55.
Lee, J.Y. and Kwak, K.H., The Optimal Allocation Model for SAM Using Multi-Heuristic Algorithm: Focused on Theater Ballistic Missile Defense, IE Interfaces, 2008, Vol. 21, No. 3, pp. 262-273.
Lee, Z.J., Su, S.F., and Lee, C.Y., Efficiently Solving General Weapon-Target Assignment Problem by Genetic Algorithms, World Wide Web, 2003, Vol. 6, No. 3, pp. 261-278.
Li, J., Wu, G., and Wang, L., A Comprehensive Survey of Weapon Target Assignment Problem: Model, Algorithm, and Application, Engineering Applications of Artificial Intelligence, 2024, Vol. 137, 109212.
Lloyd, S.P. and Witsenhausen, H.S., Weapons Allocation is NP-Complete, Proceedings of the 1986 Summer Computer Simulation Conference, 1986, pp. 1054-1058.
Lu, Y. and Chen, D.Z., A New Exact Algorithm for the Weapon-Target Assignment Problem, Omega, 2021, Vol. 98, 102138.
Manne, A.S., A Target-Assignment Problem, Operations Research, 1958, Vol. 6, No. 3, pp. 346-351.
McCormick, G.P., Computability of Global Solutions to Factorable Nonconvex Programs: Part I — Convex Underestimating Problems, Mathematical Programming, 1976, Vol. 10, No. 1, pp. 147-175.
Na, H., Ahn, J., and Moon, I.C., Weapon-Target Assignment by Reinforcement Learning with Pointer Network, Journal of Aerospace Information Systems, 2023, Vol. 20, No. 1, pp. 53-59.
Soland, R.M., Optimal Weapon Allocation with Integer Programming, Operations Research, 1973.
Wright, D. and Tracy, C., Hypersonic Weapons: Background and Issues for Congress, Congressional Research Service Report, 2023.
Zhou, Y., Li, X., Zhu, Y., and Wang, W., A Discrete Particle Swarm Optimization Algorithm Applied in Constrained Static Weapon-Target Assignment Problem, Proceedings of the 12th World Congress on Intelligent Control and Automation (WCICA), 2016, pp. 3118-3123.