가변 범위 커버링 문제에 대한 분지평가 알고리즘

1. 서론 및 기존 연구

설비위치문제(Facility Location Problem, FLP)는 통신 및 물류 분야에서 많이 발생하는 최적화 문제이다(Laporte et al., 2015). 가장 많이 알려진 FLP중에는 p-median 문제가 있다(Cooper 1963). p-median 문제는 주어진 노드들의 수요(demand)를 p개의 설비들 중에 하나의 설비에서 제공하는 것을 제약으로 하며 이 때의 비용을 최소화 하는 문제이다. 또 하나의 잘 알려진 FLP로는 용량 제약이 없는 설비 위치 문제(Uncapacitated Facility Location Problem, UFLP)이 있다(Balinski, 1965). 이 문제는 설비의 용량이 무한하다고 가정하고 하나의 고객 노드를 하나의 설비에 할당할 때 설비의 설비 비용을 최소화하는 문제이다. p-median와 UFLP는 모두 NP-hard에 속하는 문제이며 따라서 다항시간 알고리즘이 알려져 있지 않다. 이 외에도 매우 다양한 FLP 변형 문제들이 존재한다. 이러한 변형 문제들은 특정 분야의 특수한 제약을 고려하거나 또는 경로 문제 등과 결합하여 위치-경로 문제로 확장되기도 한다.

FLP와 유사한 문제 중에 커버링 위치 문제(covering location problems, CLP)가 있다. CLP와 FLP는 많은 유사점이 있으나 CLP의 가장 중요한 특징은 커버링 거리(covering distance) 혹은 유효 반경(covering radius)을 가진다는 것이다. 예를 들어 응급차량(ambulance)의 위치는 응급차량이 도달할 수 있는 한계 시간을 고려해야 하며 이는 결국 응급차량 설비의 위치가 커버하는 영역이 거리상으로 한계가 있음을 의미한다. 일반적으로 CLP는 FLP의 일반화(generalization) 문제로 볼 수 있으며 따라서 다양한 FLP 변종 문제들의 CLP 버전들이 역시 존재한다. 또한 이러한 CLP 변형 문제들은 FLP 변형 문제들과 마찬가지로 NP-hard에 속하는 문제들이다.

Berman et al.(2009)은 CLP 변형 문제로 가변 범위 커버링 문제(variable radius covering problem, VRCP)를 제시하였다. VRCP는 기존 CLP와 달리 커버링 영역 또는 유효 반경(effective covering radius)을 선택할 수 있다고 가정한다. 예를 들어 거리에 조명을 설치하고자 할 때 조명의 밝기와 높이에 따라 조명의 도달거리가 달라지게 된다. 주어진 지역을 모두 조명으로 커버링 하고자 하면 조명의 위치와 조명이 필요한 지역의 상황에 따라 다른 밝기와 높이의 조명을 설치하는 것이 더 유리할 수 있다. VRCP는 이러한 상황을 기존의 CLP가 유효 반경을 고정 값으로 주어지는 것에 비해 결정 변수화 함으로써 모형화 한다. 설비의 유효 반경을 모형화하기 위해서는 설비 유효반경에 대한 비용을 정의해야 하며, 일반적으로 VRCP에서는 유효 반경에 대해 비-감소 함수(non-decreasing cost function)을 사용한다. Berman et al. (2009)는 VRCP를 처음 정의하고 이를 풀기 위해서 유전자 알고리즘(genetic algorithm)에 기반한 휴리스틱 알고리즘을 제시하였다. 또 한 가지 VRCP의 특징으로는 기존의 CLP 연구들과 다르게 설비의 위치 후보지가 미리 주어지지 않는 것이다. 기존의 많은 CLP 또는 FLP 연구들은 설비의 후보 위치가 미리 주어지고 후보지들 중에 임의 개의 설비 위치를 선택하였다. 반면에 VRCP에서는 설비의 위치가 공간상에 어떤 위치에서도 설치 가능하다고 가정하며, 이는 VRCP의 결정 변수가 설비들의 유효 반경을 비롯해서 위치 좌표도 포함해야 함을 의미한다. 설비의 위치의 후보가 미리 주어지지 않는 현실 문제로 인공위성 편대 커버링 문제를 들 수 있다(Beste, 1978). 인공위성을 통해 지표면을 감시하거나 서비스를 제공하기 위해서는 지표면의 일정 영역을 다수의 인공위성을 이용해서 커버링 할 필요가 있다. 이런 경우에는 인공위성의 방향과 커버링 영역은 연속적인 값으로 정의되고 따라서 지표면에서 커버링 되는 영역의 중심도 특정 장소가 아닌 연속적인 값을 가진다.

본 연구는 VRCP에 대해 최적해 알고리즘을 제시한다. VRCP에 대한 연구는 휴리스틱 알고리즘에 집중하였다(Berman et al., 2010; Berman et al., 2019; Küçükaydın and Aras, 2020; Jabalameli, 2011). 반면에 최적해 알고리즘(exact algorithm)은 주로 수리모형을 이용하여 솔버(solver)로 푸는 방법들이 연구되었다(Karatas and Eriskin, 2023; Blanco et al., 2023). 본 연구는 VRCP에 대한 최적해 알고리즘으로 분지평가법(Branch-and-Price)에 기반한 새로운 알고리즘을 제시한다. 분지평가법은 열생성(column generation)과 분지한계법(branch-and-bound)을 결합한 알고리즘으로 어려운 최적화 문제들에 대해 많이 적용되었다(Barnhart et al., 1998). 분지한계법 알고리즘은 열생성을 위한 열생성 주문제(master problem)와 열생성 부문제(column generation subproblem)로 구성된다. 이때, 주문제는 집합-커버링(set covering) 또는 집합-파티셔닝(set partitioning)으로 모형화되며 열생성 부문제는 문제에 따라 다르게 정의되어야 한다. 분지평가법 알고리즘의 성능은 열생성 부문제를 얼마나 효율적으로 풀 수 있는지에 크게 의존한다.

