큰 규모의 전력망(electric power grid)을 효율적으로 유지 보수하기 위해서는 매우 복잡한 의사결정과정이 필요하다. 일반적인 전력망은 발전소(power stations), 변전소(electric substations), 전력선(power transmission lines) 등으로 구성되어 있다(Kaplan, 2014). 각각의 구성 요소들은 고유한 특징을 가지고 있어 서로 다른 의사결정 문제의 대상이 된다. 예를 들어, 발전소에서 전력을 생산하는 문제는 흔히 발전기 충당 문제(Unit Commitment Problem, Saravanan et al., 2013, 이하 UCP)라고 한다. UCP는 전력 수요를 충족하기 위해 발전기의 가동 스케줄을 정하는 문제로써 매우 많은 연구가 이루어져 왔다(Bertsimas et al., 2012; Carrión and Arroyo, 2004; Kazarlis et al., 1996). UCP의 가장 큰 특징은 발전기의 작동 특성 제약을 반영해야 한다는 것이다. 이를 위해 대부분의 경우 풀기 어려운 혼합 정수계획법(Mixed Integer Programming, MIP) 문제로 모형화된다. 전력망에 대한 또 다른 주 연구주제로는 전력망 디자인(power gird design) 문제가 있다(Keyhani, 2016). 이 문제는 실제로 전력을 공급하는 전송망을 계획하는 문제로써 네트워크 디자인(Network flows, Ahuja and Magnanti, 1998) 문제의 특수한 경우라 할 수 있다. 전력을 수송하기 위한 전력망은 교류 전기의 비선형적인 저항을 고려해야 하므로 일반적으로 매우 풀기 어려운 비선형 혼합정수계획법(Nonlinear MIP) 문제로 모형화하거나 시뮬레이션 기반의 최적화 알고리즘을 사용한다(Huang et al., 2018; Koutroulis and Blaabjerg, 2011).

최근 빠른 속도로 가속화되는 기후변화에 따라, 재사용 가능 에너지(renewable energy)를 사용하는 전력망에 관한 연구 또한 많이 이루어지고 있다(Liserre et al., 2010; Carrasco et al., 2006). Keck et al.(2019)은 배터리 전력 저장소(Battery Energy Storage 또는 Electrical Energy Storage, EES)를 이용한 전력망 효율 개선 가능성에 관해 연구하였다. 신재생에너지는 기후와 같은 외부환경에 민감하다는 단점이 있는데, 전력 저장소를 사용하여 효율을 크게 개선할 수 있음을 보였다. Schmid(2012)는 인도 사례를 중심으로 재생에너지의 정착을 위한 최적 정책 연구를 수행하였다. 전기자동차 등의 새로운 전력 수요로 인해 전력망 운영에 미치는 영향에 관한 연구 또한 이루어지고 있다(Xiao, 2014).

