A framework for the assessment of electric heating load flexibility contribution to mitigate severe wind power ramp effects

Electric Power Systems Research 142 (2017) 268–278

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

A framework for the assessment of electric heating load flexibility contribution to mitigate severe wind power ramp effects Antti Alahäivälä a,∗ , Jussi Ekström a , Juha Jokisalo b , Matti Lehtonen a a b

Department of Electrical Engineering and Automation, Aalto University, Espoo, FI-00076 Aalto, Finland Department of Energy Technology, Aalto University, Espoo, FI-00076 Aalto, Finland

a r t i c l e

i n f o

Article history: Received 5 October 2015 Received in revised form 20 September 2016 Accepted 21 September 2016 Available online 15 October 2016 Keywords: Demand response Wind generation Heating load Power ramps

a b s t r a c t A great share of wind power generation in a power system is likely to test the system flexibility with extreme power ramp events. Therefore, this paper proposes a framework that allows the assessment of electricity demand flexibility on these events. Two electric heating methods are of particular interest, since heating is an attractive source of flexibility in a cold climate. The proposed framework consists of generating wind generation scenarios, identifying severe ramps, and solving a system unit commitment (UC) problem, with and without demand response (DR), when these ramps occur. These stages are integrated into a Monte Carlo simulation technique, which allows the capturing of the electric heating load contribution to mitigate ramp effects in different system conditions. The contribution is investigated in case studies and evaluated with several measures. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Power systems with high wind power penetration are required to address the challenges of the variability and uncertainty of wind power generation. It is commonly recognized that the need to overcome these challenges increases the requirement for system reserves in different timescales [1]. Uncertainty and variability truly test the flexibility of a power system during severe power ramps, when the system needs to allocate adequate reserves while following the power ramp. Therefore, different methods to increase the system flexibility are of interest. Several methods have already been proposed in the literature to address the challenges of uncertainty and variability. Stochastic approaches are effectively utilized to tackle, for example, the problem of estimating sufficient spinning reserves [2], or the problem of clearing a market with reserves [3] in an environment with high wind power penetration. In addition to the stochastic methods, more frequent generation scheduling and higher scheduling resolution are seen as solutions [4], as well as market integration [5]. Severe ramps need special attention since they can have a great impact on power system operation. Conventionally, the supplyside maintains sufficient capacity to respond to the morning and

∗ Corresponding author. E-mail address: antti.alahaivala@aalto.fi (A. Alahäivälä). http://dx.doi.org/10.1016/j.epsr.2016.09.026 0378-7796/© 2016 Elsevier B.V. All rights reserved.

evening ramps in demand. However, wind generation fluctuation produces even steeper net load ramps. At the worst situations, they may cause loss of load or energy, and therefore uneconomic dispatch of generation [6], as the existing generation mix is not capable of following the variations. Conventionally, the system’s ability to supply the demand has been of interest in power system reliability evaluation, where system faults are a major concern [7]. The authors in [7] recently considered the ramp rates of generating units in their reliability analysis and indicated that reduced ramping capability reduces the system reliability. Any investigation of severe ramps requires that they are first properly characterized. [8] provides different definitions for ramp events and tackles the problem of detecting them. Since severe ramps occur rarely, modeling them can be beyond traditional forecasting methods. Therefore, Ganger et al. in [9] propose methods to model severe ramps. In addition to the characterization of ramp events, strategies to cope with them in power system operation can be found. In [10], the authors formulate an economic dispatch problem with ramp capability constraints. The ultimate goal in the paper is to make the system more robust against intermittent generation, and motivate suppliers and other resource owners to provide the system with flexibility. Furthermore, [11] proposes a dispatch strategy with increased dispatch resolution during ramps to deal with rapid ramp events. Many of the aforementioned strategies to advance wind power integration rely on improved operation strategies in the power

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

system. The integration can also be assisted by unleashing the flexibility of the demand-side. The flexibility of the demand-side, or demand response, has been seen to be applicable to several power system operations [12]. A typical approach to study the effect of DR on system operation is to dispatch the flexible demand with the generation. Two clearly different methods to model the DR can be distinguished. Firstly, the flexibility can be modelled by using the price elasticity of the demand, which is a particularly useful approach in market clearing problems [13,14]. However, for research aiming to study the DR potential, this method is challenged since the price elasticities are difficult to define. Another method integrates the demand flexibility into a UC problem that utilizes it to minimize the operation cost. In [15], a part of the power system demand is assumed to be shiftable, and is attached to a UC model. This is seen as a centralized approach, where the system operator has an aggregated flexible load model and access to central control of the loads. Even though this method may be difficult to implement in practice, it provides a good approach to study the full potential of DR. Thus, the integration of a flexible demand model to generation dispatch is also exploited, for example, in [16–18]. Particularly in a cold climate, electric heating load provides a promising source of flexibility. As shown in [19], electric heating load is certainly flexible, especially when it includes a storage element. With a direct approach, i.e., radiators that directly heat indoor air, the flexibility is more limited, but, thanks to the thermal inertia of a building, the heating can be controlled without instantaneous violation of customer comfort. In the literature, heating load has also been considered as a viable option to assist wind power integration. The authors in [18] couple a heat pump model with market dispatch in order to study the DR potential for wind power integration. This approach allows the authors to study how thermal storage can benefit market operation. Furthermore, the author in [20] investigates the potential of space heating to assist in wind power integration, in terms of matching the wind production and heat demand, with and without additional thermal storage. In addition to the studies on the potential of DR, there exists extensive literature covering the problem of managing residential thermostatically controlled loads, which includes applications, such as electric space heating. The methods can vary from local cost-minimization using electricity price [21] to a centralized approach with generic control signals [22,23] or direct load control [24]. As shown above, the challenges of wind generation have been comprehensively investigated in the earlier literature and various methods to tackle these challenges have been proposed. The motivation for our study arises from the fact that severe net load variations may have a undesirable effect in a power system, as a result of which they should be given special attention. Thus, this paper proposes a framework for the assessment of the demand-side contribution to mitigate the effect of severe ramps with a focus on electric heating. The framework can be used as a tool to analyze rare and severe ramps, as it captures the ramps’ effects under different system conditions. The core of the framework is a novel combination of a UC, a detailed heating load DR model, and a stochastic inclusion of severe power ramp occurrence. Due to the stochastic nature of wind power generation, the occurrence of the most severe variations vary in time, resulting in different consequences in the power system. To capture the effect of ramps on the system operation in different seasons and times of day, the proposed approach is based on Monte Carlo simulation. We divide the simulation into the following stages, which are repeated for several years. Firstly, a wind generation time series is generated for one year and this time series is subtracted from the load profile. Secondly, days with severe ramps are detected from this net load profile. Thirdly, for the detected days, the rest of the generation is first dispatched with flexible electric heating and then with fixed load by solving