본 연구의 주요 기여점은 첫째, VRCP에 대해 최초로 분지평가법에 기반한 최적해 알고리즘을 제시한다. 둘째, 열생성 부문제를 효율적으로 풀기 위해 새로운 알고리즘(스케닝 알고리즘)을 제시한다. 셋째, 제시한 최적해 알고리즘과 기존 접근방법을 전산실험을 통해 비교하여 제시한 알고리즘의 우수성을 확인한다.

본 논문의 구성은 다음과 같다. 2장은 VRCP를 정확하게 정의하고 수리모형을 제시한다. 3장은 분지평가 최적해 알고리즘을 제시하고 특히 열생성 부문제를 빠르게 풀기 위한 스케닝 알고리즘을 제시한다. 4장은 전산실험을 통해 제시한 알고리즘의 성능을 평가하고, 마지막으로 5장에서 결론을 제시한다.

2. 문제정의

본 장에서는 VRCP를 정의하고 이에 대한 수리모형을 제시한다. 커버 대상(예: 고객) 노드들의 집합 N이 주어지며 고객 노드 u∈N의 2차원 평면상의 위치가 (a_i,b_i)로 주어진다고 가정한다. 주어진 공간에서 설비(facility)의 위치는 연속된 평면상에 제한 없이 결정될 수 있다. 즉, 설비의 후보지 노드 집합이 미리 주어지지 않는다. 어떤 설비 j의 위치가 (x_j, y_j)이고 이 설비의 유효 반경(covering radius)를 r_j라고 할 때 다음 조건이 만족하면 고객 노드 i는 설비 j로 “커버”된다.

이때, $$ d_{i,j}:=\sqrt{{\left(x_j-a_i\right)}^2+{\left(y_j-b_i\right)}^2}$$로 정의된다. 즉, 고객 노드와 설비 노드 간의 거리가 설비의 유효 반경보다 크지 않으면 이 고객 노드는 설비 노드에 커버된다. 일반적인 FLP의 경우 설비의 유효 반경 r_j가 상수로 주어지지만 VRCP에서는 결정 변수로 정의되며, 또한 유효 반경은 설비마다 독립적으로 결정된다. 설비를 설치하기 위해서는 고정 비용 F와 유효 반경에 의존하는 유효 반경 비용 함수 f(r_j)를 더한 비용이 필요하다. 유효 반경 비용 함수 f(r_j)는 비-감소(non-decreasing) 함수로써 설치한 설비의 유효 반경에 따른 비용 증가를 모형화 한다. 즉,

를 항상 만족한다. 유효 반경 비용 함수 f(r_j)는 위 조건을 만족하되 꼭 볼록함수(convex function)가 아닐 수 있다. 기존 연구에서는 f(r_j)를 $$ f\left(r_j\right):=r_j^2$$ 또는 $$ f\left(r_j\right):=\mathrm{l}\mathrm{o}\mathrm{g}\left(r_j+ϵ\right)$$와 같은 비선형 함수들로 많이 정의하였다. 이때, ϵ은 함수값이 무한히 작아지는 것을 방지하는 양의 상수이다. $$ f\left(r_j\right):=r_j^2$$와 같은 비용 함수는 일반적으로 커버 거리가 증가함에 따라 비용이 빠르게 증가하는 상황에서 사용하고, $$ f\left(r_j\right):=\mathrm{l}\mathrm{o}\mathrm{g}\left(r_j+ϵ\right)$$은 일반적으로 거리가 증가함에 따라 오히려 비용의 증가가 감소하는 상황에서 사용한다. 본 연구는 어느 특정 함수의 형태를 고정하지 않는다. 단지, 비-감소 가정은 커버 거리가 늘어나되 오히려 비용이 감소하는 것이 불가능한 현실을 반영한 것이다.

VRCP의 제약은 모든 고객 노드들이 하나 이상의 설비에 커버되는 것이다. 목적함수는 총 설치 비용의 합으로 정의된다. 먼저 수리모형에 사용될 집합과 결정 변수를 정의한다. 고객 노드 집합이 N으로 주어졌을 때 설비는 최대 N개가 필요하다. 따라서 설비의 집합을 마찬가지로 N으로 정의할 수 있다. 결정 변수 z_j는 설비 j∈N이 설치되면 1 아니면 0의 값을 가진다. 상수 F는 설비가 하나 설치될 때 마다 발생하는 고정비(fixed cost)를 나타낸다. 결정 변수 w_ij는 고객 i∈N이 설비 j∈N에 커버되면 1 아니면 0의 값을 가진다. 설비 j∈N의 유효 반경은 결정 변수 r_j에 의해 결정된다. 고객 노드 i∈N와 설비 j∈N의 거리는 결정 변수 d_ij로 표현된다. 설비 j의 위치는 결정 변수 x_j와 y_j로 정해진다. VRCP에 대한 수리모형을 다음과 같이 제시한다.

목적함수 (1)은 설치할 설비의 총 고정비와 총 유효 반경 비용 함수의 합으로 정의된다. 제약식 (2)는 설비와 고객 노드 간의 거리의 하한을 제한한다. 제약식 (3)은 모든 고객 노드가 적어도 하나 이상의 설비에 할당 되어야 함을 나타낸다. 제약식 (4)는 설치된 설비만 고객을 커버할 수 있게 하며, 제약식 (5)는 설비의 유효 반경이 고객 노드와의 거리보다 작지 않을 경우에만 설비가 고객 노드를 커버 가능하게 한다. M은 매우 큰 값을 가지는 상수이다. 제약식 (6)-(9)는 각각의 결정 변수의 도메인(domain)을 나타낸다.

