Wavelet-based ESS sizing strategy to enable power peak-shaving in PV systems

This work presents an Energy Storage Systems (ESS) sizing strategy, that allows grid connected PV systems to provide peak-shaved maximum-power-Ramp-Rate (RR) compliant power while storing exceeding power. For this purpose power generation and power loss caused by complying several RR restrictions were estimated. The estimation is based on irradiance and temperature real data from a measuring station located in the Atacama desert in Chile. The presented sizing strategy is based on Discrete Wavelet Transform (DWT) and allows to obtain an energy and power ESS sizing matching different criterion. Moreover, the selection of the ad-hoc wavelet and decomposition level, required filter the estimated PV power, is also presented.


I. INTRODUCTION
T HE sustained growth in participation of renewable resources in electric market has brought grid stability concerns [1], [2]. The intermittent and variable nature of renewable resources is a major drawback when discussing further participation of renewable resources in electric markets, setting a virtual limit to high penetration rates [2]- [4].
Moreover, many countries have upgraded their grid codes addressing maximum power Ramp Rate (RR), maximum power variation per unit of time, regulations for renewables. For example, the grid-code of Hawaii limits RRs to ±2 or ±1 MW per minute, depending on the time of day [5]. In Puerto Rico RRs are limited to ±10 % rated power per minute for PV systems [6], while in Australia RRs are limited to ±16.67 % rated power per minute [7]. Other grid codes limit only positive RRs (load taking) as is the case of Ireland, Germany and Chile, where 30 MW per minute, 10% rated power per minute and 0 to 20 % rated power per minute are the corresponding positive RR restrictions [5], [8], [9].
ESSs have been vastly researched as an alternative to deal with generation-demand mismatch and provide additional services as peak shaving, base load generation, capacity firming, load shifting, etc [3], [8], [10]- [15]. All previously named additional services aim at modifying the power output of a system by storing and releasing power, according to the demand of the selected service. Peak shaving mitigates power peaks by storing (releasing) exceeding (missing) power, required to match a desired lower frequency power output, thus dealing with variability-intermittency drawbacks.
ESS sizing strategies have became a primary concern when looking for cost-effective solutions, therefore it is not surprising to find several different ESS sizing methods in the literature. In [13], a desired wind-farm-output-power profile (base load generation) is compared to a low pass (FIR) filtered version of a long-horizon wind-speed measurement with added statistical noise, generating a ESS power curve reference. The resulting wave form is fed to a cost function aiming to minimise Battery Energy Storage System (BESS) cost per service-hours; cost function factors include life-span, full charge/discharge cycles, initial investment and penalties for not meeting the agreed power curve form. Nonetheless, frequency depending phase delay introduced by low pass filtering reshapes the ESS power curve reference, hence causing over or under sizing of the ESS. In [14], an ESS sizing strategy that aims at minimising economical penalties caused by not complying day-ahead power bidding is presented. Bid power is calculated from historical wind-power data, while ESS sizing is estimated by comparing bid power with several scenarios including statistical deviations from forecasted wind power. Sizing is performed by presenting, in a histogram, the errors between bid power and the different deviation scenarios as a function of energy, the ESS energy rating is selected by choosing a certain biding compliance ratio. This method presents a high dependence on the bidding strategy and wind power forecast, which are not described in the document. Moreover, the sizing strategy considers inter-day dynamics, which are not relevant for solar daily cycles (dawn to dusk). A sizing strategy to comply with maximum power RR regulation is presented in [15]. In order to achieve such requirement, wind forecast and uncertain load behaviour are analysed. The strategy considers transforming the difference between estimated power generation and load consumption from time-domain to frequency-domain through Discrete Fourier Transform (DFT). The result is later filtered, assigning low frequencies (complying with maximum power RR) as the desired power output, medium frequencies as BESS sizing reference (with a Depth of Discharge (DoD) ≤ 80%) and high frequencies as Super Capacitor sizing reference. DFT strategy, decomposes the full signal into periodic sinusoids, loosing information regarding the time-location of frequencies, therefore causing ESS over-sizing. Short Time Fourier Transform (STFT) improves frequencies location accuracy, though resolution in time/frequency is limited due to the fixed timelength moving window used to analyse all frequencies. An strategy to limit power fluctuations in a wind farm by limiting the maximum frequency of exported power is proposed in [8]. The strategy considers the addition of a hybrid energy storage system, formed by a lithium ion BESS and a super capacitor bank. Several historical data sets are used to estimate the ESS sizing, each data set is filtered by Discrete Wavelet Transform (DWT), hence separating the data set in two, a part containing frequencies below a maximum desired frequency (desired export power profile) and a part containing higher frequencies (ESS power profile). ESS profile curves are averaged and later separated (filtered) into lower frequencies, considered BESS requirements, and higher frequencies, to be handled by the super capacitor bank. The efficiency of the converters, maximum BESS DoD and BESS charge/discharge cycles are also considered to generate the BESS and super capacitor bank respective sizing. Even though wavelet strategy improves time-location and resolution problems related to DFT and STFT, this strategy relies on averaging the results hence hiding dynamic behaviours and providing a non-accurate sizing.
This work presents a ESS sizing strategy based in DWT, the strategy provides an in depth analysis of the required energy and power ratings for the ESS and its power converter, in order to provide RR compliant peak-shaved power. An assessment of annual energy loss, caused by complying different grid-codeimposed maximum power RR regulations in an hypothetical 2 MW PV plant, is also presented. For this purpose empirical irradiance and temperature measurement, from Antofagasta (Chile were used to estimate annual exported power of the PV plant. To the best knowledge of the authors both, the assessment of annual losses due to RR compliance and the DWT-based ESS sizing strategy correspond to new knowledge.
The document is arranged in the following manner, section II provides a full description of the PV plant model and rating considered to estimate a 2 MW PV plant power generation during a full year, an analysis of annual energy loss due to RR regulation is presented in section III, the ESS sizing strategy based in DWT is presented in section IV, section V presents the conclusions of this work.

II. PV PLANT POWER ESTIMATION
A year of irradiance and temperature measurements, taken from a monitoring station located in the Atacama desert in Chile (−24.0833 • , −69.9167 • ), where used to estimate power generation of an hypothetical 2 MW PV plant, formed by 447 strings of 20 Jinko JKM 60PP PV modules each [16]. The PV array estimated power at STC has been oversized by 15% respect to the inverter rated power (2 MW), in accordance with common commercial practices to increase plant factor and correct power reduction caused by temperature.
There are several models capable to estimate power generation from a PV plant, some of them are circuital models [17], Artificial Neural Networks (ANN) [18], algebraic equation [19], etc. Circuital models emulate the behaviour of a PV module as an electrical circuit, being single and double diode circuital models the mainstream models. The set of equations describing these models are usually solved through numerical methods as Newton-Raphson or Bisection [17], [20]- [22]. These models require parameter identification, performed for example through any of the strategies described in [23]. ANNs are another modelling alternative, where inputs pass through an array of neurons, formed by bias, weights and activation (transfer) functions in order to estimate an output. The training procedure consists in adjusting the weights and bias in order to minimise the mean squared error between the estimated output and a known target [19], [24], [25]. A similar alternative to ANN corresponds to Polynomial Regression where the estimated output is obtained through a polynomial applied to a set of inputs. The polynomial coefficients are determined by applying regression to the set of inputs and target(s) [19]. However, circuital models have a highly dependence on parameter identification and the accuracy of the numerical method applied to solve the equation system. On the other hand, ANN and Polynomial Regression, depend on knowing the target (desired output), hence are not useful to estimate the power yield just from irradiance and temperature.
Alternatively, algebraic equations are also capable of estimating PV power generation from irradiance and temperature measurements, and do not depend on parameter identification, depending instead on module parameters provided by the PV module manufacturer (data sheet). This mathematical model estimates the maximum power point (P mpp in W) as a function of the current irradiance (G in W/m 2 ) and cell temperature (T in • C).
Where γ, T stc , G stc , P mpp stc , η and N are respectively respectively the Temperature Coefficient of P mpp in W / • C(or maximum power correction factor for temperature [19]), STC cell temperature in • C, STC irradiance in W/m 2 , maximum power point at STC in W, PV module efficiency and number of modules in the PV plant [19]. Cell temperature (T ) is estimated as a function of the ambient temperature (T a ) through [26]: Where G noct and T noct are the irradiance (in W/m 2 ) and cell temperature at NOCT (in • C) at Nominal Operating Cell Temperature (NOCT). This latter model was selected to estimate the annual exported power of the hypothetical 2 MW PV plant. For this purpose uniform irradiance and temperature conditions, sampled every minute, and PV plant parameters, presented in Table I, were applied.

III. MAXIMUM POWER RAMP RATE LOSSES
Chilean grid code limits positive RRs (load taking) between 0 and 20% of the PV plant nominal power per minute, while negative RRs are not restricted. The instantaneous positive RR level is stated by the Transmission System Operator ("Sistema Electrico Nacional") according to its requirements in order to guarantee grid stability [9]. Power ramps exceeding RR regulation are not allowed, forcing the system to operate in a non-optimal power point and thus causing non-estimated annual energy losses. An analysis of annual energy loss caused by complying with positive power RR limitation is shown in Fig. 2. Where RR x % corresponds to a maximum power RR of x% per minute of the nominal power of the PV plant. The top plot (Fig. 2a)) presents histograms of power exceeding positive RRs of 5, 10, 15 and 20 % of the nominal PV plant power (2 MW). While, the bottom plot (Fig. 2b)) shows the recoverable annual energy loss as function of the power loss (power exceeding each RR regulation).
From Fig. 2 applying constantly RR regulations of 5%, 10%, 15 % or 20 % will generate annual energy losses of 43, 20, 11 or 7 MWh, respectively. Which correspond to 0.087%, 0.04%, 0.022% or 0.014% of the annual energy generation, hence annual energy losses caused by RR complying are negligible. Nonetheless, RR regulation must be complied, by loosing or storing the RR % non complying. A RR 10% was selected as limit to be used during rest of the document, since this presents a good compromise between 0 and 20 %, and it matches other grid-codes RR regulations.

