2.1 전력시스템 모델링

전력시스템의 전기기계적 과도 안정성 분석에는 외란에 대한 비선형 동적 응답 계산이 포함된다. 결과 공식은 다음과 같이 나타낼 수 있는 대형 DAE 시스템으로 구성된다.

여기에서 $x$와 $y$는 각각 상태 변수와 기타 대수, non-state 변수의 벡터이다. 상미분방정식(Ordinary differential equations – ODEs)에는 동기 발전기, 관련 여자 시스템, 터빈-거버너 시스템, 동적 부하 및 기타 동적 모델이 포함되고 상태 변수는 발전기의 회전자 각도, 회전자 속도, 필드 전압 등이 있다. 대수 방정식(algebraic equations - AEs)에는 네트워크의 정적 부하 및 power balance 방정식(네트워크 방정식)이 포함되고 대수 변수는 발전기의 d 축 전류, q 축 전류, 전기적 출력 등이 있다. 본 논문에서는 발전기의 제어 시스템을 포함한 더 세부적인 모델링을 하기 위하여 일반적으로 전력시스템 동적 시뮬레이션에 널리 사용되는 $q$ 축에 두 개의 제동 권선, $d$ 축에 계자권선과 함께 하나의 제동 권선이 있는 동기 발전기(IEEE Model 2.2 ^[22]), IEEE Type 1의 여자 시스템, 1차 터빈-거버너 모델을 채택했다. 자세한 모델링은 ^[22,23]에서 찾을 수 있다. 완성된 모델은 모든 제어 장치를 포함하여 각 발전기에 대해 15개의 상태 변수를 갖는다.

2.2 배경: Parareal 알고리즘

본 절에서는 Parareal 알고리즘을 간단히 설명한다. 알고리즘의 주요 아이디어는 전체 시뮬레이션 기간을 $\triangle t$의 시간 간격을 갖는 $N$개의 독립적인 하위 간격으로 분해하고 각각의 하위 간격을 $\delta t$의 시간 간격으로 더 세분화하여 분해된 여러 하위 간격에 대한 수치 솔루션을 동시에 계산하는 것이다. 문제 (1), (2)의 하위 간격 $N$개에 대한 수치 솔루션을 병렬로 풀기 위해 모든 하위 간격의 초기 상태가 필요하며 이는 $\triangle t$의 시간 간격으로 계산하는 coarse operator(상대적으로 덜 정확하지만 계산 비용이 저렴한 방법)으로 제공된다. 이후에 $\delta t$의 시간 간격으로 계산하는 fine operator(더 정확하지만 계산 비용이 큰 수치 적분 방법)를 수행하여 실제 궤적을 찾고 각 독립 하위 간격의 해를 병렬로 수정한다. Parareal 알고리즘은 모든 coarse 솔루션이 true 솔루션에 수렴할 때까지 coarse operator와 fine operator를 반복적으로 수행한다. 다음 형태의 초기값 문제를 예로 설명하고자 한다.

전체 기간[$t_{0},\: T$]을 N개의 하위 간격으로 $[t_{0},\: T]$$=$$[t_{0},\: t_{1}]\cup$ $[t_{1},\: t_{2}]\cup\cdots\cup[t_{N-1},\: t_{N}]$, $t_{n+1}- t_{n}=\triangle t$와 같이 나눈다. 본 논문에서 fine operator와 coarse operator를 통해 얻은 시스템 상태 변수들은 각각 $x^{fine}$과 $x^{coarse}$로 표기하고 수정된 coarse 솔루션은 $x^{*}$로 표기한다. 아래 단계들은 Parareal 알고리즘의 구현을 대략적으로 요약하여 보여준다.

1) 시간 간격 $\triangle t$와 $\delta t$를 각각 사용하는 coarse ($C_{\triangle t}$) 와 fine ($F_{\delta t}$) operator를 선택한다. 시간 $t_{n-1}$에서 초기 상태 $x_{n-1}$을 가지고 선택한 operator들을 사용하여 시간 $t_{n}$에서의 문제 (3)에 대한 해를 구한다.