수리 모형 (P)는 유효 반경 비용 함수 f(r)이 볼록 함수(convex function)일 경우 정수 조건을 완화하면 볼록 최적화 문제(convex optimization problem)가 됨을 쉽게 알 수 있다. 특히, 유효 반경 비용 함수 f(r)이 선형 함수라면 정수 조건을 완화하면 다항식 시간 복잡도를 가지는 이차 콘 계획법(second-order cone programming)이 되어 일반적인 정수 조건을 완화한 분지한계법(branch-and-bound method)으로 최적해를 구할 수 있다. 그러나, 여전히 제약식 (2)의 비선형성은 선형 계획법(linear programming) 기반의 분지한계법에 비해 이차 콘 계획법 기반의 분지한계법은 (같은 다항식 시간 알고리즘이긴 하지만) 효율적이지 못할 수 있다. 더욱이, 유효 반경 비용 함수 f(r)가 비-볼록 함수(non-convex function)인 경우에는 최적해를 구하기 위해서는 더욱 복잡한 알고리즘과 특별한 솔버(solver)가 필요하다.

VRCP는 NP-hard로 알려져 있으며(Berman et al., 2009), “최소 부분집합 설비”(minimal subset facility)의 개념을 적용하면 설비의 위치 후보와 유효 반경이 주어지는 FLP로 변환할 수 있다. 유효 반경 비용 함수가 비-감소 함수라면 고객 집합 N의 부분 집합 S에 대해 가장 작은 유효 반경을 가지는 설비의 위치와 반경이 유일하게 정해진다. <Figure 1>은 같은 고객 노드 부분집합을 커버링하는 두 개의 다른 설비를 보여준다. 왼쪽 그림의 설비는 오른쪽 그림의 설비에 비해 유효 반경이 더 크고 따라서 유효 반경 비용도 더 크다. 같은 고객 노드를 커버링 하며 오른쪽 그림보다 더 작은 유효 반경을 가지는 설비는 불가능 하며 이와 같이 같은 고객 노드 집합을 커버링하며 가장 작은 유효 반경을 가지는 설비를 최소 부분집합 설비(minimal subset facility)라 부른다. 최소 부분집합 설비는 설비의 위치와 유효 반경이 유일하게 정해지며(<Figure 1> (b) 참조) 이러한 경우에는 항상 설비의 유효 범위의 경계에 두 개 또는 세 개의 고객 노드가 있음을 보일 수 있다(Berman et al., 2009). 즉, 고객 노드의 부분집합 중에 크기가 한 개, 두 개 또는 세 개인 모든 부분집합에 대해 최소 부분집합 설비가 유일하게 정해진다.

Figure 1.

Two possible facilities covering the same subset S ⊆ N

• 정리 1. Berman et al.(2009) 고객 노드가 n개 있을 때 최소 부분집합 설비는 최대 n(n² + 5)/6개 존재한다.

위 정리를 이용하면 |N|(|N|² + 5)/6개의 설비 후보지 집합을 가지는 일반적인 FLP로 변환할 수 있다. 하지만 |N|이 커지면 여전히 수리모형의 크기가 증가하여 최적해를 얻기 힘들어진다. 결과적으로 최적해를 얻기 위해서는 비선형 정수계획법 모형 (P)을 풀거나 크기가 매우 큰 선형 정수계획법 모형(FLP 로 변형한 모형)을 풀어야 한다. 이에 본 연구에서는 비교적 빠른 시간에 최적해를 얻을 수 있는 분지평가법(branch-and-price algorithm)을 제시한다.

3. 분지평가 알고리즘

분지평가 알고리즘(branch-and-price algorithm)은 열생성(column generation)과 분지한계법(branch-and-bound)을 결합한 것으로 FLP와 관련된 문제들에 많이 사용되었다(Klose andGörtz, 2007; Shu et al., 2021). 분지평가 알고리즘은 먼저 집합 커버링 문제(set covering problem)로 문제를 재모형(reformulation)하며 이때 결정 변수는 정해진 고객 노드의 부분집합을 커버하는 하나의 설비의 설치 여부를 결정하는 이진변수로 정의된다. 고객 노드의 부분집합의 수가 매우 많을 수 있기 때문에 일반적으로 집합 커버링 문제는 모든 결정 변수를 포함하지 않고 일부 결정 변수만을 포함하는 제한된 형태를 가지며 최적해를 구하는데 필요한 추가적인 결정 변수(column)를 생성하는 과정인 열생성 과정을 거치게 된다. 주어진 고객 노드 집합 N에 대한 모든 최소 부분집합 설비의 집합을 C라고 정의하자. 어떤 최소 부분집합 설비 c∈C의 유효 반경을 r_c라고 정의하고 이때 발생하는 유효 반경 비용을 f_c: f(r_c)라고 정의하자. 또한 이 설비가 고객 노드 i∈N를 커버하는지 여부를 나타내는 함수 $$ {\delta }_i^c$$를 다음과 같이 정의한다.

주어진 최소 부분집합 설비의 부분집합 $$ \widehat{C}\subseteq C$$에 대해, VRCP의 열생성 주문제(master problem)는 다음과 같이 정의된다.

목적함수 (10)은 설비의 총 설치비를 최소화함을 의미한다. 제약식 (11)은 모든 고객 노드가 하나 이상의 설치된 최소 부분집합 설비로 커버되게 한다. 제약식 (12)는 최소 부분집합 설비를 결정하는 이진변수를 정의한다. 결정 변수 u_c는 하나의 최소 부분집합 설비를 나타내며 c∈C가 커버하는 고객 노드의 집합 $$ \left\{i\in N:{\delta }_i^c=1\right\}$$과 유효 반경 r_c, 설치 위치 (x_c,y_c)는 커버하는 노드의 집합을 가장 작은 유효 반경으로 커버하는 최소 부분집합 설비에 의해 모두 미리 정해지는 고정값들이다. 문제(MP)는 모든 최소 부분집합 설비가 아닌 부분집합 $$ \widehat{C}\subseteq C$$에 대해서 정의되며, 열생성을 위한 선형완화 문제가 다음과 같이 정의된다.

