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

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.

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, 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,8,9,10,11,12,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: 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.

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 unit i on outage is given by [5]:

$$U_i=1-e^{-{\gamma }_{i}T}$$

(1)

where ${\gamma }_i$ is the failure rate of unit i and T is the lead time. Note that throughout this paper, indices i, j, and k refer to the intra-zonal conventional units, inter-zonal conventional units, and DRPs, respectively. Indices m, s, and t refer to the generation system states, net demand scenarios, and the time-periods of the scheduling horizon. Finally, capital P represents the power output of conventional and renewable energy producers, whereas the lowercase p refers 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]:

$$P_{m,i}(X)=(1-U_i)p_{m-1}+(U_i)p_{m-1}(X-P_i)$$

(2)

where $p_{m-1}\left(X\right)$ and $p_{m,i}\left(X\right)$ are, respectively, the cumulative probabilities of generation system state m before and after unit i is added, and P_i is the installed capacity of unit i. 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 of i).

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-period t, the net demand is given by:

$$D^t=L^t-P{}_W^t-P_{PV}^t$$

(3)

where $L^t$, $P_W^t$, and $P_{PV}^t$ 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]:

$$\sigma {}_L^t=\frac{Y}{100}L{}_F^t$$

(4)

where ${\sigma }_L^t$ and $L_F^t$ are the standard deviation of the load forecast error and forecasted load at time t, respectively, and Y is 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,10,11], the Weibull distribution [32,33,34,35,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]:

$$f_w(v)=(k/\gamma ){(v/\gamma )}^{k-1}\exp [-{(v/\gamma )}^k]$$

(5)

where v is the actual wind velocity, and γ and k are 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]:

$$P_W=\left\{\begin{array}{ll}0& v<v_{\min }\\ P_f(v)& v_{\min }\le v<v_r\\ P_r& v_r\le v<v_{\max }\\ 0& v\ge v_{\max }\end{array}\right\}$$

(6)

where P_r is the maximum power output of a wind turbine rated at wind velocity v_r, v_min, and v_max are the cut-in and cut-off wind velocities, respectively. Finally, $P_f\left(v\right)$ is the non-linear part of the power curve, representing the power output at wind velocity v. In this study, we utilize the power coefficient-based model given by reference [42]:

$$P_f(v)=\frac{1}{2}{\rho }_{air}\times A_{wt}\times C_{eq}\times v^3$$

(7)

where ρ_air is the air density, considered constant (1.225 kg/m³), A_wt is the wind turbine rotor swept area, and C_eq is 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:

$$f_{V,\theta }(v,\theta )=2\pi f_{\psi }(\psi )f_V(v)f_{\theta }(\theta ),0<v<v_{\max },0\le \theta \le 2\pi$$

(8)

where V and θ are the random variables representing wind speed and direction, respectively, and ψ is the random variable representing the relationship structure between V and θ, which is given by the following expression:

$$\psi =2\pi [F_V(v)-F_{\theta }(\theta )],0\le \theta \le 2\pi$$

(9)

where the terms $F_V\left(v\right)$, $F_{\theta }\left(\theta \right)$, and $F_{\psi }\left(\psi \right)$ are, respectively, the cumulative distribution functions of V, θ, and ψ.

Equation (7) for the bivariate case becomes:

$$P_f(v)=\frac{1}{2}{\rho }_{air}\cdot A_{wt}\cdot C_{eq}\cdot v^3\cdot f_{V,\theta }(v,\theta ).$$

(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)=\left\{\begin{array}{ll}\frac{\Gamma (\alpha +\beta )}{\Gamma (\alpha )\Gamma (\beta )}\times s{i}^{(\alpha -1)}\times (1-si{)}^{(\beta -1)}& for0\le si\le 1,\alpha \ge 0,\beta \ge 0\\ 0& otherwise\end{array}\right\}$$

(11)

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

$${\displaystyle \underset{0}{\overset{\infty }{\int }}{\rho }^{\lambda -1}{e}^{-\rho }d\rho }$$

(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]:

$$\beta =(1-\mu )\left(\frac{\mu (1+\mu )}{{\sigma }^2}-1\right)$$

(13)

$$\alpha =\frac{\mu \beta }{1-\mu }$$

(14)

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]:

$$P_{PV}=si\cdot A_{PV}\cdot {\eta }_{PV}\cdot {\eta }_{mppt}\cdot \cos {\theta }_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 Figure 2. These intervals are considered as scenarios with individual probabilities corresponding to the mid-point of each interval.

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:

$$EEN{S}^t={\displaystyle \sum _{m=1}^M\left({\displaystyle \sum _{s=1}^S[(D_s^t-{\displaystyle \sum _{i=1}^I{P}_{i,m}^t})\times q_s]}\right)}{q}_m$$

(16)

where I, M, and S are the total number of conventional units, generation system states, and net demand scenarios, respectively. $P_{i,m}^t$ 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 $EEN{S}_m^t$ 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].

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 \left\{{\displaystyle \sum _{j=1}^J(IR{}_j^{inst}}-IRe{}_j^t-IRr{}_{j,}^t),{\displaystyle \sum _l^L{B}_l^{\max }}\right\}$$

(17)

where $I{R}_j^{inst}$ is the installed capacity of interconnected unit j, and $IR{e}_j^t$ and $IR{r}_j^t$ are the capacity committed for energy generation and spinning reserve service provided by inter-zonal unit j at time-period t, respectively. $B_l^{\max }$ is the maximum transmission capacity of transmission line l. Finally, J and L are 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]:

$$D{R}_k^t=D{R}_k^{\min }{u}_k^t+{\displaystyle \sum _{n=1}^N{\alpha }_k^n{u}_k^n}$$

(18)

$$CD{R}_k^t=c{d}_k^{\min }D{R}_k^{\min }{u}_k+{\displaystyle \sum _{n=1}^{N}c{d}_k^n{\alpha }_k^n{u}_k^n}$$

(19)

$$\alpha {}_k^n=DR{}_k^n-DR{}_k^{n-1}$$

(20)

where $D{R}_k^t$ is the amount of demand response provided by the kth DRP at time-period t, $D{R}_k^{\min }$ is the minimum amount of demand response provided by the kth DRP, u^t_k is the binary indicator of the kth DRP’s status at time-period t (0—not providing DR, 1—providing DR), ${\alpha }_k^n$ is the kth DRP’s nth segment of the piecewise linear cost function, and uⁿ_k is the binary indicator of the kth DRP’s status during nth segment of the piecewise linear cost function (0—not providing DR, 1—providing DR). $CD{R}_k^t$ is the cost of DR provided by the kth DRP at time-period t. $c{d}_k^{\min }$ is the minimum price provided by the kth DRP, $c{d}_k^n$ is the price provided by the kth DRP during nth segment of the piecewise linear cost function, and $D{R}_k^n$ is the nth segment of the piecewise linear cost function provided by the kth 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]:

$$\begin{array}{ccc}{C}_{total}& =& {\displaystyle \sum _{t=1}^T({\displaystyle \sum _{i=1}^I{\displaystyle \sum _{n=1}^N({\lambda }_{i,n}^t{P}_{i,n}^t+C_{i\min }^t{u}_i^t+C{S}_i^t)}}}\hfill \\ & +& {\displaystyle \sum _{i=1}^I(}Cu{p}_i^{t}Ru{p}_i^t+Cd{w}_i^{t}Rd{w}_i^t)+{\displaystyle \sum _j^{J}CI{R}_j^{t}I{R}_j^t}\hfill \\ & +& {\displaystyle \sum _{k=1}^{K}CD{R}_k^t}+{\displaystyle \sum _{m=1}^{M}S{C}_m^t})\hfill \end{array}$$

(21)

where λ^t_i,n is the slope of the nth segment of the ith unit’s piecewise linear cost function at time-period t, $P_{i,n}^t$ is the ith unit’s nth segment power production (MW) during time-period t, $C_{i,\min }^t$ is the ith unit’s minimum operating cost at time-period t, $u_i^t$ is the binary indicator of the ith unit’s status at time-period t (0—not operating, 1—operating), $C{S}_i^t$ is the ith unit’s start-up cost during time-period t, $Cu{p}_i^t$ is the ith unit’s cost for provision of up-spinning reserve during time-period t, $Ru{p}_i^t$ is the amount of up-spinning reserve provided by the ith unit (MW) during time-period t, $Cd{w}_i^t$ is the ith unit’s cost for provision of down-spinning reserve during time-period t, $Rd{w}_i^t$ is the amount of down-spinning reserve provided by the ith unit (MW) during time-period t, $CI{R}_j^t$ is the cost of power generated by the jth interconnected unit during time-period t, and $I{R}_j^t$ is the amount of reserve provided by the jth interconnected unit during time-period t. The variable $S{C}_m^t$ is given by:

$$\begin{array}{ccc}SC{}_m^t& =& {\displaystyle \sum _{i=1}^I{\displaystyle \sum _{n=1}^N{\lambda }_{i,n}^t{r}_{i,n,m}^t}}+{\displaystyle \sum _{j=1}^J{\displaystyle \sum _{n=1}^N{\lambda }_{j,n}^{t}i{r}_{j,n,m}^t}}\hfill \\ & +& {\displaystyle \sum _{k=1}^K{\displaystyle \sum _{n=1}^N{\lambda }_{k,n}^{t}d{r}_{k,n,m}^t}}+VOLL\times C{E}_m^t\hfill \end{array}$$

(22)

where $r_{i,n,m}^t$ is the deployed spinning reserve from the nth block of energy offer by the ith unit in the mth scenario during time-period t, $i{r}_{j,n,m}^t$ is the deployed up-spinning reserve from the nth block of energy offer provided by the jth inter-zonal unit in the mth scenario during time-period t, $d{r}_{k,n,m}^t$ is the deployed spinning reserve from the nth block of energy offer by the kth DRP in the mth scenario during time-period t, VOLL is the value of lost load—the predefined constant that indicates the per MWh cost of interruption of electricity supply—and $C{E}_m^t$ is the amount of energy curtailed due to shortage of power generating capacity in the mth scenario during time-period t. 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:

$$\begin{array}{ccc}(D{}^t-CE{}^t)& =& {\displaystyle \sum _{i=1}^I(P_i^{t}u{}_i^t}+Ru{p}_i^{t}u{}_i^t-Rdw{}_i^{t}u{}_i^t)\hfill \\ & +& {\displaystyle \sum _{j=1}^{J}IR{}_j^{t}u{}_j^t+{\displaystyle \sum _{k=1}^{K}DR{}_k^{t}u{}_k^t}}\hfill \end{array}.$$

(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):

$$Ru{p}_i^t\le P{}_i^{inst}-P_i^t\forall t,\forall i$$

(24)

$$Rd{w}_i^t\le P_i^t-P{}_i^{\min }\forall t,\forall i$$

(25)

where $P_i^{inst}$ and $P_i^{\min }$ are the installed capacity and the minimum power output of an ith 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:

$$Ru{p}_j^t\le I{R}_j^{inst}-I{R}_j^t\forall t,\forall j$$

(26)

$$Ru{p}_k^t\le D{R}_k^{\max }-D{R}_k^{\min ,t}\forall t,\forall k$$

(27)

where $D{R}_k^{\max }$ is the maximum curtailment level of the kth 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:

$$\begin{array}{l}{q}_m{\displaystyle \sum _{i=1}^I(P_i^t+ru{p}_{i,m}^t+rd{w}_{i,m}^t)}+{\displaystyle \sum _{j=1}^{J}i{r}_{j,m}^t}\\ ={\displaystyle \sum _{k=1}^K(D_{k,fx}^t-d{r}_{k,m}^t-C{E}_{k,m}^t)}\end{array}$$

(28)

where $ru{p}_{i,m}^t$ and $rd{w}_{i,m}^t$ are the deployed up and down spinning reserves by the ith unit in the mth scenario during time-period t, respectively.

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

$$0\le L{}_{k,m}^t\le L{}_{k,\max }^t\forall t,\forall k,\forall m.$$

(29)

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

$$0\le r_{i,up,m}^t\le q_m^t{R}_{i,up}^t$$

(30)

$$0\le r_{i,dw,m}^t\le q_m{R}_{i,dw}^t.$$

(31)

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

$$0\le i{r}_{j,\mathrm{m}}^t\le q_{m}I{R}_j^t$$

(32)

$$0\le d{r}_{k,m}^t\le q_{m}D{R}_k^t.$$