본 논문이 다루는 문제는 전력망 관리(maintenance)를 위한 작업할당(job assignment) 계획 문제이다. 전력망을 안정적으로 유지하기 위해서는 주기적인 선행 관리(proactive maintenance)가 필요하다(Wu et al., 2011). 계획 수립 기간(planning horizon) 동안 수행할 작업의 집합을 생각해보자. 각각의 작업을 수행하기 위해서는 여러 가지 자원(resource)이 소모된다. 그리고 특정 기간에 사용 가능한 자원의 상한(upper bound)과 하한(lower bound)이 주어진다고 가정하자. 전력망의 목적은 전기(electric)를 안정적으로 공급하는 것이기 때문에 유지 보수과정에서 발생하는 위험(risk)에 대한 고려가 필수적이다. 따라서 목적함수는 작업을 모두 수행했을 때 발생하는 위험을 최소화하는 것이다. 실제 전력망을 관리하는 전력회사는 다양하고 고유한 위험 척도(risk measure)를 사용한다. 본 연구는 실제 프랑스의 전력망 운영 회사인 RTE(https://www.rte-france.com/)가 사용하는 위험 척도와 실제 작업할당 데이터를 사용한다. RTE는 오랜 경험을 통해 3단계 전력망 관리 체계를 수립하였다(Ruiz et al., 2020). 1단계는 작업의 위험을 많은 수의 시나리오로 표현하고 각각의 시나리오에 위험도를 부여하는 단계이다. 2단계는 위험 시나리오들을 바탕으로 작업 수행 시기를 결정한다. 3단계는 수립된 작업 수행 계획을 검증한다. 본 연구의 대상은 2단계인 작업할당 최적화 문제이며 혼합정수계획법(MIP) 문제로 수리 모형화가 가능하다. 그러나 이 MIP 수리모형을 Cplex와 같은 상용 MIP Solver로 풀기에는 매우 어렵다. 그 이유로는 첫째, 실제 RTE가 사용하는 수준의 작업의 수와 계획 수립 기간을 고려하기 위해서는 수리모형의 크기가 매우 커지게 된다. 둘째, 위험 수준을 반영하기 위해 많은 수의 시나리오를 고려해야 하며 이 때문에 수리모형이 복잡해지고 많은 수의 결정변수와 제약식이 추가된다. 셋째, 본 문제는 일반할당문제(Generalized Assignment Problem, GAP, Ross and Soland, 1977)의 일반화된 문제로 풀 수 있어 NP-hard에 속하는 문제이다.

전력공급에 대한 지속가능성은 크게 두 가지 구성요소에 좌우된다. 첫 번째는 전력발전이고 두 번째는 전력공급이다. 전력의 발전 측면에서의 지속가능성은 많은 연구가 이루어져 왔으나 전력의 안정적인 공급을 위한 연구는 상대적으로 많이 이루어지지 않고 있다. 특히, 기후변화와 탄소세(carbon tax)의 도입은 기존의 화석연료를 전기로 대체할 것을 요구하고 있다. 대표적인 예로 화석연료를 전기로 대체하는 전기 자동차가 있으며 이러한 경향은 앞으로 점점 더 가속화될 것이다. 이러한 변화는 전력망의 관리와 운영이 예전보다 더 중요한 요소로 작용하리라는 것을 의미한다. 본 논문은 전력망 점검/관리에 가장 중요한 요소인 전력망 정전(blackout)을 방지하기 위한 위험관리(risk management)를 고려한 작업할당 문제를 다룬다는 의의가 있다.

본 논문은 다음과 같은 주요 기여점을 가진다.

• - 실제 전력망 관리 회사(RTE)의 작업할당문제를 정의하고 실제 데이터를 이용하여 수리모형을 수립한다.
• - 수리계획법과 발견적 기법을 결합한 혼합 발견적 (Hybrid Heuristic) 알고리즘을 제시한다.
• - 전산실험을 통해 제시한 알고리즘이 MIP 수리모형에 비해 매우 뛰어난 성능을 보임을 확인한다.

본 논문의 나머지 구성은 다음과 같다. 제2장은 작업할당문제에 관련한 기존 연구결과를 소개한다. 제3장은 문제를 구체적으로 정의하고 MIP 수리모형을 제시한다. 수리계획-발견적 복합 알고리즘이 제4장에서 제시되고 제5장은 전산실험 결과를 보고한다, 마지막으로 제6장은 결론을 제시한다.

작업할당문제는 고전적인 조합최적화문제 중 하나인 할당문제(assignment problem, AP, Kuhn, 1955)의 일반화 문제로 볼 수 있다. AP에 대한 다항시간 알고리즘이 알려져 있지만(Kuhn, 1955), AP에서 작업기간마다 자원 소모량의 제약식이 추가된 작업할당문제는 GAP 또는 배낭문제(Knapsack Problem, KP)의 일반화문제로써 NP-hard에 속한다(Garey and Johnson, 1979). 할당문제와 관련된 다른 최적화 문제로 가중치할당문제(Weight Assignment Problem, WAP, Ross and Zoltners, 1979)가 있다. GAP과 WAP 모두 작업을 할당할 때 소모되는 여러 가지 자원들이 존재하고 소모되는 자원의 총합의 상한과 하한이 주어진다는 공통점이 있다. 따라서 이 자원 제약식들은 복합배낭문제(Multiple Knapsack Problem, MKP, Martello and Toth, 1990)로 볼 수 있다. MKP 역시 NP-hard에 속하는 문제로써 많은 연구가 이루어졌다. 특히, 배낭 제약식(Knapsack constraints)은 정수계획법 수리모형에서 매우 자주 등장하기 때문에 선형완화 한계값(linear relaxation bound)을 개선하기 위한 유효제약식(valid inequalities)에 대한 연구가 많이 이루어졌다(Nemhauser and Wolsey, 2014). 작업할당문제 또한 배낭제약식을 많이 가지고 있어 배낭문제의 유효제약식을 모두 적용할 수 있다. 이는 특히 Cplex와 같은 MIP solver로 푸는 수리모형의 성능을 크게 향상시킨다.

일반적으로 GAP에서는 작업과 작업수행자(agent)가 존재하며 각각의 작업은 정확히 한 명의 작업수행자에 할당되어야 한다. 작업수행자는 각각의 작업을 수행하면서 자원을 소모하고 소모한 자원의 총량은 주어진 한계를 넘을 수 없다. 본 논문의 작업할당문제는 GAP에서의 작업수행자가 작업 수행 시작 시점으로 대체된 형태로 볼 수 있다. 즉, 각각의 작업은 정확히 하나의 작업 수행 시작 시점에 할당되어야 하고, 어떤 작업 수행시간의 자원 소모량은 해당 시점 자원의 총량을 넘지 않아야 한다. 만약에 모든 작업이 하나의 작업 시간 단위(time period)에 종료된다고 가정하고 하나의 위험 시나리오만 있다고 가정하면 작업할당문제의 제약식은 GAP과 동일하게 표현된다.

GAP를 변형한 많은 문제들이 존재하며 다양한 알고리즘이 연구되었다. 처음 GAP를 제시한 Ross and Soland(1975)는 분지한계법(branch and bound) 알고리즘을 이용한 최적해 알고리즘(exact algorithm)을 제시하였다. Savelsbergh(1997)은 각각의 자원 제약식에 대한 열 생성(column generation)을 이용한 분지평가법(branch and price) 알고리즘을 개발하였다. GAP에 대한 최적해 알고리즘의 경우에는 선형완화 한계값을 개선하는 데 주력하였으며 이 과정에서 열 생성의 수렴속도를 개선하는 연구들도 많이 발표되었다(Wentges, 1997; Lee and Park, 2011). GAP에 대한 소개와 특히 불확실성이 존재하는 경우들에 대한 변형문제들에 대해서는 Öncan(2007)을 참조한다.

GAP의 최적해를 보장하지 못하지만, 해의 한계를 보장하는 근사(approximate) 알고리즘에 관해서도 많은 연구가 이루어졌다. Shmoys and Tardos(1993)은 주어진 파라메터 C에 대해 목적함수가 C보다 작은 해가 존재하지 않거나 목적함수가 최대 C이며 자원을 주어진 것보다 최대 두 배 사용하는 해를 찾는 다항시간(polynomial time) 알고리즘을 제시하였다. Cohen et al.(2006)은 배낭문제에 대한 근사알고리즘을 이용하여 GAP의 근사알고리즘을 개발하였다. 이 알고리즘은 배낭문제에 대한 α-근사알고리즘을 사용하여 GAP에 대한 (1+α) 근사비율(approximation ratio)을 가진다. Nutov et al.(2006)은 GAP에 약간의 가정이 추가된 Max-GAP에 대해 (1-1/e)-근사 알고리즘을 제시하여 기존의 1/2-근사 알고리즘의 성능을 크게 개선하였다.

앞서 소개된 알고리즘들은 해의 품질을 보장하기 위해 대부분 문제 크기에 대해 지수함수적인 시간 복잡도(exponential time complexity)를 가진다. 빠른 시간 내에 좋은 해를 얻기 위한 발견적 기법 또한 많이 연구되었다. Díaz and Fernández(2001)은 Tabu-탐색을 이용한 메타-발견적 알고리즘을 제시하였다. Özbakir et al.(2010)은 벌-알고리즘(Bee algorithm, BA)을 GAP에 적용하였다. BA는 다양한 인구-기반(population-based) 알고리즘 중 하나로 꿀벌이 먹이를 찾아가는 방법에 착안하여 해를 개선해 나간다. 이 밖에도 개미-군집 방법(Ant Colony Method, ACM, Lourenço and Serra, 2002; Randall, 2004), 경로-재연결(Path relinking) 알고리즘(Alfandari et al., 2003), 유전자(genetic) 알고리즘(Feltl and Raidl, 2004) 등이 GAP의 해법으로 연구되었다.

일반적인 GAP 또는 WAP에서는 문제를 정의하는 데이터가 모두 정확하게 결정되어 주어진다는 가정을 바탕으로 한다. 본 논문에서 제시하는 작업할당문제는 발생 가능한 다양한 상황에 대해 시나리오를 정의하고 각 시나리오에 발생하는 위험을 최소화하는 것을 목적으로 한다. 이렇게 불확실한 상황에 대해서 다양한 시나리오를 고려하여 위험을 회피(risk averse)하거나 미래에 대한 기댓값(expectation value)을 최적화하는 접근방법을 일반적으로 추계적 최적화(stochastic optimization)라고 한다(Birge and Louveaux, 2011). Albareda-Sambola et al.(2006)은 할당된 작업들 중 일부가 할당되지 않을 가능성이 있는 경우에 대해 다루었다. Sarin et al.(2014)는 작업 수행에 필요한 자원의 양이 불확실한 경우에 대해 분지평가법 알고리즘을 제시하였는데, 이 경우에는 자원 소모량 제약식이 추계적 배낭문제(stochastic Knapsack problem)가 되어 풀기가 어려워지며 이를 해결하기 위해 열 생성 방법을 사용하였다.

본 논문에서 다루는 작업할당문제는 위에 소개된 GAP을 대상으로 개발된 알고리즘들을 그대로 사용할 수 없다. 그 이유로는 첫째, 작업이 어떤 작업 수행 시작 시점에 할당되면 주어진 작업 지속기간 동안 계속되어야 하기 때문이다. 이는 각각의 시점에 대해 독립적으로 작업을 할당할 수 있는 것이 아니라, 시간상으로 선후 관계에 있는 작업 시간대에 의존함을 의미하기 때문에 모든 작업수행자가 독립적으로 작업을 할당되어야 하는 GAP의 가정에 위반된다. 둘째, (제3장에서 자세히 설명할) 목적함수의 위험 평가함수가 일반적으로 사용되는 기댓값 함수가 아니라 기댓값 함수와 주어진 τ로 정의되는 τ-분위(quantile) 값의 합으로 정의되기 때문이다. 이러한 목적함수는 데이터와 문제 상황을 제시한 RTE가 실제로 사용하는 형태로 기존의 추계적 GAP 또는 WAP 문제와 결정적으로 다른 형태의 수리모형과 풀이방법을 요구한다.

본 장에서는 추계적 작업할당문제(Stochastic Job Assignment Problem, 이하 SJAP)를 정의하고 수리모형을 제시한다.

계획 수립 기간 T동안 수행해야 할 작업들의 집합 I와, 각각의 작업이 수행될 때 소모되는 자원의 집합 R을 고려해보자. 각각의 작업 i∈I는 계획 수립 기간 T내에서 시작되어야 하며, 작업 수행 시작 시점 t'∈T에 따라 정해지는 작업 지속기간 d_i,t'동안 중단될 수 없다. 또한, 모든 작업은 |T| + 1 시점 이전에 종료되어야 하며 이는 작업 수행 시작 시점 t'이 |T| - d_i,t'+1보다 앞서거나 같아야 한다는 것을 의미한다. (이하 제약 유형 1).

<Figure 1>은 제약 유형 1에 대한 예시를 보여준다. 계획 수립 기간이 T = {1,2,3}, 작업 집합이 I = {i₁, i₂}로 주어졌을 때, 모든 작업의 작업 수행 시점은 4를 넘을 수 없다. 작업 i₁은 계획 수립 기간 내에 시작되고 종료되는 반면, 작업 i₂는 작업 수행 시점의 한계 시점을 넘어 배정되어 제약 유형 1을 어기게 된다.

Figure 1. 
Example of Type 1 Constraints

t'∈T에 시작되는 작업 i∈I를 t∈t', t'+1,⋯, t'+di,t'-1에 수행하기 위해서는 자원 r∈R이 bi,t'r,t만큼 소모되며, 어떤 시점 t∈T에 사용 가능한 자원 r∈R의 소모량은 상한 utr와 하한 ltr을 만족해야 한다(이하 제약 유형 2).

<Figure 2>는 제약 유형 2에 대한 예시를 보여준다. 계획 수립 기간이 T = {1,2,3}, 작업 집합이 I = {i₁, i₂}, 자원 집합 R={r}로 주어졌을 때, 모든 시점 t∈T에서의 자원 r의 총 소모량에 대한 상한과 하한이 각각 5, 0이라고 가정하자. 이때, 그림에서 각 bi,t'r,t는 해당 시점에서 사용되는 자원 r의 소모량을 의미한다. 시점 1, 3에서 자원 r의 소모량은 각각 2, 4로 해당 자원에 대한 상한과 하한을 만족하지만, 시점 2에서의 소모량은 6으로 상한을 넘게 되어 제약 유형 2를 어기게 된다.

Figure 2. 
Example of Type 2 Constraints

전력망 관리 관점에서 어떤 시점 t∈T에서 두 작업 i₁, i₂∈I의 연관성이 깊어 함께 작업할 때 전체 전력망 유지가 불가능한 경우가 있다. 이 경우에 대한 삼중 쌍(triplets) (t, i₁, i₂)의 집합 E가 주어지고 이들 삼중 쌍의 두 작업 i₁, i₂∈I는 시점 t∈T에 동시에 수행될 수 없다고 가정한다(이하 제약 유형 3).

<Figure 3>은 제약 유형 3에 대한 예시를 보여준다. 계획 수립 기간이 T = {1,2,3}, 작업 집합이 I = {i₁, i₂, i₃}, 집합 E=t,ip,iq∀t∈1,2,∀ip∈I,∀iq∈I, ip≠iq가 주어졌다고 하자. 즉, 모든 작업은 시점 1, 2에서 함께 수행될 수 없다. 시점 1의 경우 작업 i₁만 작업이 수행되며 집합 E의 조건을 만족하지만, 시점 2의 경우 작업 i₁과 작업 i₃가 동시에 수행되며 제약 유형 3을 어기게 된다.

Figure 3. 
Example of Type 3 Constraints

어떤 작업할당계획을 모든 작업 i∈I에 대해 작업 i의 작업 수행 시작 시점 ti'들의 집합 π라하고 다음과 같이 정의하자.

이때, π는 t'∈T에 시작한 작업 i∈I가 시점 t∈t', t'+1,⋯, t'+di,t'-1에서 수행된다는 의미를 갖는다. 표현 가능한 π 중 SJAP의 세 가지 제약 조건을 모두 만족하는 작업할당계획을 유효해(feasible solution)라고 부르며 π~라고 하자.

전력망 유지 보수를 위한 작업은 안정적인 전력공급 관점에서 위험을 발생시키며, 본 논문에서 다루는 SJAP의 목적함수는 작업을 모두 수행했을 때 발생하는 위험 평가함수를 최소화하는 것이다. 본 연구에서 사용하는 위험 평가함수 f(π)는 RTE에서 정의한 것으로, 주어진 해 π에 대해 일반적으로 계산되는 기댓값 함수인 기대 위험(expected value of risk) 함수 g(π)와 주어진 τ∈[0, 1]로 정의되는 τ-분위(quantile) 값의 합인 기대 초과(expected value of excess risk) 함수 h(π)로 정의된다. 위험 척도는 추계적이기 때문에 변동성을 갖는데, 전체 위험을 기댓값 측면에서만 측정할 경우 산포 등에 대한 정보를 고려하지 못할 수 있으므로 기대 초과 값을 함께 고려한다.

시점 t∈T에서 정의된 위험 시나리오의 집합 S_t를 생각해보자. 시점 t'∈T에 시작한 작업 i∈I는 작업 수행 시점 t∈t', t'+1,⋯, t'+di,t'-1마다 정의된 위험 시나리오 s∈S_t에 대해 ci,t's,t에 해당하는 위험을 발생시킨다.

정의하는 추계적 작업할당문제에서 사용되는 집합 및 파라메터를 정리하면 다음과 같다.

• 집합 및 파라메터(Sets and parameters)
• T : 계획 수립 기간의 집합
• I : 계획 수립 기간 동안 수행해야 할 작업의 집합
• R : 작업이 수행될 때 소모되는 자원의 집합
• di,t' : 작업 i∈I 의 작업 수행 시작 시점 t'∈T에 따라 정해지는 작업 지속기간
• bi,t'r,t : t'∈T에 시작되는 작업 i∈I를 t∈t', t'+1, ⋯,t'+di,t'-1에 수행하기 위해서 사용되는 자원 r∈R의 소모량
• utr : t∈T에 사용되는 자원 r∈R의 총 소모량의 상한
• ltr : t∈T에 사용되는 자원 r∈R의 총 소모량의 하한
• E : t∈T에 함께 작업할 수 없는 두 작업 i₁, i₂∈I을 나타내는 삼중 쌍 (t,i₁,i₂)의 집합
• S_t : t∈T에 정의되는 위험 시나리오의 집합
• ci,t's,t : t'∈T에 시작되는 작업 i∈I가 t∈t', t'+1,⋯, t'+di,t'-1에 수행될 때 위험 시나리오 s∈S_t에 대해 발생하는 위험

어떤 해 π가 주어질 때 모든 t∈T에 대해, It:=i∈Iti'≤t≤ti'+di,t-1, ∀ti'∈π로 정의하자. 즉, 집합 I_t는 I의 부분집합이면서 시점 t∈T에서 수행되는 모든 작업 i∈I들의 집합이다.

위험 평가함수 f(π)의 첫 번째 항인 기대 위험 함수 g(π)는 위험 시나리오, 계획 수립 기간 두 측면에서의 평균 위험을 의미하며 다음과 같이 정의된다.

위험 평가함수 f(π)의 두 번째 항인 기대 초과 함수 h(π)는 계획 수립 기간의 시점별 위험 시나리오에 대한 평균 위험 대비 분위 수 측면에서의 평균적인 초과 위험을 의미하며, 이는 τ-분위 값에 의존한다.

공집합이 아닌 유한실수들의 집합 W⊆R이 주어졌다고 할 때, W에서의 τ-분위 값 Q_τ(W)는 W의 원소 중 ⌈τㆍ|W|⌉번째 큰 수를 의미하며 다음과 같이 정의된다.

그렇다면, 기대 초과 함수 h(π)는 다음과 같이 정의된다.

정의된 두 함수 g(π), h(π)에는 의사결정권자의 위험 수용 정책에 따라 가중치 α∈[0, 1]가 부여될 수 있다. 이렇게 정의되는 전체 위험평가함수 f(π)는 다음과 같다. 예제를 통한 자세한 계산과정이 부록 1 (Appendix 1)에 수록되어 있다.

fπ=αgπ+1-αhπ.

따라서 본 논문에서 다루는 SJAP는 정의된 제약을 모두 만족시키며 모든 작업을 수행했을 때 위험평가함수 f(π)를 최소화하는 문제이다.

본 장에서는 정의된 추계적 작업할당 문제를 해결하기 위한 수리계획-발견적 복합 알고리즘의 각 단계를 자세히 설명한다. 복합 발견적 기법은 수리계획법과 기존의 발견적 기법을 통합하여 점차적인 해의 개선을 이뤄낸다. 이를 위해서 수리모형을 간소화하고 수리모형의 해를 이용한 국소탐색(local search)을 메타-휴리스틱(Meta-heuristic)과 접목한다.

<Figure 4>는 전체 알고리즘의 흐름을 보여준다. 그림에서 확인할 수 있듯이 본 논문에서 제시하는 알고리즘은 크게 두 단계로 구성된다. 첫 번째 단계에서는 근사된(approximated) 수리모형을 통해 좋은 유효해를 얻는다. SJAP는 지켜야 할 자원의 종류가 많고 (제약 유형 2), 상호 배타적인 작업들을 고려해야 하기 때문에 (제약 유형 3) 간단한 방법으로 유효해를 얻기가 힘들다. 또한, RTE가 사용하는 위험평가함수는 3.2절에서 제시한 수리모형(P)을 매우 풀기 어렵게 만든다. 따라서 질 좋은 초기 유효해를 얻기 위해서 수리모형을 더 풀기 쉬운 형태로 근사한 후 초기 유효해를 얻는다. 초기 유효해는 원래 문제와 다른 목적함수를 가지는 수리모형에서 얻기 때문에 원 문제의 최적해임을 보장할 수 없다. 이를 개선하기 위해 국소 이웃해 탐색(local neighborhood search)을 수행한다. 이 과정을 수리-발견적 국소탐색 단계(Matheuristic Local Search Phase)로 부르며 4.1절에 자세한 설명을 제시한다. 두 번째 단계에서는 앞 단계에서 얻은 유효해를 이웃해(neighborhood solution) 탐색을 통해 지속해서 개선하는 과정을 메타-휴리스틱 알고리즘인 담금질 기법(Simulated Annealing, 이하 SA)과 함께 적용해서 개선한다. 이 과정을 메타-휴리스틱 전역탐색 단계(Meta-heuristic Global Search Phase)라고 부르며 이에 대한 자세한 내용은 제4.2절에서 설명한다.

Figure 4. 
Flowchart of Proposed Algorithm

4.1 수리-발견적 국소탐색

4.1.1 근사 수리모형을 이용한 유효해 탐색

제3.2절에서 제시한 수리모형(P)의 목적함수 (1)은 두 항으로 이루어져 있다. 첫 번째 항은 위험 시나리오, 계획 수립 기간 두 측면에서의 평균 위험을 의미한다. 두 번째 항은 계획 수립 기간의 시점별 위험 시나리오에 대한 평균 위험 대비 분위 수 측면에서의 평균적인 초과 위험을 의미한다. 목적함수의 두 번째 항을 구하기 위해서는 (9) - (11), (13)의 많은 제약식과 (20) - (21), (23)의 이진 변수가 추가로 필요하게 되며, 이는 문제를 풀기 어렵게 만든다. 그러나 이 제약식들은 해의 유효성(feasibility)에 영향을 미치지 않고 단지 목적함수 값을 계산하는 데 사용된다. 따라서 빠른 시간 안에 실용적인 유효해(feasible solution)를 얻을 수 있도록 수리모형 (P)의 목적함수의 두 번째 항을 좀 더 다루기 쉽게 근사한 새로운 수리모형을 제시한다.

수리모형 (P)에서 목적함수 (1)의 두 번째 항은 3장에서 정의된 주어진 해 π에 대한 기대 초과 함수 h(π)를 표현하는 항이며, 기대 초과 함수 h(π)는 주어진 해 π에 대해 다음과 같이 정의된다.

hπ=1T∑t∈Tmax0, Qτ∑i∈Itci,ti's,ts∈St-1St∑s∈St∑i∈Itci,ti's,t.

편의상 ∑i∈Itcs,ti,ti'을 Λst, ∑s∈StΛst를 Ω_t라고 하자. 어떤 시점 t∈T에 대해, Λst는 위험 시나리오 s∈S_t에 대해서 수행되는 작업 집합 I_t 기준 위험 합, Ω_t는 위험 시나리오 집합 S_t 기준 위험 합을 의미한다. 그렇다면, h(π)를 다음과 같이 다시 표현할 수 있다.

hπ=1T∑t∈Tmax0, QτΛtss∈St-1StΩt.

h(π)에서 어떤 시점 t∈T에 대해 QτΛsts∈St와 Ω_t를 구하기 위해서는 해당 시점 t에 수행되는 작업 집합 I_t에 대한 각각의 위험 값이 모든 위험 시나리오에 대해 연산되어야 한다. 이때, I_t는 모든 작업 i∈I의 작업 수행 시작 시점 ti'이 결정되어야 정의되는 집합이다. 즉, QτΛsts∈St와 Ω_t의 연산에서 각 시점 t∈T에서 수행되는 작업 집합 I_t와 작업 i∈I_t의 작업 수행 시작 시점 ti'는 서로 종속적이다. 이는 수리모형으로 표현했을 때 많은 이진 결정변수와 제약식을 필요로 하므로 문제를 풀기 어렵게 만드는 가장 큰 이유가 된다. 해당 연산에서 각 시점 t∈T에서 수행되는 작업 집합 I_t와 각 작업 i∈I_t의 작업 수행 시작 시점 ti'을 독립적으로 고려해서 연산할 수 있다면, 기대 초과 함수 h(π)는 각 시점에서 수행되는 작업 집합과 각 작업의 작업 수행 시작 시점에 대해서 위험 시나리오 집합에 대한 초과 위험 값을 단순히 더하면 되기 때문에 훨씬 다루기 쉬운 형태가 된다.

QτΛsts∈St와 Ω_t의 연산에 대해 각 시점 t∈T에서 수행되는 작업 집합 I_t와 각 작업 i∈I_t의 작업 수행 시작 시점을 독립적으로 고려하면서, h(π)와 비슷한 값을 가지는 함수 j(π)를 다음과 같이 정의한다.

jπ=1T∑t∈T∑i∈It∑t'∈Ttimax0, Qτci,t's,ts∈St-1St∑s∈Stci,t's,t.

이때, Tti:=t'∈Tt'≤t≤t'+di,t'-1로 정의되며, 집합 T^i_t는 시점 t∈T에서 수행되는 작업 i∈I_t에 대해, 작업 i를 시점 t에 작업하기 위해 해당 작업을 수행 시작해야 하는 시점 t'들의 집합이다.

위 식에서 알 수 있듯이 j(π)의 연산에서 각 시점 t∈T에서 수행되는 작업 집합 I_t와 작업 i∈I_t의 작업 수행 시작 시점 t'∈T^i_t은 독립적이다. 이러한 성질을 이용해서 문제를 풀기 전 모든 t∈T에 대해 정의되는 각각의 i∈I_t와 t'∈T^i_t에 대해 계산한 값들의 합을 미리 구한다면, 이를 통해 j(π) 다음과 같이 다시 표현할 수 있다.

jπ=1T∑t∈T∑i∈It∑t'∈Ttiξi,t't.

이때, ξti,t'는 다음과 같이 정의되며, t∈T에 따른 각각의 i∈I_t와 t'∈T^i_t의 위험 시나리오 집합 S_t에 대한 초과 위험을 의미한다.

ξi,t't=max0, Qτci,t's,ts∈St-1St∑s∈Stci,t's,t, if i∈It and t'∈Tti,0                 ,otherwise.

이렇게 정의된 j(π)는 t'∈T, t∈t', t'+1,⋯, t'+di,t'-1, I_t, S_t 에 따라 주어지는 cs,ti,ti'의 분포가 t'에 대해 비슷한 경향을 가질수록 h(π)와 가까운 값을 가지게 된다. 위험평가함수 f(π)의 두 번째 항을 h(π) 대신 j(π)로 고려하는 수리모형의 정의를 위해 결정변수 J를 다음과 같이 정의한다.

• J : 계획 수립 기간의 시점별 위험 시나리오에 대한 평균 위험 대비 분위 수 측면에서의 평균 초과 위험의 근사 값.

그렇다면, SJAP의 새로운 수리모형 (P2)를 다음과 같이 정의한다.

P2 min αG+1-αJ,

(26)

s.t. J≥1T∑t∈T∑i∈I∑t'∈Ttiξi,t'tyi,t't,

(27)

J≥0,

(28)

(2) - (8), (12), (14), (16) - (19), (22), (24).

그리고, 수리모형 (P2)를 풀어서 얻은 최적해에 대한 집합 θ를 다음과 같이 정의한다.

θ:=i,t'xi,t'=1, ∀i∈I, ∀t'∈T.

SJAP의 이웃해 탐색을 위한 초기 유효해 π~0을 다음과 같이 정의한다.

π~0=ti'ti'=t', ∀i,t'∈θ.

4.1.2 목적함수 섭동 알고리즘

4.1.1항에서 설명한 방법으로 얻은 초기 유효해 π~0를 기반으로 이웃해를 탐색하여 해의 개선을 시도한다. 수리모형 (P2)는 근사한 목적함수를 사용하기 때문에 어느 정도 좋은 해를 구할 수 있으나, 여전히 최적해와는 차이가 있을 수 있다. 이때 목적함수 계수에 변동을 주면 목적함수의 기울기가 바뀌어 다른 해가 도출될 수 있고, 이 해의 실제 목적함수 값은 개선될 수 있다. 이러한 과정은 목적함수 섭동(perturbation)이라고 알려져 있다.

어떤 해 π의 정확한 목적함수 값은 f(π)이며, 이는 해당 해를 얻은 후에 수리모형 (P)의 목적함수 (1)을 계산하면 쉽게 얻을 수 있다. 어떤 해 π의 어떤 작업 i∈I의 작업 수행 시작 시점을 ti'에서 t^∈T로 옮긴 해를 π i:ti'→t^라고 정의하자. 이때 전체 알고리즘에서 계산시간을 많이 소모할 수 있으므로 변화된 해의 유효성(feasibility)은 판단하지 않는다. π i:ti'→t^의 정확한 목적함수 값의 변화량 δi,t^π를 다음과 같이 정의한다.

δi,t^π:=fπ i:ti'→t^-fπ.

어떤 해 π에 대해, 모든 작업의 시작 시점을 다른 모든 시점으로 이동할 때의 목적함수 변화량을 모두 계산하면 |I×|T|| 크기의 변화량 행렬(change matrix) Δ_π를 얻는다. 변화량 행렬 Δ_π의 i번째 행(row), t^번째 열(column)은 δi,t^π의 값을 가진다. 변화량 행렬 Δ_π은 어떤 해 π를 약간 변화시킬 때 정확한 목적함수가 어떻게 변화할지 알려주는 역할을 한다. 즉, 변화량 행렬 Δ_π은 어떤 해 π에 대해, 작업의 시작 시점을 변경하는 경우에 대한 목적함수 민감도 정보를 준다. 본 연구에서는 목적함수 섭동을 위해 변화량 행렬 Δ를 사용한다.

목적함수 섭동 알고리즘은 수리모형 (P2)의 최적해를 기반으로 유효해 π~를 구하고 π~에 대한 변화량 행렬 Δπ~을 이용해 목적함수에 변동을 주어 수리모형 (P2)를 반복적으로 풀면서 새로운 유효해를 구한다. 새롭게 구해지는 유효해를 π^라고 하자. 목적함수 섭동 알고리즘은 fπ^가 fπ~보다 클 경우, 다시 말해 더 이상의 해의 개선이 이뤄지지 않을 때까지 실행된다. 유효해 π~의 변화량 행렬 Δπ~를 이용한 목적함수 섭동의 결과, 수리모형 (P2)의 목적함수 (26)은 다음과 같이 변동이 더해진 목적함수 (29)로 대체되어 수리모형을 풀게 된다.

min αG+1-αJ+∑i∈I ∑t^∈Tδi,t^π~xi,t^.

(29)

<Table 1>은 부록 1 (Appendix 1)의 Figure 9의 예시에서 사용되는 해 π에 대한 변화량 행렬 Δ_π의 예시를 나타내며, 이 예시의 경우 (P2)의 목적함수 (26)에 변동이 더해진 목적함수 (30)이 새로운 목적함수로 사용된다.

Table 1. 
change matrix Δ_π of example solution π

I	t^=1	t^=2	t^=3
T
i1	0	-3	2
i2	0	-2	-1
i3	-1	0	1

수리-발견적 국소탐색 알고리즘의 의사 코드는 <Figure 5>와 같다.

min αG+1-αJ-3x1,2+2x1,3-2x2,2-x2,3-x3,1+x3,3.

(30)

Figure 5. 
Pseudocode of Matheuristic Local Search Algorithm

4.2 메타-휴리스틱 전역탐색 (Meta-heuristic Global Search Algorithm)

제4.1절에서 제시된 수리-발견적 국소탐색 알고리즘은 비교적 빠른 시간 내에 좋은 유효해를 도출할 수 있다. 하지만 수리모형 (P2)는 목적함수로서 어떤 해 π에 대한 정확한 위험평가함수인 αg(π)+(1-α)h(π)를 고려하는 대신αg(π)+(1-α)j(π) 를 고려하기 때문에 최적의 해를 구하지 못할 수 있다. 또한, 수리모형 (P2)의 목적함수에 변동을 더해서 수리모형을 반복적으로 푸는 목적함수 섭동 알고리즘의 특성상, 새롭게 구해지는 유효해 π^는 경우에 따라 국소 최적해(local optima)에 머무를 수 있다. 본 절에서는 국소 최적해 상태에서 벗어나 전역 최적해가 존재할 수 있는 탐색 공간(search space)을 고려하는 SA 알고리즘을 활용한 메타-휴리스틱 전역탐색 알고리즘(이하 MHGSA)을 제시한다. 다양한 메타-휴리스틱 알고리즘 중에 특별히 SA를 사용하는 가장 큰 이유는 SA가 해 탐색과정에서 유효하지 않는 해(infeasible solutions)로의 이동을 확률적으로 허용하기 때문이다. 본 알고리즘은 특성상 이웃해 생성 단계에서 유효하지 않은 해를 생성할 가능성이 매우 높다. 그러나 더 좋은 해를 얻기 위해서는 유효하지 않은 해의 이웃해를 탐색하여 유효해를 얻을 필요가 있고, SA는 이러한 과정을 간단한 파라메터를 통해 조절할 수 있는 방법을 제공한다.

SA 알고리즘은 현재해에 대한 국소탐색 알고리즘을 통해 이웃해를 찾고, 현재해에서 이웃해로 이동할 것인지를 메트로폴리스 규칙(metropolis rule)에 따라 확률적으로 결정하는 과정을 반복하며 해를 개선한다. MHGSA에서는 현재해에 대한 여러 개의 이웃해를 찾고, 각 이웃해들에 대해 병렬적으로 메트로폴리스 규칙을 적용함으로써 해의 개선속도를 향상시킨다.

4.2.1 국소탐색을 통한 이웃해 생성

SJAP의 모든 해의 집합을 Π라고 정의하자. 국소탐색 알고리즘은 어떤 해 π∈Π로부터 이웃해 π'∈Π로의 이동을 반복함으로써 이뤄지며 이를 위해 4.1.1항에서 정의한 변화량 행렬 Δ_π를 사용한다.

어떤 해 π∈Π에 대해, 모든 i∈I와 t^∈T에 대한 π i:ti'→t^들 중 δi,t^π가 0보다 작거나 같으면서 제약 유형 3을 어기지 않는 π i:ti'→t^의 집합 Π'을 다음과 같이 정의하자.

Π':π i:ti'→t^δi,t^π≤0, ∀i∈I, ∀t^∈T \Tπ,iE.

이때, Tπ,i1E:=∪i2∈I \i1:t,i1,i2∈Etti2'≤t≤ti2'+di2,ti2'-1로 정의되며 TEπ,i1는 어떤 해 π에서 작업 i₁이 모든 작업 i₂∈I＼{i₁}에 대해 제약 유형 3을 만족하게 할 수 있는 i₁의 작업 수행 시점의 집합을 의미한다. 그러면 Π'는 π에 대해 어떤 작업 i∈I의 작업 수행 시작 시점을 ti'에서 t^∈T \Tπ,iE로 옮긴 제약 유형 3을 만족하는 이웃해 π i:ti'→t^ 중, 실제 위험평가함수 값이 개선될 수 있는 이웃해들의 집합을 의미한다. 또한 π'∈Π'는 t'∈T에 시작한 작업 i∈I가 시점 t∈{t', t'+1,⋯, t'+di,t'-1}에서 수행된다는 의미를 갖기 때문에 제약 유형 1을 만족함을 알 수 있다. 하지만 제약 유형 2에 대한 유효성(feasibility)은 여전히 보장되지 않기 때문에 해의 개선 여부를 판단하기 위해 f(π)와 f(π')을 직접 비교한다면 유효해를 보장하지 못할 수 있다. 따라서 제약 유형 2를 고려할 수 있도록 적합도 함수(fitness function) f'(π')을 다음과 같이 정의한다.