제약식 (15)는 제약식 (12)에서 정수조건을 완화한 것이다. 설치 고정비 F와 변동비 f_c는 모두 음의 값을 가지지 않는다고 가정한다. 결정 변수 u_c의 상한제약(u_c ≦ 1)은 선형완화 문제 (MP-R)에서 사용하지 않았는데, 상한 제약식이 없어도 F ≧ 0과 f_c ≧ 0의 조건을 만족하면 열생성을 통해 (MP-R)이 상한 제약식의 유무에 상관없이 같은 최적해를 가진다는 것을 쉽게 보일 수 있다. 제약식 (14)에 해당하는 쌍대 변수(dual variables)를 π_i로 정의하자. (MP-R)을 풀었을 때 쌍대 변수 값을 $$ {\left({\pi }_i^{\mathrm{*}}\right)}_{i\in N}$$하고 하면, 그때 얻은 해가 최적해일 조건은 모든 결정 변수의 할인가(reduced cost)가 비음(nonnegative)이어야 함이며 다음과 같이 표현된다.

3.1 열생성 부문제(column generation subproblem)

열생성 부문제는 모든 최소 부분집합 설비 중에 할인가(reduced cost)가 음(negative)인 것을 찾는 것이며 다음과 같은 최적화 문제로 모형화 된다.

$\[ \left(\text{SP}\right) \xi ={\text{min}}_{c\in C}\left\{F+f_c-\sum _{i\in N}{\pi }_i^{\mathrm{*}}{\delta }_i^c\right\} \]$

만약에 ξ < 0이면 을 만족하는 $$ c^{\mathrm{*}}\in C$$을 찾을 수 있다. 새로운 열(column) $$ c^{\mathrm{*}}$$을 (MP-R)에 추가하고($$ \widehat{C}\leftarrow \widehat{C}\cup \left\{c^{\mathrm{*}}\right\}$$) (MP-R)을 다시 푼 다음 위 과정을 반복한다. 열생성은 ξ ≧ 0을 만족하면 종료한다.

열생성 부문제 (SP)를 풀기 위한 가장 간단한 방법은 모든 최소 부분집합 설비들에 대해 할인가를 모두 계산하는 것이다. 2장에서 언급했듯이 최소 부분집합 설비의 개수는 O(n(n²+5)/6)로 주어지며 (n := |N|), 이는 다항 시간 복잡도 함수(polynomial time complexity function)이다. 따라서 모든 최소 부분집합 설비에 대해 할인가 $$ f_c-\sum _{i\in N}{\pi }_i^{\mathrm{*}}{\delta }_i^c$$를 계산하여 가장 작은 값을 가지는 최소 부분집합 설비 $$ c^{\mathrm{*}}$$를 찾을 수 있다. 하나의 최소 부분집합 설비로는 최대 n개의 고객 노드가 커버될 수 있다. 설비 c에 의해 커버되는 고객 노드들의 집합을 N_c이라고 정의하자. 열생성 부문제를 풀기 위해서는 커버된 고객 노드들에 대한 쌍대 값의 합 $$ \sum _{i\in N_c}{\pi }_i^{\mathrm{*}}$$를 모두 더하는 과정이 필요하다. 결과적으로 열생성에는 O(n(n²+5)/6)의 다항 시간 복잡도 알고리즘이 필요하다. 참고로 이 과정은 모든 c∈C에 대해 커버되는 고객 노드들에 해당하는 쌍대 변수 값의 합을 구하는 것으로 이해할 수 있다. 이때, 집합 C는 변하지 않으나 쌍대 변수값은 열생성 과정에서 매번 변하기 때문에 쌍대 변수 값을 더하는 과정은 매번 반복해서 실행해야 한다.

(1) 스케닝 방법(scanning method)

단순하게 모든 최소 부분집합 설비에 대해 쌍대값을 구하는 방법은 비록 다항 시간 안에 풀이가 가능하나 여전히 열생성 과정에서 여러 번 실행해야 하기 때문에 실행 속도가 만족스럽지 않을 수 있다. 이 방법의 가장 큰 단점은 하나의 최소 부분집합 설비 c∈C의 할인가를 계산하는 과정에서 내부의 커버되는 같은 고객 노드들의 쌍대 값을 더하는 반복 계산이 많이 발생한다는 것이다. 먼저 최소 부분집합 설비 c∈C가 커버하는 고객 노드들의 집합을 N_c라고 정의하자(즉, $$ N_c:=\left\{i\in N:{\delta }_i^c=1\right\}$$). 예를 들어, 두 개의 인접한 최소 부분집합 설비 c₁, c₂∈C를 생각해 보자. c₁과 c₂는 커버하는 고객 노드의 집합이 비슷하여 $$ \left|N_1\setminus N_2\right|=1$$, $$ \left|N_2\setminus N_1\right|=1$$, 그리고 $$ \left|N_1\cap N_2\right|\gg 1$$를 만족한다고 가정하자. c₁과 c₂에 대한 할인가를 독립적으로 별도로 계산하면 $$ \left|N_1\cap N_2\right|$$번의 중복 계산이 발생한다. 본 절에서 제안하는 방법은 이러한 중복 계산을 줄여 열생성 부문제를 푸는 시간을 줄이는 것이 목표이며 이를 스케닝(scanning) 방법이라 부르기로 한다. 스케닝 방법의 기본 아이디어는 비슷한 최소 부분집합 설비인 c₁과 c₂의 할인가를 독립적으로 계산하지 않고 먼저 c₁의 할인가를 계산한 후 이 값을 이용해 c₂의 할인가를 쉽게 계산하는 것이다.