IV. ESS SIZING STRATEGY
The ESS sizing strategy consists in filtering the estimated PV power (P mpp ) through DWT (see Appendix V-A), separating the PV power into a peak-shaved RR-10%-compliant power curve (P out ) and a non-peak-shaved RR 10% noncompliant power curve (P ess ). This latter curve is obtain from the subtraction P mpp − P out , and corresponds to the power reference to be consider for the ESS sizing process. Daubechies 8 decomposition level 4 was applied to perform the filtering process. The selection of this Wavelet and the required decomposition level is presented in subsection IV-A.
ESS power analysis is presented in Fig.3. This analysis is performed over the ESS maximum accumulable daily energy the ESS is capable to store if limiting bi-directional power to a certain power level. For this purpose daily ESS charge depletion was considered. ESS maximum accumulable daily energy will be addressed as EMADE and accordingly the maximum annual value of EMADE, the annual average of EMADE and annual standard deviation of EMADE will be identified asÊ made , E made and σ emade . Curves depicted in Fig.3 correspond to C m =Ê made , C a3s = E made + 3 · σ emade , C a2s = E made + 2 · σ emade , C a1s = E made +1·σ emade and C a = E made . The bottom plot in Fig.3 presents a histogram of positive (store) and negative (release) power. The decrease in EMADE curves after reaching the maximum, is due to limiting ESS charge to positive values, hence negative power surpassing positive power implies that more power is being drawn from the ESS thus reducing the ESS energy requirements.
Finally, the ESS sizing strategy consists in selecting a bi-directional power level and choosing an ESS maximum accumulated daily energy curves from Fig. 3. As example we have selected a bidirectional power rating of 200 kW and selected curve C a3s (E made + 3 · σ emade ) as the compliance criteria, hence ESS energy rating is 120 kWh. It must be noted that most ESS power reference (P ess ) is concentrated between 0 and 200 kW.