2) 처음 coarse operator를 통해 구해진 해들은 순차적으로 구해진다.

위 첨자는 반복 횟수를 나타내고 $x_{0}^{*,\: 0}$은 시간 $t_{0}$의 초기값 을 나타낸다.

3) $k = 1$으로 반복이 시작된다. $N$개의 하위 간격에 대한 해를 fine operator를 사용하여 병렬로 동시에 구한다.

$x_{n}^{f\in e,\: k}$은 시간 $t_{n}$에서의 해를 나타낸다.

4) coarse 솔루션을 순차적으로 최신화한다.

for $n = k : N$

end

5) 단계 3으로 가서 $n=[1,\: \cdots ,\: N]$에 대해 $x_{n}^{*,\: k}-x_{n}^{*,\: k-1}\le tol$를 만족할 때까지 반복적으로 coarse 솔루션을 최신화한다.

알고리즘의 동작을 이해하기 위해 첫 번째 반복 ($k = 1$)이후 $n = 1$에서 최신화된 coarse 솔루션을 고려하고자 한다; 즉, $x_{1}^{*,\: 1}= x_{1}^{coarse,\: 1}+ x_{1}^{f\in e,\: 1}- x_{1}^{coarse,\: 0}$. $x_{1}^{coarse,\: 1}= x_{1}^{coarse,\: 0}$ 때문에 $n=1$에서 coarse 솔루션은 fine 솔루션으로 최신화된다 ($x^{*,\: 1}=x_{1}^{f\in e,\: 1}$). 그러므로 모든 coarse 변수들은 $N$ 번의 반복에서 fine 솔루션(true 솔루션)으로 최신화되고 이는 모든 간격 $N$ 에 대해 fine operator를 순차적으로 사용한 솔루션과 동등하다. 따라서 속도 향상은 Parareal 반복 횟수인 $k$가 $N$ 보다 작은 경우에만 일어난다. 이것은 coarse operator의 솔루션 시간과 병렬계산과 연관된 다른 요소(예: 통신 시간)들을 무시할 수 있다는 전제하에 Parareal 알고리즘의 이론적 속도 향상은 $\dfrac{N}{k}$임을 의미한다. 아래의 그림 1은 Parareal 알고리즘을 그래픽으로 표현한 것이다.

그림 1. Parareal 알고리즘의 그래픽 표현

Fig. 1. The graphical illustration for the coarse and fine operator

2.3 Parareal 알고리즘 구현

전력시스템 시뮬레이션에서 DAEs를 푸는 방법으로는 크게 두 가지 유형이 있다. 첫 번째 유형은 ODE의 연속적인 해와 AE를 따로 푸는 분할 명시적 방법 (partitioned-explict - PE)이다. 두 번째 방법은 시스템의 DAE를 동시에 푸는 방법이다. 명시적 적분을 사용한 PE 방법은 상용 소프트웨어 프로그램에서 주로 사용된다. 따라서, 본 논문은 PE 방법으로 Parareal 알고리즘을 구현하고자 한다. 기존의 연구 ^[11]에 기반하여 아래의 Midpoint-Trapezoidal 예측-수정(Trap) 방법이 DAEs를 풀기 위한 Parareal 알고리즘의 coarse operator로 선택되었다.

Trap 방법에서는 시뮬레이션의 근사해를 오직 한 번의 반복 횟수 ($j = 1$)로 구하였다. Fine operator는 전력시스템 동적 시뮬레이션에서 널리 쓰이는 아래의 Runge-Kutta 4th order (RK4) 방법이 사용되었다.