최소 부분집합 설비는 고객 노드 1개, 2개 또는 3개로 유일하게 정의된다. 노드 1개일 경우에는 그 노드만을 커버링 하는 노드와 같은 위치에 있는 설비가 최소 부분집합 설비가 되고 노드가 2개인 경우는 두 노드를 양쪽 끝으로 하는 설비가 최소 부분집합 설비가 된다. 또한 노드가 3개인 경우는 3개 노드가 같은 직선 위에 위치하지 않는다면 3개 노드를 모두 지나는 설비의 위치가 유일하게 하나로 정해진다. 1개 또는 2개로 정의되는 경우는 각각 n개와 n(n - 1)/2개의 최소 부분집합 설비가 존재하고, 이들의 중심점과 유효 거리는 쉽게 미리 계산할 수 있다. 따라서 이들이 커버하는 고객 노드 집합도 미리 계산할 수 있고 이를 이용해 할인가도 쉽게 계산할 수 있다. 이제 3개의 고객 노드로 정의되는 최소 부분집합 설비를 생각해 보자. 어떤 두 개의 고객 노드들에 대해서, 이들 노드를 연결하는 가상의 선을 그리고 이 선의 중심에서 이 선과 수직인 가상의 직선 l을 그린다. 이미 선택된 두 개의 노드와 다른 어떤 노드로 정의되는 최소 부분집합 설비의 중심점은 직선 l위에 있으며 이때 중심점의 좌표도 쉽게 계산할 수 있다.

<Figure 2>는 두 개의 고객 노드 4와 5를 연결하는 직선과 수직인 직선 l_4,5 위에 3개의 노드 4,5, 그리고 1로 정의되는 최소 부분집합 설비의 중심점 c₁이 위치함을 보여준다. 이때, 중심점은 중심점에서 3개의 노드 간의 거리가 동일해야 함을 이용해서 연립방정식으로 간단하게 구할 수 있다. 3개 노드의 좌표를 각각 (x₁, y₁), (x₂, y₂), 그리고 (x₃, y₃)로 정의하자. 이들 3개의 노드로 정의되는 최소 부분집합 설비의 중심점 좌표 (x,y)는 다음과 같이 주어진다.

$\[ \begin{array}{c}x=\frac{\left(x_1^2+y_1^2\right)\left(y_2-y_3\right)+\left(x_2^2+y_2^2\right)\left(y_3-y_1\right)+\left(x_3^2+y_3^2\right)\left(y_1-y_2\right)}{D}\hfill \\ y=\frac{\left(x_1^2+y_1^2\right)\left(x_3-x_2\right)+\left(x_2^2+y_2^2\right)\left(x_1-x_3\right)+\left(x_3^2+y_3^2\right)\left(x_2-x_1\right)}{D}\hfill \end{array} \]$

Figure 2.

The Center of Minimal Subset Facility Determined by Three Nodes 4,5, and 1

이때, $$ D=2\left(x_1\left(y_2-y_3\right)+x_2\left(y_3-y_1\right)+x_3\left(y_1-y_2\right)\right)$$이다. 만약 D = 0이면 3개의 노드는 하나의 직선상에 있음을 의미하며 이 경우에는 거리가 먼 두 개의 노드로 정의되는 최소 부분집합 설비에 이미 포함되어 여기서는 무시할 수 있다. 두 개 노드를 고정하고 이 두 노드를 제외한 다른 노드들 각각에 대해 위의 과정은 반복하면 최소 부분집합 설비들의 중심점들이 모두 직선 k위에 위치하게 된다.

<Figure 3>은 총 7개의 고객 노드가 있을 때 두 개의 노드 4와 5를 기준으로 생성되는 모든 최소 부분집합 설비들을 보여준다. 직선 l_4,5 위에 있는 중심점들을 왼쪽에서부터 오른쪽으로 스케닝(scanning)하면 $$ c_1\to c_2\to c_3\to c_4\to c_5$$의 순서가 되며 이 순서는 노드 4와 5를 연결하는 직선에서부터의 거리와 직선의 왼쪽/오른쪽 위치하는 여부를 이용해서 쉽게 구할 수 있다. 스케닝 알고리즘은 모든 노드 쌍에 대해 위 예제와 같은 중심점의 순서를 미리 계산하고 이를 이용해 열생성 과정에서 이용한다.

Figure 3.

The centers of minimal subset facilities determined by three nodes 4,5, and 1 are on the line l4,5

스케닝의 장점은 스케닝을 통해 최소 부분집합 설비의 할인가 계산을 빠르게 할 수 있다는 것이다. <Table 1>은 <Figure 2>의 예제에 대한 스케닝 과정을 보여준다. 직선 l_4,5 왼쪽에서부터 미리 계산된 중심점을 차례로 방문한다. <Table 1>에서 스텝(“step”) 열은 방문 순서를 나타낸다. 중심점들은 미리 순서가 정해지고 다음 스텝의 중심점으로 정의되는 최소 부분집합 설비가 커버하는 고객 노드 집합 대비 추가되는 노드(“in”)와 제외되는 노드(“out”)를 미리 쉽게 계산할 수 있다. 이때 “touching node”는 최소 부분집합 설비가 정의되는 세 번째 노드를 나타낸다. 매 스케닝 스텝마다 최대 두 개(“in” 또는 “out”)의 노드들이 발생할 수 있다. 이는 결과적으로 두 번째 스케닝 스텝에서부터는 해당하는 최소 부분집합 설비의 할인가를 최대 두 번의 연산으로 구할 수 있음을 의미한다. 예를 들어 스텝 1에서 5개의 커버되는 노드들의 할인가 총합을 계산한다. 다음 스텝에서는 커버되는 노드들 {2,3,4,5}의 할인가 합을 계산하기 위해서는 스텝 1에서 구한 할인가에서 “out” 노드인 1의 할인가를 빼주는 한 번의 연산만 필요하다. 즉, 각각의 스텝에서의 할인가는 이전 스텝의 할인가에서 “in” 노드의 할인가를 더해주고 “out” 노드의 할인가를 빼주면 구할 수 있다.

Table 1.

Scanning Example

• 정리 2. 스케닝 알고리즘은 열생성 부문제를 O(n³) 다항시간에 푼다.