A. Wavelet selection
There are several wavelet families capable of filtering the estimated PV power (P mpp ) and producing a peak-shaved RR 10 % compliant power output, though each wavelet will produce a different ESS requirement. In order to select the best alternative the filtering capabilities of Wavelets Haar, Daubechies 2 to 10, Symlets 2 to 8 and Coiflets 1 to 5 were compared. ESS daily depletion was considered since solar cycle (dawndusk) allows to inject the stored power within a day, reducing energy losses associated to longer storage periods caused by ESS self discharge. Therefore a day-oriented analysis was performed, for this purpose a total of 1440 pairs of irradiance and temperature measurements were taken every day, since each wavelet level halves (down samples) the amount of data samples, hence only 5 DWT levels are possible without modifying the length of the data set.
As means to illustrate the effects of Wavelet filtering, to obtain peak-shaved RR 10% compliant power, Figs. 4 and 5  show respectively filtering of P mpp through Wavelets Symlet 2 (decomposition level 5) and Daubechies 8 (decomposition level 4), where the top plots present the wavelet-filtered output power (P out ) compared to the estimated PV power (P mpp ) and the bottom plot shows the required ESS power (P ess = P mpp − P out ). A summary of the decomposition level each wavelet requires to comply with RR 10 % is presented in Table II. Note that Haar, Daubechies 2 and Symlets 2 do not comply with RR 10 %, despite of applying a 5 level decomposition. Fig. 6 presents the annual maximum positive and negative RRs obtained by filtering the PV power through the 3 worst performing wavelets (Haar, Symlet 2 anc Coiflet 1) and best performing Wavelets (Daubechies 8, Symlet 6 and Coiflet 2) as function of the decomposition level (1 to 5).
Daubechies 8 level 4 was selected as the Wavelet to filter the estimated PV powerb (P mpp ), since it complies with RR 10 %, presents even and low positive and negative RRs. Lower wavelet decomposition levels were preferred, since it requires less computational effort to perform the decomposition.