269

Fig. 1. The overview of the framework.

a UC problem. Thus, the results capture the contribution of electric heating load under different base load, heat demand, and wind generation conditions. The following summarizes the core contributions of the paper: • A methodology enabling the systematic identification of severe power ramps and an evaluation of DR impact on their effects is introduced. • Performance measures to assess the influence of DR on the ramp effects are proposed. This paper first provides an overview of the framework in Section 2, after which it presents the wind production and heating load modeling in Sections 3 and 4, respectively. In Section 5, the actual simulation algorithm is described, together with ramp detection and performance measures. Case studies are conducted in Section 6 and conclusions are drawn in Section 7. 2. Overview of the framework The proposed framework with its main components is summarized in Fig. 1. The core is the actual simulation algorithm, which generates wind production scenarios, detects severe ramps from a net load profile, solves a UC with and without DR, and calculates and stores the investigated measures. The implementation naturally requires the specification of different components and the definition of a severe ramp, i.e., the ramp rate threshold describing the minimum steepness of the studied ramps. Since severe ramps are mainly caused by wind production, it is modeled with a time series model capable of capturing the variable and stochastic nature of wind production. Furthermore, before the simulation, the studied power system is modeled with desired accuracy by combining the system data with a UC formulation and a DR model. The simulation itself produces a significant amount of output data, from which the DR influences can be analyzed. The following sections describe the framework and its components in more detail. 3. Generation of wind production scenarios In order to study the ramps, a plausible method to generate wind power generation scenarios is required. The reader should note that the utilized methodology is only one possible option but it provides an effective way to generate aggregated wind power time series according to desired specifications through Monte Carlo simulations (e.g. the spatial distribution and installed capacity of the wind power generation). The method consists of a validated transformed autoregressive (AR) model with Cholesky decomposition (ARC) to generate wind speed time series and a turbine model to convert the wind speed time series to wind power time series. This approach allows us to simulate wind generation data for several wind farms located near each other. More formally, ARC stands for an approach where the simulated univariate autoregressive time

270

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

the off-peak period at night time. In addition to the heat demand, the sizing is dictated by several factors, such as the electricity tariff and the length of the off-peak period [27]. In this paper, a fully mixed storage model is utilized as a buffer between the electricity use and heat demand. The fully mixed storage assumes that all the water is perfectly mixed in the tank, and can thus be modeled with a unique temperature [28]: V ·  · Cp

Fig. 2. Flowchart of the wind generation modeling.

series have been transformed to multivariate cross-correlated time series using the Cholesky decomposition of the correlation matrix. The simulation of a wind generation scenario is divided into six consecutive steps, presented in the flowchart in Fig. 2. The approach used is similar to those introduced in [25,26]. The complete estimation and simulation process utilized to create the wind power generation time series is described in detail in [26]. 4. Modeling the electric heating The electric heating load types considered are direct electric space heating (DESH) and electric storage space heating (ESSH). The former utilizes thermostatically controlled radiators to maintain the indoor temperature within acceptable limits in each room of a building. Since the radiators directly heat the indoor air, the available flexibility is relatively limited and any performed control is typically followed by a rebound. In order to properly capture these dynamics, we exploit a two-capacity building thermal model that couples together the temperatures of the indoor air and building mass [19,23]. Mathematically, this is given by Ca

C

dT a dt

m dT

m

dt

= H ae · (T e − T a ) + H am · (T m − T a ) + H ag · (T g − T a ) +H as

· (T v

am

a

=H

− T a) +  m

· (T − T ) + H

me

e

m

· (T − T )

(1)

(2)

t: time; C a : thermal capacity of the indoor air; C m : thermal inertia of the building structures; T a : indoor air temperature; T e : external outdoor air temperatures; T m : building fabric and mass temperature; T g : ground temperature; T v : ventilation supply air temperature; Hae : sum of the infiltration heat capacity flow and window heat conductance; H am : sum of heat conductance in the solid wall and convection on the surface; H ag : floor heat conductance; H as : ventilation air heat capacity flow; H me : sum of heat conductance in the solid wall and convection on the surface; : total heating power (including solar radiation, heat gains, and electric power converted to heating power). Storage heating, on the other hand, conventionally stores the heat in a central water tank, from which the energy is further distributed around the building. This effectively decouples the electricity usage and heating, providing significant load shifting potential. The ESSH systems are typically designed to be able to store the greatest daily heat demand either partially or fully during

dT s = −U · A · (T s − T set ) + qs − qd dt

(3)

