Optimal Allocation of Spinning Reserves in Interconnected Energy Systems with Demand Response Using a Bivariate Wind Prediction Model

Article

Optimal Allocation of Spinning Reserves in Interconnected Energy Systems with Demand

Response Using a Bivariate Wind Prediction Model

Yerzhigit Bapin^1,* , Mehdi Bagheri^1,2 and Vasilios Zarikas¹

1 School of Engineering and Digital Sciences, Nazarbayev University, 53 Kabanbay Batyr ave,

010000 Nur-Sultan, Kazakhstan; [email protected] (M.B.); [email protected] (V.Z.)

2 National Laboratory Astana, Center for Energy and Advanced Material Science, Nazarbayev University, 010000 Nur-Sultan, Kazakhstan

* Correspondence: [email protected]; Tel.:+7-701-444-8411

Received: 23 July 2019; Accepted: 3 September 2019; Published: 9 October 2019 Abstract:The proposed study presents a novel probabilistic method for optimal allocation of spinning reserves taking into consideration load, wind and solar forecast errors, inter-zonal spinning reserve trading, and demand response provided by consumers as a single framework. The model considers the system contingencies due to random generator outages as well as the uncertainties caused by load and renewable energy forecast errors. The study utilizes a novel approach to model wind speed and its direction using the bivariate parametric model. The proposed model is applied to the IEEE two-area reliability test system (RTS) to analyze the influence of inter-zonal power generation and demand response (DR) on the adequacy and economic efficiency of energy systems. In addition, the study analyzed the effect of the bivariate wind prediction model on obtained results. The results demonstrate that the presence of inter-zonal capacity in ancillary service markets reduce the total expected energy not supplied (EENS) by 81.66%, while inclusion of a DR program results in an additional 1.76% reduction of EENS. Finally, the proposed bivariate wind prediction model showed a 0.27% reduction in spinning reserve requirements, compared to the univariate model.

Keywords:bivariate probability density function; demand response; equivalent assisting unit method;

interconnected power systems; spinning reserves; renewable energy

1. Introduction

Nowadays, the world faces serious challenges, such as global warming, air pollution, and shortage of natural resources, that must be addressed in order to ensure sustainable development of the human race. Since renewable energy sources do not pollute the atmosphere, nor do they depend on the availability of fossil fuels, their usage may help to alleviate (or even eliminate) these problems.

Within the last 10 years, the worldwide renewable capacity grew from 1.061 TW to 2.179 TW [1].

According to the projections made by the International Energy Agency (IEA), renewable power will amount to 40% of the total energy production by 2040 [2]. Nonetheless, in order to achieve the renewable energy goals successfully, many serious problems must be resolved in the near future.

Due to the unpredictable and variable nature of renewables, integration of the high shares of renewable power will require having a great deal of flexibility from the energy systems.

The balance between supply and demand in the grid is traditionally achieved by adjusting the power output of conventional generating capacity or by regulating the power consumption on the demand side. The demand side regulation, or demand response (DR), can be obtained through establishment of economic incentives and installment of proper infrastructure [3]. Nowadays,

Energies2019,12, 3816; doi:10.3390/en12203816 www.mdpi.com/journal/energies

various types of DR programs are utilized by independent system operators (ISOs) to offer incentives stimulating reduction in electricity consumption during contingency periods or when electricity prices are high. In this regard, the utilization of DR programs can benefit both the electricity consumers and ISOs.

Proper introduction of renewable power and DR into the energy market and grid system infrastructure demands revision of traditional operational approaches with considerable emphasis on reliability and adequacy [4]. Conventional adequacy assessment methods used by the grid system operators can be divided into deterministic and probabilistic approaches. The deterministic approach sets the power system reliability criteria such that the system is capable of withstanding a single contingency (N-1) or a k number of contingencies happening at the same time (N-k).

Purely deterministic security criteria do not consider random processes taking place during the grid operation. Implementation of deterministic security criteria is relatively simple and straightforward;

thus, it is quite popular among ISOs all over the world [5]. As opposed to deterministic criteria, the probabilistic security criteria are more complex, requiring thorough statistical data, such as availability of power generating units and transmission lines, load and renewable forecast errors, etc.

The advantage of probabilistic security criteria in comparison to deterministic criteria lies in its ability to quantify the likelihood and significance of existing uncertainties relative to the adequacy of a power system [6].

The probabilistic power system adequacy assessment is gaining more attention with the growth of renewable generation. Various adequacy assessment studies that take into account renewable power have been presented in the last few years [7–13]. Most of these studies consider random events either by imposing an upper bound to reliability indices determining the loss of load or loss of energy expectation, or by adding a financial penalty into the cost function. For example, reference [9]

proposes a hybrid reliability criterion by setting an upper limit to the loss of load probability (LOLP) and expected load not served (ELNS). The authors calculate LOLP and ELNS by considering only the single- and double-outage events. Although setting a target value to risk indices can lead to reduction in computational intensity [10], models with bounded reliability metrics may fail to find a globally optimal solution. References [11] and [12] propose methodologies, where estimation of spinning reserve requirements is conducted using cost–benefit analysis. The methodologies find the optimal spinning reserve requirements by minimizing the total system costs, i.e., the unit operating expenses and the cost of load curtailment. However, both studies model the energy systems as an isolated entity; thus, cannot be applied to the multi-area power grids. A similar approach based on the cost–benefit analysis has been proposed in [13]. The methodology considers uncertainties due to the load and wind forecast errors and employs the equivalent assisting unit method to account for inter-zonal power flows, yet this approach considers only the wind power as a source of renewable energy. In reference [14], the authors propose a method to calculate the spinning reserves in systems with renewable energy sources. The computational efficiency of the method is achieved by application of the cross entropy (CE) concept. The main idea of CE is to modify the outage replacement rates (ORR) of the system components based on their importance in terms of the system reliability. The uncertainty associated with renewable generation is taken into account by a multi-state short-term Markov model and a capacity factor. The Markov model captures random outages of renewable generators, whereas the capacity factor captures intermittency of a primary energy source. The method mainly focuses on reliability aspects rather than spinning reserve optimization. Reference [15] presents a method for quantification of spinning reserves in interconnected power systems based on the cost–benefit analysis.