ODE와 AE의 번갈아가는 해에 따른 오차를 완화시키기 위해 네트워크 방정식은 fine과 coarse operator 모두에서 풀어진다. 본 논문에서는 상기의 일반적인 수치해석 기법 기반의coarse operator를 사용하였고 근해석적 기법 기반의 다양한coarse operator를 고려하여 Parareal 알고리즘의 성능을 분석한 내용은 ^[11]에서 확인할 수 있다.

3.1 확률미분대수방정식

Stochastic differential equation(SDE)의 일반적인 구조는 다음과 같이 나타낼 수 있다.

여기에서 $\eta$는 무작위 변수를 나타낸다;$\alpha$는 SDE의 드리프트 계수($\eta$이 $dt$의 시간 동안 어떻게 움직이는지)를 나타낸다;$\beta$는 SDE의 확산 계수를 나타낸다; 그리고 $B$는 브라운 운동(무작위성의 원천으로써 표준 Wiener 과정과 같다)을 나타낸다.

Wiener 과정은 여러 형식의 공식적인 정의가 있다. 가장 실용적인 정의는 다음과 같은 특징을 갖는 실수 값의 연속적인 Stochastic process w(t)이다: w(0)은 0이고, w(t)는 모든 구간에서 무한한 변동을 갖는다.

시간 영역 시뮬레이션에서 브라운 운동은 다음과 같은 특성을 갖는다: $B(0)= 0$이고, 모든 시간 $0\le t_{1}\le t_{2}\le\cdots\le t_{N}$ 에서 독립적인 증분을 갖는다. 즉, $B(t_{2})-B(t_{1}),\: \cdots ,\: B(t_{N})-B(t_{N-1})$은 모두 독립적인 무작위 변수이다. 또한, 증분 $B(t+h)- B(t)$는 평균이 0이고 분산이 $h$인 정규분포를 따른다. 즉, 모든 $t,\: h > 0$에 대해$B(t+h)-B(t)\sim N(0,\: h)$이다. 특히, 브라운 운동의 미분형은 white noise로 정의된다.

^[24]의 Ito의 해석과 Stratonovich의 해석을 사용하여 (11)은 적분형으로 표현할 수 있고 이와 (12)을 사용하여 (11)은 다음과 같이 일반적인 미분방정식으로 나타낼 수 있다.

시간 $t$는 표기를 최소화하기 위해 생략되었다.

본 논문에서는 ^[21]에 근거해 SDAEs로 구성되어있는 확률적 전력시스템 동역학적 모델에 대해 논의한다. 이전에 논의한 바와 같이, 네트워크 이상으로 인한 비선형 동적 응답을 얻기 위해 전력시스템의 전기기계적 거동과 안정성에 대한 해석은 수치적인 시간 영역 시뮬레이션이 필요하다. 해결해야 하는 해당 수식은 DAE의 큰 집합으로 보편적으로 (1)과 (2)로 표현된다. (1)과 (2)는 널리 받아들여졌으며, 관련된 시간 영역 시뮬레이션 도구들은 전력시스템 네트워크의 DSA에서 흔히 사용되고 있다. 그러나 대규모 DER의 통합은 전력시스템 네트워크에서 동적 거동의 복잡성과 불확실성을 증가시킨다. 따라서 이러한 결정론적 DAE들은 특별한 주의를 요구한다. 특히 (1), (2)에서는 시스템 내의 불확실성과 확률적 변동이 전력시스템의 과도 현상에 미치는 영향을 무시한다. 더욱 현실적인 DSA를 수행하기 위해서는 부하와 신재생에너지의 랜덤 변화 및 제어 장치의 측정 오차 등으로 인한 확률적 변동을 고려하는 것이 필요하다.

계통에 존재하는 확률적 변동의 영향은 회전자의 위상 각, 전압 크기 등과 같은 주요 상태 및 대수 변수에 확률적인 거동을 일으킬 수 있다. 전력시스템 변수와 파라미터는 상호작용하여 이와 같은 확률적 변동에 영향을 미칠 수 있다. 이러한 상호작용과 (13)을 토대로 하여 전력시스템 변수와 파라미터를 포함한 확률적 변동은 다음과 같이 모델링될 수 있다.