where V is the storage volume,  and Cp are the density and specific heat capacity of water, Ts is the water temperature, U is the coefficient of heat loss, A is the surface area of the storage, Tset is fixed indoor temperature, qs is the electric charging power heating the water, and qd is the space heating demand. The models represented by (1)–(3) are later integrated in a UC model to allow their utilization in a cost-minimizing manner. The formulation is given in Appendix A. One should note that the integration requires discretization of the continuous system, which is a source of error. Regardless of this error, the model can better capture the thermal dynamics of a building than a 1-capacity model, which is more commonly used in the literature [18,21]. The main difference between the model structures is the lack of a building mass temperature node in the 1-capacity model. Thus, we can get the 1-capacity model from the 2-capacity model by removing (2) and substituting Tm with Te and Ham with Hame = 1/(1/Ham + 1/Hme ) in (1). Note that the thermal capacitance, Ca , also needs to be reestimated. The models are briefly compared in Section 6.1. Both the above described models represent a single load, or alternatively, they describe the behavior of several identical loads. In practical systems, however, the load characteristic vary across a load population but modeling the dynamics separately in the UC is likely to make the problem intractable. Furthermore, a population of thermostatically controlled loads is prone to power oscillations if the characteristics of the loads are similar [29] or if their control is not properly designed [23]. Therefore, using the above models to represent the dynamics of a load population includes an assumption that the loads are centrally controlled, for example, by the strategy presented in [23] or [30]. Clustering techniques can be employed to reduce the number of loads with different model parameters, as suggested in [23], which allows modeling the population with a few representative load dynamics. 5. Simulation algorithm This section presents the simulation algorithm for the assessment of the electric heating load contribution to mitigate the effect of ramps. The section is divided into three parts. Firstly, the severe ramp events are defined and their detection is described. The second part shows the actual simulation algorithm, whereas the last part presents the performance measures of interest. 5.1. Ramp detection In this paper, the ramp severity is defined by its rate (MW/hour). Since the conventional power generation and DR need to respond not only to the wind production but also demand variations, it is more relevant to study net load (load minus available wind generation) ramps. From hourly resolution data, a straightforward way to detect ramps is to set a ramp rate threshold and compare it with the net load variations. For an up-ramp to be severe enough, the condition is lk − lk−k ≥ Rth

(4)

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

271

Fig. 3. Example of a two-hour up-ramp. The change in the net load exceeds the selected ramp rate threshold.

where lk is the net load, k is ramp length in hours, and Rth is the ramp rate threshold. A similar definition is used for the down-ramp detection. Fig. 3 illustrates a detected two-hour up-ramp that starts at time step k − 2 and ends at k. In this paper, the load profile is assumed to remain the same in each year, but with an appropriate load model, the year-to-year variations of the demand can also be included in the analysis. The question regarding which ramps can be considered to be severe is power system related and depends, for example, on the system’s generation mix and interconnections. General guidelines to define the ramp rate threshold are not treated in this paper. 5.2. The algorithm The wind power generation as well as fixed and flexible load are combined in a UC problem that seeks the cost optimal mix of generators to supply the demand for the coming day. It is argued that solving the UC problem is required to truly address the effect of ramps. The minimum on and off times of the generators and their ramping constraints are the main factors affecting the capability of a power system to follow the net load. The formulation for the UC is provided in Appendix A, and is based on the methodology presented in [31]. The UC problem is to minimize the objective function, which is a sum of electricity production cost, generator start-up cost, the cost of DR, and the penalty for unserved load. The objective is solved subject to the supply-demand balance, spinning reserve requirement, and the technical constraints of the generation units, which are generation limits, ramping constraints, and minimum up and down times. The aforementioned general UC formulation is further extended with the DR constraints. The reader should note that the UC model can be modified to better represent any given power system, while the employed formulation is a trade-off between computation time and the solution accuracy. As mentioned, the objective function of the UC considers the cost of utilizing the demand-side flexibility that is defined based on the assumed control strategy. For the centrally controlled load population, the cost of upward and downward flexibility can be given by a predetermined cost per megawatt hour [15]. This is to say, the customers participating in the DR program are provided with a monetary compensation if they are controlled. In this paper, we study the full potential of the heating load flexibility, i.e., assuming the flexibility to be free of charge. There also exist other approaches to consider the cost of utilizing the demand-side flexibility. For example, the cost of flexibility can be defined by price elasticities, i.e., the loads make decisions locally by reacting to the electricity price level [15]. However, these approaches are out of the scope of this paper.

Fig. 4. Flowchart of the simulation algorithm.

The simulation algorithm for the assessment of the electric heating load contribution to mitigate the effect of ramps is illustrated in Fig. 4. At the top of the flowchart, before the actual Monte Carlo simulation, a severe ramp is defined by setting the desired ramp rate threshold level. During the simulation, the program generates yearly wind generation profiles and uses them to form net load profiles. If a severe ramp is detected from the net load, the UC problem is solved. It is solved only for the actual day with the severe ramp and the preceding day, to find the initial state of the different generators. This is repeated with and without the flexible heating load. Since a large number of years need to be studied and solving the UC problem can be computationally slow, limiting the number of studied days is crucial. After solving the UC with and without the DR, the performance measures of interest are calculated and the parameter convergence is checked. 5.3. Measuring DR effect Several performance measures can be used to evaluate the flexible electric heating load contribution. Large net load variations may cause situations with excess electricity supply or alternatively there can be a deficit of energy, leading to loss of load. Thus, loss of load expectation (LOLE) is exploited to measure the system reliability [32], whereas expected wind curtailment (EWC) indicates the system’s ability to utilize wind generation during the extreme ramps. These indices can also be interpreted to measure the ability or benefit of the system to achieve supply-demand balance during an extreme ramp. The following gives the mathematical representation of the indices in the case of an up-ramp: 1  c EWC = up wj N N up

(5)

j=1

1  LOLE = up LOLj N N up

(6)

j=1

where wjc is the curtailed wind generation during ramp j and Nup is the total number of studied up-ramps. Loss of load, LOLj , gets a value of one if there is a deficit of electricity production during the jth ramp and zero otherwise. Regarding jth ramp, the definition

272

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

for a single severe ramp was earlier given in Section 5.1 and an illustrative example was provided in Fig. 3. The calculation of the performance measures presented in this section is based on the given definition. This is to say that only the values observed during the jth ramp are considered in the calculation and j ∈ {1, 2, . . ., Nup }. For example, in the case of the two-hour ramp in Fig. 3, the detected ramp includes only the values observed at time steps k − 2, k − 1, k. In the absence of DR, all the system flexibility to cope with extreme ramps needs to be provided by the generation-side. This is to say that a more flexible but costlier generation mix is typically committed or the generators are required to run in non-optimal operating points. Thus, the average benefit of demand response (AB) is employed to indicate the change in the electricity production cost. For up-ramps, it is defined by N NoDR − pc DR 1  pcj j N up pcjNoDR up

AB =

(7)

j=1