The proposed methodology utilizes the radial-equivalent-independent approach to reduce the total number of buses in the grid system. The optimization of spinning reserve requirements is conducted using either security constrained unit commitment (SCUC) or security constrained economic dispatch (SCED). Although the proposed methodology introduces an effective way to quantify the spinning reserves in large interconnected systems, it lacks the ability to account for renewable and load forecast errors. Reference [16] proposes a methodology for energy and spinning reserve market clearing

(ESRMC) involving conventional units and wind power generators. The multi-objective optimization problem is solved by a genetic algorithm with the goal to minimize the total system operating cost and risk level. The methodology uses parametric assumptions to model the uncertainties associated with load and wind power, yet it is unclear how these uncertainties were incorporated into the optimization problem. Reference [17] proposes a model for co-optimized energy and spinning reserve market, including the provision of demand response from industrial and commercial customers. The proposed model utilizes a standard unit commitment to optimize the energy generation schedules, but the reserve allocation is carried out based on the deterministic security criterion equal to 10% of the net demand.

The models described above provide a more adequate estimation of operating reserve requirements compared to deterministic reliability and adequacy assessment methods. However, due to the emergence of new technologies and market mechanisms, there is a need for development of reserve estimation approaches capable of considering these novelties.

Although DR services are relatively new in the electric power industry, substantial research works have been carried out to incorporate DR into different market structures. Reference [18] describes the demand response reserve project (DRRP) developed by the New England independent system operator (ISO-NE). The goal of the DRRP is to integrate DR into the ancillary service market by providing operating reserve and forward reserve services. However, this study does not explain how the system reliability is evaluated by the ISO-NE during the implementation of the DRRP. In reference [19], a methodology for enhancing the spinning reserve capacity by application of the direct load control (DLC) program is proposed. The DLC is a load control mechanism where the contracted customers are disconnected from the grid by the grid utility in contingency period. The methodology has been tested on the IEEE-14 bus system to find out how implementation of the DLC program affected the amount of energy not supplied. However, the study examines only 3 out of 10 possible contingency scenarios and neglects simultaneous contingencies. The DLC participation scenarios were set to be equal to 15%, 30%, and 50%. Reference [20] present a probabilistic model considering provision of DR reserves in the ancillary service demand response (ASDR) market. The proposed model considers only uncertainties due to failure of conventional capacity, transmission line outages, and load forecast error, and neglects the uncertainty caused by renewable generation. The optimization of SCUC problem is accomplished in two stages using the mixed-integer linear program (MILP). The model was tested separately on the IEEE one-area reliability test system (RTS) and on the customized six-bus grid system, yet both systems were assumed to have no interconnection between each other, or any other system whatsoever.

A similar approach has been presented in reference [21]; however, this model neglects to consider the uncertainty due to load forecast error, and the contingencies are presented only in terms of random outages of generating units. Reference [22] presents a security-constrained joint forward energy and ancillary service market clearing mechanism, which incorporates DR reserves; the method utilizes SCUC formulation developed in their previous studies [23,24]. In this SCUC formulation, the reliability criteria require procurement of spinning reserves provided by conventional capacity and demand side reserves to withstand a pre-determined contingency. The method is applicable to isolated systems and does not consider uncertainties due to load and renewable forecast errors. In [25], the authors introduce a complete DR model in which demand side flexibility is exploited in energy and reserve markets.

The DR model is integrated into the spinning reserve estimation model described in reference [12].

The optimal spinning reserve, including DR reserve, is determined using a two-stage SCUC in which the market clearing of energy production and reserve service is done simultaneously. However, the methodology neglects multiple-order contingencies in order to reduce computational complexity.

Reference [26] presents a DR model integrated in conventional SCUC with deterministic security criteria, yet the model assumes constant demand elasticity and linear correlation between electricity price and demand. The model is tested on the IEEE one-area RTS to determine reduction in total system costs. In reference [27], the authors incorporate incentive-based and time-based DR programs into a stochastic optimization scheduling model, considering only generating unit failure as stochastic input. In this study, the electricity demand and prices are modeled in the same fashion as it was

performed in reference [26] (i.e., inelastic demand and linear price-demand correlation). Reference [28]

proposes a model incorporating the incentive-based DR program with the dynamic economic dispatch (DED) problem. The study develops non-linear models of responsive load for the incentive-based DR programs. The combined DED and DR model is solved by the application of the random drift particle swarm optimization (DPSO) algorithm. However, the study utilized conventional N-1 security criterion to determine spinning reserve requirements. Reference [29] proposed a reliability assessment based on the truncated state space (TSS) method combined with demand response. The truncation of the state space is done to reduce the computational complexity of the problem and is defined as elimination of system states with low probability. The model analyzes the power systems at generation and transmission levels and was implemented on the IEEE one-area RTS. Although the objective of this approach was to evaluate reliability of the large and complex power grids in a fast and efficient manner, it can only be applied to isolated systems due to the lack of the ways to account for interconnection power flows.

The main contribution of this study is to present a day-ahead spinning reserve estimation model for interconnected power systems based on a two-stage probabilistic SCUC. The model takes into consideration risks associated with capacity deficit due to random generator outage and renewable forecast error, spinning reserve provided by conventional generators (CG), as well as up-spinning reserve provided by demand response providers (DRPs) and interconnected capacity. Renewable power production is modeled using parametric univariate and bivariate models. The provision of spinning reserve and DR reserve is carried out by ISO through ancillary service markets. The influence of interconnected capacity on the adequacy of the test system is accounted for via the equivalent assisting unit approach.

The rest of this paper is organized as follows: Section2describes the proposed market structure and characterizes the key components of the model. Section3describes the case study implemented on the IEEE two-area RTS, as well as presenting the numerical results and discussion. Section4presents the conclusion based on the results of this study.

2. Model Description

2.1. Market Structure

The procurement of energy and ancillary services from the market participants is executed by ISO on a day-ahead basis. The flowchart of the multi-period market clearing process is depicted in Figure1.

Reference [28] proposes a model incorporating the incentive-based DR program with the dynamic economic dispatch (DED) problem. The study develops non-linear models of responsive load for the incentive-based DR programs. The combined DED and DR model is solved by the application of the random drift particle swarm optimization (DPSO) algorithm. However, the study utilized conventional N-1 security criterion to determine spinning reserve requirements. Reference [29]

