지금까지의 선형계획 관련 포스팅들에서는, 모든 변수들을 ‘음이 아닌 실수’ 라고 가정했다. 그러나, 만약 규격 용량이 정해진 발전기를 도입한다면, ‘이 발전기를 3.5대 도입하는 것이 최적이다’ 라고 보고하는 것은 비현실적이다. 발전기 대수는 3대 또는 4대이기 때문이다.

이런 발전기의 ‘대수’는 정수 (integer) 이기 때문에, 현실적인 결과를 원한다면 최적해에서도 정수로 도출되어야 한다. 이번 포스팅에서는 특정 변수가 정수로 강제되는 선형계획 문제에 대해 설명한다.


정수 변수 예시: 풍력 발전기의 ‘대수’

정수 변수 사용 예시로, 필자가 Applied Energy에 제1저자로써 게재했던 논문의 최적화 문제를 들고자 한다. 해당 논문에 대해 필자가 이미 설명해놓은 포스팅이 있으니, 참고 바란다.

대상 에너지시스템은 태양광 + 풍력 + 디젤 + ESS로 전기부하를 공급하며 중앙계통과는 분리되어 있다 (녹색섬이라 생각하자).

혹시 재생발전량이 0이더라도, 디젤발전기만으로 부하 충족 자체는 가능하다. 반대로 재생발전량이 너무 많아 부하 충족은 물론 ESS 완충까지 하더라도 전기가 남을 경우 이는 특별한 용도 없이 dump load에서 소비된다, 즉 출력제한(curtail)된다. 이 때 녹색섬이라 불리려면, 전체 전기에너지 수요의 일정 비율은 재생발전으로 충당해야 한다.

islandsystem 태양광/ 풍력/ 디젤/ ESS로 이루어진 섬 전기 공급 시스템 도식도. 잉여전력을 버릴 수 있음.

특정 기간의 시간별 전기부하, 일사량, 풍속을 input 자료로 사용하여, 재생발전기 초기설치비용과 운영기간 내 디젤 연료비의 총 합을 최소화하는 태양광/ 풍력/ ESS 각각의 용량을 도출하는 최적화 문제를 풀어서 결정할 수 있다. 수식은 아래와 같다.

Minimize $c_{\text{pv}} s_{\text{pv}} + c_{\text{wt}} s_{\text{wt}} + c_{\text{batt}} s_{\text{batt}} + c_{\text{diesel}} \sum_{t} p_{\text{diesel}}[t] + \epsilon \sum_t ( p_{\text{disch}}[t] + p_{\text{ch}}[t] )$

Subject to

$p_{\text{load}}[t] = p_{\text{diesel}}[t] + s_{\text{pv}} p_{\text{pv}}[t] + s_{\text{wt}} p_{\text{wt}}[t] + p_{\text{disch}}[t] – p_{\text{ch}}[t] – p_{\text{dump}}[t] \quad \forall t$

$e_{\text{batt}}[t+1] = e_{\text{batt}}[t] + \mu p_{\text{ch}}[t] – p_{\text{disch}} /\mu \quad \forall t$

$ e_{\text{batt}}[t] \leq s_{\text{batt}} \quad \forall t$

$ p_{\text{ch}}[t] \leq \gamma s_{\text{batt}} \quad \forall t$

$ p_{\text{disch}}[t] \leq \gamma s_{\text{batt}} \quad \forall t$

$ \sum_{t} p_{\text{diesel}}[t] \leq (1-\alpha_{\text{renewable}}) \sum_{t} p_{\text{load}}[t] $

위에서 $c$는 비용, $s$는 용량, $p$는 전력, $e$는 저장된 에너지, $\mu$는 효율, $\alpha$는 규정 비율이며, pv는 태양광, wt는 풍력, batt는 배터리, load는 부하, diesel은 디젤엔진, ch는 충전, disch는 방전, dump는 출력제한이다. $s_{\text{wt}}$는 풍력 발전기의 대수이므로 음이 아닌 정수(nonnegative integer)이며, 나머지는 음이 아닌 실수들이다. 비용은 태양광/ 풍력/ 배터리에 대해서는 단위 용량당 초기투자비, 디젤에 대해서는 연료비이다.

$\epsilon$은 ‘같은 시간에는 충전과 방전 둘 중 하나 이상은 0’이 되게 만드는 데 필요한 penalty term의 계수인데, 이에 대해서는 필자의 논문 설명 포스팅을 참고하기 바란다.