where pcjDR and pcjNoDR are the average costs of producing a one MWh of electricity during the jth ramp, with and without flexible heating load, respectively. The definition of AB is to say that it considers only the electricity production cost during a ramp and possible start-up costs are not included. Note that the production cost, pc, is only one possible money related figure that can be used. For example, many electricity markets utilize the marginal cost of generation as an electricity price [33], which could be an alternative to pc. However, the marginal cost is not solved in this paper. Finally, the electric heating load effect on the ramp rate is monitored. If the steepness is reduced, the system can be seen to have more ramping capacity to cope with possible prediction errors and even more extreme ramp events. The ramp rate reduction (RR) for up-ramp is N NoDR − RDR 1  Rj j RR = up N RjNoDR up

(8)

j=1

where RjDR is the net load ramp rate with the flexible heating load and RjNoDR without it. All the aforementioned performance measures are similarly defined for down-ramps. After termination, the simulation provides expected values for the aforementioned performance measures and their distributions. The expected values indicate the possible benefit of DR, whereas the distribution can be used to identify separate events (e.g. extreme cases) and analyze possible patterns (e.g. when curtailment occurs). 6. Case studies In this paper, case studies are conducted to demonstrate the contribution of heating load to mitigate the effect of ramps. The main goals are to identify the effect of the heating method and study how the heating load penetration affects the results. It is assumed that the activation of DR and wind generation are free of charge q (ck (q − qd ) = 0 and ckw (w) = 0 in (A.1)), i.e., the simulation provides an upper bound for the DR benefit and wind generation utilization. This section first presents the simulation setup, after which the results of the case studies are shown. 6.1. Simulation setup For the generation of the wind production scenarios, the correlation structure between 12 wind farms is calculated from 19 wind speed measurement locations in Finland. The monthly diurnal variations are estimated from two high altitude measurement locations in Finland. The Weibull parameters defining the wind conditions in

Fig. 5. Power curve used in converting wind speed to wind power. The dashed vertical lines indicate the cut-in (1 m/s) and cut-out (27 m/s) speeds.

Table 1 Wind power generation parameter values [26,35]. Parameter

Value

Unit

Weibull shape parameter k Weibull scale parameter  Turbine nominal power Turbine hub height Turbine yearly availability Turbine cut-in speed Turbine nominal speed Turbine cut-out speed

8.3 2.4 5 140 8322 1 12 27

– – MW m h m/s m/s m/s

each wind farm are identical for all locations for simplicity. The chosen Weibull parameters are averages of the parameters obtained from Finnish Wind Atlas [34] database for 20 operating or planned wind farms in Finland, and represent good wind speed conditions for wind power generation. To transform the wind speed time series to wind power time series, a piecewise power curve presented in Fig. 5 is utilized. The turbine specific parameters include yearly availability, the nominal power, cut-in, nominal and cut-out speeds. In this paper, Gamesa 5 MW G128 [35] turbines are considered for all wind farms. The parameter values for the wind generation model are presented in Table 1. The building model utilized in the simulations is as presented in Section 4 and the parameter values for the model are given in Table 2. The building model represents a 180 m2 medium massive structure two-floor single-family house, in which the thermal insulation is based on the Finnish 2010 building regulations [36,37]. The parameter estimation procedure and the visual validation of Table 2 Heating load model parameter values. The notation corresponds to the notation presented in Section 4 and Appendix A. Parameter DESH: Ca Cm Tg Tv Hae Ham Hag Has Hme Q¯ T¯ a T

a

Value 2344 20,183 7 18 52.33 933.3 28.8 87.43 60.48 6000 23 19

Unit kJ/K kJ/K ◦ C ◦ C W/K W/K W/K W/K W/K W ◦ C ◦ C

Parameter ESSH: V  Cp U A Q¯ s T¯ s Ts Common: Tset t

Value

Unit

1 1000 4.2 0.83 5.5 10,000 90 50

m3 kg/m3 kJ/kg,K W/m2 ,K m2 W ◦ C ◦ C

21 3600

◦

C s

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

Fig. 6. The indoor temperature solved using (1)–(2) with a discretization of one second and one hour. The temperature is compared with a corresponding 1-capacity building model. A thermal capacity value of 11,918 kJ/K is used for the 1-capacity building model. In the figure, t is the discretization time step.

the model are presented briefly in [19]. The effect of the discretization on the model accuracy is illustrated in Fig. 6, which shows the indoor temperature during a four-hour interruption in the heating (the outdoor temperature is set to the yearly mean, 6.5 ◦ C in the illustration). The same heating power sequence is fed to each model (the 2-capacity model with one hour and one second discretization, and the 1-capacity model with one second discretization) and the resulting indoor temperatures are compared in the figure. The solar radiation heating power for the model originates from test year weather data that has been developed to reflect the current climate conditions in southern Finland [38,39]. The outdoor temperature profile used by the model was measured in the Helsinki region in 2012. The heat gains are assumed to consist of the heating power produced by inside lighting, ventilation, and various household appliances. For the ESSH, a fully mixed storage model with losses is utilized (see (3)) and the heat demand drawn from the storage is solved by the described building model. The selection of storage size is based on the assumption that the storage is able to partially store the greatest daily heat demand [27]. The parameter values are shown in Table 2 and they correspond to a storage size of 47 kWh. Area A of the IEEE Reliability Test System (RTS-96) provides the power system model employed to demonstrate the DR effect during extreme ramp events. This system is selected, as it is welldocumented in the literature in terms of the technical constraints of its power plants. The power generation mix is given in [40], whereas the technical constraints and the quadratic cost curves can be found in [41]. The following modifications are made to the system: the hydro plants are omitted, nuclear plants are must-run units, and unit group U12 as well as U20 is aggregated to form a larger unit. Thus, the generation mix comprises 19 units with a total capacity of 3105 MW. In addition to the above modifications, the fixed electricity load corresponds to the hourly load of the Finnish system area in 2012 [42]. The peak load of this profile is scaled in such a way that it equals the total generation capacity minus the spinning reserve requirements (SR = 400 MW), leading to a peak load of 2705 MW. The amount of flexible heat demand is defined from this peak load by using the following scaling factor G = x · (P − SR)/ˆqd

(9)