proposed a reliability assessment based on the truncated state space (TSS) method combined with demand response. The truncation of the state space is done to reduce the computational complexity of the problem and is defined as elimination of system states with low probability. The model analyzes the power systems at generation and transmission levels and was implemented on the IEEE one-area RTS. Although the objective of this approach was to evaluate reliability of the large and complex power grids in a fast and efficient manner, it can only be applied to isolated systems due to the lack of the ways to account for interconnection power flows.

The main contribution of this study is to present a day-ahead spinning reserve estimation model for interconnected power systems based on a two-stage probabilistic SCUC. The model takes into consideration risks associated with capacity deficit due to random generator outage and renewable forecast error, spinning reserve provided by conventional generators (CG), as well as up-spinning reserve provided by demand response providers (DRPs) and interconnected capacity.

Renewable power production is modeled using parametric univariate and bivariate models. The provision of spinning reserve and DR reserve is carried out by ISO through ancillary service markets. The influence of interconnected capacity on the adequacy of the test system is accounted for via the equivalent assisting unit approach.

The rest of this paper is organized as follows: Section 2 describes the proposed market structure and characterizes the key components of the model. Section 3 describes the case study implemented on the IEEE two-area RTS, as well as presenting the numerical results and discussion. Section 4 presents the conclusion based on the results of this study.

2. Model Description

2.1. Market Structure

The procurement of energy and ancillary services from the market participants is executed by ISO on a day-ahead basis. The flowchart of the multi-period market clearing process is depicted in Figure 1.

Figure 1. Flowchart of the electricity market clearing process.

The entire process begins with submission of the day-ahead generation forecasts by renewable energy producers (REP) and load forecasts by electricity distribution companies (DisCo), and

Load and Renewable Forecasts

Net Demand Calculation

Start of Trading

Submission of Bids

Evaluation of Bids

Winning Bids Committed to Energy Production and Reserve Provision

Figure 1.Flowchart of the electricity market clearing process.

The entire process begins with submission of the day-ahead generation forecasts by renewable energy producers (REP) and load forecasts by electricity distribution companies (DisCo), and wholesale industrial consumers (WIC) to ISO. In this market model, we assume that REPs are not involved in competitive bidding; instead, to enhance their deployment, renewables are supported by incentives, such as priority dispatch and feed-in tariffs.

At the second stage, ISO initiates trading by submitting the aggregated net demand to the trading platform. The market structure considered in our model includes energy and ancillary service procurement. During the third stage, the market participants submit their offers. The proposed model takes into consideration inter-zonal power trading; therefore, interconnected conventional generators (ICG) can also participate in the ancillary service market. According to the proposed market structure, conventional generators (CG) can bid in energy and ancillary service markets, while DRPs and ICGs can only provide up-spinning reserve. It is important to note that, although, the market designs with both inter-zonal energy and reserve trading is more common, the proposed model focuses only on inter-zonal spinning reserve procurement, since the described market structure is based on the energy market of Kazakhstan with minor modifications applied to the optimization of the energy and ancillary service market. In Kazakhstan, the electricity demand is fully covered by domestic power generating units, whereas the cross-zonal capacity is used for balancing, and for the short- and mid-term reserve purposes.

At the fourth stage of the market clearing process, the winning bids are committed for energy production and ancillary services using two-stage stochastic SCUC. In this model, the energy production and ancillary service are considered as competitive commodities; therefore, the market clearing of these commodities is exercised simultaneously. The main objective of unit commitment and economic dispatch optimization is to minimize the system operating costs.

2.2. Generation System

In this study, the generation system is modeled using the capacity outage probability table (COPT).

COPT is created using the Markovian representation of the generation system states. In a two-state Markov process, a power generating unit can be either up and generating the required amount of power, or down due to a technical issue.

The probability of capacity deficit due to the unit failure, also known as the outage replacement rate, is dependent on the rate at which a unit transitions from operating state to a failure state and the rate at which the unit is repaired or replaced. Assuming the exponential distribution of the time to failure of each unit and neglecting the repair process, the probability of finding unition outage is given by [5]:

U_i=1−e^−γiT (1)

whereγiis the failure rate of unitiandTis the lead time. Note that throughout this paper, indicesi,j, andkrefer to the intra-zonal conventional units, inter-zonal conventional units, and DRPs, respectively.

Indicesm,s, andtrefer to the generation system states, net demand scenarios, and the time-periods of the scheduling horizon. Finally, capitalPrepresents the power output of conventional and renewable energy producers, whereas the lowercaseprefers to probability.

The calculation of generation system states is carried out using the recursive technique described in reference [5] and includes information on available capacity and corresponding probabilities of the system states. The recursive equation is given by reference [5]:

Pm,i(X) = (1−Ui)pm−1+ (Ui)pm−1(X−Pi) (2) wherepm−1(X)andpm,i(X)are, respectively, the cumulative probabilities of generation system statem before and after unitiis added, andPiis the installed capacity of uniti. Note that the total number of possible system states is equal to 2ⁿ, where n is the total number of units included in the generation system (i.e., the maximum value ofi).

2.3. Net Demand

Alongside the generator failure, the uncertainty due to load and renewable forecast error may significantly affect the reliability of a grid system [30,31]. In this study, only wind and solar photovoltaic power is considered as renewable energy. The incentives aimed to enhance the deployment of renewable energy sources described in the previous section, as well as the uncorrelated nature of load and renewable forecast errors, allow one to consider renewable energy as negative load [9]; therefore, for any given time-periodt, the net demand is given by:

D^t=L^t−P^t_W−P^t_PV (3)

where L^t, P^t_W, and P^t_PV are the system load, wind, and solar energy production at time-period t, respectively.

2.3.1. Load Forecast Error

Load forecasting has been evolving since the outset of electricity, and to this day, great abundance of sophisticated techniques has been developed for accurate prediction of the consumer load. For this reason, as well as the repeating nature of the load, it is reasonable to assume that the standard deviation of the forecast error equates to a portion of the actual load [9]:

σ^t_L= ^Y

100L^t_F (4)

whereσ^t_LandL^t_Fare the standard deviation of the load forecast error and forecasted load at timet, respectively, andYis the function depending on the accuracy of the forecasting framework.

2.3.2. Wind Power Output