이 문제에서 도출한 태양광/ 풍력/ 배터리 용량은, 디젤 발전기를 통해 공급되는 전기에너지의 비중이 전체의 $1-\alpha_{\text{renewable}}$ 이하가 되게 만드는 용량이다. 이를테면 $\alpha_{\text{renewable}}=0.6$인 경우, 연간 총 전기에너지 수요 중 최대 40%까지만 디젤 발전기로 충당될 수 있다는 제약 하에 시스템 구성을 도출한다.

이 때, 풍력 대수 변수 $s_{\text{wt}}$는 ‘정수로 강제’한다.

(태양광과 배터리의 경우 구성 단위인 모듈 당 용량이 작으므로, $s_{\text{pv}}$와 $s_{\text{batt}}$를 연속변수로 봐도 실제 설치 가능 용량과 큰 차이가 나지 않을 것이라 보았다.)

참고로 풍력 대수 변수 $s_{\text{wt}}$의 계수 $p_{\text{wt}}[t]$는 풍력 ‘1대당’ 시간별 발전량이며, 이는 풍속과 풍력발전의 성능곡선을 통해 계산된 전력에 변환기/ 변압기 효율을 곱한, 즉 부하에 최종적으로 전달되는 교류전력 값이다.

시간별 풍력 발전량 데이터가 이상하지 않은지 확인하기 위해서는, 이전 포스팅에서 언급한 ‘이용률’을 체크해 보아야 한다. 풍력의 이용률은 장소 및 터빈 높이에 따라 그 편차가 크지만, 국내에서는 대략 15~40% 정도이다.


풍력 대수 변수 (정수) 를 포함하는 문제 코드

위 문제를 Python cvxopt로 푸는 코드는 아래와 같다. (GitHub Repo 링크)

from cvxopt import matrix, spmatrix, sparse, glpk
import numpy as np
import time

penetration = 0.6
eff_batt = 0.94
cost_diesel= 300 # 디젤발전기로 발전된 전기 1kWh당 투입연료비
c_rate_batt = 0.5
initialenergy_batt = 0
penaltycoeff = 0.01 # 동시 충/방전 오류를 막기 위한 penalty term의 계수 (너무 낮으면 오류 생김)

r = 0.05 # 화폐가치 하락에 대한 할인율

initialcost_pv = 3000000 # 태양광 1kW당 초기투자비
initialcost_wt = 80000000 # 풍력 1대당 초기투자비
initialcost_batt = 900000 # 배터리 1kWh당 초기투자비 (SOC 고려한 실용량 가정)
mtncost_pv = initialcost_pv * 0.02 # 태양광 1kW당 유지보수비 (매년 발생) 
mtncost_wt = initialcost_wt * 0.02 # 풍력 1대당 유지보수비 (매년 발생)
mtncost_batt = initialcost_batt * 0.01 # 배터리 1kWh당 유지보수비 (매년 발생)


t=8760

demand_island = np.loadtxt("load_hourly.txt") 
renewable_output = np.loadtxt("renewables_hourly.txt")  
totaldemand = matrix(np.sum(demand_island))

pvoutput = matrix(renewable_output[0:t,0]) # 1kW 태양광 패널의 시간별 전기 출력
wtoutput = matrix(renewable_output[0:t,2]) # 대상 규격 풍력터빈 1대의 시간별 전기 출력

def block_eye(size):
    return spmatrix(1,range(size),range(size)) 

def block_zeros(row,col):
    return sparse(matrix(np.zeros((row,col))))

def block_ones(row,col):
    return sparse(matrix(np.ones((row,col))))

def block_batt(size):
    return sparse([[block_zeros(size,1)],[block_eye(size)]]) - sparse([[block_eye(size)],[block_zeros(size,1)]])


Aeq_balance = sparse([[block_eye(t)], [-block_eye(t)], [block_eye(t)], [block_zeros(t,t+1)], [-block_eye(t)], [pvoutput],[wtoutput],[block_zeros(t,1)]])
beq_balance = matrix(demand_island[0:t],tc='d')

Aeq_batterydynamic = sparse([[block_zeros(t,t)],[-eff_batt*block_eye(t)],[block_eye(t)/eff_batt],[block_batt(t)],[block_zeros(t,t)],[block_zeros(t,1)],[block_zeros(t,1)],[block_zeros(t,1)]])
beq_batterydynamic = block_zeros(t,1)

A_maximumstored = sparse([[block_zeros(t,t)],[block_zeros(t,t)],[block_zeros(t,t)],[block_zeros(t,1)],[block_eye(t)],[block_zeros(t,t)],[block_zeros(t,1)],[block_zeros(t,1)],[-block_ones(t,1)]])
b_maximumstored = block_zeros(t,1)

