1. 서론
풍력터빈을 설계하거나 설치할 때 바람의 정보는 필수적이다. 바람의 등급을 4가지로 나누어서 하중해석을 통해 구조물을 설계하며, 풍력단지 개발 시에는 최소 1년의 바람 계측 자료를 필요로 한다. 계측지역의 바람에 대한 통계는 평균풍속, 난류도, 바람장미, Weibull 분포 등을 예로 들 수 있다 [1, 2]. 풍력터빈에 대한 디지털 트윈 기술에서는 풍력터빈 혹은 단지를 모델링 한 후 가상의 시계열 바람자료를 통해 운전 상태를 파악할 수 있다. 이와 같이 계측을 통해서든 아니면 수치적인 절차를 통해서든 풍력터빈 시스템을 연구하기 위해서는 바람자료를 요구한다. 그러나 계측데이터는 비용과 지속 시간 측면에서 한계가 있어서 무한 생성이 가능한 고품질 인공 바람자료(synthetic wind data)가 풍력터빈 로터의 기계적 출력이 그리드의 전기적 출력으로 연결되는 통합해석 과정에서 필수적인 역할을 하게 되었다. 따라서 실측 데이터를 대신하여 통계적으로 바람을 잘 모사하는 대리 시계열(surrogated time series) 생성을 위한 풍속 모델링의 필요성이 대두되었다 [3].
풍력터빈의 작동 상태를 실시간으로 파악하기 위해서는 시계열 바람자료를 필요로 한다. 이 시계열 바람자료가 어떤 지역의 통계적인 특성을 반영할 수 있다면 더 훌륭한 바람자료라고 볼 수 있다. 이러한 바람자료를 생성하는데 마르코프 사슬 (Markov chain) 방법이 적용된다 [4-6]. 이 방법은 선형회귀법이나 데이터 간의 표준편차의 비를 이용하여 시계열 바람자료를 생성하는 MCP (measure-correlate-predict) 방법과는 달리 바람의 통계자료를 기반으로 하기 때문에 예측 정확도를 더 향상시키는 것으로 알려져 있다 [7]. 이러한 마르코프 사슬 방법은 풍속 뿐 만 아니라 태양복사 등의 자연적인 과정을 모델링하는 데에도 적용되어왔다 [8].
마르코프 사슬에서는 현재의 정보가 미래의 상태를 결정하는 데 영향을 미친다. 한 상태에서 다음 상태로 넘어가는 과정이 사슬과 같이 연결된다. 또한 마르코프 사슬에는 차수(order)가 존재한다. 1차 마르코프 사슬은 계측자료의 정보로부터 구한 2차원 배열의 누적상태 전이확률과 2개의 난수 발생을 통해서 현재의 정보가 미래의 자료를 생성하는데 관여하도록 한다. 2차 마르코프 사슬은 계측자료의 정보로부터 구한 3차원 배열의 누적상태 전이확률과 2개의 난수 발생을 통해서 과거와 현재의 정보가 미래의 자료를 생성하는데 관여한다. 많은 연구에서 1차 마르코프 사슬이 시계열 자료를 생성하는데 충분하다고 추천한다. 반면에 2차 마르코프 사슬이 정확도를 높이는데 기여한다는 주장도 존재한다. 본 연구에서는 시계열 바람자료를 생성하는 방법으로 1차 마르코프 사술과 2차 마르코프 사슬을 적용하여 두 방법의 차이점이나 개선점이 무엇인지를 비교 분석하고자 한다.
2. 마르코프 사슬
2.1 1차 마르코프 사슬
참조지역의 관측자료가 없이 처음부터 마르코프 사슬을 통해 인위적인 바람 자료의 생성이 불가능한 것은 아니지만, 마르코프 사슬을 통해 시계열 바람자료를 손쉽게 얻기 위해서는 특정 지역의 시계열 풍속자료를 참고한다. 풍속의 최소값과 최대값이 정해지면 n개의 구간으로 풍속을 나누어서 각 구간을 상태(state)라고 정의한다. 먼저 1차 마르코프 사슬에 대해서 알아보기로 한다. 1차 마르코프 사슬은 시간간격이 Δt일때 시간 t + Δt에서의 풍속은 시간 t의 풍속에 의존한다. 상태 i에서 상태 j로 전이되는 확률을 구해야 하는데 이를 상태전이확률(state transition probability) pij로 정의한다. 다시 설명하면 시계열 풍속데이터를 순차적으로 모든 시간 t에서의 풍속이 어느 상태에 있는지를 정하는 1차원 배열을 구한 다음, 그 배열이 상태 i에서 상태 j로 천이된 경우의 수를 처음부터 모두 구해나간다. 물론 i와 j는 1에서 n까지의 자연수이다. 이를 통해 구한 n × n의 2차원 배열이 상태전이확률 행렬이다.
상태전이확률 행렬(state transition probability matrix) P를 수식으로 표현하면 다음과 같다.
\(\begin{align}P=\left[\begin{array}{cccc}p_{11} & p_{12} & \ldots & p_{1 n} \\ p_{21} & p_{22} & \ldots & p_{2 n} \\ \mid & \mid & \mid & \mid \\ p_{n 1} p_{n 2} & \cdots & p_{n n}\end{array}\right]\end{align}\) (1)
식 (1)에서 각 행의 요소 n개에 대한 합은 1로 정규화 한다. 따라서 pij는 0과 1 사이의 값이다. 인공 바람 자료를 생성하기 위해서는 식 (1)로부터 각 행에대해 누적 상태전이확률 행렬 Pc (cumulative state transition probability matrix)를 구성한다. 이 행렬의요소 pcik 는 다음과 같다.
\(\begin{align}p_{i k}^{c}=\sum_{j=1}^{k} p_{i j}\end{align}\) (2)
어느 지역의 바람자료에 대한 누적 상태전이확률을 알게 되면 다음과 같은 절차를 통해 비로소 마르코프사슬을 통한 시계열 바람자료의 생성이 가능해진다.
- 시간 t에서 상태 i를 구한 후 0과 1 사이의 난수 r1을 발생시켜 식 (2)의 pci,k-1 ≤ r1 < pci,k인 k를 찾아낸다.
- 또 다른 독립된 0과 1 사이의 난수 r2의 발생을 통해 다음 시간에서의 풍속 V = (k - 1 + r2)ΔV를 구한다.
- 위 순서를 순차적으로 반복하여 시계열 풍속자료를 생성한다.
2.2 2차 마르코프 사슬
2차 마르코프 사슬에서 X를 확률적인 과정이라 할때, 상태전이확률 행렬 P의 요소 pijk는 다음과 같다.
pijk = p{X(t) = si |X(t - 1) = sj∧X(t - 2) = sk} (3)
여기서 i, j, k∈{1, 2, …, n}이며, 현재의 상태가 j이고 그 이전의 상태가 k일 때 다음의 상태 i에서의 확률을 표현한다. 다시 설명하면 시계열 풍속데이터를 순차적으로 모든 시간 t에서의 풍속이 어느 상태에 있는지를 정하는 1차원 배열을 구한 다음, 그 배열이 상태 i에서 상태 j로 변한 다음, 다시 상태 k로 전이된 경우의 수를 처음부터 모두 구해나간다. 물론 i, j, k는 1에서 n까지의 자연수이다. 이를 통해 구한 n × n × n의 3차원 배열이 2차 마르코프 사슬을 위한 상태전이확률 행렬이다. 따라서 pijk는 다음과 같은 상태전이확률 행렬 P의 한 요소가 된다. 편의상 3차원 배열을 2차원 행렬형식으로 재배열 한 모양을 갖는다.
\(\begin{align}P=\left[\begin{array}{cccc}p_{111} & p_{112} & \ldots & p_{11 n} \\ p_{121} & p_{122} & \ldots & p_{12 n} \\ \mid & \mid & \mid & \mid \\ p_{1 n 1} & p_{1 n 2} & \ldots & p_{1 n n} \\ p_{211} & p_{212} & \ldots & p_{21 n} \\ p_{221} & p_{222} & \ldots . & p_{22 n} \\ \mid & \mid & \mid & \mid \\ p_{n n 1} & p_{n n 2} & \ldots . & p_{n n n}\end{array}\right]\end{align}\) (4)
식 (4)에서 각 행의 요소(상태) n개에 대한 합은 1로 정규화 한다. 따라서 pijk는 0과 1 사이의 값이다. 인공 바람 자료를 생성하기 위해서는 식 (4)로부터 각 행에 대해 누적 상태전이 확률 행렬 Pc를 구성한다.이 행렬의 요소 pcijk는 다음과 같다 [9].
\(\begin{align}p_{i j k}^{c}=\sum_{l=1}^{k} p_{i j l}\end{align}\) (5)
2차 마르코프 사슬은 시간 t - Δt와 시간 t에서의 풍속이 정해진 상태에서 시간 t - Δt에서의 풍속을 구한다는 것 외에 절차는 1차 마르코프 사슬의 절차와 동일하다.
3. 결과 및 토의
본 연구를 위해 사용된 바람 계측 자료는 기상청의 AWS로 관측된 전라북도 부안군 위도면의 2020년 1월 1일부터 1년 동안의 풍황자료이다. Fig. 1은 관측 위치의 Weibull 변수로 측정된 풍속을 Weibull 확률분포함수의 2가지 변수라 할 수 있는 척도계수와 형상계수는 각각 C = 2.71, k = 1.40이다 [10].
Fig. 1 Weibull distribution at Wido in 2020 [10]
풍속은 관측 자료를 참고하여 최고 14 m/s로 선택하였으며 구간은 0.5 m/s 등간격으로 총 28구간으로 정하였다. 1년 동안 위도에서 관측된 풍속 자료를 바탕으로 전이확률 행렬과 누적 전이확률 행렬을 구하면 Table 1 및 Table 2와 같다. 식 (4)와 식 (5)로부터 각표의 열 방향으로는 28 × 28 = 784개의 요소가 존재하며, 행 방향으로는 28개의 요소가 존재한다.
Table 1 A 28-state transition probability matrix for annual hourly mean wind speed at Wido island
Table 2 A 28-state cumulative transition probability matrix for annual hourly mean wind speed at Wido island
누적 전이확률이 구해지면 난수 발생기를 통해 0에서 1 사이의 난수를 발생시켜 현재 상태에서 다음 상태로의 풍속을 순차적으로 구해나가면서 시계열 자료를 생성한다. 시계열 자료는 총 5,000〜50,000개를 생성시킨 후 1년 동안의 계측자료(길이는 8,784)와 비교하도록 하였다. Fig. 2(a)는 위도에서 AWS로 계측된 1년 동안의 풍속 자료 중 1,000 시간까지의 시계열 풍속 자료를 보여준다. Fig. 2(b)는 1차 마르코프 사슬과 2차 마르코프 사슬을 이용한 인공 시계열 바람자료를 보여준다. 시계열 풍속자료를 분석하는 것은 쉽지 않기 때문에 확률밀도 함수로 나타낸 것이 Fig. 3이다. Fig. 3에 1차 및 2차 마르코프 사슬 결과와 계측자료의 확률밀도 함수를 상호 비교하였다.
Fig. 2 Time series wind speed
Fig. 3 Probability distribution for measured and synthetic wind speed (N = 5,000, 계측자료의 길이는 8,784)
Fig. 4는 시간에 따른 자기 상관도를 계측 자료 및 1차와 2차 마르코프 사슬 결과를 상호 비교한 것이다. 이틀의 시간동안에 대해서는 2차 마르코프 사슬을 이용한 방법이 확실한 개선 효과를 보여 주고 있음을 확인할 수 있다. 다만 시간이 지날수록 자기 상관도는 0에 수렴해야 할 것으로 생각되는데 이러한 현상이 보이지 않고 진동하는 현상을 보이고 있다. 이러한 현상은 데이터 샘플링의 길이가 짧아서 생기는 문제로 간주되어 그 양을 2 배에서 10 배까지로 증가시키면서 확률통계치의 결과를 관찰하기로 하였다.
Fig. 4 Auto-correlation for measured and synthetic wind speed (N = 5,000, 계측자료의 길이는 8,784)
본 연구에서 시계열 자료의 시간에 따른 변동의 양상을 나타내는 변수로 자기상관도 (autocorrelation)를 고려한다. 두 속도의 시간간격이 kΔt일 때, 두 시계열값에 대한 자기상관도 ck는 다음과 같다 [10].
\(\begin{align}c_{k}=\frac{\frac{1}{N-k} \sum_{i=1}^{N-k}\left(v_{i}-\bar{v}\right)\left(v_{i+k}-\bar{v}\right)}{\frac{1}{N} \sum_{i=1}^{N}\left(v_{i}-\bar{v}\right)\left(v_{i}-\bar{v}\right)}\end{align}\) (6)
여기서 i = 1, 2, 3, …, N이고, k = 0, 1, 2, …, N - 1이며, \(\begin{align}\overline {v}\end{align}\) 는 평균 풍속이다. 이론적으로 풍속간의 시간차가 없으면 자기상관도는 1이 되지만 풍속 간의 시간 간격이 커질수록 자기 상관도는 0에 근접하는 것이 특성이다.
마르코프 사슬을 이용한 인공바람자료 생성은 무한한 길이의 시계열 바람자료의 생성이 가능하다. 데이터 샘플링의 개수 N을 5,000, 8,784, 20,000, 30,000, 40,000, 50,000로 증가시켜가면서 인공바람자료의 통계분석을 하였으며 N = 40,000 이상에서 자기상관도가 시간 증가에 따라서 안정적으로 0에 수렴하게 됨을 확인하였다. 이 경우에 대해서 2차 마르코프 사슬의 돗수분포 결과가 1차 마르코프 사슬의 돗수분포 결과에 비해 계측치와 더 잘 일치하는 결과를 보여줌이 확인되었다. 시계열 자료의 길이 N을 50,000개로 하였을 때 풍속의 돗수분포 즉 확률밀도 함수를 Fig. 5에 나타내었다. Fig. 5에서 보여주듯이 풍속의 돗수분포는 1년 동안의 계측치와 1차 및 2차 마르코프를 이용한 결과가 서로 높은 정확도를 가지고 일치하고 있다.
Fig. 5 Probability distribution for measured and synthetic wind speed (N = 50,000, 계측자료의 길이는 8,784)
Fig. 6은 계측치와 2차 마르코프의 인공 바람자료에 대한 자기상관도를 보여주고 있다. 대략적으로 이틀동안의 바람자료는 매우 높은 정확도를 가지고 실험치와 서로 잘 일치하며, 그 이후는 자기 상관도가 0의 값으로 수렴하는 것을 잘 보여준다. Fig. 4와 마찬가지로 1차 마르코프 사슬은 2차 마르코프 사슬에 비해 자기 상관도가 떨어지는 결과를 보여주고 있다. Fig. 6으로부터 계측치가 시간이 지나도 계속 진동하는 현상을 보이는 이유는 데이터 샘플링의 길이 N이 영향을 주는 것으로 판단할 수 있다.
Fig. 6 Auto-correlation for measured and synthetic wind speed (N = 50,000, 계측자료의 길이는 8,784)
Table 3은 2차 마르코프 사슬을 이용하였을 때 풍속 자료의 최소값/최대값 및 평균과 표준 편차를 보여주고 있다. Table 3으로부터 통계적인 값들을 볼 때 실측치와 2차 마르코프 사슬을 이용한 인공 바람자료간의 주요 통계 변수들의 상대오차가 5% 범위 이내에서 모두 잘 일치하고 있음을 알 수 있다. Table 3에서 보여주지는 않지만 1차 마르코프 사슬을 이용한 경우에도 유사한 정확도를 나타낸다 [10].
Table 3. Comparison of statistical parameters between measured and synthetic data
4. 결론
본 연구에서는 1차 및 2차 마르코프 사슬을 이용하여 인공 바람자료를 생성한 후 통계적인 주요 변수들을 상호 비교하였다. 여기서 통계적인 주요 변수는 확률밀도함수와 자기상관도이며, 추가적으로 평균 및 표준편차를 실측치와 비교하였다. 풍속의 평균 및 표준편차, 그리고 확률밀도함수 분포는 2차 마르코프 사슬이 1차 마르코프 사슬에 비해 나아지기는 하지만 뚜렷한 개선을 보인다고 보기는 어렵다. 자기 상관도는 확실히 2차 마르코프 사슬이 1차 마르코프 사슬에 비해서 초기 대략 이틀 분량의 길이에 대한 일치도가 매우 높은 결과를 보여주었다. 데이터 샘플링의 길이가 자기 상관도에 영향을 미치는 것이 확인되었으며, 데이터의 길이를 조절하여 시간이 지날수록 자기 상관도가 0에 수렴하도록 하는 데이터 샘플링의 길이를 찾아낸 후 평균 및 표준편차에 대한 통계분석을 시도하는 것이 의미가 있는 것으로 확인된다.
본 연구에서는 실측 바람자료를 참조하여 마르코프 사슬을 이용한 인공바람자료를 생성한 후 그 통계적인 특성이 실측자료의 통계 특성을 잘 모사하고 있음을 확인하였다. 인공바람 자료는 계측자료가 갖고 있는 시간의 한계를 벗어나 무한한 길이의 바람자료 생성이 가능하다. 마르코프 사슬을 이용한 시계열 인공바람자료의 생성 기법은 가상의 풍력터빈 시스템의 성능을 평가하는데 매우 유용할 것으로 사료된다.
후기
본 연구는 한국전력공사의 2021년 착수 연구개발 과제 연구비에 의해 지원되었음. (과제번호: R21XO01-6)
References
- Oh, K. Y, Kim, J. Y, Lee, J. K., Ryu, M. S., and Lee, J. S., 2011, "Wind resources assessment and performance evaluation of wind turbines for the test-bed of offshore wind farm," J. Wind Energy, Vol. 2, No.. 1, pp. 15~20 (in Korean).
- Km, S. K, and Km, B. S., 2022, 'Load comparison evaluation and meteorological mast wind data analysis for wind turbine site-specific load analysis," J. Wind Energy, Vol. 13, No.. 3, pp. 43~52 (in Korean).
- DzAmbrosio, D., Schoukens, J., De Troyer, T., Zivanovic, M., Runacres, M.C., 2022, "Data-driven generation of synthetic wind speeds: A comparative study," IET Renewable Power Generation, Vol. 16, pp. 922~932
- Balouktsis, A., Tsanakas, D., and Vachtsevanos, G., 1986, "Stochastic simulation of hourly and daily average wind speed sequences," Wind Engineering, Vol. 10, No. 1, pp. 1~11.
- Lei, M., Shiyan, L., Chuanwen, J., Hongling, L., and Yan, Z., 2009, "A review on the forecasting of wind speed and generated power," Renewable and Sustainable Energy, Vol. 13, No. 4, pp. 915~920.
- Pesch, T., Schroders, S., Allelein, H J., and Hake, J. F., 2015, "A new Markov-chain-related statistical approach for modelling synthetic wind power time series," New Journal of Physics, Vol. 17, pp. 1~15.
- Kang, M., Son, E., Lee, J., and Kang, S., 2022, "Accurate wind speed prediction using effective Markov transition matrix and comparison with other MCP models," New and Renewable Energy, Vol. 18, No. 1, pp. 17~28 (in Korean).
- Poggi, P., Notton, G., Muselli, M., and Louche, A., 2000, "Stochastic study of hourly total solar radiation in Corsica using a Markov model," Int. J. of Climatology, Vol. 20, pp. 1843~1860.
- Shamshad, A., Bawadi, M.A., Wan Hussin, W. M. A., Majid, T. A., and Sanusi, S. A. M., 2005, "First and second order Markov chain models for synthetic generation of wind speed time series," Energy, Vol. 30, pp. 693~708.
- Ryu, K. W., 2021, "Use of Markov chains for synthetic wind data generation and its statistical verification," J. Wind Energy, Vol. 12, No.. 3, pp. 13~18 (in Korean).