(14)를 사용하여 (1), (2)의 DAEs는 아래의 일반적인 SDAEs로 변형할 수 있다.

여기에서 확률적 변동 $\eta$의 영향을 포함하도록 기존의 $f$와 $g$는 수정되었다. 초기값이 $\eta(t_{0})=\eta_{0}$으로 주어진다고 가정한다. 즉, 15(c)에서 $\dot{\eta}= 0$ 그리고 $\eta_{0}= 0$을 가져 확률적 변동의 영향을 무시한다면 (15)은 일반적인 결정론적 초기값 문제의 형태를 갖는 것을 확인할 수 있다.

3.2 Ornstein-Uhlenbeck Process

시간 영역 시뮬레이션을 위해 본 논문은 일반적인 stochastic process인 (11)을 사용하는 대신 mean-reverting process로도 불리는 Ornstein -Uhlenbeck’s process를 사용한다. 일부 연구(예: ^[25])에서는 전력시스템의 확률적 부하를 모델링하기 위해 이 process를 사용했으며, Ornstein-Uhelnbeck’s process를 설명하는 일반적인 SDE공식은 다음과 같이 정의된다.

여기에서 드리프트 계수 $\alpha$와 확산 계수 $\beta$는 상수이다; $\mu$는 미리 지정된 평균값이다. 이 process에서 $\eta(t)$는 평균값 $\mu$로 이동하려 하고 시간이 무한히 지남에 따라 정규분포 $N(\mu ,\: \dfrac{\beta^{2}}{\alpha})$를 나타낸다. $\alpha$는 평균회귀 속도를 결정하고 $\beta$는 편차를 조절한다. 다른 함수로 드리프트 계수와 확산 계수를 사용할 수 있지만 Ornstein-Uhlenbeck’s process는 이 함수들을 크게 간소화하고, 분산이 경계가 있으므로^[21] 물리적 시스템에서 확률적 변동을 모델링할 때 적절한 선택으로 여겨진다. 결론적으로, 본 논문은 기존의 확률미분대수방정식 (15)c을 (16)으로 대체하였다.

3.3 확률론적 부하 모델링

전력 시스템의 시간 영역 시뮬레이션에서 부하의 불확실성은 무시되어왔지만 그들의 확률론적인 행동은 분산 에너지 자원의 성장으로 더욱 증가할 것으로 예상된다. 이를 위하여 다음과 같은 ZIP 모델을 고려한다:

$P_{L},\: Q_{L},\: $과 $v_{m}$은 각각 유효 전력 부하, 무효 전력 부하 그리고 전압의 크기를 나타낸다; $P_{L0},\: Q_{L0}$과 $v_{m0}$은 정상상태에서의 유효 전력 부하, 무효 전력 부하와 전압의 크기를 나타낸다; $a_{1},\: a_{2}$와 $a_{3}$는 각각 유효 전력 부하에서 일정한 전력, 전압 및 임피던스 성분의 비율을 나타내고 이들의 합은 1이 되어야 한다; 비슷하게 $b_{1},\: b_{2}$와 $b_{3}$는 각각 무효 전력 부하에서의 비율을 나타내고 이들의 합은 1이 되어야 한다.

이러한 결정론적 궤적에 대한 확률론적 부하 모델은 식 (14), (15)를 기반으로 다음과 같이 확률적인 변동을 도입함으로써 모델링된다.

다음으로는 식 (16)에 의하여, 일반적으로 표현된 확률론적 부하 모델(18)은 다음과 같은 확률미분대수방정식의 집합으로 만들 수 있다:

위의 확률적 모델은 원래의 미분대수방정식에서 추가적인 미분 방정식을 간단히 포함한다. 본 논문은 확률적 부하 모델(즉, 시스템 매개변수)에 중점을 두고 있지만 전력 시스템의 상태 및 대수 변수의 궤적에 대한 확률적 모델도 어려움 없이 유도될 수 있다^[21].