A_maximumcharge = sparse([[block_zeros(t,t)],[block_eye(t)],[block_zeros(t,t)],[block_zeros(t,t+1)],[block_zeros(t,t)],[block_zeros(t,1)],[block_zeros(t,1)],[-c_rate_batt*block_ones(t,1)]])
b_maximumcharge = block_zeros(t,1)

A_maximumdischarge = sparse([[block_zeros(t,t)],[block_zeros(t,t)],[block_eye(t)],[block_zeros(t,t+1)],[block_zeros(t,t)],[block_zeros(t,1)],[block_zeros(t,1)],[-c_rate_batt*block_ones(t,1)]])
b_maximumdischarge = block_zeros(t,1)

A_penetration = sparse([[block_ones(1,t)],[block_zeros(1,t)],[block_zeros(1,t)],[block_zeros(1,t+1)],[block_zeros(1,t)],[block_zeros(1,3)]])
b_penetration = matrix([(1-penetration)*totaldemand])

lowerbounds = matrix([[block_zeros(1,t)],[block_zeros(1,t)],[block_zeros(1,t)],[matrix([initialenergy_batt])],[block_zeros(1,t)],[block_zeros(1,t)],[block_zeros(1,3)]]) 
upperbounds = matrix([[np.inf*block_ones(1,t)],[np.inf*block_ones(1,t)],[np.inf*block_ones(1,t)],[matrix([initialenergy_batt])],[np.inf*block_ones(1,t)],[np.inf*block_ones(1,t)],[np.inf*block_ones(1,3)]]) 

A_lowerbound = -block_eye(5*t+4)
b_lowerbound = -lowerbounds.T

A_upperbound = block_eye(5*t+4)
b_upperbound = upperbounds.T

Aeq = sparse([Aeq_balance,Aeq_batterydynamic])
beq = matrix([beq_balance,beq_batterydynamic])

A = sparse([A_maximumstored, A_maximumcharge, A_maximumdischarge, A_penetration, A_lowerbound, A_upperbound])
b = matrix([b_maximumstored, b_maximumcharge, b_maximumdischarge, b_penetration, b_lowerbound, b_upperbound])

npvfactor = (1/(1+r))*(1-(1/(1+r))**20)/(1-(1/(1+r))) # 초항이 (1+r)^{-1}, 마지막항이 (1+r)^{20}인 등비수열의 합
costfun_pv = initialcost_pv + mtncost_pv * npvfactor # 초기투자비 + 매년 유지보수비(npvfactor를 곱해 현가화)
costfun_wt = initialcost_wt + mtncost_wt * npvfactor # 초기투자비 + 매년 유지보수비(npvfactor를 곱해 현가화)
costfun_batt = initialcost_batt * (1 + (1/(1+r))**10) + mtncost_batt * npvfactor # 초기투자비 + 10년 후 재투자비(현가화) + 매년 유지보수비 (현가화)
costfun_diesel = cost_diesel * block_ones(1,t) * npvfactor  # 전기 사용량 요금 (현가화)


c = matrix([[costfun_diesel],[penaltycoeff*block_ones(1,t)],[penaltycoeff*block_ones(1,t)],[block_zeros(1,t+1)],[block_zeros(1,t)],[matrix([costfun_pv])],[matrix([costfun_wt])],[matrix([costfun_batt])]]) # 이거 sparse로 하면 셧다운되는데 matrix로 하면 잘 됨...왜지..

idx_pv = (t+t+t+(t+1)+t+1)-1
idx_wt = (t+t+t+(t+1)+t+2)-1 # 변수 벡터 x에서 풍력터빈 대수 변수의 인덱스
idx_batt = (t+t+t+(t+1)+t+3)-1

start_time = time.time()
(status,x)=glpk.ilp(c,A,b,Aeq,beq,I=set([idx_wt])) # 풍력터빈 대수 변수를 정수로 제한
elapsed_time = time.time() - start_time

p_ch = np.array(x[(t):(2*t)])
p_disch = np.array(x[(2*t):(3*t)])
indicator_error = (p_ch.reshape(1,-1) @ p_disch.reshape(-1,1))[0][0]

if indicator_error > 1:
    print('경고: 같은 시간에 충전과 방전이 동시에 양수가 되는 오류가 발생했습니다.')
print('태양광 {}kW, 풍력 {}기, 배터리 {}kWh입니다'.format(x[idx_pv],x[idx_wt],x[idx_batt]))
print('{}초 걸렸습니다'.format(round(elapsed_time,2)))
print((c*x)[0] - penaltycoeff*np.sum(p_ch+p_disch)) # 총 비용 (penalty term의 값 제함)