A relatively straightforward and effective way to model renewable power output is to utilize a parametric probability density function (PDF) of choice. To date, a large number of studies have been conducted to determine the PDFs that provide for accurate renewable power generation prediction models. The variety of distributions used to model wind and solar forecast error range from normal distribution [9–11], the Weibull distribution [32–36], mixed distribution based on Laplace and normal distributions [37], the Beta distribution [38], the hyperbolic distribution [39], and the Levyα-stable distribution [40].

Commonly, wind speed and wind forecast error are modeled using Weibull and normal distributions, respectively, described by marginal (univariate) probability density functions. In our work, we perform an analysis using a traditional approach, which utilizes a Weibull distribution, and for the first time in the power system reliability assessment, bivariate probability density function, which represents joint distribution of wind speed and direction.

Univariate Wind Prediction Model

The univariate parametric wind speed PDF is given by the following equation [41]:

fw(v) = (_k/γ)(_v/γ)^k−1exp[⁻(_v/γ)^k] (5) wherevis the actual wind velocity, andγandkare the scale and shape factors of the PDF. The calculation of wind turbine power output should be accomplished based on the wind turbine power curve provided

by the manufacturer. In this study, the wind turbine power output calculation is conducted using the following relation [42]:

PW=



















0 v<v_min P_f(v) v_min≤v<vr

Pr vr ≤v<vmax

0 v≥vmax



















(6)

wherePr is the maximum power output of a wind turbine rated at wind velocityvr,v_min,andvmax

are the cut-in and cut-offwind velocities, respectively. Finally,Pf(_v)is the non-linear part of the power curve, representing the power output at wind velocityv. In this study, we utilize the power coefficient-based model given by reference [42]:

P_f(v) = ¹

2ρair×Awt×Ceq×v³ (7)

whereρairis the air density, considered constant (1.225 kg/m³),Awtis the wind turbine rotor swept area, andCeqis the dimensionless power coefficient equivalent, assumed to be equal to 0.4 [41].

Bivariate Wind Prediction Model

The univariate wind prediction model lacks the ability to consider other parameters that may have a significant impact on the prediction results. This issue can be addressed through utilization of the bivariate models consisting of a mixture of two marginal probability distributions of the parameters of interest. Reference [43] suggests that bivariate models show superior results as compared to the univariate counterparts in terms of goodness-of-fit criteria. In this study, we propose a bivariate model that considers the wind speed and wind direction through utilization of the Johnson–Wehrly bivariate PDF with a mixture of two Weibull distribution functions. The Johnson–Wehrly bivariate PDF is expressed as follows:

fV,θ(v,θ) =2πf_ψ(ψ)fV(v)f_θ(θ), 0<v<vmax, 0≤θ≤2π (8) whereVandθare the random variables representing wind speed and direction, respectively, andψis the random variable representing the relationship structure betweenVandθ, which is given by the following expression:

ψ=2π[FV(v)⁻F_θ(θ)], 0≤θ≤2π (9) where the termsFV(v),F_θ(_θ), andF_ψ(_ψ)are, respectively, the cumulative distribution functions ofV, θ, andψ.

Equation (7) for the bivariate case becomes:

Pf(v) = ¹

2ρair·Awt·Ceq·v³·fV,θ(v,θ). (10) It should be noted that in this study, the wind velocity is expressed in m/s, while wind direction is expressed in terms of the following cardinal directions: North (N), South (S), East (E), West (W), Northeast (NE), Northwest (NW), Southeast (SE), and Southwest (SW).

2.3.3. Solar Power Output

The rate at which a PV module generates electricity is affected by the intensity of the solar radiation, PV module temperature, and technical characteristics of the module [44]. Similar to wind power production, solar generation can be modeled using a parametric probability density function. In this

study, we employ the Beta distribution that has recently been utilized by other research groups [32,33].

The Beta PDF for this case is given by the following expression [32]:

f_b(si) =











Γ(α+β)

Γ(α)Γ(β)×si^(α−1)×(1−si)^(β−1) f or0≤si≤1,α≥0,β≥0

0 otherwise











(11)

where f_b(si) represents the probability density function of solar irradiance and Γ represents the following Gamma function [33]:

∞

Z

0

ρ^λ−1e^−ρdρ (12)

whereρis an integration variable.

The shape parameters of the Beta distribution are denoted byαandβand are given by the following formulas [32]:

β = (1 − µ) ^µ(1 + µ) σ² − 1

!

(13) α = ^{µ β}

1 − µ (14)

whereµandσare the mean and standard deviation of solar irradiance. Given the solar irradiance, total area of PV modules—Apv, the PV module conversion efficiency—ηpv, the maximum power point tracking efficiency—ηmppt, and the solar irradiance angel—θs, the total power output of a PV farm can be calculated as follows [45]:

PPV=si·APV·ηPV·ηmppt·cosθs. (15) The resulting net demand dataset must be discretized into an odd number of equal intervals for further calculations. The discretization is performed by dividing the net demand probability distribution into equal intervals as indicated in Figure2. These intervals are considered as scenarios with individual probabilities corresponding to the mid-point of each interval.

Energies 2019, 12, x FOR PEER REVIEW 8 of 25

2

(1 )

(1 )   1

 



 + 

= −   −  

1

  

= 

−

where µ and σ are the mean and standard deviation of solar irradiance. Given the solar irradiance, total area of PV modules—A^pv, the PV module conversion efficiency—η^pv, the maximum power point tracking efficiency—η^mppt, and the solar irradiance angel—θ^s, the total power output of a PV farm can be calculated as follows [45]:

PV PV PV mppt cos s

P = si A     . (15)

The resulting net demand dataset must be discretized into an odd number of equal intervals for further calculations. The discretization is performed by dividing the net demand probability distribution into equal intervals as indicated in Figure 2. These intervals are considered as scenarios with individual probabilities corresponding to the mid-point of each interval.

Figure 2. Seven-interval approximation of probability density function (PDF).

Further calculations are required to determine the level of reliability of the grid system. The reliability metric used in this study, the expected energy not supplied (EENS), represents the amount of unsupplied energy for a given time-period, and is given by the following expression:

,

1 1 1

[( ) ]

M S I

t t t

s m

m

i

s i

m s

EENS D P q q

= = =

 

=  −  

 

  