4.1 Parareal 알고리즘 검증

첫째로, Parareal 알고리즘의 수렴성을 검증하였다. 개인용 컴퓨터의 허용 가능한 프로세서를 고려하여, fine solver의 병렬계산을 위한 프로세서 수는 50개로 정하였다. 시뮬레이션 기간 10초는 50개의 하위 간격으로 나누어지고 각 하위 간격에서의 fine operator는 다시 100개의 하위 간격으로 나누어져 계산하였다.(즉, fine operator의 시간 간격은 ($\delta t$$=\dfrac{10}{50\times 100}$$= 0.002s$). 분자의 10은 전체 시뮬레이션 시간이고 분모의 50은 coarse operator에 사용되는 하위 간격의 개수이다. 이 시뮬레이션에서 각 프로세서의 coarse operator는 10개의 하위 간격을 사용하고, true 솔루션은 시간 간격 0.002s로 계산하는 RK4 방법을 사용하여 구했다.

그림 2는 Parareal 알고리즘을 적용하여 구한 New England 시스템에 4주기의 1 모선에서 3상 단락 시뮬레이션 결과 중 발전기 1의 위상 각과 슬립 속도의 궤적을 나타낸다. 그림 2(a)는 변수 중 발전기의 위상각을 관찰하였으며 파란색 선은 true 솔루션, 초록색 선은 parareal 알고리즘을 사용하여 나온 결과이다. 각 선의 궤적이 일치하는 것을 통해 parareal 알고리즘이 true 솔루션에 수렴하는 것을 확인할 수 있었고 이와 같은 수렴을 그림 2(b)에서는 다른 변수인 슬립속도를 통해 확인하였다. 다른 외란에 대해서 해당 내용을 확인함으로써 Parareal 알고리즘의 수렴성을 검증하였다.

또한, Parareal 반복 횟수에 대한 coarse operator의 시간 간격의 영향을 분석하기 위해 fine operator의 시간 간격은 0.002s로 고정하고 coarse operator의 시간 간격을 0.05s부터 0.005s까지 감소시켰다. 마찬가지로, fine operator의 시간 간격의 영향을 분석하기 위해 coarse operator의 시간 간격은 0.008s로 고정하고 fine operator의 시간 간격을 0.002s부터 0.0001s까지 감소시켰으며 해당 시뮬레이션 결과는 표 1에 정리하였다.

그림 2. (a) 4주기의 1모선에서 3상 단락 시뮬레이션의 발전기 1의 위상각에 대한 Parareal 시뮬레이션 결과

Fig. 2. (a) Parareal simulation results for rotor angle of generator 1 followed by the 3-phase fault at bus 1 with 4 cycles

표 1에 나온 결과를 통해, coarse operator의 시간 간격이 감소하면 Parareal 반복 횟수가 감소하는 경향과 일정 시간 간격 후에는 더 이상 감소하지 않는 것을 확인할 수 있다. Fine operator의 경우, 모두 동일한 Parareal 반복 횟수를 보여준다. 이유는 fine operator가 0.002s에서 이미 높은 정확도로 true 솔루션을 제공하였기 때문이고, 시간 간격을 더 줄이면 솔루션의 정확도가 약간 향상되는 것을 의미한다. 반대로, fine operator의 시간 간격을 증가하는 시뮬레이션도 수행하였고 Parareal 알고리즘의 반복 횟수가 약간 증가하는 것을 확인할 수 있었다.

4.2 Parareal 알고리즘 속도 향상

Parareal 알고리즘의 실용적 가치와 속도 향상을 검증하기 위해 Compute and Data Environment for Science (CADES) HPC 컴퓨팅 플랫폼 ^[26]을 활용하여 시뮬레이션을 수행하고 프로세서의 수를 변경하며 계산 시간을 기록하였다. 초대규모 네트워크인 Eastern Interconnection (EI) 5,617-generator 70,285–bus 시스템을 사용하여 4주기의 모선 1에서 3상 단락을 시뮬레이션했다. 또한, 이 테스트에서는 EI시스템의 DAEs를 단순화하지 않고 이전 절에서 사용된 것과 동일한 동적 모델을 사용하였다. 아래의 표 2는 다양한 프로세서에서 HPC 플랫폼 을 사용한 Parareal 알고리즘의 계산 시간 (초)를 나타낸다.

그림 2. (b) 4주기의 1모선에서 3상 단락 시뮬레이션의 발전기 1의 슬립속도에 대한 Parareal 시뮬레이션 결과

Fig. 2. (b) Parareal simulation results for slip speed of generator 1 followed by the 3-phase fault at bus 1 with 4 cycles

표 2를 보면, 프로세서 수가 증가함에 따라 계산성능이 크게 향상되는 것을 확인할 수 있었다. 500개의 프로세서를 사용한 Parareal 알고리즘은 10개 프로세서를 사용한 것보다 약 6배 빠르게 완료된다. 이는 프로세서로의 수를 확장할 때, 알고리즘이 각 프로세서에서 더 많은(더 작은) 하위 간격을 갖게 하므로 coarse 및 fine operator에서 더 많은 병렬계산을 허용하여 각 operator의 계산 시간을 크게 줄이기 때문이다. 더 많은 프로세서로 시뮬레이션하는 것은 통신 시간을 증가시키지만, 더 많은 병렬계산으로 인한 감소가 더 크기 때문에 총 시뮬레이션 시간을 줄어드는 것을 확인할 수 있다. 이러한 결과는 전력시스템 동적 거동을 예측하기 위한 신흥 솔루션으로써 Parareal 알고리즘의 큰 잠재력을 보여준다.

4.3 확률론적 전력시스템 시뮬레이션

일반적으로, SDAEs는 해석적 솔루션을 거부하기 때문에 수치적인 방법을 활용해야 한다. 본 절에서는 가장 일반적으로 사용되는 Euler와 Milstein schemes ^[27]과 같은 분할 명시적 방법을 사용하였다. 확률적 부하 모델(즉, 시스템 매개변수)을 사용하였다(전력시스템의 상태, 대수 변수를 포함한 확률적 모델의 자세한 유도방식은 ^[21]에서 확인할 수 있다). 확률적 ZIP 부하 모델을 사용하여 모든 부하에 대해 확률적 변동을 포함하였고 확률적 부하 모델에 대한 Wiener process의 표본 경로는 fine operator와 같은 시간 간격(즉, $\delta t =\triangle t_{B}= 0.002s$)을 갖는다. 전력시스템의 불확실성을 반영하기 위해 각 고장 시나리오에 대해 100개의 다른 Wiener process을 시뮬레이션하였다.

그림 3은 모선 7의 전압 크기에 대한 확률적 거동을 보여준다. 그림에서처럼 고장이 해결된 후 전압 크기의 결정론적 값(즉, 평균값)은 점진적으로 nominal 값으로 회복된다. 전압 크기의 평균값은 고장 후 약 1초 정도의 시간이 지나면 최소 전압 한계값 (값 0.9의 회색 점선)보다 커진다. 그러나 SDAE 시스템의 많은 솔루션은 1초 이상 동안 한계값보다 작으며, 결정론적 모델이 포착할 수 없는 전압 회복이 지연되는 정보(예: 전압 크기가 한계값보다 커지기 위해서는 5초 이상의 시간이 필요)를 보여준다. 이러한 정보는 DSA 연구에 도움이 되며, 보호 계전기의 운영에 대한 더 정확한 지침을 제공할 수 있을 것이다.

그림 3. 확률론적 부하로 인한 모선 7의 전압 크기

Fig. 3. Bus voltage magnitude at bus 7 with stochastic loads