where P is the conventional generation capacity (3105 MW), qˆ d is the peak heat demand (5.2 kW) solved from the building model, and x is the penetration of the flexible heating load. It is assumed that the fixed heating load, qd , is included into the electricity demand, d. The DR is defined as the deviation between the fixed and flexible heating load (q − qd ), which is also the flexibility utilized in the UC (see (A.2)). The fixed heating load is solved from the same building thermal model that is used in the UC problem. By using this definition, the heat demand with both heating methods can be kept equal, allowing us to compare the dynamics of DESH and ESSH. The

273

Fig. 7. Probability distribution of one, two, and three hour ramp rates. The vertical lines represent the up- and down-ramps of different lengths.

reader can also note that the simulation case without DR can be obtained by setting the penetration, x, to be zero. Installed wind generation capacity is set to 932 MW, corresponding to 30% of the conventional generation capacity, and the geographical spread of the farms is chosen so that steep ramps can be achieved. This is to say, the wind farms are located close to each other, resulting in high spatial correlations between the farms (i.e. between the wind speeds in each farm). This directly causes fast changes in wind speeds, and therefore in output power, to appear contemporaneously, which consequently results in more severe ramps. No conventional generation is replaced by wind generation, which makes the ramping capability of the power system slightly optimistic, as extra capacity is available and can be committed if a severe ramp occurs. The unit commitment is solved for a 24 hour period (K = 24) using one hour resolution. The penalty for unserved load is set high enough to keep loss of load uneconomical all the time. The quadratic cost curves are approximated with ten pieces. The simulations are run in Matlab and the UC problem is solved with Matlab’s built-in mixed-integer linear programming solver intlinprog.

6.2. Results Firstly, 3500 year-long wind production time series are generated, from which net load profiles are formed. From these profiles, one, two, and three hour up- and down-ramps that occur with a probability of 0.0023%, approximately once in 5 years, are the ramps that warrant studying. This corresponds to the first box, define threshold, of the flowchart in Fig. 4. Fig. 7 illustrates the probability distributions of one, two, and three hour net load ramps as well as the ramp rates that separate 0.0023% of the probability mass from the distribution tails. The corresponding ramp rates of the one, two, and three hour up-ramps in MW are 423, 617, and 768, respectively, while the down-ramps are 376, 511, and 673. As these rates separate the mass from the distribution tail, they defined the gentlest slope for the studied case. For comparison, if wind production is not installed in the system, the greatest observed rates in MW for the corresponding cases are 231, 346, and 417 for up-ramps, and 139, 217, and 282 for down-ramps. Thus, the wind production increases the expected ramp rates considerably. In this paper, we estimate the desired performance measures by simulating the severe ramps observed from the 3500 net load profiles. 3500 is basically the number of loops executed, see Fig. 4, so that approximately every five years the outer loop is executed, i.e., a severe ramp is detected and the UC is solved. The number is selected as a trade-off between simulation time and the convergence of different parameters. The ramps are studied with three different flexible load penetrations (3, 6, and 9%, as defined in Section 6.1) and with both heating methods. These penetrations equal the maximum power of 94, 187, and 281 MW in the case of DESH, and 156, 312, and 468 MW in the case of ESSH.

274

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

Fig. 8. Ramp rate reduction with each ramp length and penetration. (a) DESH, up-ramp, (b) ESSH, up-ramp, (c) DESH, down-ramp, and (d) ESSH, down-ramp.

The ramp rate reduction for each of the studied cases is presented in Fig. 8. Fig. 8a and b shows the results of DESH and ESSH for up-ramps, respectively, and Fig. 8c and b for down-ramps, with each of the penetrations and ramp lengths. As can be seen in the figure, DR clearly reduces the ramp rate in all of the cases. ESSH shows a greater reduction but in the case of the down-ramps, the difference is smaller. Heating load flexibility also seems to reduce the ramp steepness less during down-ramps, which is particularly true with ESSH. In such a case, the system rather reduces the conventional generation (ramps and shuts down plants) than exploits the DR flexibility to maintain the generation levels. As an extreme down-ramp typically results from increasing wind generation and decreasing demand in the late evening, the period of low net load may last several hours. In Fig. 8, when the penetration of flexible load increases, the reduction is greater but the differences between reduction increments decrease as a function of penetration. Furthermore, the figures reveal how the system reduces the ramp steepness most during one hour ramps, which suggests that the DR is the most valuable in balancing short variations. Regardless of the greater demand response potential of ESSH, the ramp length has a somewhat similar influence on the ramp rate reduction of both heating methods. The average benefit of flexible heating load is illustrated in Fig. 9. Again, the upper row (Fig. 9a and b) presents the results of DESH and ESSH for up-ramps, respectively, and Fig. 9c and d for down-ramps, for each penetration and ramp length. Considering the significant influence of DR on the ramp rate in Fig. 8, the obtained average benefit is relatively minor. This is understandable, as the UC problem finds a different optimal generation mix in each case, and is also able to curtail wind in order to allow more cost-efficient mixes during the ramps. Furthermore, the reader should note that a part

of the DR benefit comes from the start-up costs of the generators. The figure, and AB, considers only the cost of producing electricity during a ramp. The start-up costs are not included in AB but they are considered in the UC problem. From the figure, it is clear that the production costs decrease as a function of the DR penetration and that ESSH again provides greater reduction, particularly during up-ramps. In the case of down-ramps, the longer ramp lengths reduce the average benefit but such a clear trend is not visible in the case of up-ramps. The expected wind curtailment, EWC, is presented in Table 3. The curtailment is considerably greater during up-ramps than down-ramps. This is somewhat intuitive, as the ramp rate of an up-ramp can typically be reduced by curtailing wind at the beginning of a ramp when there is plenty of wind generation and the load is low. Furthermore, the utilized conventional generation mix has a greater ramp-down capability than ramp-up, as a result of which it is able to cope with the down-ramps more easily. From the table, it is also possible to observe how EWC decreases from one-hour up-

Table 3 EWC [MWh] for each ramp length, DR penetration, and ramp direction. Case

Up 1 h

Up 2 h

Up 3 h

Dn 1 h

Dn 2 h