where I, M, and S are the total number of conventional units, generation system states, and net demand scenarios, respectively. 𝑃_𝑖,𝑚^𝑡 is the power that is available to the system from a conventional unit i during time-period t at system state m, and q^s, and q^m are the probability of net demand and generation system availability scenarios, respectively. It is worth noting that although variable 𝐸𝐸𝑁𝑆_𝑚^𝑡 highly depends on the level of capacity forced out of service, the probability of this outage may have an even stronger impact on the loss of energy expectation. For instance, a simultaneous failure of two or more units may cause a significant disruption of the electricity supply. However, the likelihood of this event is very low; thus, the overall loss of energy expectation would be lower as compared to a single-unit outage event.

2.4. Interconnected Capacity and Transmission Lines

Nowadays, it is quite rare to see a power system operating in island mode (i.e., in complete isolation from neighboring power systems), since grid interconnection greatly improves reliability and adequacy of interconnected systems. Interconnection can help to keep the balance between generation and load during peak load hours more efficiently, utilize generation capacity in a cost-effective manner, and reduce capacity reserve requirements [46].

Figure 2.Seven-interval approximation of probability density function (PDF).

Further calculations are required to determine the level of reliability of the grid system.

The reliability metric used in this study, the expected energy not supplied (EENS), represents the amount of unsupplied energy for a given time-period, and is given by the following expression:

EENS^t= XM

m=1







 XS

s=1

[(D^t_s− XI

i=1

P^t_i,m)^×qs]









qm (16)

whereI,M, andSare the total number of conventional units, generation system states, and net demand scenarios, respectively. P^t_i,mis the power that is available to the system from a conventional uniti during time-periodtat system statem, andqs, andqmare the probability of net demand and generation system availability scenarios, respectively. It is worth noting that although variableEENS^t_mhighly depends on the level of capacity forced out of service, the probability of this outage may have an even

stronger impact on the loss of energy expectation. For instance, a simultaneous failure of two or more units may cause a significant disruption of the electricity supply. However, the likelihood of this event is very low; thus, the overall loss of energy expectation would be lower as compared to a single-unit outage event.

2.4. Interconnected Capacity and Transmission Lines

Nowadays, it is quite rare to see a power system operating in island mode (i.e., in complete isolation from neighboring power systems), since grid interconnection greatly improves reliability and adequacy of interconnected systems. Interconnection can help to keep the balance between generation and load during peak load hours more efficiently, utilize generation capacity in a cost-effective manner, and reduce capacity reserve requirements [46].

This study uses the equivalent assisting unit approach described in references [5] and [11] to account for interconnected power flows. The maximum power assistance that can be provided by the interconnected system is the minimum of free inter-zonal capacity and the tie-line capacity [47]:

IR^max,t = min











J

X

j=1

(IR^inst_j −IR e^t_j−IR r^t_j,),

L

X

l

B^max_l











(17)

where IR^inst_j is the installed capacity of interconnected unit j, and IRe^t_j and IRr^t_j are the capacity committed for energy generation and spinning reserve service provided by inter-zonal unit j at time-periodt, respectively.B^max_l is the maximum transmission capacity of transmission linel. Finally,J andLare the total number of interconnected reserve units and transmission lines.

The maximum capacity assistance level is utilized to create a capacity model in the same way as it was described in the previous subsection. The resulting COPT is regarded as an equivalent multi-stage generating unit, which can be integrated in an existing capacity model of an assisted system.

A similar approach is employed to generate the transmission system model, where all possible scenarios are determined given the availability rate of each transmission line. The combined generation–transmission system model represents the table with state probabilities and corresponding available capacities.

2.5. Demand Response

Along with interconnected capacity, electricity consumers can significantly improve a system’s ability to withstand sudden capacity outages by participating in DR programs. Various types of DR programs allowing consumer participation in electricity markets as a service provider have been developed by ISOs. The complete description of the most common DR program types can be found in reference [48]. The main goal of the DR program employed by the proposed model is to provide up-spinning reserve service in the ancillary service market. In the proposed market structure, the DRPs serve as aggregators of DR offers provided by retail customers. A bid submitted by the DRP includes the energy and capacity costs of reserve, as well as the minimum and maximum load reduction duration. The energy cost of reserve is paid to the DRP only in cases of actual reserve deployment, whereas the capacity cost is covered in any circumstance. The capacity cost of reserve can be submitted in the form of a price–quantity pair, which denotes the amount of remuneration for reduction of certain amounts of load. For instance, the pair 15 $/MW–10 MW implies that the DRP requires 15 $/MW for reducing its load by up to 10 MW [48].

Note that the amount of DR reserve offered by the DRP must exceed the minimum amount specified by ISO. In the ancillary service market, DRPs can be treated in the same way as conventional generators, thus the DRP price–quantity functions are given by the following expressions [20]:

DR^t_k =_DR^min_{k u}^t_k+ XN

n=1

αⁿ_kuⁿ_k (18)

CDR^t_k = cd^min_k DR^min_k u_k+ XN

n=1

cd_kⁿα_kⁿu_kⁿ (19)

αⁿ_k =DRⁿ_k−DRⁿ_k⁻¹ (20)

whereDR^t_kis the amount of demand response provided by thekth DRP at time-periodt,DR^min_k is the minimum amount of demand response provided by thekth DRP,u^tkis the binary indicator of thekth DRP’s status at time-periodt(0—not providing DR, 1—providing DR),αⁿ_kis thekth DRP’snth segment of the piecewise linear cost function, anduⁿ_kis the binary indicator of thekth DRP’s status duringnth segment of the piecewise linear cost function (0—not providing DR, 1—providing DR).CDR^t_kis the cost of DR provided by thekth DRP at time-periodt.cd^min_k is the minimum price provided by thekth DRP,cdⁿ_kis the price provided by thekth DRP duringnth segment of the piecewise linear cost function, andDRⁿ_k is thenth segment of the piecewise linear cost function provided by thekth DRP.

2.6. Stochastic SCUC 2.6.1. Objective Function

As it was noted previously, the unit commitment problem of this model is expressed as a two-stage stochastic MILP. The first-stage calculations are performed for the most favorable scenario, which is called the base case. The base case denotes the generation system state when all units are available for energy production, resulting in the most efficient unit commitment. The objective of the first stage is to find the unit commitment configuration resulting in the lowest system operating cost, where the system operating cost is given by [46]:

Ctotal = PT t=1

PI i=1

PN n=1

(λ^t_i,nP^t_i,n+C^t_iminu^t_i+CS^t_i) + ^PI

i=1

(Cup^t_iRup^t_i+Cdw^t_iRdw^t_i) +

J

P

j

CIR^t_jIR^t_j +

PK k=1

CDR^t_k+ PM m=1

SC^t_m

!

(21)

whereλ^ti,nis the slope of thenth segment of theith unit’s piecewise linear cost function at time-period t,P^t_i,n is theith unit’snth segment power production (MW) during time-periodt, C^t_i,minis theith unit’s minimum operating cost at time-periodt,u^t_i is the binary indicator of theith unit’s status at time-periodt(0—not operating, 1—operating),CS^t_iis theith unit’s start-up cost during time-periodt, Cup^t_iis theith unit’s cost for provision of up-spinning reserve during time-periodt,Rup^t_iis the amount of up-spinning reserve provided by theith unit (MW) during time-periodt,Cdw^t_iis theith unit’s cost for provision of down-spinning reserve during time-periodt,Rdw^t_iis the amount of down-spinning reserve provided by theith unit (MW) during time-periodt,CIR^t_jis the cost of power generated by

thejth interconnected unit during time-periodt, andIR^t_jis the amount of reserve provided by thejth interconnected unit during time-periodt. The variableSC^t_mis given by:

SC^t_m = ^PI

i=1

PN n=1

λ^t_i,nr^t_i,n,m+

J

P

j=1

PN n=1

λ^t_j,nir^t_j,n,m + ^PK

k=1

PN

n=1λ^t_k,ndr^t_k,n,m + VOLL×CE^t_m

(22)

wherer^t_i,n,mis the deployed spinning reserve from thenth block of energy offer by theith unit in the mth scenario during time-periodt,ir^t_j,n,mis the deployed up-spinning reserve from thenth block of energy offer provided by thejth inter-zonal unit in themth scenario during time-periodt,dr^t_k,n,mis the deployed spinning reserve from thenth block of energy offer by thekth DRP in themth scenario during time-periodt,VOLLis the value of lost load—the predefined constant that indicates the per MWh cost of interruption of electricity supply—andCE^t_mis the amount of energy curtailed due to shortage of power generating capacity in themth scenario during time-periodt. The objective of the second stage is to find the most optimal reserve schedule by comparing different scenarios.

2.6.2. First-Stage Constraints

The objective function (21) must be minimized by enforcing the set of constraints presented by Equations (23)–(33). The first-stage constraints are given by the Equations (23)–(27).

The balance between power production and load is specified by the following equation:

(D ^t − CE^t) = ^PI

i=1

(P_i^tu^t_i+Rup_i^tu^t_i−Rdw^t_i u^t_i) +

J

P

j=1

IR^t_ju^t_j+ PK k=1

DR_k^tu^t_k

. (23)

Power generating units are subject to the operating constraints, such as ramping rate, maximum and minimum operating capacity, and minimum up and down time.

Additional unit operating constraints that should be specified separately are given by Equations (24) and (25):

Rup_i^t≤P^inst_i −P_i^t∀t,∀i (24) Rdw_i^t≤P^t_i−P^min_i ∀t,∀i (25) whereP^inst_i andP^min_i are the installed capacity and the minimum power output of anith unit.

As it was noted above, this model assumes that the interconnected capacity and DRPs provide only up-spinning reserve; thus, their operating constraints are given by:

Rup_j^t≤IR_j^inst−IR_j^t∀t,∀j (26)

Rup_k^t≤DR_k^max−DR_k^min,t∀t,∀k (27)

whereDR^max_k is the maximum curtailment level of thekth demand response provider.

2.6.3. Second-Stage Constraints

The second-stage scenarios are generated using the Monte-Carlo method; the constraints associated with the second-stage scenarios are presented below. For all time periods and scenarios, the power balance equation is given by:

qm I

P

i=1

(P_i^t+rup_i,m^t +rdw_i,m^t ) +

J

P

j=1

ir_j,m^t

= ^PK

k=1

(D_k,^t_{f x} − dr_k,m^t −CE_k,m^t )

(28)

whererup^t_i,mandrdw^t_i,mare the deployed up and down spinning reserves by theith unit in themth scenario during time-periodt, respectively.

The maximum involuntary load curtailment for each DRP is specified by the following equation:

0≤L^t_k,m≤L^t_k,max∀t,∀k,∀m. (29)

Inequalities (30) and (31) specify the limits that are enforced to deployed spinning reserves:

0≤r^t_i,up,m≤q^t_mR^t_i,up (30)

0≤r^t_i,dw,m≤qmR^t_i,dw. (31)

Deployed interconnected capacity and DRP reserve constraints are specified by inequalities (32) and (33):

0≤ir^t_j,m≤qmIR_j^t (32)

0≤dr^t_k,m≤qmDR^t_k. (33)

This formulation of reliability criteria allows for the ability to find the optimal spinning reserve requirement by balancing the spinning reserve operating costs and socioeconomic costs associated with load curtailment. On one hand, the reduced spinning reserve requirement may have negative effects on the adequacy of a power system; on the other hand, this reduction may be reasonable if the risk associated with capacity deficit is insignificant, or the social value of the curtailed load is very low.

3. Case Study

3.1. Test System Description

This section presents a case study performed on the IEEE two-area RTS [49]. The two-area RTS is formed by merging two single-area RTSs labeled as “Area A” and “Area B” [49] connected through three interconnection lines. The original test system presented in reference [49] was slightly modified by replacing the hydro unit located at bus 122 with one 600 MW wind farm and a 200 MW solar farm.

In this test case, Area A and Area B are considered as the assisted and assisting systems.

The assisting system provides spinning reserves only when its intra-zonal needs for energy generation and reserve capacity are fully covered. The intra-zonal capacity requirements for the assisting system are determined using the same approach as that for the assisted system, but without considering inter-zonal power flows and DR.

The cost function coefficients and unit operational constraints were taken from reference [50].

The analysis was conducted for a 24-h operating horizon. In this test system, DRPs are represented by load busses. We assume that the maximum amount of offered spinning reserve amounts to 30% of DRPs power consumption. The DR capacity fee is set to 10 $/MW with 20 $/MWh marginal cost and VOLL is set to 4000 $/MWh.