증명. 두 개의 노드 쌍마다 한 번의 스케닝을 수행하기 때문에 스케닝은 총 n(n - 1)/2회 실행된다. 하나의 스케닝에서는 최대 n - 2개의 스텝이 있을 수 있고 하나의 스텝에는 최대 두 번의 연산이 필요하다. 이때, 중심점에서 같은 거리에 또 다른 고객 노드가 존재하는 경우에는 하나의 스텝에서 둘 보다 많은 연산이 필요할 수도 있다. 그러나 이 경우에는 다른 스텝에서 연산이 감소하게 되고 선택된 노드 쌍이 아닌 모든 다른 노드는 최대 한 번씩 들어오고(in) 나갈(out) 수 있기 때문에 한 번의 스케닝에서 발생하는 연산은 스텝의 수인 n - 2의 두 배가 된다. 결과적으로 $$ n\left(n-1\right)/2\cdot 2\left(n-2\right)=n\left(n-1\right)\left(n-2\right)=O\left(n^3\right)$$를 얻는다. □

스케닝 방법을 사용하기 위해서는 먼저 모든 두 개의 노드 쌍에 대해 스텝마다의 in 노드와 out 노드를 구해야 한다. 이 과정은 다음과 같이 수행될 수 있다.

• 과정 1. 선택된 두 개의 노드 쌍에 대해 나머지 노드들로 정의되는 최소 부분집합 설비의 중심점을 구하고 이 중심점으로 정의되는 최소 부분집합 설비의 유효 반경을 구한다. 이때, 중심점이 두 노드 쌍의 왼쪽에 있으면 유효 반경에 –1를 곱해준다.
• 과정 2. 유효 반경들을 오름 차순으로 정렬(sorting)한다. 정렬한 최소 부분집합 설비들의 순서가 스케닝 알고리즘의 스케닝 순서가 된다.
• 과정 3. 처음 스케닝 스텝의 최소 부분집합 설비가 커버하는 노드 집합을 구하고 다음 스케닝 스텝의 최소 부분집합 설비가 커버하는 노드 집합과 비교하여 in 노드와 out 노드를 구한다. 이 과정을 마지막 스케닝 스텝까지 반복한다.
• 과정 4. 과정 1로 돌아가서 다른 노드 쌍에 대해 같은 과정을 반복한다.

위 과정을 통해 모든 스케닝 스텝의 in 노드와 out 노드를 한번 구하면 이후에 열생성 과정에서는 스케닝 과정에서 최소 부분집합 설비의 할인가를 쉽게 계산할 수 있다.

4. 계산 실험

본 연구에서 제시한 알고리즘의 성능을 평가하기 위해 전산실험을 수행하였다. 모든 실험은 Intel i9-13900K CPU와 64GB RAM을 가진 Linux 머신에서 수행했으며 모든 알고리즘은 Python 3.10으로 구현되었다. 열생성 스케닝 알고리즘은 Python상에서 numba 패키지를 이용해서 구현하였다(Numba 2018). MIP Solver로는 Cplex 22.1을 사용하였다. 모든 실험에서 최대 8개의 스레드를 사용하였다.

4.1 문제 생성 방법

전산 실험을 위해 Berman et al.(2009)과 같은 방법으로 문제를 생성하였다. 먼저 1000 * 1000 평면에 주어진 값인 n개의 노드를 임의의 위치로 생성한다. 같은 크기의 문제(n이 같은 문제들) 10개의 다른 문제를 렌덤 시드(random seed)를 다르게 설정하여 생성한다. 이들 문제를 n숫자-s숫자로 표기한다. 예를 들어 문제 이름 n100-s1은 고객 노드 100개를 렌덤 시드 1을 이용해서 생성한 문제이다. 설비의 설치 비용 함수는 다음과 같이 정의한다.

$\[ F+r^{\alpha } \]$

이 때, F와 α는 주어지는 파라메터이고 r은 설비의 유효 커버 반경이다. 예를 들어 F = 100이고 α = 2라면 설비의 고정비는 1000이고 유효 반경의 제곱으로 변동비가 증가한다.

문제의 크기는 노드 개수 최소 50개에서 최대 350개까지 25개 차이로 생성하였다. 즉, 총 130개 문제를 생성하였다.

(1) 알고리즘 성능 비교

제안한 분지평가법 알고리즘의 성능을 기존 알고리즘과 비교하였다. 아래는 본 절에서 비교한 알고리즘이다.

• - BnP: 본 연구에서 제시한 분지한계법 알고리즘
• - BF: 부르트-포스(brute-force) 알고리즘. 문제(MP)를 모든 최소 부분집합 설비 C에 대해 정의하고 결정 변수를 이진 변수로 변경

BF는 주어진 문제에 대해 최소 부분집합 설비를 모두 미리 생성하고 이것들을 모두 한꺼번에 수리모형에 포함해 문제를 푸는 방법이다(Berman et al., 2009). 따라서 모든 부분집합 설비를 찾아내는 과정이 필요하고 이때 많은 시간이 소요될 수 있다. 본 비교 실험에서는 BF를 위해 모든 부분집합 설비를 생성하는 시간도 문제 푸는 시간에 포함했다. BF는 이미 모든 변수가 포함되어 있기 때문에 열생성 과정이 필요 없는 일반적인 이진 변수 최적화 문제(Binary Programming Problem)이다. BF는 이 모형을 Cplex로 푼 것을 의미한다. 이와는 별도로 수리모형 (P)를 비선형 정수계획법 solver인 couenne(Couenne, 2006)로 성능을 평가하였다. 그러나 couenne는 노드 개수가 10개인 경우에도 몇 시간 이내에 최적해를 구하지 못하였다. 따라서 couenne는 이후의 전산 실험 비교 대상에서 제외하였다. 모든 알고리즘은 제한시간 1800초(30분)를 사용하였으며 cplex는 스레드 수 8개를 제외하고는 모두 디폴트 설정값을 사용하였다. 분지한계법에서 열생성 과정은 노드 쌍에 대해서 독립적이고 병렬적으로 실행할 수 있다. 따라서 numba를 사용해서 8개의 스레드를 이용해서 열생성 과정의 성능을 개선하였다. 그리고 이 과정에서 각 스레드는 최대 20개까지의 음수 할인가를 가지는 열(column)을 찾아 (MP-R)에 추가하였다. 이 절의 모든 실험은 설비 비용 파라메터 F = 2000과 α = 2를 사용하였다.