Dn 3 h

No DR DESH: 3% 6% 9% ESSH: 3% 6% 9%

200

185

180

13

8

12

152 116 87

154 131 112

159 130 112

9 8 9

7 5 6

13 7 8

89 36 27

100 76 62

108 96 74

5 5 4

6 4 3

10 9 7

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

275

Fig. 9. Average benefit with each ramp length and penetration. (a) DESH, up-ramp, (b) ESSH, up-ramp, (c) DESH, down-ramp, and (d) ESSH, down-ramp.

ramps to three-hour without DR, while with the DR the behavior is the opposite. One explanation is the gentler ramp rate of longer ramps, which the system is able to follow even without the heating load flexibility. On the other hand, although the heat demand is low in summer and sometimes in spring and fall, there is some DR available, enough to assist the system with the short ramps but not with the longer ramps. Nevertheless, one should note that in practice other reasons affect the results as well. For example, with a great amount of heating load flexibility available, the system is able to commit more economical but less flexible units compared to the case without DR, which results in increased wind curtailment but reduced system operation cost. This explains the fact that doubling the DR may not halve the EWC. In other words, the system has more freedom to find an economical generation mix instead of tackling the wind generation variability. This is to say that the severity and exceptional nature of the steep ramps are reduced if we increase the system flexibility. As the loss of load is significantly penalized in the UC problem, the LOLE is not closely considered here. The highest observed value, 0.019, is obtained without the DR, whereas with the flexible heating load the values are only a fraction of this or zero. Finally, the probability distributions of the studied measures are illustrated and compared in Fig. 10. Only the case with one hour ramps and a DR penetration of 6% is shown. Fig. 10a and d reveals how the shape of the distributions of ramp rate reduction depends on the ramp direction. For up-ramps, the reduction is either relatively great or zero, whereas the reduction is more evenly distributed in the case of down-ramps. Furthermore, in Fig. 10a and d, the maximum observed values are clearly visible. In each of the distributions, zero values are strongly present. This is to say that DR is not required, or more likely, that there is no

flexibility available due to low or really high heat demand. In the case of wind curtailment (Fig. 10c and f), the utilization of flexibility reduces the curtailment, which conversely emphasizes the value zero in the distributions. The distributions of average benefit presented in Fig. 10b and e shows slightly more probability mass on the positive side, indicating benefit on average. However, there clearly exists variation in the obtained benefit, which is understandable as the average benefit index considers only the generation cost during a ramp. This is to say that the flexibility is also of value outside the one hour lasting ramp incidence. Furthermore, it should be noted that the studied system accommodates a great share of wind generation. The generation does not only produce steep ramps but it is also highly variable in general due to the spatial correlation, thus testing the system’s ability to follow the variations. Generally, in the case with a DR penetration of 6%, DESH was found to reduce the total operation cost of the system with 5% and ESSH with 4% compared to the case without the DR. These numbers also consider the start-up costs. The reader should note that ESSH increases electricity demand due to its thermal losses, which also increases the total operation cost. The presented results strongly indicate the beneficial influence of heating load flexibility during rare events such as the studied extreme ramps. The power system can clearly utilize the flexibility to reduce the net load ramp rate, while dispatching more economical generation with more wind power. Regardless of the shown benefits, the results cannot yet be used to analyze the total benefit from DR and to justify any investments in the control architecture required to unleash the flexibility. For such a justification, the market behavior should be analyzed over a long time period. The results, however, support this kind of analysis by separating

276

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

Fig. 10. Probability distributions of the performance measures in the case of one hour ramps and 6% DR penetration. (a) RR, up-ramp, (b) AB, up-ramp, (c) EWC, up-ramp, (d) RR, down-ramp, (e) AB, down-ramp, and (f) EWC, down-ramp.

and identifying the DR effect on the system reliability and market behavior, despite the fact that steep ramps only rarely occur. In the simulations, the optimization was solved only for a day with a detected severe ramp, while initial values for the daily optimization were obtained by solving the preceding day. This procedure may not necessarily capture the market behavior that would emerge when analyzed over a long time period. The solution methodology for initial values depends on the studied system, and is not profoundly considered in this study. One approach is to include a long enough look-ahead horizon for each daily optimization and simulate several days before a detected ramp. The presented simulations and results only aimed to demonstrate the contribution of heating load to mitigate the effect of rare and severe ramps. Therefore, any continuous market behavior over a long time period was not modeled. 7. Conclusion This paper provides a framework that allows studying the influence of flexible heating load on the mitigation of severe wind power ramp effects. The framework couples together a heating load model, a tool to generate wind production time series, and a unit commitment problem, which are further integrated into a Monte Carlo simulation. Furthermore, the study suggests several performance measures to assist the evaluation of heating load influence. The ultimate purpose of the Monte Carlo simulation is to provide expected values for the measures in such a way that they capture the different system conditions. As suggested by the case studies, the power system clearly utilizes the flexible heating load to reduce the ramp rates during severe ramps, as well as using it to reduce wind curtailment and operation cost. The results also show how the greatest benefit is derived in the case of short ramps. Naturally, increasing the DR penetration improves the results, but already minor penetration is

advantageous and further increment may not improve the outcome linearly. Furthermore, including a storage element in the heating system seems to have a beneficial effect on the results in many cases. For future work, we are interested in studying the effect of the spatial distribution of the Finnish wind generation on the emerging ramps and their influence on the operation of the Finnish power system. To mitigate the possible undesired effects, the role of electric heating will be essential.

Acknowledgements This work was supported by Aalto University as part of the project Smart Control Architectures for Smart Grids (SAGA) under the Aalto Energy Efficiency Research Program (AEF). The authors would like to thank Merkebu Zenebe Degefa and the SGEM project for providing the heat gain data that was used in the heat demand modeling. We would also like to thank the Finnish Meteorological Institute for providing the wind speed and outdoor temperature data.

Appendix A. Unit commitment The unit commitment problem with and without the flexible electric heating load is formulated in this appendix. When the problem is solved without DR, only the mixed-integer linear program provided in Appendix A.1 is considered. With DESH, the constraints presented in Appendix A.2 are added to the general formulation. Similarly, the constraints given in Appendix A.3 are added when ESSH load is studied. The general problem is based on the one presented in [31] and therefore, only the parts deviating from the original formulation are given. All the variables are non-negative if not stated otherwise.

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