위 코드에서 주목할 부분은, $\texttt{glpk.ilp}$ 이다. 정수 변수가 없을 때의 메서드는 $\texttt{glpk.lp}$ 였는데, 마침표 뒤에 $\texttt{i}$가 추가되었다. I는 integer를 의미하며, 이 메서드로 ‘정수’선형계획 (Mixed-integer Linear Programming, MILP) 문제를 푼다.

또한 메서드의 마지막 인수로 $\texttt{ I=set([idx_wt]) }$ 가 추가되었다. 이는 최적화 문제의 모든 변수들을 담는 벡터 $\textbf{x}$ 내에서 정수 변수들의 위치를 담는 리스트이다. 이 문제에서는 $\texttt{idx_wt}$의 값이 43082인데, 이는 $s_{\text{wt}}$가 $\textbf{x}$의 43803번째 원소임을 의미한다 (Python에서는 배열의 인덱스가 0부터 시작하므로).


‘이진수’ 도입을 통한 모델 현실 설명력 증가: 용량에 비례하지 않는 비용

정수 변수 활용 사례는 이 외에도 많다. 여기서는 최적화 문제의 현실 설명력을 높이기 위한 ‘이진수’ 변수 활용 사례를 설명한다.

지금까지의 분석에서는 초기투자비가 설비 용량에 정비례한다고 암묵적으로 가정했다. 그러나 실제로는 용량에 정비례하지 않는 비용이 존재할 수 있다. 경우에 따라선 이 비용을 경제성 분석에 반영해야 할 수 있다.

이를 테면 위 사례에서, 풍력을 설치하지 않으면 전혀 들지 않지만, 풍력을 한 대라도 설치하면 고정적으로 필요한 유관 비용 $c_{\text{misc}}$가 존재한다고 가정하자. 그리고 일단 풍력을 도입한다면, 이 유관 비용은 풍력 대수가 몇 대이든 비슷하다고 가정하자.

(이는 풍력발전 도입 사업 준비 과정에서의 풍황계측/ 타당성조사/ 인허가취득 관련 soft cost일 수도 있고, 계측/ 통신 및 안전 관련 인프라 구축 비용일 수도, 부지 관련 공사 비용일 수도 있다. 필자도 사실 정확히 아는 건 아니다만… 아무튼 일반적으로 설비 규모보다는 설비 도입 여부 자체에 크게 좌우되는 비용이 존재할 수 있다.)

즉 $s_{\text{wt}}=0$이면 목적함수에 $c_{\text{misc}}$가 더해지면 안 되지만, $s_{\text{wt}}>0$이면 목적함수에 $c_{\text{misc}}$가 더해져야 한다.


이를 반영하기 위해서는, ‘이진수’ 변수 $b_{\text{wt}} \in \lbrace 0,1 \rbrace$와 ‘매우 큰 양수 계수 $M$’을 도입해 아래 제약을 추가한다.

$ s_{\text{wt}} \leq M b_{\text{wt}} $

만약 $b_{\text{wt}}=0$이면, $s_{\text{wt}} \leq 0$이 되므로 반드시 $s_{\text{wt}} = 0$이다. 반대로 $b_{\text{wt}}=1$이면, $s_{\text{wt}} \leq M$이고 $M$은 매우 큰 양수이므로 사실상 $s_{\text{wt}}$에 대한 추가 제약이 없는 셈이다.

그리고 기존의 목적함수에는 $c_{\text{misc}} b_{\text{wt}}$ 를 더한다. 그러면 풍력발전기가 한 대도 도입되지 않을 경우 목적함수 최소화를 위해 $b_{\text{wt}}=0$이 되어 $c_{\text{misc}}$가 비용에 반영되지 않고, 한 대라도 도입되면 반드시 $b_{\text{wt}}=1$이므로 $c_{\text{misc}}$가 비용에 반영된다.

마지막으로, 변수 $b_{\text{wt}}$를 ‘이진수’로 강제해야 하는 코드를 추가한다. 이는 $\texttt{glpk.ilp}$의 메서드에서, 위 풍력 대수 구하는 문제에서 언급한 $\texttt{ I=set([]) }$를 $\texttt{ B=set([]) }$로, 즉 $\texttt{I}$를 $\texttt{B}$로 바꾸면 된다 (B는 binary를 의미한다).


‘이진수’ 도입을 통한 모델 현실 설명력 증가: 용량당 단가가 다를 경우

이진수 변수 도입은, 용량에 따라 단위용량당 비용을 다르게 적용해야 하는 경우에도 쓸 수 있다.