f'π'=fπ'+∑r∈R∑t∈Tmax0,∑i∈Itbi,ti'r,t-utr.               +∑r∈R∑t∈Tmax0,ltr-∑i∈Itbi,ti'r,t

적합도 함수 f'(π')는 기존 위험평가함수인 f(π')와 부적합도(penalty)의 합으로 이뤄진다. 이때, 부적합도는 주어진 해 π'가 모든 자원 r∈R에 대해 어긴 소모량의 합을 의미한다.

4.2.2 메트로폴리스 규칙을 통한 현재해 이동

SA 알고리즘에서는 현재해 π를 통해 이웃해 π'를 생성했을 때 허용비율(accept ratio)에 따라 이웃해로의 이동이 이뤄지며, 허용비율은 π, π'의 적합도 함수 값과 온도(temperature) 파라메터 Γ에 따라 아래와 같은 확률로 정의된다.

Pr.=1, if f'π'≤f'π,exp-f'π'-f'πΓ, otherwise..

메트로폴리스 규칙은 이웃해 집합 Π'에 대해 병렬적으로 적용되며, 허용비율에 의해 이동하기로 결정된 이웃해 중 가장 좋은 적합도 함수 f'(π')값을 가지는 해 π'으로 이동한다.

한 번의 국소탐색을 통한 이웃해 생성과 메트로폴리스 규칙을 통한 현재해 이동이 시행되는 것을 하나의 반복(iteration)이라고 한다.

4.2.3 초기 온도(Initial Temperature) 및 동적 냉각 계획(Dynamic Cooling Schedule)

SA 알고리즘에서는 허용비율을 조절하는 온도 파라메터 Γ가 작아질수록 전역 최적해로 수렴하며, 이를 위해 초기온도와 냉각속도를 결정해야 한다. 일반적으로 SA 알고리즘에서 초기온도는 대부분의 해에 대해 이동을 허용할 수 있도록, 즉 허용비율을 1에 가깝게 초기온도를 높게 설정한다. 이후 반복할 때마다 냉각비율 β에 따라 온도 Γ를 감소시킨다.

본 연구에서는 초기온도 Γ₀를 0에 가까운 값을 가지도록 1fπ~0로 설정한다. 이는 허용비율을 매우 낮게 해서 수리-발견적 국소탐색 알고리즘을 통해 구한 좋은 품질의 유효해에 대한 해 연마(solution polishing)를 하기 위해서다. 이때, 초기온도를 낮게 설정했기 때문에 Γ가 빠르게 0으로 수렴하게 된다. 이는 지역 최적해로의 수렴으로 이어질 수 있으므로, 지정된 임계 값 값 ϵ보다 Γ가 작아질 경우 온도를 다시 올려주어 허용비율을 높일 수 있는 동적 냉각 계획을 적용한다. 본 연구에서는 수리-발견적 국소탐색 알고리즘을 통해 도출된 유효해가 π~0 일 때, 초기온도 Γ0=1fπ~0, 냉각비율 값 β = 0.95, 임계 값 ϵ=1105로 설정한다. 온도 Γ가 임계 값 ϵ보다 작아질 경우, 해당 시점의 해 π'에 대해 Γ∼U0,f'π'에 따라서 온도를 조정한다.

4.2.4 중단 조건(Termination Criteria)

중단 조건은 전체 알고리즘이 끝나는 조건으로, 이론적인 SA 알고리즘의 중단 조건은 온도 Γ가 0에서 수렴할 때이다. 하지만 본 연구의 알고리즘에서는 동적 냉각 계획을 사용하기 때문에 매 반복이 끝나는 시점의 적합도 함수 값의 개선 여부를 비교하여 일정 시간 동안 해의 개선이 없으면 전체 알고리즘을 중단한다. 본 논문의 실험결과에서는 200초를 기준으로 한다.

SA 알고리즘을 활용한 메타-휴리스틱 전역탐색 알고리즘의 의사코드는 다음 <Figure 6>과 같다.

Figure 6. 
Pseudocode of Meta-heuristic Global Search Algorithm

본 장에서는 제안한 수리계획-발견적 복합 알고리즘에 대한 계산 실험결과를 소개한다. 모든 실험은 iMac 3.5GHz 4-Core Intel Core i5 프로세서, 16GB 메모리 환경에서 Python 3.8을 이용해서 수행되었으며, 병렬 연산의 경우 최대 4개의 thread를 이용했다.

5.2 알고리즘별 성능 비교

본 절에서는 제3.2절에서 제시한 MIP 수리모형 (P)를 이용한 최적해 알고리즘(exact algorithm)과 4장에서 제시한 수리계획-발견적 복합 알고리즘의 결과를 비교한다. MIP 수리모형을 이용한 최적해 알고리즘의 실험과 수리계획-발견적 복합 알고리즘의 근사 수리모형에 대해 MIP solver인 Cplex 20.10의 Python 라이브러리를 사용하여 각각 (P)와 (P2)를 풀었다. 수리계획-발견적 복합 알고리즘에서 사용되는 변화량 행렬의 계산과 메타-휴리스틱 전역탐색 알고리즘의 경우 JIT(Just-In-Time) 컴파일러를 사용하는 numba(http://numba.pydata.org/)를 이용하여 Python으로 구현되었다. 두 알고리즘 모두에 대해 제한시간은 7200초로 설정했다.

<Table 5>, <Table 6>, <Table 7>은 앞 절에서 제시된 문제들에 대한 계산 실험결과를 보여준다. 각각의 표에서 최적해 알고리즘에 대한 결과는 MIP 열에, 수리계획-발견적 복합 알고리즘의 결과는 HMPHA열에 나타난다. 알고리즘별 전체 수행시간인 Tot_time 열은 초 단위로 측정되었으며, 알고리즘별 최종 목적함수 값은 Obj. 열에 나타난다. 수리계획-발견적 복합 알고리즘에서 수리-발견적 국소탑색 알고리즘을 이용해 초기 유효해 구하기까지의 시간이 Init_time 열에, 해당 유효해의 목적함수 값이 Obj. 열에 나타난다. 해당 실험에서는 위험 수용 정책 가중치 α를 0.5로 고정하여 사용했으며, 이는 목적함수인 위험평가함수 f에서 기대 위험 g와 기대 초과 h를 동일한 비율로 고려함을 의미한다.

Table 5. 
Experimental Results for Exact Algorithm and Proposed Heuristic Algorithm on A Set

instance	MIP		HMPHA
	Tot_time(sec)	Obj.	Init_time(sec)	Init_Obj.	Tot_time(sec)	Obj.
A_01	3.9	1767.82	18.8	1767.82	218.8	1767.82
A_02	>7200	-	9.2	4727.13	3830.7	4680.57
A_03	0.8	848.18	7.0	848.18	207.0	848.18
A_04	639.8	2085.88	3222.5	2085.88	3422.5	2085.88
A_05	>7200	-	71.6	646.42	3945.9	636.03
A_06	9.4	590.62	64.4	590.62	264.4	590.62
A_07	0.1	2272.78	0.2	2279.99	432.5	2272.78
A_08	>7200	756.16	0.3	749.95	250.0	746.93
A_09	0.1	1507.28	0.1	1507.32	250.0	1507.28
A_10	5.9	2994.85	2.5	3015.74	2192.6	2995.08
A_11	>7200	-	2.8	504.37	1117.2	495.48
A_12	4.1	789.63	0.9	790.07	570.0	789.63
A_13	>7200	-	17.6	2009.50	3854.2	2000.15
A_14	>7200	-	4.3	2309.90	2112.6	2269.54
A_15	>7200	-	5.5	2315.81	2143.7	2281.68

Table 6. 
Experimental Results for Exact Algorithm and Proposed Heuristic Algorithm on B Set

instance	MIP		HMPHA
	Tot_time(sec)	Obj.	Init_time(sec)	Init_Obj.	Tot_time(sec)	Obj.
B_01	>7200	-	20.7	4023.63	377.1	3996.30
B_02	>7200	-	14.3	4326.85	3831.1	4299.64
B_03	>7200	-	572.5	35294.10	3864.5	35293.85
B_04	>7200	-	267.6	34852.78	4312.9	34829.29
B_05	>7200	45329.99	75.9	2426.48	3342.4	2399.25
B_06	>7200	-	16.1	4309.37	3954.9	4289.68
B_07	>7200	-	47.7	7601.88	5847.2	7554.61
B_08	>7200	-	12.2	7437.39	1522.9	7435.72
B_09	>7200	-	14.6	7512.74	3808.9	7496.88
B_10	>7200	-	59.1	10772.47	4039.4	10658.11
B_11	>7200	-	21.1	3651.17	3809.9	3626.66
B_12	>7200	-	241.6	38479.87	7200	37883.86
B_13	>7200	-	32.5	5208.43	4039.8	5052.90
B_14	>7200	-	338.3	12044.07	7200	11963.19
B_15	>7200	-	1093.1	23172.46	7200	22868.60

Table 7. 
Experimental Results for Exact Algorithm and Proposed Heuristic Algorithm on C Set

instance	MIP		HMPHA
	Tot_time(sec)	Obj.	Init_time(sec)	Init_Obj.	Tot_time(sec)	Obj.
C_01	>7200	-	19.9	8519.93	688.2	8518.33
C_02	>7200	-	21.6	3587.73	3849.9	3553.08
C_03	>7200	-	468.5	33561.40	3990.8	33547.98
C_04	>7200	-	272.6	37620.84	3285.9	37595.91
C_05	>7200	-	71.7	3201.21	2448.5	3177.68
C_06	>7200	-	95.2	8504.24	4303.2	8432.07
C_07	>7200	-	21.3	6093.66	3961.7	6086.94
C_08	>7200	-	103.5	11380.17	4607.6	11242.23
C_09	>7200	-	26.0	5648.32	3986.4	5609.23
C_10	>7200	-	262.6	44652.37	7200	44005.87
C_11	>7200	-	35.0	6003.97	3968.1	5821.38
C_12	>7200	-	333.4	13001.96	7200	12871.71
C_13	>7200	-	1589.2	43334.10	7200	43165.61
C_14	>7200	-	1128.1	27079.45	7200	26925.44
C_15	>7200	-	1994.1	43484.35	7200	42630.78

Cplex의 경우, 문제 그룹 A에 속한 크기가 작은 문제들에 대해서 빠른 시간 내에 최적해를 찾았다. 특히, 최적해를 구한 대부분의 경우에는 휴리스틱 알고리즘보다 빠르게 최적해를 찾을 수 있었다. 하지만 최적해 알고리즘은 문제 크기가 커질수록 기하급수적으로 시간이 증가하여 문제 그룹 B, C의 문제에 대해서는 한 문제를 제외하고는 정수 해(integer solution)를 찾지 못했다. 반면 수리계획-발견적 복합 알고리즘은 모든 문제에 대한 유효해(feasible solution)를 찾았으며, 최적해 알고리즘이 최적해를 찾은 문제 중 A_10을 제외하고 동일한 최적해를 찾을 수 있었다. 또한 수리계획-발견적 국소탐색 알고리즘을 통해 구한 초기 유효해의 목적함수 값은 최적해 알고리즘이 해를 구한 경우 최적해의 목적함수와 근접한 값을 가지며 질 좋은 초기 유효해를 빠른 시간 내에 구할 수 있음을 보여준다.

<Table 8>, <Table 9>, <Table 10>은 수리계획-발견적 복합 알고리즘의 단계별 소요 시간의 비율과 초기 유효해에 대한 목적함수 개선율을 비교해서 보여준다. 각각의 표에서 수리-발견적 국소탐색 알고리즘 종료 시점에 대한 값은 MHLSA(Phase 1) 열에 나타나고, 메타-휴리스틱 전역탐색 알고리즘 종료 시점에 대한 값은 MHGSA(Phase 2) 열에 나타난다.

Table 8. 
Comparison of experimental results for MHLSA and MHGSA on A set

instance	Time (ratio)		Obj. (improvement)
	MHLSA (Phase 1)	MHGSA (Phase 2)	MHLSA (Phase 1)	MHGSA (Phase 2)	Total improv.
A_01	8.58%	91.42%	0.00%	0.00%	0.00%
A_02	20.49%	79.51%	0.63%	0.36%	0.98%
A_03	3.39%	96.61%	0.00%	0.00%	0.00%
A_04	94.16%	5.84%	0.00%	0.00%	0.00%
A_05	24.58%	75.42%	1.41%	0.20%	1.61%
A_06	24.35%	75.65%	0.00%	0.00%	0.00%
A_07	16.14%	83.86%	0.08%	0.24%	0.32%
A_08	0.11%	99.89%	0.33%	0.07%	0.40%
A_09	0.03%	99.97%	0.00%	0.00%	0.00%
A_10	0.18%	99.82%	0.63%	0.06%	0.69%
A_11	0.56%	99.44%	1.60%	0.17%	1.76%
A_12	64.91%	35.09%	0.05%	0.00%	0.05%
A_13	48.23%	51.77%	0.45%	0.02%	0.47%
A_14	0.36%	99.64%	1.62%	0.13%	1.75%
A_15	0.39%	99.61%	1.16%	0.32%	1.47%
Avg.	20.43%	79.57%	0.53%	0.10%	0.63%

Table 9. 
Comparison of Experimental Results for MHLSA and MHGSA on B Set

instance	Time (ratio)		Obj. (improvement)
	MHLSA (Phase 1)	MHGSA (Phase 2)	MHLSA (Phase 1)	MHGSA (Phase 2)	Total improv.
B_01	46.97%	53.03%	0.68%	0.00%	0.68%
B_02	64.31%	35.69%	0.62%	0.00%	0.63%
B_03	19.40%	80.60%	0.00%	0.00%	0.00%
B_04	33.93%	66.07%	0.05%	0.02%	0.07%
B_05	24.47%	75.53%	1.08%	0.04%	1.12%
B_06	78.04%	21.96%	0.37%	0.09%	0.46%
B_07	7.50%	92.50%	0.30%	0.33%	0.62%
B_08	86.87%	13.13%	0.02%	0.00%	0.02%
B_09	32.18%	67.82%	0.21%	0.00%	0.21%
B_10	4.33%	95.67%	0.74%	0.32%	1.06%
B_11	93.69%	6.31%	0.64%	0.03%	0.67%
B_12	18.46%	81.54%	1.13%	0.43%	1.55%
B_13	24.69%	75.31%	2.42%	0.58%	2.99%
B_14	18.36%	81.64%	0.32%	0.35%	0.67%
B_15	15.19%	84.81%	0.00%	1.31%	1.31%
Avg.	37.89%	62.11%	0.57%	0.23%	0.80%

Table 10. 
Comparison of Experimental Results for MHLSA and MHGSA on C Set

instance	Time (ratio)		Obj. (improvement)
	MHLSA (Phase 1)	MHGSA (Phase 2)	MHLSA (Phase 1)	MHGSA (Phase 2)	Total improv.
C_01	70.94%	29.06%	0.02%	0.00%	0.02%
C_02	3.33%	96.67%	0.86%	0.11%	0.97%
C_03	30.27%	69.73%	0.04%	0.00%	0.04%
C_04	12.28%	87.72%	0.03%	0.04%	0.07%
C_05	11.43%	88.57%	0.72%	0.02%	0.73%
C_06	46.72%	53.28%	0.52%	0.33%	0.85%
C_07	81.79%	18.21%	0.09%	0.02%	0.11%
C_08	38.65%	61.35%	0.82%	0.40%	1.21%
C_09	64.58%	35.42%	0.65%	0.05%	0.69%
C_10	24.99%	75.01%	1.28%	0.17%	1.45%
C_11	87.95%	12.05%	0.59%	2.47%	3.04%
C_12	6.80%	93.20%	0.20%	0.80%	1.00%
C_13	22.10%	77.90%	0.00%	0.39%	0.39%
C_14	15.68%	84.32%	0.00%	0.57%	0.57%
C_15	27.72%	72.28%	0.00%	1.96%	1.96%
Avg.	36.35%	63.65%	0.39%	0.49%	0.87%

Time(ratio) 열은 단계별 알고리즘의 수행시간이 전체 알고리즘의 수행시간에서 차지하는 비율을 나타낸다. Obj.(improvement) 열은 이전 단계 종료 시점의 목적함수 값 대비 감소율을 나타낸다. 이때, 수리-발견적 국소탐색 알고리즘의 해의 비교 대상은 초기 유효해이다. Total improv. 열은 전체 알고리즘이 종료되었을 때 초기 유효해의 목적함수 값을 얼마나 감소시켰는지를 감소율로 나타낸다.

표에서 확인할 수 있듯이 전체 알고리즘을 구성하는 세부 알고리즘별 차지하는 시간 비율과 목적함수의 감소율은 문제의 그룹에 따라 다른 것을 확인할 수 있다. 문제의 크기가 비교적 작은 A그룹의 경우, Phase 1이 전체 알고리즘의 수행시간에서 차지하는 비율은 평균 약 20%이며, Phase 2의 경우 수행시간은 약 80%를 차지한다. 반면 상대적으로 문제의 크기가 큰 B, C그룹의 경우 수리-발견적 국소탐색 알고리즘이 전체 알고리즘의 수행시간에서 차지하는 비율은 각각 평균 35% 이상 차지하며 A그룹에 비해 해당 비율이 늘어남을 확인할 수 있다.

목적함수 감소율 측면에서, A그룹의 경우 Phase 1에서 초기 유효해를 평균 0.53% 개선하고 Phase 2에서 Phase 1의 해를 평균 0.10% 개선하는 것을 확인할 수 있다. 반면 B그룹의 경우 목적함수는 Phase 1과 Phase 2에서 각각 평균 0.57%, 0.23% 감소하고, C그룹의 경우 각각 평균 0.39%, 0.49% 감소하는 것을 확인할 수 있다. 이러한 결과는 문제의 크기에 따라 세부 알고리즘별 해의 개선에 차지하는 비율이 유동적으로 변하면서 해를 효과적으로 개선할 수 있음을 보여준다.

해당 실험결과에 대한 세부사항을 정리한 <Table 12>가 부록 2 (Appendix 2)에 수록되어 있다.

<Figure 7>은 수리계획-발견적 복합 알고리즘이 해를 개선하는 추이를 보여준다. 그림에서 수리-발견적 국소탐색 알고리즘(MHLSA)에 따른 목적함수 값에 해당하는 부분은 점선으로, 메타-휴리스틱 전역탐색 알고리즘(MHGSA)에 따른 목적함수 값에 해당하는 부분은 실선으로 나타냈으며 알고리즘의 전환이 이루어지는 시점의 목적함수 값은 빨간색 점으로 표현했다.

Figure 7. 
Performance of the Hybrid Mathematical Programming-Heuristic Algorithm

그림에서 확인할 수 있듯이 수리계획-발견적 복합 알고리즘은 수리-발견적 국소탐색 알고리즘을 통해 빠른 시간내에 유효해를 구하고 개선하며, 이후 메타-휴리스틱 전역탐색 알고리즘을 통해 추가적으로 해를 개선한다. 또한 문제에 따라 각 알고리즘이 해를 개선하는 정도가 다른 것을 확인할 수 있다. Figure 7의 (b), (c), (e)에 해당하는 문제에 대해서는 수리-발견적 국소탐색 알고리즘이 해를 개선하는 데 있어서 큰 비중을 차지했다. 반면 (a), (d), (f)에 해당하는 문제에 대해서는 메타-휴리스틱 전역탐색 알고리즘이 큰 비중을 차지한다.

이러한 결과는 제시하는 수리계획-발견적 복합 알고리즘이 추계적 위험 목적함수를 고려하는 SJAP에 대해 많은 시나리오 개수를 갖는 경우에도 효율적임을 보여준다.

5.3 위험 수용 정책 가중치의 영향에 따른 영향 비교

제3장에서 설명했듯이 위험평가함수는 두 항의 조합으로 정의된다. 이는 두 개의 독립적인 목적함수를 가지는 경우로 볼 수 있고 두 목적함수의 중요도는 의사결정 파라메터 α에 의해 결정된다. 본 절에서는 의사결정자의 위험 수용 정책에 따라서 해가 어떻게 변화하는지 확인하기 위해 같은 문제를 의사결정 파라메터 α를 변화시켜 가면서 반복적으로 푼 결과를 제시한다.

<Figure 8>은 각각의 문제에 대해 위험 수용 정책 가중치 α가 0.1, 0.2, ⋯, 0.9로 변화할 때, 각 α에 따라 구해지는 해의 기대 위험 값 g와 기대 초과 값 h의 변화를 보여준다. α가 커질수록 g는 작은 값을, h는 큰 값을 가지는 것을 확인할 수 있다. 이는 평균적인 측면에서 위험 척도를 더 많이 고려하는 정책일 경우 기대 위험 g가 작은 해를 도출하는 것이 더 이익임을 의미하고, 반대로 위험 척도를 변동성의 측면에서 더 많이 고려하는 정책일 경우 기대 초과 h가 작은 해를 도출하는 것이 더 이익임을 의미한다. 이러한 결과는 목적함수로 사용되는 위험평가함수 f가 αg+(1-α)h로 정의되었기에 자명하다. Figure 8은 일반적으로 다-목적함수 최적화 문제(Multi-criteria optimization problems)에서 많이 사용되는 효율해(efficient solutions 또는 Pareto solutions)의 경계(frontier)를 표시한다. 그래프의 효율해를 연결한 선을 기준으로 오른쪽 위의 해들은 모두 효율해보다 두 가지 목적함수 모두 좋지 않은 해들에 해당한다. 따라서 정책 수립자는 효율해의 경계(efficient frontier) 위에서 해를 선택해야 하며 어떤 해를 선택할지는 상황에 따라 다르다. 예를 들어 파라메터 α가 큰 경우는 평균적인 위험 함수를 더 중요시하겠다는 의미이고, 반대로 α가 작은 경우에는 τ-분위 위험 함수에 더 큰 가중치를 부여하여 평균적인 위험보다 위험의 변화를 줄이는 경향의 해를 선택하겠다는 의미이다.

Figure 8. 
Distribution of g, h for different α

본 논문에서는 일정 기간 진행되는 작업들의 수행시기를 결정하는 추계적 작업할당문제를 풀기 위한 효율적인 알고리즘을 제시하였다. 각 작업은 시점마다 위험을 발생시키며, 위험은 시점마다 정해진 시나리오에 따라 측정되기 때문에, 이 문제를 추계적 작업할당문제라고 한다. 특히 목적함수는 작업 수행에 따라 발생하는 두 가지 위험요소(평균 위험과 τ-분위 위험)를 동시에 고려하며, 이는 실제 전력망 운영회사(RTE)가 제시하는 위험평가함수이다.

이 문제는 NP-hard에 속하는 문제이며, 문제의 크기가 커질수록 선형 최적화 기반 분지한계법을 통해 최적해를 얻는 시간이 지수적으로 증가한다. 이를 해결하기 위해 본 논문에서는 수리계획-발견적 복합 알고리즘을 제시한다. 제안하는 알고리즘은 추계적 작업할당문제의 특성을 이용하여 유효해를 얻는 간소화된 수리계획문제와 이후 해를 개선하는 발견적 기법의 두 단계로 구성되어 있다. 특히, 두 번째 단계인 해 개선 발견적 기법 알고리즘은 국지해의 개선 방향을 구하는 방법과 국지 최적해를 탈출하기 위한 방법을 조합하여 개발되었다. 실험결과, 제시하는 수리계획-발견적 복합 알고리즘이 작은 크기의 문제들에 대해 대부분의 경우 최적해를 구할 수 있었으며, 큰 크기의 문제들에 대해서는 MIP 수리모형을 이용한 Cplex의 성능을 능가함을 확인할 수 있었다.

본 연구가 제시하는 알고리즘은 최적해를 구하기는 어렵지만, 수리모형을 이용해서 유효해를 비교적 쉽게 구할 수 있는 경우에 적용할 수 있다. 예를 들어, 최적화 문제의 목적함수가 다루기 어려운 경우(비선형, 추계적 목적함수 등)에는 이 목적함수를 표현하기 위해 많은 결정변수와 제약식이 추가될 수 있다. 이런 경우에 목적함수를 간략화하고 유효해만을 표현하는 수리모형을 수립한 후 이를 이용해 정확한 목적함수 값을 개선하는 발견적 기법을 적용할 수 있다.