A.1. General formulation

On the other hand, if the heating is on, the temperature is required to be limited, in order to preserve the comfort level inside the building. Thus, we set M = Q¯ , where Q¯  T¯ a and Q¯ is the rated heating power of the space heating system. The daily mean indoor temperature is also bounded, which prevents the optimization from keeping the temperature at the lower level over a day. This is given by

The objective is to minimize the total operation cost. I K  

min

p

s (ci,k (v, p) + ci,k (v) + ckw (w))

k=1 i=1

 K

+

(A.1)

(A.9)

k=1

where K is the length of the optimization period, i is the generation p unit index, I is the number of units, ci,k (v, p) represents electricity s production cost, ci,k (v) is the start-up cost of a unit, p is the electricity production, and v is a binary-variable indicating the on/off status of a unit. ckw (w) is the cost of dispatched wind generation, q

ck (q − qd ) represents the cost model of the flexible demand, PEd is the penalty for unserved load, whereas dkns is the unserved load. ¯ k − wkc , where w ¯ k is the The dispatched wind is defined as wk = w available and wkc the curtailed wind generation. The optimization is subject to the following constraints: Supply-demand balance:

The parameters in (A.5) and (A.6) are as follows: A1 = 1 +

t ae (H + H am + H ag + H as ), Ca

A3 =

t , Ca

B1 =

t am H , Cm

A4k =

A2 =

t am H , Ca

t ae e g (H · Tk + H ag · T g + H as · T v + qk + qsr ) k Ca (A.10)

B2 = 1 +

t am (H + H me ), Cm

B3k =

t me e H · Tk Cm (A.11)

g

pi,k + wk − wkc = dk + G · (qk − qdk ) − dkns

∀k

(A.2)

i=1

Spinning reserve requirement:



1 a Tk ≥ T set K K

q

(ck (q − qd ) + PE d · dkns )

k=1

I 

277

where qk and qsr are the heating power of the internal heat gains k and solar radiation, respectively, and t is the discretization time step. A.3. ESSH model

I

¯ k − wkc ≥ dk + G · (qk − qdk ) − dkns + SR p¯ i,k + w

∀k

(A.3)

S1 · T sk − T sk−1 = S2 · T set + S3 · (q k − q dk ) ∀k

i=1

Upper boundaries for the unserved load and the curtailed wind: dkns

The evolution of storage temperature is modelled with

≤

dk − G · qdk

∀k,

wkc

¯k ≤w

∀k

(A.4)

where dk is the electricity demand, qk is the flexible heating load, qdk is the fixed heating load, p¯ i,k indicates the maximum available power output of unit i, and SR is the system’s spinning reserve requirement. Other constraints are generation limit and ramping constraints, and minimum up and down times, for which formulations can be found in [31]. Furthermore, the same reference provides mathematical representations for the piecewise linear p approximation of the quadratic production cost, ci,k (v, p), and for s the start-up cost, ci,k (v). Shut-down cost is not considered in our study. A.2. DESH model The thermal dynamics of a building are obtained by discretizing the differential equations (1) and (2). The evolution of indoor temperature is defined by a A1 · Tka − Tk−1 − A2 · Tkm = A3 · qk + A4k

∀k

(A.5)

Building mass temperature: m −B1 · Tka + B2 · Tkm − Tk−1 = B3k

∀k

(A.6)

Allowed band for the indoor temperature variation: T a ≤ Tka ≤ T¯ a + M · uk

∀k

(A.7)

The upper boundary of heating power: qk ≤ M · (1 − uk ) ∀k a

(A.8)

where T and T¯ a are the lower and upper indoor temperature boundaries, respectively. The binary variable, uk , is needed to relax the upper boundary of the indoor temperature if the temperature level starts to rise due to the internal heat gains and solar radiation.

(A.12)

where water temperature and charging power are bounded by T s ≤ T sk ≤ T¯

s

∀k,

q k ≤ Q¯ s

∀k

(A.13)

and the parameters of the storage model are S1 = 1 + S2,

S2 = U · A · S3,

S3 =

t V ·  · Cp

(A.14)

where T s and T¯ s are the lower and upper storage temperature boundaries, respectively, and Q¯ s is the rated charging power of the thermal storage. References [1] H. Holttinen, et al., Methodologies to determine operating reserves due to increased wind power, IEEE Trans. Sustain. Energy 3 (4) (2012) 713–723. [2] M. Ortega-Vazquez, D. Kirschen, Estimating the spinning reserve requirements in systems with significant wind power generation penetration, IEEE Trans. Power Syst. 24 (1) (2009) 114–124. [3] J. Morales, A. Conejo, J. Perez-Ruiz, Economic valuation of reserves in power systems with high penetration of wind power, IEEE Trans. Power Syst. 24 (2) (2009) 900–910. [4] J. Kiviluoma, et al., Short-term energy balancing with increasing levels of wind energy, IEEE Trans. Sustain. Energy 3 (4) (2012) 769–776. [5] T. Aigner, S. Jaehnert, G. Doorman, T. Gjengedal, The effect of large-scale wind power on system balancing in Northern Europe, IEEE Trans. Sustain. Energy 3 (4) (2012) 751–759. [6] M. Milligan, et al., Operational analysis and methods for wind integration studies, IEEE Trans. Sustain. Energy 3 (4) (2012) 612–619. [7] P. Wang, Z. Gao, L. Bertling, Operational adequacy studies of power systems with wind farms and energy storages, IEEE Trans. Power Syst. 27 (4) (2012) 2377–2384. [8] R. Sevlian, R. Rajagopal, Detection and statistics of wind power ramps, IEEE Trans. Power Syst. 28 (4) (2013) 3610–3620. [9] D. Ganger, J. Zhang, V. Vittal, Statistical characterization of wind power ramps via extreme value analysis, IEEE Trans. Power Syst. 29 (6) (2014) 3118–3119. [10] N. Navid, G. Rosenwald, Market solutions for managing ramp flexibility with high penetration of renewable resource, IEEE Trans. Sustain. Energy 3 (4) (2012) 784–790. [11] S. Wang, D. Yu, J. Yu, A coordinated dispatching strategy for wind power rapid ramp events in power systems with high wind power penetration, Int. J. Electr. Power Energy Syst. 64 (0) (2015) 986–995.