필자가 이전에 수행했던 프로젝트 중에, ‘태양광 용량이 1000kW 이하면 kW당 단가를 $c_{\text{pv}.1}$로, 1000kW 이상이면 kW당 단가를 $c_{\text{pv}.2}$로 적용’할 것을 요구받은 적이 있다. 1000kW 이상의 대량 발주 시 단가를 상대적으로 낮출 수 있다는 이유였던 것으로 기억한다. 즉 $c_{\text{pv}.1}>c_{\text{pv}.2}$이다.

(이는 용량이 1500kW인 경우 1000kW에 대해서는 단가 $c_{\text{pv}.1}$을, 나머지 500kW에 대해서는 단가 $c_{\text{pv}.2}$를 적용하는 것이 아님에 주의. 1500kW 전부에 대해 단가 $c_{\text{pv}.2}$를 적용하는 것이다.)

이 경우, 우선 태양광 용량 변수도 ‘1000kW 이하인 경우에 대한 변수’ $s_{\text{pv}.1}$과, ‘1000kW 이상인 경우에 대한 변수’ $s_{\text{pv}.2}$ 두 개를 정의한다. 그리고 목적함수에는 $c_{\text{pv}.1} s_{\text{pv}.1} + c_{\text{pv}.2} s_{\text{pv}.2} $ 를 추가한다.

이제 두 용량 변수의 범위를 표현하는 제약조건들을 추가해야 한다. 아래 두 제약조건을 추가하면 될까?

$s_{\text{pv}.1} \leq 1000$, $s_{\text{pv}.2} \geq 1000$

그렇지 않다. 제약조건 $s_{\text{pv}.2} \geq 1000$ 때문에 $s_{\text{pv}.2}$가 반드시 1000kW 이상이 되어야 하기 때문이다. 실제로는 태양광을 아예 도입하지 않거나 1000kW 미만으로 도입하는 경우도 있을 텐데, 그런 경우를 표현할 수 없다.

이 경우에도 이진수 $b$를 도입해야 한다. 이진수 $b$가 0이면 1000kW 이하 도입에 대응하고 (도입 안 하는 경우도 포함), $b$가 1이면 1000kW 이상 도입에 대응한다고 하자. 그러면 ‘매우 큰 양수 계수’를 $M$이라 할 때, 아래 제약조건들을 추가하면 된다.

$s_{\text{pv}.1} \leq 1000 \times (1-b)$, $ 1000 b \leq s_{\text{pv}.2} \leq M b$

$b=0$이면 $s_{\text{pv}.1} \leq 1000$이고 $s_{\text{pv}.2} \leq 0$이 되어, $s_{\text{pv}.1}$만 양수가 되고 그 값은 1000 이하로 제한된다. 반대로 $b=1$이면 $s_{\text{pv}.1} \leq 0$이고 $1000 \leq s_{\text{pv}.2} $가 되어, $s_{\text{pv}.2}$만 양수가 되고 그 값은 1000 이상으로 제한된다.


(이런 식의 이진수 활용 사례는 몇 가지가 더 있다. 특히 추후에 ‘주택의 누진요금제’ 관련 포스팅 및 ‘화석연료를 소비하는 중대형 발전기’의 운전 조건 관련 포스팅 등에서 다시 등장할 예정이다.)

선형계획법 기반 분산에너지시스템 최적화

01) 최소비용 시스템과 시간별 자료의 중요성
02) 배터리의 충/방전 스케줄 결정: 수식
03) 배터리의 충/방전 스케줄 결정: Python 코드 및 결과
04) 태양광과 배터리의 '용량' 결정: 목적함수 ('현재가치' 비용) 및 수식
05) 태양광과 배터리의 '용량' 결정: 코드, 결과, 투자회수기간 및 절감량 계산
06) 정수 (integer) 변수 도입으로 현실 설명력 증대
07) 공동주택의 '누진제' 전기요금 (단일계약) 수식
08) 전기 부하와 냉/난방 부하를 동시에 고려 (섹터커플링)
09) '부분'부하 성능 관련 제약들
10) 출력 조정 관련 제약들


(선형계획법 시리즈의 지식들은, 필자가 2012년부터 (주)블루이코노미전략연구원 (대표: 오시덕 박사)과 함께 선형계획법 기반으로 신재생 및 열병합 에너지 시스템 경제성분석/ 최적 에너지시스템 도출 Tool 개발/ 에너지 정책 효과 분석 관련 프로젝트들을 여러 건 수행하며 축적한 ‘기본 지식’에 해당하는 부분임을 밝힌다.)