(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 Section 2.4 and Section 2.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 (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.

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 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.

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. Total EENS (MWh) calculated for different system configurations.

Configuration 1	Configuration 2	Configuration 3
9.05	1.66	1.5

. 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 Figure 4 and Figure 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.

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.

3.2.3. Total Costs of Scheduled Reserves

The total costs of scheduled reserves are presented in

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

Configuration 1	Configuration 2	Configuration 3
1,203,800	969,740	958,790

. 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 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.

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].

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

As it can be seen, the reserve requirements obtained by implementation of the bivariate wind prediction model are somewhat less as compared to the univariate case. In fact, the average reserve requirement decline was 0.27%. The main reason for this decline is reduction in wind forecast error and relatively large share of wind power production in the test system (5.84% in terms of generated energy). The same can be observed in EENS (Figure 11).

3.4. Computation Time

The simulations presented in the previous sections were performed on a desktop PC with an Intel Core i5-7200 CPU running at 2.50 GHz and 4.00 GB RAM on Windows 10 operating system. The test system consisted of 53 units and 48 buses. The total time required to run the model averages to 19 min, which is mostly spent carrying out the unit commitment and spinning reserve allocation optimization. The computation time required to run the simulations can be reduced by decreasing the COPT resolution, or by neglecting the system states with low probability. However, it is important to note that this reduction will negatively affect the accuracy of the model.

To evaluate the performance of the proposed model on a larger scale, similar tests were conducted on a set of test systems with 2, 3, and 5 times as many generators and buses as in the original two-area RTS. The results show that the time required to run the simulations grow exponentially, with an average of almost 80 min required to carry out the calculations for a system consisting of 265 generators and 240 buses. Clearly, this kind of performance is unacceptable for a real grid system with much more complex topology, and a larger number of generators and buses. Nevertheless, the computational time can be easily reduced to an acceptable range by implementing the model in a higher-order programming language along with running it on a high-performance computer.

4. Conclusions

A probabilistic model for estimation of spinning reserves in interconnected systems with demand side reserve provision was discussed in this work. The proposed model accounts for random outages of conventional units as well as uncertainties caused by the load, solar, and wind forecast errors. Random unavailability of generating capacity was modeled as a two-state Markov process. The load, wind, and solar forecast errors were modeled using normal, Weibull, and Beta distributions. In cases of wind power production, two prediction models (namely univariate and bivariate) were considered. Additionally, the model considered the interconnected capacity through utilization of the equivalent assisting multi-state unit approach. In this study, interconnected capacity and DRPs were assumed to provide only up-spinning reserve service. The model was developed according to the specific market structure in which the intra-zonal conventional units participate in the energy and ancillary service markets, whereas the inter-zonal capacity and DRPs participate only in the ancillary service market. The two-stage unit commitment problem was formulated such that MILP could be applied to conduct the optimization. In this formulation, the optimal unit commitment for energy production was determined in the first stage, while quantification of spinning reserves was conducted in the second stage of the program by recursion through the most probable scenarios. Furthermore, to minimize the cost function, the cost of reserve along with the cost of possible load shedding were taken as objectives and the spinning reserve optimization was conducted.

The proposed model was evaluated on the modified version of IEEE two-area RTS. The goal of the case study was to determine how the adequacy and economic performance of Area A system would be affected by participation of inter-zonal capacity and DRPs in the ancillary service market. In addition, the case study aimed at evaluation of the influence of the bivariate wind prediction model on the overall reliability of the test system. Among three considered system configurations, the configuration with inter-zonal capacity and DR has shown the best results with desirable reliability and system costs. According to the analysis, based on the given cost of DR, the presence of DRPs did not affect the reserve requirements during off-peak hours because of their inability to compete with low-cost conventional units. However, during peak-time, DR positively affected the system reliability by reducing the loss of load expectation. Further analysis showed that reduction of DR marginal costs positively affected DRPs ability to compete with conventional units, resulting in inclusion of DR reserves in the off-peak reserve schedule. Finally, the analysis showed that the bivariate wind prediction model reduces the overall reserve requirements and EENS in all three configurations by reducing uncertainty associated with wind forecasting.

This model has demonstrated that diversity in spinning reserve capacity positively affects not only reliability of energy systems, but also the overall social welfare. Future work will focus on extending this model to the multi-connected energy systems with distributed renewable energy sources, as well as on reduction of the uncertainty related to the consumer behavior by implementation of Bayesian network theory.