A cyclic time-dependent Markov process to model daily patterns in wind turbine power production

A cyclic time-dep enden t Mark o v pro cess to mo del daily patterns in wind turbine p o w er pro duction T eresa Sc holz a , Vitor V. Lop es a, ∗ , Ana Estanqueiro a a LNEG, National L ab or atory for Ener gy and Ge olo gy, Estr ada do Pa¸ co do Lumiar, 22, 1649-038, Lisb o a, Portugal Abstract Wind energy is b ecoming a top con tributor to the renew able energy mix, whic h raises p oten tial reliability issues for the grid due to the fluctuating nature of its source. T o ac hieve adequate reserv e commitment and to pro- mote mark et participation, it is necessary to pro vide mo dels that can capture daily patterns in wind p o w er pro duction. This pap er presents a cyclic in- homogeneous Mark ov pro cess, which is based on a three-dimensional state- space (wind p o wer, speed and direction). Eac h time-dep enden t transition probabilit y is expressed as a Bernstein p olynomial. The mo del parameters are estimated by solving a constrained optimization problem: The ob jectiv e function combines tw o maximum likelihoo d estimators, one to ensure that the Mark ov pro cess long-term behavior repro duces the data accurately and another to capture daily fluctuations. A conv ex form ulation for the o v erall optimization problem is presen ted and its applicability demonstrated through the analysis of a case-study . The prop osed mo del is capable of repro ducing the diurnal patterns of a three-y ear dataset collected from a wind turbine lo cated in a moun tainous region in P ortugal. In addition, it is shown how to compute p ersistence statistics directly from the Marko v pro cess transition matrices. Based on the case-study , the p o wer pro duction persistence through the daily cycle is analysed and discussed. Keywor ds: Cyclic Mark ov pro cess, wind p o wer, p ersistence, diurnal pattern ∗ Corresp onding author Email addr esses: teresa.scholz@lneg.pt (T eresa Scholz), vitor.lopes@lneg.pt (Vitor V. Lopes), ana.estanqueiro@lneg.pt (Ana Estanqueiro) Pr eprint submitte d to Ener gy Octob er 14, 2013 1. In tro duction 1 The EC Europ ean Parliamen t ob jectiv e to achiev e 20% of the consumed 2 energy from the renewable energy sector by 2020 in tro duced a serious chal- 3 lenge to the planning and op erating of p o wer systems. Wind energy is b e- 4 coming a top contributor to the renewable energy mix due to rather high 5 capacities and generation costs that are b ecoming comp etitiv e with conv en- 6 tional energy sources [28]. Ho wev er, wind energy systems suffer from a ma jor 7 dra wback, the fluctuating nature of their source, which affects the grid secu- 8 rit y , the p o w er system op eration and market economics. There are several 9 to ols to deal with these issues, suc h as the knowledge of wind p o wer p ersis- 10 tence and wind sp eed or p o wer simulation. Persistence is related to stabilit y 11 prop erties and can pro vide useful information for bidding on the electricity 12 mark et or to maintain reliability , e.g. by setting reserve capacity . 13 Wind p o wer or sp eed sim ulation can b e used to study the impact of wind 14 generation on the p o wer system. F or this task, a sufficien tly long time series 15 of the p ow er output from the wind plants should b e used. How ever, real 16 data records are commonly of short length and th us synthetic time series 17 are generated b y sto c hastic simulation techniques to mo del wind activit y 18 [16]. Shamshad et al. [23] used first and second-order Mark ov c hain mo d- 19 els for the generation of hourly wind sp eed time series. They found that a 20 mo del with 12 wind sp eed states (1 m/s size) can capture the shap e of the 21 probabilit y density function and preserve the prop erties of the observ ed time 22 series. Additionally , they concluded that a second-order Marko v chain pro- 23 duces b etter results. Nfaoui et al. [15] compared the limiting b eha vior o f their 24 Mark ov c hain mo del with the data histograms gotten from hourly a veraged 25 wind sp eed and sho w ed that the statistical c haracteristics w ere faithfully re- 26 pro duced. Sahin and Sen [22] rep orted the use of a first-order Mark ov c hain 27 approac h to simulate the wind sp eed, where: a) both transitions b et w een 28 consecutiv e times and within state wind sp eeds are sampled using an uni- 29 form distribution; and, b) extreme states are sampled with an exp onential 30 distribution. They show ed that statistical parameters w ere preserved to a 31 significan t extent; how ever, second-order Marko v c hain mo dels could yield 32 impro ved results. 33 Although wind p o w er can b e computed from syn thetic wind sp eed time 34 series, Papaefth ymiou and Kl¨ oc kl [16] show that a sto chastic mo del using 35 wind p o wer leads to a reduced n umber of states and a lo wer Mark ov c hain 36 mo del order. They compared a Marko v chain based metho d for the direct 37 2 generation of wind p o w er time series with the transformed generated wind 38 sp eed. Both the auto correlation and the probability densit y function of the 39 sim ulated data show ed a go o d fit. Thus, they concluded that it is b etter to 40 generate wind p o wer time series. Chen et al. [7] also mo deled wind p o wer b y 41 using differen t discrete Mark ov chain mo dels: the basic Marko v mo del; the 42 Ba yesian Mark ov mo del, whic h considers the transition matrix uncertaint y; 43 and, the birth-and-death Marko v mo del, which only allows state transitions 44 b et w een immediately adjacen t states. After comparing the wind p o wer au- 45 to correlation function, the authors find the Ba yesian Mark ov model b est. 46 Lop es et al. [13] prop osed a Marko v chain mo del using states that combine 47 information ab out wind sp eed, direction and p o w er. F rom the transition 48 matrix, they compute statistics, such as the stationary p o wer distribution 49 and p ersistence of p o wer pro duction, which show a close agreemen t with 50 their empirical analogues. The mo del was then used for the t wo-dimensional 51 sto c hastic mo deling of wind dynamics b y Raisc hel et al. [21]. They aim at 52 studying the in teractions b et ween wind velocity , turbine aero dynamics and 53 con troller action using a system of coupled sto chastic equations describing 54 the co-evolution of wind p o w er and sp eed. They show ed that b oth the de- 55 terministic and sto c hastic terms of the equations can b e extracted directly 56 from the Mark ov chain mo del. 57 The knowledge of wind p o w er pro duction p ersistence pro vides useful in- 58 formation to run a wind park and to bid on the electricity mark et, since it 59 pro vides information ab out the exp ected p o wer steadiness. It can b e seen 60 as the av erage time that a system remains in a giv en state or a subset of 61 states. Existen t literature fo cuses mainly on wind sp eed p ersistence, which 62 is used for assessing the wind p o w er p oten tial of a region. Persistence can b e 63 determined directly from the data [20, 19]; how ever, the presence of missing 64 data leads to an underestimate of actual p ersistence. Alternative metho ds 65 are based on wind sp eed duration curv es [14, 10], the auto correlation func- 66 tion or conditional probabilities. Ko¸ cak [11] and Cancino-Sol´ orzano et al. [5] 67 compare these tec hniques, and both conclude that wind speed duration curve 68 yields the b est results, i.e. results that follow the geographical and climatic 69 conditions of the analyzed sites. Moreov er, Cancino-Sol´ orzano et al. [5] an- 70 alyze the concept of “useful p ersistence”, which is the time schedule series 71 where the wind sp eed is b etw een the turbine cut-in and cut-out sp eed. The 72 results gotten from this analysis coincide with the p ersistence classification 73 obtained using the sp eed duration curv es. In addition, Ko¸ cak [12] suggests a 74 detrended fluctuation analysis to detect long-term correlations and analyze 75 3 the persistence properties of wind sp eed records. Sigl et al. [24], Corotis et al. 76 [8] and Po je [19] prop osed an approac h based on the use of a p o wer law or 77 exp onen tial probabilit y distributions for the p ersistence of wind sp eed ab ov e 78 and b elow a reference v alue. A Mark ov c hain based metho d to derive the 79 distribution of persistence is introduced by Anastasiou and Tsek os [1], who 80 sho w its capability on wind sp eed data. 81 Most metho ds in literature of wind sp eed and p o w er syn thesis fail to rep- 82 resen t diurnal patterns in the artificial data. How ever, these are relev ant for 83 energy system mo deling and design, since their kno wledge allo ws to plan and 84 sc hedule b etter. F or instance, a p o wer pro duction b eha vior tha t b est matc hes 85 demand needs smaller reserve capacity . Recently , Suomalainen et al. [26, 25] 86 in tro duced a metho d for syn thetic generation of wind sp eed scenarios that 87 include daily wind patterns by sampling a probabilit y distribution matrix 88 based on five selected daily patterns and the mean sp eed of eac h da y . Cara- 89 p ellucci and Giordano [6] adopt a physical-statistical approac h to syn thesize 90 wind sp eed data and ev aluate the influence of the diurnal wind sp eed profile 91 on the cross-correlation b et w een pro duced energy and electrical loads. The 92 parameters of their mo del, such as diurnal pattern strength or p eak hour of 93 wind sp eed are determined through a m ulti-ob jectiv e optimization, carried 94 out using a genetic algorithm. 95 This pap er in tro duces a cyclic time-v arian t Marko v mo del of wind p o wer, 96 sp eed and direction designed to consider the daily patterns observed in the 97 data. The mo del can b e used to synthesize data for the three v ariables 98 and is capable of repro ducing the daily patterns. Moreo ver, it allo ws to 99 compute p ersistence statistics depending on the time of the da y . The pap er is 100 organized as follows: Section 2 in tro duces the prop osed mo del as an extension 101 of the “regular” Mark o v chain mo del, whic h is then used for comparison. 102 F urthermore it is sho wn, how to compute the time-of-the-day dep enden t 103 p ersistence statistics directly from the Marko v mo del transition matrices. 104 In section 3 the constrained con vex optimization problem to get the mo del 105 parameters is introduced and explained. It is applied to the analysis of a 106 case-study based on real dataset, section 4. Since the mo del describ es the 107 join t statistics for wind p o w er, sp eed and direction, Section 5 explains ho w 108 to create synthetic time-series for these v ariables. Section 6 compares the 109 syn thesized data of b oth the time-v arian t and the time-inv ariant versions of 110 the mo del. Moreov er, it is sho wn ho w the p ersistence of p o wer pro duction 111 v aries through the daily cycle. 112 4 113 Nomenclature 114 α 0 Initial state distribution at time step t = 0 115 β i,j µ Co efficien ts of the Bernstein p olynomial mo deling the transition proba- 116 bilit y p i,j ( t ) 117 1 A unit column v ector of the same size as subset A 118 P P 0 · . . . · P T − 1 119 A Subset of the state space, con taining the states of interest for p ersistence 120 S Set of observ ed state transitions 121 S z Set of transitions observ ed in the data together with the scaled time of 122 the da y z at whic h they are observ ed 123 ω W eight of the extra transitions added to the ob jectiv e function 124 π Stationary distribution of a time-inv ariant Marko v chain 125 π ∗ lim t →∞ P t 126 π r Stationary distribution at time r of a time-v arian t cyclic Marko v pro cess 127 π r ( j ) Stationary probability , of state j at time of the day r 128 π r ( A ) V ector whose elements are the stationary probabilities of the states in 129 the set A at time of the da y r 130 τ Persistence 131 τ r Time-dep enden t p ersistence in a cyclic Marko v pro cess 132 b µ,k ( z ) µ -th Bernstein basis p olynomial of order k 133 E [ ] Exp ected v alue op erator 134 P t t -th step transition matrix of a Mark o v pro cess 135 p i,j ( t ) t -th step transition probabilit y of a Mark ov pro cess 136 5 p av g i,j Daily a verage probability of transition from state s i to s j 137 r t Remainder of time step t mo dulo T 138 S Mark o v pro cess state space 139 s i i -th state of a Mark ov pro cess 140 T P erio d of a cyclic Marko v pro cess 141 t Time step of a Mark o v pro cess 142 X t Mark ov pro cess 143 z Scaled time of the day 144 π A Stationary probabilit y distribution of the states in subset A 145 r time of the da y 146 6 2. Time-inhomogeneous Mark o v mo del 147 2.1. Definition 148 A discrete finite Mark ov pro cess { X t ∈ S, t ≥ 0 } is a sto c hastic pro cess 149 on a discrete finite state space S = { s 0 , s 1 , ..., s n } , n ∈ N , whose future 150 ev olution dep ends only on its current state [9]. This Marko v prop ert y is 151 expressed mathematically b y 152 P r { X t +1 = s j | X t = s i ∧ X l ∈ S ∀ l = 0 , ..., t − 1 } = P r { X t +1 = s j | X t = s i } . 153 P r { X t +1 = s j | X t = s i } describ es the probability of the Marko v pro cess 154 mo ving to state s j at time step t + 1 giv en that it is in state s i and is called 155 the t -th step transition probabilit y , denoted as p i,j ( t ). Th us, for eac h time 156 step t the Mark ov pro cess has an asso ciated transition probabilit y matrix 157 P t , a n b y n matrix with entries [ P t ] i,j = p i,j ( t ) for all i, j ∈ { 0 , . . . , n } . 158 Eac h P t satisfies the follo wing prop erties: p i,j ( t ) ≥ 0 and P j p i,j ( t ) = 1 159 ∀ i, j ∈ { 0 , . . . , n } , ∀ t . A Marko v pro cess is called cyclic with p erio d T ∈ N , 160 if T is the smallest num b er, such that p i,j ( mT + r ) = p i,j ( r ) for all m in 161 N , 0 ≤ r < T [18]. Thus, a cyclic Marko v pro cess is describ ed by T transi- 162 tion matrices P r , r = 0 , ..., T − 1. The remainder of time step t mo dulo T 163 will b e denoted as r t and th us r t = r t + mT ∀ t, m ∈ N . 164 If the transition probabilities are time-independent, i.e. p i,j ( t ) = p i,j , the 165 pro cess is called a (time-homogeneous) Mark ov c hain and its probabilit y ma- 166 trix P ∈ R n +1 × n +1 is giv en b y [ P ] i,j = p i,j . By analogy to the time-dep enden t 167 case it holds that p i,j ≥ 0 and P j p i,j = 1 ∀ i, j ∈ { 0 , . . . , n } . 168 169 2.2. Communic ation classes and irr e ducibility 170 2.2.1. Time-invariant Markov chain 171 The probability of reaching a state s j from a state s i in l time steps is given 172 b y P l ( i, j ), i.e. the l -th p o wer of the transition matrix P . If a state s j can b e 173 reac hed from a state s i in a finite n umber of time steps and vice v ersa, i.e. 174 ∃ l ∈ N P l ( i, j ) ≥ 0 ∧ P l ( j, i ) ≥ 0, the states s i and s j comm unicate. All states 175 that communicate with eac h other are said to b e in the same communication 176 class. If all states of a state space are in the same communication class, i.e. 177 if it is p ossible to reac h ev ery state from an y other state in a finite n umber 178 of time steps, the corresp onding transition matrix P is called irreducible. 179 7 2.2.2. Cyclic time-variant Markov pr o c ess 180 A cyclic Mark ov pro cess with p erio d T is describ ed by T transition ma- 181 trices P r , one for each time of the day r = 0 , ..., T − 1. The probabilit y of the 182 pro cess reac hing state s j from state s i in l time steps at time t = 0 is giv en as 183 ( P 0 · . . . · P T − 1 ) m · P 0 · . . . · P r l ( i, j ) with l = mT + r l . F or an arbitrary time-step 184 t , the form ula must b e m ultiplied from the left with the term P r t · . . . · P T − 1 . 185 Th us, the Marko v pro cess is irreducible, if the matrix P = P 0 · . . . · P T − 1 is 186 irreducible, i.e. if ∃ l ∈ N P l ( i, j ) ≥ 0 ∀ i, j . 187 2.3. The stationary distribution 188 If a Marko v c hain is irreducible and ap erio dic then the long-term statistics of a Marko v c hain are describ ed b y the stationary probabilit y distribution: π = lim t →∞ α 0 P t . The distribution is indep enden t of the initial distribution α 0 and satisfies the balance equation π = π P . By the Perron-F robenius theorem it can b e computed as the normalized eigenv ector corresp onding to the eigen v alue 1 of the transition matrix [17]. In the case of the cyclic time-inhomogeneous Marko v pro cess there is also a stationary distribution π r , for all r < T . It can b e in terpreted as the limiting distribution of the Mark ov pro cess considering only the datap oin ts sampled at time of the day r . If the matrix P is irreducible, i.e. if ∃ π ∗ , such that π ∗ = lim t →∞ α 0 · P t and the process is ap erio dic, the stationary distribution π r exists and is given b y π ∗ · P 0 · . . . · P r − 1 , since it satisfies the balance equation 1: π r = π ∗ · P 0 · . . . · P r − 1 π r = π ∗ · P · P 0 · . . . · P r − 1 π r = π ∗ · P 0 · . . . · P r − 1 · P r · P r +1 · . . . · P T − 1 · P 0 · . . . · P r − 1 π r = π r · ( P r · P r +1 · . . . · P T − 1 · P 0 · P 1 · . . . · P r − 1 ) . (1) 2.4. Persistenc e 189 The p ersistence of a given state s i is related with the num b er of steps 190 the system consecutiv ely remains in this state. In the time-homogeneous 191 case, it follo ws a geometric distribution with exp ected v alue (1 − p i,i ) − 1 and 192 is denoted b y τ . Anastasiou and Tsekos [1] sho wed that it is p ossible to 193 determine the exp ected time that a Mark ov chain sta ys consecutively inside 194 a given subset of states using a simple closed-form expression. F or example, in 195 wind p o wer applications, a typical subset of interest could contain all states 196 8 corresp onding to p o w er pro duction ab o ve a giv en threshold. T o compute 197 this estimate, the states are ren umbered, s.th. they can be partitioned in to 198 t wo disjoin t subsets: A = { s ν , ..., s n } con taining the states of interest; and 199 A = { s 0 , ..., s ν − 1 } , its complemen t. Then, the transition matrix is rearranged 200 in to the following blo c k structure: 201 P =  A B C D  =          p 0 , 0 · · · p 0 ,ν − 1 p 0 ,ν · · · p 0 ,n . . . . . . . . . . . . p ν − 1 , 0 · · · p ν − 1 ,ν − 1 p ν − 1 ,ν · · · p ν − 1 ,n p ν, 0 · · · p ν,ν − 1 p ν,ν · · · p ν,n . . . . . . . . . . . . p n, 0 · · · p n,ν − 1 p n,ν · · · p n,n          , where the first and last blo c k of rows and columns corresp ond to the 202 states in subset A and A , resp ectiv ely . The exp ected v alue of p ersistence, 203 i.e. the exp ected n um b er of time steps the Marko v pro cess consecutiv ely 204 remains in the subset A once it is en tered, is giv en by: 205 E { τ } = π A 1 A π A C 1 B , where π A is the stationary probability distribution of the states in subset 206 A and 1 A is the unit column v ector of size ( n − ν + 1) × 1 [1]. 207 F or the time-inhomogeneous case, p ersistence τ t is defined as the num b er of 208 time steps the Marko v pro cess is expected to remain in a state (set of states), 209 once it is en tered at time t . F or a cyclic Mark o v pro cess, the p ersistence τ t 210 is equal for all t that are congruen t mo dulo T . Thus, it is only necessary 211 to compute the persistence for τ r , r = 0 , ..., T − 1. This can b e ac hiev ed b y 212 adapting the deriv ation of equation 2.4, pro vided by Anastasiou and Tsek os 213 [1], to time-v arian t cyclic Marko v pro cesses. 214 After renaming, s.th. the subset of interest is A = { s ν , ..., s n } , the states of 215 eac h transition matrix P r are rearranged as in equation 2.4. 216 9 P r =  A r B r C r D r  =            p 0 , 0 ( r ) · · · p 0 ,ν − 1 ( r ) p 0 ,ν ( r ) · · · p 0 ,n ( r ) . . . . . . . . . . . . p ν − 1 , 0 ( r ) · · · p ν − 1 ,ν − 1 ( r ) p ν − 1 ,ν ( r ) · · · p ν − 1 ,n ( r ) p ν, 0 ( r ) · · · p ν,ν − 1 ( r ) p ν,ν ( r ) · · · p ν,n ( r ) . . . . . . . . . . . . p n, 0 ( r ) · · · p n,ν − 1 ( r ) p n,ν ( r ) · · · p n,n ( r )            , The probabilit y of τ r to b e equal to l is given as: 217 P r ( τ r = l ) = P r ( X r ∈ A , ..., X r + l ∈ A , X r + l +1 / ∈ A | X r ∈ A , X r − 1 / ∈ A ) = X i ∈A P r ( X r = i | X r ∈ A , X r − 1 / ∈ A ) · P r ( X s ∈ A , r < s ≤ l , X r + l +1 / ∈ A | X r = i ) = X i ∈A X k / ∈A X j ∈A ˜ π r ( i ) · P i,j ( r , l − 1 , A ) · p j,k ( r + l ) (2) with 218 P i,j ( r , l, A ) = P r ( X r + l = j, X k ∈ A , 0 < k < l | X r = i ) = D r · . . . · D r + l − 1 · C r + l · 1 A and 219 ˜ π r ( i ) = P r ( X r = i | X r ∈ A , X r − 1 / ∈ A ) = P j / ∈A P r ( X r = i | X r − 1 = j ) · P r ( X r − 1 = j ) P i ∈A P j / ∈A P r ( X r = i | X r − 1 = j ) · P r ( X r − 1 = j ) = P j / ∈A π r − 1 ( j ) p j,i ( r − 1) P k ∈A P j / ∈A π r − 1 ( j ) p j,k ( r − 1) , i ∈ A , (3) where π r ( j ) is the long term probabilit y of o ccurrence (stationary proba- 220 bilit y) of state j at time of the day r ; also note that π t ( j ) = π r ( j ) for t = mr , 221 ∀ m ∈ N . Equation 3 can b e rewritten in the matrix form to include all states 222 in the subset A : 223 10 ˜ π r ( A ) = π r − 1 ( A ) · B r − 1 π r − 1 ( A ) · B r − 1 · 1 A , where 1 A is a unit v ector of dimension ( n − ν + 1) × 1 and π r − 1 ( A ) is a 224 v ector of dimensions 1 × ν , whose elements are the stationary probabilities 225 of the states in the set A at time of the da y r − 1. 226 Th us, equation 2 can b e rewritten as: 227 P r ( τ r = l ) = ˜ π r ( A ) · D r · . . . · D r + l − 1 · C r + l · 1 A = π r − 1 ( A ) · B r − 1 π r − 1 ( A ) · B r − 1 · 1 A · D r · . . . · D r + l − 1 · C r + l · 1 A . The exp ected v alue of p ersistence at time r can then b e derived as: 228 E ( τ r ) = ∞ X l =1 l · π r − 1 ( A ) · B r − 1 π r − 1 ( A ) · B r − 1 · 1 A · D r · . . . · D r + l − 1 · C r + l · 1 A Making use of the cyclicit y of the Marko v pro cess, this can b e expressed 229 as: 230 E ( τ r ) = ∞ X l =1 l · π r − 1 ( A ) · B r − 1 π r − 1 ( A ) · B r − 1 · 1 A · D m · D r · . . . · D r + r l − 1 · C r + r l · 1 A where D = D r · . . . · D T · D T +1 · . . . · D r − 1 and l = mT + r l . 231 232 It can b e seen that the sum conv erges after splitting it into T partial 233 sums, one for each time of the day r . F or each partial sum, the only term not 234 constan t is the matrix p o wer D m , which conv erges b ecause all D eigen v alues 235 are smaller than 1. The infinite sum for the exp ected v alue of p ersistence at 236 time r can b e appro ximated to an arbitrary degree of accuracy  b y defining 237 f l = l · π r − 1 ( A ) · B r − 1 π r − 1 ( A ) · B r − 1 · 1 A · D m · D r · . . . · D r + r l − 1 · C r + r l · 1 A . and successively adding f l , l = 0 , 1 , ..., L until the difference b etw een t w o 238 consecutiv e sums is smaller than  , i.e. un til | f L | <  . 239 11 3. P arameter estimation 240 3.1. Time-homo gene ous Markov chain 241 The common approac h to estimate the Marko v c hain transition matrix 242 P is through the optimization of a constrained maxim um lik eliho od func- 243 tion, which describ es the realization probabilit y of a given dataset [2]. F or 244 a sequence of M states { X 0 = s i 0 , ..., X M = s i M } with s i 0 , ..., s i M ∈ S 245 and i 0 , ..., i M ∈ { 0 , ..., n } , its probabilit y can b e computed as P r { X 0 = 246 s i 0 } p i 0 ,i 1 p i 1 ,i 2 · . . . · p i M − 1 ,i M . Since the term P r { X 0 = s i 0 } is constant, giv en a 247 set of observed state transitions S , it is p ossible to estimate P by maximizing 248 the lik eliho o d 249 O F 1 = Y ( i,j ) ∈S p i,j , (4) where a transition is described b y an ordered pair ( i, j ) indicating the ori- 250 gin and the destination of the transition. In practice, instead of maximizing 251 O F 1 with resp ect to the p i,j v ariables it is preferable to minimize the nega- 252 tiv e log-lik eliho o d function, i.e. − log( O F 1 ), since it transforms the original 253 mathematical programming problem in to an equiv alent one that is conv ex 254 and, thus, has a unique solution [4]. The ov erall optimization problem is 255 form ulated as follows: 256 min − X ( i,j ) ∈S log( p i,j ) sub ject to p i,j ≥ 0 ∀ i, j = 0 , ..., n X j p i,j = 1 ∀ i = 0 , ..., n The constrain ts ensure non-negativit y of the transition probabilities and that 257 they sum up to 1 for eac h ro w of the transition matrix. 258 3.2. Cyclic time-variant Markov pr o c ess 259 The goal of this time-v ariant Marko v pro cess is to get a mo del that accu- 260 rately repro duces the long-term b eha vior while considering the daily patterns 261 observ ed in the data. Th us, the prop osed ob jective function com bines tw o 262 maxim um likelihoo d estimators: the first term maximizes the lik eliho o d of 263 12 the cycle-av erage probabilit y; and, the second term maximizes the likeli- 264 ho od of the time-dep enden t probability . The final optimization problem is 265 transformed in to a conv ex one using the negative logarithm of the ob jective 266 function. This section pro vides a detailed description of the ob jective func- 267 tion, the parametrization of the time-v ariant probability functions, and the 268 constrain ts that must b e added to the optimization problem to ensure its 269 Mark ov prop erties. 270 3.2.1. Obje ctive function 271 The transition probabilities are considered to b e time-v arian t and cyclic 272 with a p erio d of T , i.e. for eac h time of the day r (= 0 , ..., T − 1) there is a 273 differen t transition matrix P r . In this pap er, the time-dep enden t transition 274 probabilities p i,j ( r ) are mo deled b y Bernstein p olynomials. This has sev eral 275 adv an tages: a) a p olynomial representation of the transition probabilities 276 leads to a con vex ob jectiv e function and constraints, i.e. the optimization 277 problem has a unique solution; b) a p olynomial represen tation allows to de- 278 crease the n um b er of v ariables in the optimization problem: for each transi- 279 tion, instead of T v ariables only k + 1 are needed for a k order p olynomial; c) 280 Bernstein p olynomials are non-negativ e, whic h simplifies probability mo del- 281 ing, when compared to other p olynomial bases; and d) they ha ve the conv ex 282 h ull property , which, com bined with de Casteljau algorithm, allows to easily 283 write probabilit y b oundary conditions. 284 Bernstein p olynomials are linear combinations of Bernstein basis p olynomi- 285 als b µ,k ( z ), z ∈ [0 , 1]. The k + 1 Bernstein basis p olynomials of order k are 286 defined as: 287 b µ,k ( z ) =  k µ  z µ (1 − z ) k − µ with µ = 0 , ..., k and  k µ  the binomial co efficien t. Th us, the transition 288 probabilities p i,j ( z ) are describ ed b y 289 p i,j ( z ) = k X µ =0 β i,j µ b µ,k ( z ) , with β i,j µ ∈ R and z = r T , since the p olynomial v ariable has to b e scaled, 290 s.th. it is b et ween 0 and 1. 291 292 13 Hence, to maximize the lik eliho o d of the time-dep endent transition prob- 293 abilities given the data, the ob jective function must consider the time of the 294 da y z when the transition happ ens. Therefore, the ob jectiv e function in tro- 295 duced in (4) b ecomes P ( i,j ) z ∈S z log( p i,j ( z )), where S z is the set of observ ed 296 transitions together with the time z when they happ ens. This ob jectiv e func- 297 tion allows to compute the in tra-cycle transition probability functions, and 298 th us to represent the daily patterns present in the data. 299 A second term is added to this function, namely P ( i,j ) ∈S log( p av g i,j ), where 300 S is the set of transitions observ ed in the data as defined in section 3.1 and 301 p av g i,j is the cycle-av erage (daily) probabilit y of transition from state s i to s j . 302 It can b e computed as follows: 303 p av g i,j = 1 1 − 0 Z 1 0 p ij ( z ) dz = Z 1 0 k X µ =0 β i,j k b µ,k ( z ) dz = k X µ =0 β i,j k Z 1 0 b µ,k ( z ) dz = 1 k + 1 k X µ =0 β i,j k This second term is the maxim um lik eliho o d estimator for the daily a v- 304 erage probability and its addition to the ob jective function increases the 305 consistency of the long-term b ehavior of the Marko v pro cess with the data. 306 Therefore, the o verall ob jectiv e function O F 2 is giv en by: 307 O F 2 = − X ( i,j ) ∈S log( p av g i,j ) − X ( i,j ) z ∈S z log( p i,j ( z )) = − X ( i,j ) ∈S log( 1 k + 1 k X µ =0 β i,j k ) − X ( i,j ) z ∈S z log( k X µ =0 β i,j µ b µ,k ( z )) and minimization is p erformed with resp ect to the co efficien ts β i,j µ (mo del 308 parameters). 309 3.2.2. Constr aints 310 The estimation of the mo del parameters requires the transition probabil- 311 it y functions to comply with several constraints, to ensure: 312 • C 0 - and C 1 -con tinuit y at z = 0, 313 14 • ro w-sto c hasticity of the transition matrices at ev ery time of the da y z 314 and 315 • that the transition probabilit y functions are non-negative and bounded 316 b y 1. 317 Th us, to complete the sp ecification of the optimization problem this section 318 explains all the necessary constraints required for the mo del parameters to 319 describ e a cyclic Marko v pro cess. 320 Perio dicity 321 The transition probability functions are mo deled using Bernstein p olyno- 322 mials, whic h are smo oth, i.e. C ∞ -con tinuous functions. In general, the v alues 323 at b oth ends of their domain (0 and 1) need not b e equal. Thus, to av oid 324 sudden changes in the v alue and slope of eac h probabilit y function b et w een 325 the cycles, t wo constrain ts are added to ensure C 0 and C 1 -con tinuit y . An- 326 other reason is the arbitrariness of the cycle starting p osition, whic h affects 327 the p osition of the discontin uity if these conditions are not used. 328 The first constrain t is p i,j (0) = p i,j (1). Since b µ,k (0) = δ µ, 0 and b µ,k (1) = 329 δ µ,k the constrain t can b e reform ulated as β i,j 0 = β i,j k , where δ is the Kroneck er 330 delta. The second constrain t is added to ensure C 1 -con tinuit y , i.e. dp i,j dz (0) = 331 dp i,j dz (1). The first deriv ative of a Bernstein basis p olynomial can b e written 332 as a com bination of tw o p olynomials of lo wer degree: 333 db µ,k dz ( z ) = k ( b µ − 1 ,k − 1 ( z ) − b µ,k − 1 ( z )) Th us, the first deriv ative of a transition probability p i,j ( z ) is giv en by: 334 dp i,j dz ( z ) = k ( k X µ =1 ( β i,j µ − β i,j µ − 1 ) b µ − 1 ,k − 1 ( z ) − β i,j k b k,k − 1 ( z )) Hence, using b µ,k (0) = δ µ, 0 and b µ,k (1) = δ µ,k as w ell as the first constraint 335 β i,j 0 = β i,j k ∀ i, j = 0 , ..., n , the constrain t dp i,j dz (0) = dp i,j dz (1) reduces to the 336 follo wing linear constraint: 337 β i,j k = β i,j 0 = 0 . 5( β i,j 1 + β i,j k − 1 ) R ow sto chasticity of tr ansition matric es 338 15 T o ensure ro w sto chasticit y of the time-v arian t transition matrices, it is 339 necessary to ensure that P j p i,j ( z ) = 1 for all i and z . Since the Bernstein 340 basis p olynomials of order k form a partition of unity , i.e. 341 k X µ =0 b µ,k ( z ) = 1 the constrain t can b e re-written as a linear combination of the p olynomial 342 co efficien ts: 343 X j p i,j ( z ) = 1 ⇔ X j k X µ =0 β i,j µ b µ,k ( z ) = 1 ⇔ X j β i,j µ = 1 Non-ne gative tr ansition pr ob abilities ar e b ounde d by 1 344 The most straightforw ard wa y to implemen t this constrain t is to add t wo 345 inequalities for eac h time of the da y and each transition probability p i,j , i.e. 346 p i,j ( z ) ≥ 0 ∀ i, j, z p i,j ( z ) ≤ 1 ∀ i, j, z (5) Ho wev er, this constrain t significan tly increases the problem size, since it 347 requires 2 · T · n 2 inequalities. An alternative constrain t can b e form ulated by 348 using the conv ex hull property of the Bernstein p olynomials. This constraint 349 mak es the o verall optimization problem size smaller, but is more restrictive. 350 Ev ery Bernstein p olynomial P k µ =0 β µ b µ,k ( z ) alwa ys lies in the con vex h ull 351 defined b y its control p oin ts ( k µ , β µ ), µ = 0 , ..., k . Thus the constraint 352 0 ≤ p i,j ( z ) ≤ 1 ∀ i, j = 0 , ..., n can b e reformulated in terms of the p olynomial co efficien ts as 353 0 ≤ β i,j µ ≤ 1 ∀ µ = 0 , ..., k (6) Since constraint 6 is a sufficien t but not necessary condition for con- 354 strain t 5, the reformulation leads to a more restrictiv e o verall minimization 355 16 problem, i.e. the optimum ob jective function v alue is alwa ys higher or equal 356 when compared with the problem with original constrain t 5. The conv ex 357 h ull b ound of Berstein polynomials can b e tigh tened by sub division, i.e. 358 b y sub dividing the domain in tw o regions and finding new con trol points 359 β i,j 0 (1) , ..., β i,j k (1) and β i,j k +1 (1) , ..., β i,j 2 k (1) such that the function output re- 360 mains unc hanged. With each sub division, the control p oints form a tighter 361 b ound around the p olynomial and th us the p olynomial co efficien ts can as- 362 sume v alues in a wider range. The new control p oin ts represent the polyno- 363 mial restricted to the tw o sub-interv als [0 , z ∗ ] and [ z ∗ , 1], where z ∗ ∈ [0 , 1] 364 is the cutting p oin t of the division. F or simplicity , z ∗ is fixed to 0.5 in all 365 transition probabilities. The new control p oin ts can b e determined b y linear 366 com binations of the original con trol p oin ts β i,j 0 , ..., β i,j k . This can b e p erformed 367 efficien tly using de Casteljau algorithm, which in matrix form is giv en as: 368    β i,j 0 (1) . . . β i,j k (1)    =      b 0 , 0 ( z ∗ ) 0 · · · 0 b 0 , 1 ( z ∗ ) b 1 , 1 ( z ∗ ) · · · 0 . . . . . . . . . . . . b 0 ,k ( z ∗ ) b 1 ,k ( z ∗ ) · · · b k,k ( z ∗ )         β i,j 0 . . . β i,j k    = C l    β i,j 0 . . . β i,j k    = C l · β ij (7) and 369    β i,j k +1 (1) . . . β i,j 2 k (1)    =      b 0 ,k ( z ∗ ) b 1 ,k ( z ∗ ) · · · b k,k ( z ∗ ) 0 b 0 ,k − 1 ( z ∗ ) · · · b k − 1 ,k − 1 ( z ∗ ) . . . . . . . . . . . . 0 · · · 0 b 0 , 0 ( z ∗ )         β i,j 0 . . . β i,j k    = C r    β i,j 0 . . . β i,j k    = C r · β ij (8) The sub division can b e applied recursively to further impro ve the con vex 370 b ound around the p olynomial. The corresp onding coefficients are computed 371 b y applying equations 7 and 8 to the left and righ t co efficien t ve ctors. Defin- 372 ing C = ( C l , C r ) T and I z as the iden tity matrix of dimension z × z , the 373 co efficien ts β i,j ( w ) = ( β i,j 0 ( w ) , ..., β ij 2 w k ( w )) after w sub divisions can b e ob- 374 tained b y: 375 β i,j ( w ) = ( C ⊗ I 2 w − 1 ) · ( C ⊗ I 2 w − 2 ) · . . . · ( C ⊗ I 2 ) · ( C ⊗ I 1 ) · β i,j where ⊗ denotes the Kronec ker pro duct. The num b er of inequalities 376 needed for the implemen tation of this constrain t is ( k + 1) · 2 ω +1 · n 2 . Th us, 377 17 its use only mak es sense if it decreases the problem size, i.e. for a num b er of 378 sub divisions ω suc h that ( k + 1) · 2 ω +1 ≤ T . 379 3.2.3. Pr oblem formulation 380 The ov erall optimization problem to b e solv ed for the estimation of the 381 transition probabilit y co efficien ts β i,j µ can b e written as: 382 min β i,j µ − X ( i,j ) ∈S log( 1 k + 1 k X µ =0 β i,j µ ) − X ( i,j ) z ∈S z log( k X µ =0 β i,j µ b µ,k ( z )) (9) sub ject to β i,j 0 = β i,j k ∀ i, j = 1 , ..., n (10) β i,j 0 = 0 . 5( β i,j 1 + β i,j k − 1 ) ∀ i, j = 0 , ..., n (11) X j β i,j µ = 1 ∀ i = 0 , ..., n ; ∀ µ = 0 , ..., k (12) β i,j ( w ) ≤ 1 ∀ i, j = 0 , ..., n (13) 0 ≤ β i,j ( w ) ∀ i, j = 0 , ..., n (14) where w is the num b er of sub divisions and k is the order of the Bern- 383 stein p olynomials, whic h ha ve to b e sp ecified. The ob jective function (9) is 384 a combination of tw o negativ e log-lik eliho o d functions to ensure the Mark o v 385 pro cess captures both the daily patterns and the long-term b eha vior of the 386 original data. The optimization is p erformed with respect to sev eral con- 387 strain ts: constraints (10) and (11) ensure C 0 - and C 1 -con tinuit y at z = 0. 388 The ro w-sto c hasticit y of the transition matrix is ensured b y constrain t (12). 389 The last tw o constraints (13) and (14) b ound the transition probabilities 390 b et w een 0 and 1. 391 It is exp ected that the ob jective function decreases with the polynomial 392 order and the num b er of sub divisions. The parameters of the Mark ov c hain 393 mo del β i,j µ are estimated b y solving the optimization problem using a rigor- 394 ous n umerical solver. The mo del w as formulated making use of the casadi 395 computation framew ork [3] and the optimization w as p erformed by ip opt, a 396 nonlinear in terior-p oin t solv er [27], which ensures con vergence to the global 397 optim um in the case of conv ex optimization problems. 398 18 4. Application of the cyclic Mark o v pro cess to wind turbine mo d- 399 eling 400 4.1. The data 401 The data for this study w as obtained from a wind p o wer turbine in a 402 wind park lo cated in a mountainous region in Portugal. The time series 403 consists of a three-year p erio d (2009-2011) of historical data gotten from the 404 turbine data logger. The sampling time of 10 min utes leads to 144 samples 405 eac h da y . The data-set comprises three v ariables, wind p o wer, sp eed and 406 direction (nacelle orien tation). The wind speed information w as collected 407 from the anemometer placed in the wind turbine h ub. Due to confiden tialit y , 408 wind p ow er and sp eed data v alues are reported as a fraction of the rated 409 p o w er and the cut-out sp eed, resp ectiv ely . 410 4.2. Markov chain state definition 411 Discrete Marko v chain mo dels require the definition of the states when 412 applied to describ e contin uous v ariables. This w ork prop oses to c haracterize 413 the wind turbine states using three differen t v ariables: wind p o wer, sp eed and 414 direction. As such, each state is defined by all the p oin ts inside a p olyhedron 415 in three-dimensional space. 416 Figure 1: Representation of all data p oin ts pro jected into the: a) wind direc- tion and sp eed plane (left); and, b) wind p o wer and speed plane (righ t). Each rectangle is the pro jection of a state polyhedron into the t wo planes. Overall, they define the final state partition for the three-dimensional v ariable space. Fig. 1 presen ts all data observ ations and the state partitions pro jected 417 in to: a) the wind direction and speed plane; and, b) the wind p ow er and 418 19 sp eed plane. As exp ected, the observ ations pro jected into the wind p ow er 419 and sp eed plane define the characteristic p ow er curv e of the wind turbine. 420 It sho ws the three operational regions of a wind turbine: a) below the cut-in 421 sp eed no p o wer is pro duced; b) b et ween cut-in and rated wind sp eed the 422 p o w er increases prop ortionally to the cub e of wind sp eed; c) at wind sp eeds 423 b et w een the rated and the cut-out wind sp eed, the turbine control system 424 limits the p o wer output to a constant v alue. In the wind direction and sp eed 425 plane, data is widely scattered and shows the dominan t wind patterns at the 426 site. Three accumulation regions can b e identified: one for low wind sp eeds, 427 cen tered on 0 . 25, whic h is the mo de of the wind sp eed, and tw o defining the 428 dominan t wind directions around 100 ◦ and 300 ◦ . 429 The data space is discretized unev enly to get a goo d resolution of the high- 430 slop e region of the p o wer curve. In a previous work [13], this partition was 431 used in a time-homogeneous Mark ov c hain and pro v ed to lead to an accurate 432 represen tation of the original data. The wind direction and p o wer are divided 433 b y an equally spaced grid leading to 12 ( { d 1 , ..., d 12 } ) and 20 ( { p 1 , ..., p 20 } ) 434 classes, resp ectiv ely . The wind sp eed is divided as follows: v alues b elo w the 435 cut-in sp eed define one class sp 1 ; b et ween the cut-in and rated wind sp eed 436 the discretization is narro w ed by selecting 10 classes ( { sp 2 , ..., sp 11 } ); and b e- 437 t ween the rated and cut-out wind sp eed discretization is widened and 4 classes 438 ( { sp 12 , ..., sp 15 } ) are defined. Data p oints with wind sp eed ab ov e the cut-out 439 wind sp eed are discarded. The complete state set is constructed b y listing 440 all p ossible combinations of the classes of eac h v ariable. Due to ph ysical 441 constrain ts b et ween the v ariables, most of the states are empty (fig. 1(left)) 442 and can are discarded. This reduces the num b er of states from 3840 to 778, 443 for this turbine. 444 4.3. A dditional tr ansitions to pr omote a single c ommunic ation class 445 The solution of the optimization problem describ ed in section 3 comprises 446 a set of transition matrices P r , r ∈ { 0 , ..., 143 } . Ho wev er, the constrain ts in 447 the optimization problem definition do not force the matrix P = P 0 · . . . · P 143 448 to b e irreducible and thus the Marko v pro cess to ha ve a single comm unication 449 class. So, during data synthesis, the Marko v pro cess can get “trapp ed” 450 within a comm unication class. T o induce the Mark ov pro cess to hav e a 451 single comm unication class, additional transition coun ts are introduced into 452 the data. The goal is to add a small set of transitions to promote state 453 connectivit y without distorting the original data. Thus, the set is comp osed 454 of transitions that connect neighboring states in the state space, since those 455 20 are the ones most lik ely to o ccur. 456 F or a state s i = ( p l , sp p , d q ), its neigh b orho od V is defined as 457 V ( s i ) := { ( p ll , sp pp , d q q ) : l l ∈ { l − 1 , l , l + 1 } , pp ∈ { p − 1 , p, p + 1 } , q q ∈ { q − 1 , q , q + 1 }} \ s i . (15) 458 It should b e noted that, unlike p o wer and sp eed, direction is a circular 459 v ariable, e.g. states d 0 and d 11 are considered neigh b ors. If a neighbor state 460 s j ∈ V ( s i ) is presen t in the dataset, a transition ( i, j ) is added to the set 461 of extra transitions S E . F or this dataset, originally consisting of 150601 462 transitions, 13610 transitions are added. 463 The extra transitions m ust be considered to happen at an unkno wn time of 464 the da y z . Thus, they can only b e accounted for in the ob jectiv e function 465 term without time information, i.e. only in the time-v ariant part of the 466 ob jective function. This directly affects the v alues for p av g i,j and, indirectly , 467 the mo del parameters. Since the aim is to cause a minimal impact on the 468 transition probabilities, the additional term is w eighed b y a factor ω < 1 469 to directly con trol its influence. In this work it is fixed to 0.05. Th us, the 470 follo wing term is added to the ob jectiv e function: 471 − ω · X ( i,j ) ∈S E log( p av g i,j ) = − ω · X ( i,j ) ∈S E log( 1 k + 1 k X µ =0 β i,j k ) (16) Although the use of the extra transition set do es not ensure the time- 472 v arian t Mark o v pro cess to ha ve a single comm unication class, results show a 473 decrease of the n umber of comm unication classes from 13 to 1 in this dataset. 474 5. Sim ulation of wind p o wer, sp eed and direction time series 475 T o simulate wind pow er, sp eed and direction time series the metho d de- 476 scrib ed by Sahin and Sen [22] is adapted to the cyclic time-v ariant Mark o v 477 mo del as follo ws. First, the cumulativ e probability transition matrices P cum r 478 with P cum r ( i, j ) = P j k =0 p i,k ( r ) are computed. Then an initial state s i , i.e. 479 X 0 = s i , is randomly selected. A new datap oin t X t +1 is generated by uni- 480 formly selecting a random num b er  betw een zero and one. The corresp ond- 481 ing state s new ( X t +1 = s new ) is chosen such that the probability of reac hing 482 it from the curren t state s i is bigger than  , i.e. such that P cum r t ( i, new ) ≥  . 483 Based on this discrete state sequence, a real v alue for the wind p o wer/speed/direction 484 v ariables is generated b y sampling each state partition uniformly . 485 21 6. Results and discussion 486 6.1. Daily p atterns in the data 487 The wind p ow er, sp eed and direction time-series clearly show a daily 488 time-dep enden t b eha vior. 489 Figure 2: Tw o-dimensional histograms of the original time series data: sp eed- time (top left), direction-time (top right) and p o wer-time (b ottom left). The subfigure on the b ottom right shows the p-v alue of the Kolmogorov-Smirno v test used to compare the wind sp eed distribution on the different times of the da y . Figure 2 shows that, on av erage, the turbine does not pro duce p o wer 490 b et w een 10 am and 3 pm. In this time in terv al, lo w wind sp eeds (0.1 - 0.25) 491 are the most lik ely even ts. There are tw o dominan t wind directions: around 492 100 ◦ and 300 ◦ . Moreo ver, they o ccur at sp ecific times of the da y; b et w een 5 493 and 10am, the wind typically blo ws from the 100 ◦ direction, the rest of the 494 da y from 300 ◦ . 495 22 T o assess whether these tw o dominant directions migh t b e due to sum- 496 mer/win ter seasonality , the dataset was divided in t wo subsets, one cov ering 497 the p erio d from April to September and the other from Octob er to March. 498 The histogram analysis sho ws that b oth, summer and winter subset, ha v e the 499 same tw o dominan t directions (figures not shown). Th us, it was concluded 500 that the time-dep enden t pattern is not induced by this seasonality . 501 Figure 2 b ottom-righ t plot shows the p-v alues of the Kolmogoro v-Smirnov 502 test applied to the wind sp eed distributions at different times of the da y . The 503 Kolmogoro v-Smirnov test is a nonparametric test for the equalit y of con- 504 tin uous one-dimensional probability distributions. Th us, the high p-v alues 505 around the diagonal illustrates that wind sp eed distributions for consecutive 506 times of the da y are similar. The same holds for wind sp eeds in the morning 507 and evening, i.e. b efore 9am and after 4:30pm. The wind sp eed distribution 508 b et w een 10am and 3pm is clearly differen t. 509 6.2. Choic e of p olynomial or der and numb er of sub divisions 510 The mo del in tro duced in section 3.2 has tw o parameters that need to 511 b e defined: k , the order of the Bernstein p olynomials used to mo del the 512 transition probabilities; and w , the n umber of sub divisions used to tigh ten 513 the conv ex h ull that b ounds the p olynomials. T o c ho ose prop er v alues for 514 these parameters, different mo dels w ere computed b y v arying k = 4 ... 10 515 and w = 0 ... 3. F or each model, synthetic data was generated follo wing the 516 pro cedure describ ed in section 5 and compared with the real dataset. 517 23 Figure 3: Bar plots comparing the ob jective function v alue (top left), the daily a v erage Jensen-Shannon distance b et w een original and syn thesized data (top right), the num b er of inequalities in the problem form ulation (b ottom left) and the CPU time sp en t in IPOPT solving the optimization problem (b ottom right) of the tested mo dels. Figure 3 sho ws bar plots comparing the different mo dels using four crite- 518 ria: the ob jective function v alue, the daily a verage Jensen-Shannon distance 519 b et w een original and synthesized wind direction data, the num b er of inequal- 520 ities in the problem formulation and the CPU time sp en t in IPOPT solving 521 the optimization problem. The Jensen-Shannon distance is the square ro ot 522 of the Jensen-Shannon div ergence d j s , whic h, for tw o discrete probability 523 distributions q 1 and q 2 is defined as: 524 d j s = 1 2 X i q 1 ( i ) q 2 ( i ) · q 1 ( i ) + 1 2 X i q 2 ( i ) q 1 ( i ) · q 1 ( i ) . Comparing the mo dels using the ob jectiv e function v alue (figure 3 top 525 left) shows a decrease of the ob jective function as the mo del order and the 526 24 n umber of subdivisions increase. It can b e seen that the impact of the num b er 527 of sub divisions is higher for mo dels with higher p olynomial order. Moreo ver, 528 the first sub division has the highest impact since it leads to the highest 529 decrease of the ob jective function v alue. The daily av erage of the Jensen- 530 Shannon distance (figure 3 top righ t) decreases with the p olynomial order 531 un til sixth order. The same b eha vior can be observ ed for the num b er of 532 sub divisions: un til the sixth order, the Jensen-Shannon distance decreases 533 with the n umber of sub divisions. The num b er of inequality constraints in 534 the optimization problem as well as the n umber of CPU seconds sp en t in the 535 solv er show the exp ected b ehavior (figure 3 b ottom). They increase linearly 536 with the p olynomial order and exp onen tially with the n um b er of sub divisions. 537 Based on these observ ations, a basis o rder of 6 with 2 subdivisions was chosen 538 as the b est trade-off betw een an accurate representation of the av erage daily 539 patterns and computational costs. 540 6.3. Capturing long-term statistics 541 This section compares the main statistical prop erties deriv ed from the 542 original data with the ones deriv ed from the data generated by the time- 543 v arian t Marko v mo del. 544 Figure 4: Comparison of the probability distribution of wind p o wer (left), wind sp eed (middle) and wind direction (righ t) of the original with the syn- thesized data. Figure 4 compares the wind p o w er (left), sp eed (middle) and direction 545 (righ t) distribution of the original with the synthetic data generated using 546 the Mark ov mo del. In general, the distributions are in close agreemen t. The 547 wind p ow er distribution is bimo dal, with the mo des lo cated at the minim um 548 and maximum p o wer. It shows that the intermediary p o wer lev els are rather 549 rare, for instance, the states corresp onding to a p o w er pro duction b et ween 550 0.4 and 0.9 ha ve a low probability . The wind sp eed distribution follo ws the 551 25 exp ected b eha vior, a single mo de distribution with a long tail for the high 552 wind sp eeds (W eibull distribution). The wind direction distribution is bi- 553 mo dal with the t wo mo des at 100 and 300 degrees, which are the prev ailing 554 wind directions at the turbine site (figure 2). 555 556 Figure 5: Tw o dimensional p o w er-time histograms of the original (left) and syn thetic (right) time-series data. Figure 5 shows tw o plots: on the left, the empirical 2D distribution of 557 the wind p o wer and direction computed from the data and, on the righ t, 558 the same distribution computed using the data generated by the Mark ov 559 mo del. Its comparison sho ws that the mo del captures the joint statistics for 560 the wind p ow er and direction from the data. It is p ossible to see the tw o 561 dominan t directions asso ciated with high wind p o w er pro duction, namely 562 the sectors from 100 to 120 and from 290 to 320 degree. Figures 2 and 5 563 clearly demonstrate the capabilit y of this Mark ov mo del to capture the com- 564 bined characteristics of the wind p o w er, sp eed and direction. The long-term 565 b eha vior of the mo del is close to what is observed in the dataset. 566 6.4. Capturing time-dep endent b ehavior 567 As shown in section 6.1, the original data clearly exhibits a time-dep enden t 568 b eha vior. T o test, if the time-dep enden t Marko v mo del can capture it, syn- 569 thetic data was generated and the histograms compared to the ones of the 570 original data. Moreov er, to obtain a comparison with the “regular” w ay of 571 data syn thesis with Mark ov mo dels, data w as also generated from the time- 572 in v ariant Marko v chain. 573 26 Figure 6: Two dimensional histograms of the syntheti c time-series data, generated with the time-v ariant Marko v mo del (left) and the time-inv ariant Mark ov c hain (right): p o wer-time (top), sp eed-time (middle) and direction- time (b ottom). The comparison of figures 2 and 6(first column) shows, that the time- 574 v arian t Marko v mo del is capable of repro ducing the time-v ariant b eha vior 575 of the data. Figure 6(second column) presen ts the results of using a time- 576 in v ariant Marko v c hain mo del, i.e. b y using constan t transition probability 577 27 functions. As exp ected, each v ariable statistic distribution remains constant 578 during the daily cycle. 579 6.5. Time-dep endent p ersistenc e of pr o duction 580 The time-dep enden t Mark ov model allo ws to compute the p ersistence of 581 p o w er-pro duction dep ending on the time of the day . Figure 6 sho ws the 582 time-dep enden t p ersistence of p ow er pro duction for differen t p o wer levels 583 ( i · 0 . 05 · p max , for i = 1 , ..., 19). The persistence analysis is presen ted for tw o 584 p o w er pro duction levels: a) p ersistence of useful p o wer pro duction (PUPP) 585 defined as ab o v e 0 . 15 · p max , i.e. the p o wer level corresp onding to the wind 586 sp eed mo de at the turbine site; and, p ersistence for high p o wer pro duction 587 (PHPP), i.e. abov e > 0 . 7 · p max . It can b e seen, that the higher the p o w er 588 lev el, the low er the p ersistence. Moreov er, for all p o wer lev els, p ersistence is 589 minimal b et ween 5 and 10 am. PHPP is fairly constant throughout the day 590 (dark line), the maximal differences are b et w een 10 and 30 min utes, whereas 591 PUPP reac hes a maximum at around 5 pm (white line). 592 28 Figure 7: Time-dep enden t p ersistence of p o wer pro duction ab ov e i · 0 . 05 · p max , for i = 1 , ..., 19. The lines highligh t the time-dep enden t p ersistence, for t wo conditions: a) in white, for p o wer pro duction ab o ve 0 . 15 · p max (PUPP); and, b) in blac k, the p o wer pro duction ab o ve 0 . 7 · p max (PHPP). Since the data sho ws tw o different dominan t directions (figure 4), figure 8 593 presen ts the p ersistence of p o wer production conditioned to each dominan t 594 direction, i.e. for the direction sectors from 90 ◦ to 180 ◦ and 270 ◦ to 360 ◦ . 595 29 Figure 8: Time-dep enden t p ersistence of p o wer pro duction ab ov e i · 0 . 05 · p max , for i = 1 , ..., 19 for direction sectors 90 ◦ -180 ◦ (left) and 270 ◦ -360 ◦ (righ t). As exp ected, the p ersistence in b oth direction sectors is low er than the 596 unconditioned p ersistence. F or wind directions in the sector 90 ◦ -180 ◦ , all 597 lev els of p o wer pro duction hav e a minim um p ersistence b et ween 80 and 100 598 min utes at around 1 pm. Maximum p ersistence is around midnigh t v arying 599 b et w een 220 minutes (PUPP) and 140 minutes (PHPP). It can be seen that 600 for p ow er lev els b elo w 50% of maxim um pro duction the time of da y dep en- 601 dency of p ersistence is very similar. F or p o wer lev els ab o ve 50% persistence 602 decreases as p o wer pro duction increases. How ever, the p ersistence v ariability 603 with the p o wer level is rather low, for example, maximum p ersistence at a 604 lev el of 75% is almost 180 minutes whereas for a level ab o v e 0.05% is 200 605 min utes. 606 F or wind directions in the sector 270 ◦ -360 ◦ , it shows that, for all p o wer lev- 607 els, the curv es for b oth PUPP and PHPP are similar, i.e. their minima and 608 maxima are lo cated around the same time of the day . F or instance, maximal 609 p ersistence of pro duction is reac hed at around midday . How ever, for this 610 direction sector, the higher the p ow er pro duction level, the lo wer the p ersis- 611 tence. F or p o w er pro duction abov e 0.05% of maxim um p o w er the p ersistence 612 is 250 minutes, p ersistence of pro duction ab o ve 75% of maxim um p o wer is 613 only 100 min utes. 614 Comparing with the other dominant direction, it can b e seen, that they hav e 615 v ery different p ersistence b eha vior. The maxima and minima are at different 616 times of the day for ev ery p o wer level. The p ersistence increases as p o wer 617 pro duction decreases, for all p o wer levels in the case of the 270 ◦ -360 ◦ sec- 618 tor. F or the 90 ◦ -180 ◦ sector it decreases only un til 50% of maxim um p o w er 619 30 pro duction. Below that, it remains approximately constant. 620 7. Conclusions 621 This pap er presents an inhomogeneous Mark ov pro cess to model wind 622 p o w er pro duction. It is developed using states, whic h combine information 623 ab out the wind sp eed, direction and p o w er v ariables, using real data recorded 624 b y a wind turbine in P ortugal. The join t partition of the three-dimensional 625 v ariable space allo ws to decrease the num b er of the mo del states and, si- 626 m ultaneously , enco des the wind p o wer curv e into the Marko v c hain mo del. 627 The transition probabilities are considered to b e functions that dep end on the 628 time of the da y and mo deled as Bernstein polynomials. The estimation of the 629 transition matrices is p erformed by solving a constrained conv ex optimiza- 630 tion problem. Its ob jectiv e function combines tw o log-likelihoo d functions 631 with the purp ose to accurately represent b oth the long-term b eha vior and 632 the daily fluctuations seen in the original data. T o ev aluate the statistical 633 prop erties of the estimated Marko v mo del, syn thetic time-series are gener- 634 ated and compared with the original data statistics. Results demonstrate 635 that the proposed Marko v mo del can reproduce the diurnal patterns in the 636 data. Moreov er it is demonstrated ho w the p ersistence of p o wer pro duction 637 throughout the time of the day can b e estimated from the Mark ov pro cess 638 transition matrices. 639 Ac kno wledgments 640 The authors thank the F unda¸ c˜ ao para a Ci ˆ encia e a T ecnologia for fi- 641 nancial supp ort (SFRH/BD/86934/2012, F COMP-01-0124-FEDER-016080 642 (PTDC/SENENR/1141718/2009)) and GENER G, SA. 643 31 References 644 [1] Anastasiou, K., Tsekos, C., 1996. Persistence statistics of marine envi- 645 ronmen tal parameters from marko v theory , part 1: analysis in discrete 646 time. Applied Ocean Researc h 18 (4), 187 – 199. 647 [2] Anderson, T. W., Go o dman, L. A., 1957. Statistical inference ab out 648 mark ov chains. The Annals of Mathematical Statistics 28 (1), 89–110. 649 [3] Andersson, J., Housk a, B., 2010. T o wards a Computer Algebra System 650 with Automatic Differen tiation for use with Ob ject-Orien ted mo delling 651 languages. Ob ject-Oriented Mo deling Languages. 652 [4] Bo yd, S. P ., V anden b erghe, L., 2004. Conv ex optimization. Cam bridge 653 Univ Pr. 654 [5] Cancino-Sol´ orzano, Y., Guti´ errez-T rashorras, A. J., Xib erta-Bernat, J., 655 2010. Analytical metho ds for wind p ersistence: Their application in 656 assessing the b est site for a wind farm in the state of V eracruz, Mexico. 657 Renew able Energy 35 (12), 2844 – 2852. 658 [6] Carap ellucci, R., Giordano, L., 2013. The effect of diurnal profile and 659 seasonal wind regime on sizing grid-connected and off-grid wind p o wer 660 plan ts. Applied Energy 107 (0), 364 – 376. 661 [7] Chen, P ., Berthelsen, K., Bak-Jensen, B., Chen, Z., 2009. Mark ov mo del 662 of wind p o wer time series using Bay esian inference of transition matrix. 663 In: Industrial Electronics, 2009. IECON ’09. 35th Annual Conference of 664 IEEE. pp. 627 –632. 665 [8] Corotis, R. B., Sigl, A. B., Klein, J., 1978. Probability mo dels of wind 666 v elo cit y magnitude and p ersistence. Solar Energy 20 (6), 483 – 493. 667 [9] Kemen y , J. G., Snell, J. L., 1976. Finite Marko v Chains. New Y ork : 668 Springer-V erlag. 669 [10] Ko¸ cak, K., 2002. A metho d for determination of wind sp eed p ersistence 670 and its application. Energy 27 (10), 967 – 973. 671 [11] Ko¸ cak, K., 2008. Practical w ays of ev aluating wind sp eed p ersistence. 672 Energy 33 (1), 65 – 70. 673 32 [12] Ko¸ cak, K., 2009. Examination of p ersistence prop erties of wind sp eed 674 records using detrended fluctuation analysis. Energy 34 (11), 1980 – 675 1985. 676 [13] Lop es, V. V., Sc holz, T., Estanqueiro, A., No v ais, A. Q., 2012. On 677 the use of Mark ov c hain models for the analysis of wind p o wer time- 678 series. In: En vironment and Electrical Engineering (EEEIC), 2012 11th 679 In ternational Conference on. pp. 770 –775. 680 [14] Masseran, N., Razali, A., Ibrahim, K., Zin, W. W., 2012. Ev aluating the 681 wind sp eed p ersistence for sev eral wind stations in p eninsular Malaysia. 682 Energy 37 (1), 649 – 656. 683 [15] Nfaoui, H., Essiarab, H., Sa yigh, A., 2004. A sto chastic Mark ov chain 684 mo del for sim ulating wind sp eed time series at T angiers, Moro cco. Re- 685 new able Energy 29 (8), 1407 – 1418. 686 [16] P apaefthymiou, G., Kl¨ oc kl, B., 2008. MCMC for wind pow er sim ulation. 687 Energy Con version, IEEE T ransactions on 23 (1), 234 –240. 688 [17] Pillai, S., Suel, T., Cha, S., 2005. The p erron-frobenius theorem: some 689 of its applications. Signal Pro cessing Magazine, IEEE 22 (2), 62–75. 690 [18] Platis, A., Limnios, N., Le Du, M., 1998. Dep endabilit y analysis of 691 systems mo deled b y non-homogeneous Mark ov c hains. Reliability Engi- 692 neering and System Safet y 61 (3), 235 – 249. 693 [19] P o je, D., 1992. Wind p ersistence in Croatia. International Journal of 694 Climatology 12 (6), 569–586. 695 [20] Pry or, S., Barthelmie, R., 2002. Statistical analysis of flo w c haracter- 696 istics in the coastal zone. Journal of Wind Engineering and Industrial 697 Aero dynamics 90 (3), 201 – 221. 698 [21] Raisc hel, F., Sc holz, T., Lop es, V. V., Lind, P . G., 2012. Uncov ering 699 wind turbine prop erties through t wo-dimensional stochastic mo deling 700 of wind dynamics. Unpublished results, arXiv:1210.7161. 701 [22] Sahin, A. D., Sen, Z., 2001. First-order Mark ov c hain approach to wind 702 sp eed mo delling. Journal of Wind Engineering and Industrial Aero dy- 703 namics 89 (34), 263 – 269. 704 33 [23] Shamshad, A., Baw adi, M., Hussin, W. W., Ma jid, T., Sanusi, S., 2005. 705 First and second order mark o v c hain mo dels for synthetic generation of 706 wind sp eed time series. Energy 30 (5), 693 – 708. 707 [24] Sigl, A. B., Corotis, R. B., W on, D. J., 1978. Run duration analysis 708 of surface wind sp eeds for wind energy application. Journal of Applied 709 Meteorology 18, 156–166. 710 [25] Suomalainen, K., Silv a, C., ao, P . F., Connors, S., 2013. Wind p o w er de- 711 sign in isolated energy systems: Impacts of daily wind patterns. Applied 712 Energy 101 (0), 533 – 540. 713 [26] Suomalainen, K., Silv a, C., F erro, P ., Connors, S., 2012. Syn thetic wind 714 sp eed scenarios including diurnal effects: Implications for wind p o w er 715 dimensioning. Energy 37 (1), 41 – 50. 716 [27] W¨ ac hter, A., Biegler, L. T., 2006. On the implementation of an interior- 717 p oin t filter line-searc h algorithm for large-scale nonlinear programming. 718 Mathematical Programming 106, 25–57. 719 [28] W en, J., Zheng, Y., Donghan, F., 2009. A review on reliabilit y assess- 720 men t for wind p o w er. Renew able and Sustainable Energy Reviews 13 (9), 721 2485 – 2494. 722 34