The simulation was performed in MATLAB R2018a (Version: 9.3.0.713579 (R2017b), Manufacturer:

MathWorks Inc., Natick, MA, USA) [51]. The MILP optimization was done in IBM ILOG CPLEX Optimization Studio (Version: 12.7.1, Manufacturer: IBM Corp., Armonk, NY, USA) [52] using YALMIP (Version: R20180413, Manufacturer: Johan Löfberg, Taipei, Taiwan, China) [53]. The computational efficiency of the model is achieved by considering the system state probabilities below 10⁵.

3.2. Comparison of System Configurations

The proposed model was applied to three different configurations of the modified IEEE two-area RTS to demonstrate how the presence of interconnected capacity and DR affects the adequacy of Area A system and the total costs associated with reserves schedules. The three analyzed system configurations are as follows:

Configuration 1: Area A system is assumed to operate in complete isolation (i.e., no interconnections with Area B system), and DR program is not in use;

Configuration 2: Area A system is connected to Area B system through interconnection lines, but DR program is not in use;

Configuration 3: Area A system is connected to Area B through interconnection lines, and DRPs can provide ancillary services as described in Sections2.4and2.5.

The objective of this analysis is to compare the following main indicators: Spinning reserve requirements (in MW), EENS, and the total costs of scheduled reserves (including only the cost of reserve procurement and the cost of load curtailment). Concurrently, these indicators are compared with respect to different DR marginal costs.

3.2.1. Spinning Reserve Allocation

First of all, to clarify some points and conclusions that will be made further, it is important to get acquainted with the load profile of the assisted system (Figure3). As can be seen, the load peak-time occurs between 7th and 21st hours, whereas the off-peak load occurs during the nighttime.

Energies 2019, 12, x FOR PEER REVIEW 13 of 25

3.2.1. Spinning Reserve Allocation

First of all, to clarify some points and conclusions that will be made further, it is important to get acquainted with the load profile of the assisted system (Figure 3). As can be seen, the load peak-time occurs between 7th and 21st hours, whereas the off-peak load occurs during the nighttime.

Figure 3. Load profile of the assisted system.

The spinning reserve allocation under different system configurations are presented in Figure 4.

The results show that the least amount of reserve capacity for all time instances is allocated under the system configuration 1, since under this configuration the availability of cheap capacity is reduced as compared to configurations 2 and 3. However, identical amounts of reserves are allocated during off-peak time under configurations 2 and 3. It should be noted that system configurations 2 and 3 differ from each other only in terms of the presence of DR. During the off-peak time, the abundance of cheap conventional capacity limits DRPs’ ability to compete in the reserve market. Thus, the off-peak time reserve schedule includes only conventional inter- and intra-zonal units, which both system configurations have in the same amount. However, during peak hours, DRPs become more competitive, as most of the cheap units are dedicated for energy production. Therefore, the peak time reserve requirement includes conventional intra- and inter-zonal units and DRPs under system configuration 3; whereas under system configuration 2, only intra- and inter-zonal units make up the reserve.

Figure 3.Load profile of the assisted system.

The spinning reserve allocation under different system configurations are presented in Figure4.

The results show that the least amount of reserve capacity for all time instances is allocated under the system configuration 1, since under this configuration the availability of cheap capacity is reduced as compared to configurations 2 and 3. However, identical amounts of reserves are allocated during off-peak time under configurations 2 and 3. It should be noted that system configurations 2 and 3 differ from each other only in terms of the presence of DR. During the off-peak time, the abundance of cheap conventional capacity limits DRPs’ ability to compete in the reserve market. Thus, the off-peak time reserve schedule includes only conventional inter- and intra-zonal units, which both system configurations have in the same amount. However, during peak hours, DRPs become more competitive,

as most of the cheap units are dedicated for energy production. Therefore, the peak time reserve requirement includes conventional intra- and inter-zonal units and DRPs under system configuration 3; whereas under system configuration 2, only intra- and inter-zonal units make up the reserve.

Energies 2019, 12, x FOR PEER REVIEW 14 of 25

Figure 4. Day-ahead spinning reserve allocation under different system configurations.

Figure 5 presents the spinning reserve allocation under system configuration 3, which was calculated using DR marginal costs of 10 $/MWh, 20 $/MWh, and 30 $/MWh. In this case study, the marginal cost of 20 $/MWh is considered as the base case. The results indicate that with a marginal cost of 10 $/MWh, the reserve schedules include DR reserve—not only during peak-time but also during off-peak time. According to the calculation results, the lowest marginal cost provided by conventional generators is approximately equal to 9 $/MWh. With DR marginal cost equal to 10

$/MWh, DRPs are capable of competing with cheap conventional units, resulting in the DRPs’

inclusion in the off-peak reserve schedule.

Figure 4.Day-ahead spinning reserve allocation under different system configurations.

Figure5presents the spinning reserve allocation under system configuration 3, which was calculated using DR marginal costs of 10 $/MWh, 20 $/MWh, and 30 $/MWh. In this case study, the marginal cost of 20 $/MWh is considered as the base case. The results indicate that with a marginal cost of 10 $/MWh, the reserve schedules include DR reserve—not only during peak-time but also during off-peak time. According to the calculation results, the lowest marginal cost provided by conventional generators is approximately equal to 9 $/MWh. With DR marginal cost equal to 10 $/MWh, DRPs are capable of competing with cheap conventional units, resulting in the DRPs’ inclusion in the off-peak reserve schedule.

Energies 2019, 12, x FOR PEER REVIEW 15 of 25

Figure 5. Day-ahead spinning reserve requirements under different demand response (DR) marginal costs.

3.2.2. EENS

The total EENS for each system configuration is calculated by integrating EENS over the scheduling horizon and presented in Table 1. A relatively high availability of cheap reserves results in the lowest EENS under system configuration 3 (Figure 6). On the other hand, insufficient capacity leads to increasing the load curtailment risk due to random outage of intra-zonal units. Figure 6 shows that a relatively low amount of reserves under system configuration 1 results in high loss of energy expectations. It should be noted that a direct correlation between scheduled reserve capacity and EENS can be observed by comparing Figures 4 and 6. Similarly, as in the case with reserve allocation, EENS values under system configurations 2 and 3 are identical for off-peak hours, whereas peak time EENS under system configuration 3 is lower.