<Table 2>는 BnP와 BF 알고리즘의 성능을 비교한다. 주어진 크기 10개 문제들에 대해 문제를 푸는 시간(초)의 평균(mean)을 비교한다. 마지막 열의 BnP/BF (%)는 BF 대비 BnP 알고리즘의 풀이 시간의 비율(%)을 보여준다. 문제의 크기가 작은 경우(예: 50)에는 두 알고리즘 모두 최적해를 빠른 시간에 구할 수 있었으나 문제의 크기가 커질수록 BF 풀이 시간이 빠르게 증가하였다. 모든 경우에서 BnP는 BF 대비 매우 빠른 시간에 최적해를 구하였는데, BnP 알고리즘은 BF 대비 4~8%의 시간을 사용하여 최적해를 구하였다. 특히, BF 알고리즘은 문제 크기 250개 이상은 알고리즘 실행 시 메모리 부족 에러가 발생하여 최적해를 구하지 못하였다. 이러한 결과의 이유로는 노드 개수가 250개라면 BF 알고리즘에서 필요한 이진 결정변수의 수가 최대 2,604,375에 이르기 때문이다. <Figure 4>는 알고리즘의 문제 풀이 시간을 비교하는 그래프이다. BnP 대비 BF는 매우 빠른 속도로 문제 풀이 시간이 증가하는 것을 명확하게 보여준다. 결과적으로 BnP는 BF 대비 매우 뛰어난 문제 풀이 성능 보임을 확인하였다.

Table 2.

Mean Solving Times of BnP and BF.

Figure 4.

Comparison of mean solving times between BnP and BF

<Table 3>은 알고리즘에서 사용한 열(column)의 개수를 보여준다. BF 알고리즘에서도 최소 부분집합 설비 C의 크기는 이론상 최댓값보다 작을 수 있다. 이유는 일직선 상에 있는 세 개의 노드는 최소 부분집합 설비를 정의하지 못하기 때문이며 <Table 3>의 # columns은 실제로 문제를 푸는데 사용된 열의 개수를 나타낸다. BnP의 경우에는 열생성 과정에서 추가된 모든 열의 개수를 표시한다. BnP/BF (%)는 BF 대비 BnP에서 사용된 열 개수의 비율(%)을 보여준다. 표에서 확인할 수 있듯이 BnP는 BF 대비 매우 작은 수의 열만을 사용하였다(최대 ~1%). 특히, 이러한 비율은 문제의 크기가 커질수록 작아지며 이는 실제로 문제를 필요한 최소 부분집합 설비의 개수는 가능한 설비의 집합의 개수보다 훨씬 작을 수 있음을 보여준다.

Table 3.

Mean Solving Times of Algorithms BnP and BF

(2) 설비 비용 파라메터에 따른 최적해 비교

본 연구의 대상인 VRCP의 특징은 설비의 비용함수를 다양하게 사용할 수 있다는 점이다. 본 절에서는 설비 비용 파라메터인 F와 α에 따른 최적해의 변화에 대해 설명한다. 직관적으로 설비의 고정비 F는 설비의 개수에 그리고 변동비 파라메터 α는 설비의 유효 반경의 크기에 영향을 줄 것이다. 구체적으로 F가 커지면 설비의 개수는 감소하고 α가 커지면 설비의 유효반경이 작은 최적해를 얻을 것으로 기대된다. <Figure 5>는 본 절의 실험 대상인 문제 n200-s9를 보여준다. 그림에서 회색 점은 고객 노드를 나타내며 이 문제에는 200개의 고객 노드가 존재한다.

Figure 5.

Problem n200-s9

<Figure 6>은 같은 문제(n200-s9)에 대해 같은 α = 2값을 가지면서 다른 F를 가지는 문제들의 최적해를 비교한다. 그림에서 붉은색 점은 설치된 설비의 중심(또는 설치 위치)을 의미한다. 설비의 유효 반경은 회색 원(circle)으로 표시되어 있다. 각 최적해의 목적함수 값은 ( )안에 표시되어 있다. 예상한 대로 F값이 커질수록 설비의 개수는 줄어들고 반면에 설치된 설비들의 유효 반경은 커지는 것을 알 수 있다.

Figure 6.

Optimal solutions for n200-s9 problem with α = 2 and different F values

<Figure 7>은 같은 F값(2000)을 가지며 다른 α값(1.9, 2.0, 2.1)을 가지는 경우의 최적해 변화를 보여준다. 유효 반경의 비용함수의 증가율을 나타내는 α가 커질수록 최적해의 설비들의 유효 반경이 줄어듬을 확인할 수 있다.

Figure 7.

Optimal solutions for n200-s9 problem with F = 2000 and different α values

5. 결 론

본 논문은 가변 범위 커버링 문제(variable radius covering problem, VRCP)에 대한 최적해 알고리즘을 제시한다. VRCP는 기존의 고정된 크기를 가지는 설비 위치 문제(FLP)를 일반화한 것으로 설비의 유효 반경에 따른 비용 변화를 고려할 수 있다는 점에서 보다 현실적인 문제라 할 수 있다. VRCP는 Berman et al.(2009)에 의해 처음 정의되고 휴리스틱 알고리즘이 제안되었다. 본 연구는 VRCP에 대해 분지평가법에 기반한 최적해 알고리즘을 최초로 제시하였으며 이는 기존의 수리모형을 이용한 방법(브루트-포스 모형)에 비해 훨씬 뛰어난 성능을 보임을 다양한 전산실험을 통해 확인하였다. 또한 문제의 파라메터에 따른 최적해의 변화를 확인하였고 이는 이 문제를 현실 상황에 적용하는데 계획 수립자에게 유용한 판단 정보를 제공할 수 있다. 전산 실험에서 사용한 변동비 파라메터 α는 실험의 편의를 위해 임의로 도입한 것으로 실제 변동비 함수는 훨씬 더 복잡한 구조를 가질 수 있다. 그럼에도 불구하고 이 함수가 비-감소 가정을 만족한다면 본 연구가 제안한 알고리즘을 적용할 수 있으며 이는 본 연구의 장점 중의 하나이다.