278

A. Alahäivälä et al. / Electric Power Systems Research 142 (2017) 268–278

[12] D. Callaway, I. Hiskens, Achieving controllability of electric loads, Proc. IEEE 99 (1) (2011) 184–199. [13] A. Yousefi, H.H.C. Iu, T. Fernando, H. Trinh, An approach for wind power integration using demand side resources, IEEE Trans. Sustain. Energy 4 (4) (2013) 917–924. [14] C. De Jonghe, B. Hobbs, R. Belmans, Value of price responsive load for wind integration in unit commitment, IEEE Trans. Power Syst. 29 (2) (2014) 675–685. [15] K. Dietrich, J. Latorre, L. Olmos, A. Ramos, Demand response in an isolated system with high wind integration, IEEE Trans. Power Syst. 27 (1) (2012) 20–29. [16] A. Papavasiliou, S. Oren, Large-scale integration of deferrable demand and renewable energy sources, IEEE Trans. Power Syst. 29 (1) (2014) 489–499. [17] B. Dupont, K. Dietrich, C.D. Jonghe, A. Ramos, R. Belmans, Impact of residential demand response on power system operation: a Belgian case study, Appl. Energy 122 (0) (2014) 1–10. [18] G. Papaefthymiou, B. Hasche, C. Nabe, Potential of heat pumps for demand side management and wind power integration in the German electricity market, IEEE Trans. Sustain. Energy 3 (4) (2012) 636–642. [19] M. Ali, A. Safdarian, M. Lehtonen, Demand response potential of residential HVAC loads considering users preferences, in: IEEE PES Innov Smart Grid Technol Eur, 2014, pp. 1–6. [20] L. Hughes, Meeting residential space heating demand with wind-generated electricity, Renew. Energy 35 (8) (2010) 1765–1772. [21] M. Ali, J. Jokisalo, K. Siren, M. Lehtonen, Combining the demand response of direct electric space heating and partial thermal storage using LP optimization, Elect. Power Syst. Res. 106 (0) (2014) 160–167. [22] S. Pourmousavi, S. Patrick, M. Nehrir, Real-time demand response through aggregate electric water heaters for load shifting and balancing wind generation, IEEE Trans. Smart Grid 5 (2) (2014) 769–778. [23] W. Zhang, J. Lian, C.Y. Chang, K. Kalsi, Aggregated modeling and control of air conditioning loads for demand response, IEEE Trans. Power Syst. 28 (4) (2013) 4655–4664. [24] N. Ruiz, I. Cobelo, J. Oyarzabal, A direct load control model for virtual power plant management, IEEE Trans. Power Syst. 24 (2) (2009) 959–966. [25] M. Koivisto, et al., Statistical analysis of large scale wind power generation using Monte Carlo simulations, in: Power Syst Comput Conf (PSCC), 2014, pp. 1–7. [26] J. Ekström, et al., Assessment of large scale wind power generation with new generation locations without measurement data, Renew. Energy 83 (2015) 362–374.

[27] R. Kara, The Handbook of Electrical Heating, The Electrical Contractors’ Association of Finland, Espoo, Finland, 1994. [28] A.C. Celador, M. Odriozola, J. Sala, Implications of the modelling of stratified hot water storage tanks in the simulation of CHP plants, Energy Convers. Manag. 52 (8-9) (2011) 3018–3026. [29] M. Nehrir, R. Jia, D. Pierre, D. Hammerstrom, Power management of aggregate electric water heater loads by voltage control, in: IEEE Power Eng Soc Gen Meet, 2007, pp. 492–497. [30] N. Lu, Y. Zhang, Design considerations of a centralized load controller using thermostatically controlled appliances for continuous regulation reserves, IEEE Trans. Smart Grid 4 (2) (2013) 914–921. [31] M. Carrión, J. Arroyo, A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem, IEEE Trans. Power Syst. 21 (3) (2006) 1371–1378. [32] R. Billinton, R.N. Allan, Reliability Evaluation of Power Systems, Pitman, Boston, MA, USA, 1984. [33] E. Ela, M. Milligan, A. Bloom, A. Botterud, A. Townsend, T. Levin, Evolution of wholesale electricity market design with increasing levels of renewable generation, 2014, Available at: www.nrel.gov/publications. [34] B. Tammelin, et al., Production of the Finnish Wind Atlas, Wind Energy 16 (1) (2013) 19–35. [35] Gamesa, Gamesa G128 5 MW wind turbine specifications, 2014, Available: http://www.gamesacorp.com/recursos/doc/productos-servicios/ aerogeneradores/catalogo-g10x-45mw-eng.pdf. [36] Ministry of the Environment, C3: Finnish code of building regulations. (Thermal insulation in a building) [In Finnish], 2010. [37] Ministry of the Environment, D3: Finnish code of building regulations. (Energy efficiency of buildings) [In Finnish], 2012. [38] Finnish Meteorological Institute, Energy reference year in current climate, 2012, Available: http://ilmatieteenlaitos.fi/energialaskennan-testivuodetnyky [In Finnish]. [39] T. Kalamees, et al., Development of weighting factors for climate variables for selecting the energy reference year according to the EN ISO 15927-4 standard, Energy Build. 47 (0) (2012) 53–60. [40] C. Grigg, et al., The IEEE reliability test system-1996. a report prepared by the reliability test system task force of the application of probability methods subcommittee, IEEE Trans. Power Syst. 14 (3) (1999) 1010–1020. [41] C. Wang, S. Shahidehpour, Effects of ramp-rate limits on unit commitment and economic dispatch, IEEE Trans. Power Syst. 8 (3) (1993) 1341–1350. [42] Fingrid Oyj, Load and generation, 2015, Available: http://www.fingrid.fi/en/ electricity-market/load-and-generation/Pages/default.aspx.