V. CONCLUSION
This document presents an ESS sizing strategy based in DWT. The strategy aims at providing peak-shaved PV power, while complying with a Ramp Rate restriction of 10 % of the PV plant nominal power. For this purpose annual power generation was estimated based on empiric irradiance and temperature measurements from the Atacama desert in Chile.
Ramp Rate restrictions limit the maximum power point tracking of solar systems, thus forcing the PV plant to perform in a non-optimal power point and causing non-estimated annual energy losses. This work presents an annual energy loss estimation due to complying several RR restrictions. From the results complying with RR restrictions does not generate meaningful annual losses (0.087% for RR 5%), though its compliance is mandatory.
The presented ESS sizing strategy generates a peak-shaved RR 10% compliant power export. The strategy is based in DWT filtering, which generates a ESS power reference. ESS power reference is analysed presenting a relation between ESS bi-directional power and ESS maximum accumulated daily energy, this latter result enables to size the ESS energy rating from a selected bidirectional power rating. In this line, the methodology applied to select a proper Wavelet and decomposition level are presented. Moreover, the presented sizing technique can be extrapolated to other renewable sources, provided the proper modelling and relevant modifications, as resource cyclic behaviours.

A. Wavelet Transform
Wavelet Transform is a time-frequency signal analysis tool that allows to decompose simultaneously a time signal (x(t)) in time and frequency domains [27]- [29]. In order to perform the decomposition a wavelet function (Ψ j,n (t)) and a scaling function (Φ j,n (t)) are required. The inner product between the scaling function with the signal generates an approximation (L x (j = 1, n), lower frequencies) of the signal x(t), while the inner product between the wavelet function with the signal x(t) keeps the details (W x (j = 1, n), higher frequencies); thus working respectively as a low-pass and high-pass filters. This decomposition process can be repeated separating the obtained approximation (L x (j = 1, n)) into a following level of details (W x (j = 2, n)) and approximation (L x (j = 2, n)), and so on (as shown in eq. (3) and (4)).
x(t) · Ψ * j,n (t)dt (4) The algorithm applied in this work corresponds to a DWT calculated through Fast Wavelet Transform (FWT) algorithm introduced in [27]. FWT mathematical procedure is shown in eqs. (5) and (6) [29], where a j [n] is decomposed into an approximation (a j+1 [n]) and details (d j+1 [n]) by convoluting a j [n] with a low-pass and high-pass conjugate mirror filters (h and g). This procedure can be cascaded generating further decomposition of the approximations obtained. The sub-index j indicates the level of decomposition.
The re-composition of the signal is perform through the Inverse Fast Wavelet Transform (IFWT) [29].