본 연구는 다양한 방법으로 확장될 수 있다. 첫째, 설비 위치 문제에서 발생하는 다양한 추가 제약을 반영할 수 있다. 예를 들어, 위성을 통한 통신 구역 설정 문제에서는 설비의 위치가 위성의 통신 커버리지 영역 문제가 된다. 이 경우에는 통신 영역이 커버링 하는 지역의 중복(overlap)은 위성 통신의 장애를 유발할 수 있으며, 따라서 중복 영역을 금지하거나 최소화하는 추가적인 제약이 필요하다. 커버링 지역의 중복을 금지하기 위해서는 특정 조건을 만족하는 설비들이 동시에 선택되지 않아야 하며 이런 제약은 열생성 과정에서 기존 알고리즘을 적용하기 어렵게 만든다. 둘째, 위치 선정 문제를 경로 문제와 결합하여 새로운 위치-경로 문제를 정의할 수 있다. 예를 들어 전쟁 상황에서 레이더를 무력화 하는 제밍(jamming)은 제밍의 범위에 따라 소모하는 에너지와 장비가 달라진다. 즉, 제밍의 범위와 장비를 탑재한 차량의 경로를 동시에 결정해야 하며 이는 새로운 변동 가능한 범위를 가지는 위치-경로 문제로 정의된다. 이 새로운 문제는 기존의 최적해 알고리즘을 적용하기 어렵고 새로운 접근 방법이 필요하면 본 연구가 제시한 알고리즘은 그 기반을 제공하는 것이 가능하다.

Acknowledgments

본 연구는 2025학년도 한국외국어대학교 교내학술 연구비의 지원에 의하여 이루어졌음.

References

• Balinski, M. L. (1965), Integer Programming: Methods, Uses, Computations, Management Science, 12(3), 253-313. [https://doi.org/10.1287/mnsc.12.3.253]
• Barnhart, C., Johnson, E. L., Nemhauser, G. L., Savelsbergh, M. W., and Vance, P. H. (1998), Branch-and-price: Column generation for solving huge integer programs, Operations Research, 46(3), 316-329. [https://doi.org/10.1287/opre.46.3.316]
• Berman, O., Drezner, Z., Krass, D., and Wesolowsky, G. O. (2009), The variable radius covering problem, European Journal of Operational Research, 196(2), 516-525. [https://doi.org/10.1016/j.ejor.2008.03.046]
• Berman, O., Drezner, Z., and Krass, D. (2010), Generalized coverage: New developments in covering location models, Computers & Operations Research, 37(10), 1675-1687. [https://doi.org/10.1016/j.cor.2009.11.003]
• Berman, O., Drezner, Z., and Krass, D. (2019), The multiple gradual cover location problem, Journal of the Operational Research Society, 70(6), 931-940. [https://doi.org/10.1080/01605682.2018.1471376]
• Beste, D. (1978), Design of Satellite Constellations for Optimal Continuous Coverage, IEEE Transactions on Aerospace and Electronic Systems, AES-14(3), 466-473. [https://doi.org/10.1109/TAES.1978.308608]
• Blanco, V., Gázquez, R., and Saldanha-da-Gama, F. (2023), Multi-type maximal covering location problems: Hybridizing discrete and continuous problems, European Journal of Operational Research, 307(3), 1040-1054. [https://doi.org/10.1016/j.ejor.2022.10.037]
• Cooper, L. (1963), Location-allocation problems, Operations Research, 11, 331-343 [https://doi.org/10.1287/opre.11.3.331]
• Cornuejols, G., Nemhauser, G. L., and Wolsey, L. A. (1990), The uncapacitated Facility Location Problem, Discrete Location Theory, John Wiley and Sons, New York, 119-171.
• Couenne (2006), https://github.com/coin-or/Couenne, .
• Jabalameli, M. S., Tabrizi, B. B., and Javadi, M. M. (2011), A Simulated Annealing method to solve a generalized maximal covering location problem, International Journal of Industrial Engineering & Computations, 2, 439-448. [https://doi.org/10.5267/j.ijiec.2011.01.003]
• Karatas, M. and Eriskin, L. (2023), Linear and piecewise linear formulations for a hierarchical facility location and sizing problem, Omega, 118, 102850. [https://doi.org/10.1016/j.omega.2023.102850]
• Klose, A. and Görtz, S. (2007), A branch-and-price algorithm for the capacitated facility location problem, European Journal of Operational Research, 179(3), 1109-1125. [https://doi.org/10.1016/j.ejor.2005.03.078]
• Küçükaydın, H. and Aras, N. (2020), Gradual covering location problem with multi-type facilities considering customer preferences, Computers & Industrial Engineering, 147, 106577. [https://doi.org/10.1016/j.cie.2020.106577]
• Laporte, G, Nickel, S., and Saldanha, F. (2015), Location science, Springer Cham, Switzerland. [https://doi.org/10.1007/978-3-319-13111-5]
• Ni, W., Shu, J., Song, M., Xu, D., and Zhang, K. (2021), A branch-and-price algorithm for facility location with general facility cost functions, INFORMS Journal on Computing, 33(1), 86-104. [https://doi.org/10.1287/ijoc.2019.0921]
• Numba (2018), https://numba.pydata.org/