Figure 5. Day-ahead spinning reserve requirements under different demand response (DR) marginal costs.

3.2.2. EENS

The total EENS for each system configuration is calculated by integrating EENS over the scheduling horizon and presented in Table1. A relatively high availability of cheap reserves results in the lowest EENS under system configuration 3 (Figure6). On the other hand, insufficient capacity leads to increasing the load curtailment risk due to random outage of intra-zonal units. Figure6shows that a relatively low amount of reserves under system configuration 1 results in high loss of energy

expectations. It should be noted that a direct correlation between scheduled reserve capacity and EENS can be observed by comparing Figures4and6. Similarly, as in the case with reserve allocation, EENS values under system configurations 2 and 3 are identical for off-peak hours, whereas peak time EENS under system configuration 3 is lower.

Energies 2019, 12, x FOR PEER REVIEW 16 of 25

Figure 6. Expected energy not supplied (EENS) under different system configurations.

Table 1. Total EENS (MWh) calculated for different system configurations.

Configuration 1 Configuration 2 Configuration 3

9.05 1.66 1.5

The comparison of EENS under different DR marginal costs showed that availability of low-cost spinning reserve positively affected the system reliability. As it can be observed from Figure 7, a low DR marginal cost resulted in low EENS. The spinning reserves scheduled under the 10$/MWh-scenario led to the lowest EENS as compared to the base case and the case with 30

$/MWh.

Figure 6.Expected energy not supplied (EENS) under different system configurations.

Table 1.Total EENS (MWh) calculated for different system configurations.

Configuration 1 Configuration 2 Configuration 3

9.05 1.66 1.5

The comparison of EENS under different DR marginal costs showed that availability of low-cost spinning reserve positively affected the system reliability. As it can be observed from Figure7, a low DR marginal cost resulted in low EENS. The spinning reserves scheduled under the 10$/MWh-scenario led to the lowest EENS as compared to the base case and the case with 30 $/MWh.

Energies 2019, 12, x FOR PEER REVIEW 17 of 25

Figure 7. EENS under different DR marginal costs.

3.2.3. Total Costs of Scheduled Reserves

The total costs of scheduled reserves are presented in Table 2. These costs are calculated for the reserve schedules determined during the second stage of unit commitment optimization and include the operating costs of reserve and the costs of possible load curtailment. In Figure 8, both the highest and lowest costs can be observed under the system configurations 1 and 3, respectively. Each system configuration resulted in lower costs during off-peak time relative to peak time, due to low EENS and excess of cheap reserve capacity during off-peak hours. During peak time, increased EENS and lack of cheap reserves result in relatively high costs. Similarly, as in the previous section, reserve schedules under system configurations 2 and 3 resulted in identical costs during off-peak hours, whereas peak hours show slightly lower costs for system configuration 3. Thus, for this particular case, it can be concluded that the main reason for the lowest total cost under system configuration 3 was a relatively low cost of DR.

Figure 7.EENS under different DR marginal costs.

3.2.3. Total Costs of Scheduled Reserves

The total costs of scheduled reserves are presented in Table2. These costs are calculated for the reserve schedules determined during the second stage of unit commitment optimization and include the operating costs of reserve and the costs of possible load curtailment. In Figure8, both the highest and lowest costs can be observed under the system configurations 1 and 3, respectively. Each system

configuration resulted in lower costs during off-peak time relative to peak time, due to low EENS and excess of cheap reserve capacity during off-peak hours. During peak time, increased EENS and lack of cheap reserves result in relatively high costs. Similarly, as in the previous section, reserve schedules under system configurations 2 and 3 resulted in identical costs during off-peak hours, whereas peak hours show slightly lower costs for system configuration 3. Thus, for this particular case, it can be concluded that the main reason for the lowest total cost under system configuration 3 was a relatively

low cost of DR. Energies 2019, 12, x FOR PEER REVIEW 18 of 25

Figure 8. Total costs of scheduled reserves under different system configurations.

Table 2. Total operational costs ($) calculated for different system configurations.

Configuration 1 Configuration 2 Configuration 3

1,203,800 969,740 958,790

Figure 9 presents the total costs of scheduled reserves calculated for DR marginal costs of 10

$/MWh, 20 $/MWh, and 30 $/MWh. The costs are almost identical for all scenarios during off-peak hours, whereas the peak-time costs are the lowest for the 10 $/MWh-scenario, and the highest for the 30 $/MWh-scenario. The reduction of the total cost is achieved due to two main reasons: First, although the total amount of scheduled spinning reserves is higher in the 10 $/MWh-scenario, the DR reserve operating costs are almost similar due to the low marginal cost of this scenario; second, the increased reserve schedule resulted in reduced EENS, which in turn, reduced the costs of load curtailment.

Figure 8.Total costs of scheduled reserves under different system configurations.

Table 2.Total operational costs ($) calculated for different system configurations.

Configuration 1 Configuration 2 Configuration 3

1,203,800 969,740 958,790

Figure9 presents the total costs of scheduled reserves calculated for DR marginal costs of 10 $/MWh, 20 $/MWh, and 30 $/MWh. The costs are almost identical for all scenarios during off-peak hours, whereas the peak-time costs are the lowest for the 10 $/MWh-scenario, and the highest for the 30 $/MWh-scenario. The reduction of the total cost is achieved due to two main reasons: First, although the total amount of scheduled spinning reserves is higher in the 10 $/MWh-scenario, the DR reserve operating costs are almost similar due to the low marginal cost of this scenario; second, the increased reserve schedule resulted in reduced EENS, which in turn, reduced the costs of load curtailment.

Energies 2019, 12, x FOR PEER REVIEW 19 of 25

Figure 9. Total costs of scheduled reserves under different DR marginal costs.

3.3. Comparison of Wind Prediction Models

First of all, it should be noted that the purpose of this section is not to compare the accuracy of univariate and bivariate wind prediction models, but rather to demonstrate the difference between two models in terms of the grid system adequacy. The readers who are interested in evaluation of these models in terms of precision are suggested to see reference [43].

Figures 10 and 11 display the reserve requirements and EENS of different system configurations obtained using univariate and bivariate wind prediction models.

Figure 9.Total costs of scheduled reserves under different DR marginal costs.