Two-Stage Coordinated Operational Strategy for Distributed Energy Resources Considering Wind Power Curtailment Penalty Cost

Abstract

The concept of virtual power plant (VPP) has been proposed to facilitate the integration of distributed renewable energy. VPP behaves similar to a single entity that aggregates a collection of distributed energy resources (DERs) such as distributed generators, storage devices, flexible loads, etc., so that the aggregated power outputs can be flexibly dispatched and traded in electricity markets. This paper presents an optimal scheduling model for VPP participating in day-ahead (DA) and real-time (RT) markets. In the DA market, VPP aims to maximize the expected profit and reduce the risk in relation to uncertainties. The risk is measured by a risk factor based on the mean-variance Markowitz theory. In the RT market, VPP aims to minimize the imbalance cost and wind power curtailment by adjusting the scheduling of DERs in its portfolio. In case studies, the benefits (e.g., surplus profit and reduced wind power curtailment) of aggregated VPP operation are assessed. Moreover, we have investigated how these benefits are affected by different risk-aversion levels and uncertainty levels. According to the simulation results, the aggregated VPP scheduling approach can effectively help the integration of wind power, mitigate the impact of uncertainties, and reduce the cost of risk-aversion.

1. Introduction

Virtual power plant (VPP) is termed as a single entity that integrates the portfolio of different types of distributed energy resources (DERs) [1,2]. The integration can make the summed capacity of many diverse small-scale DERs be large enough to trade energy or provide network support services in pool markets. Meanwhile, VPP creates a single operating profile from a composite of the parameters characterizing each DER [3,4]. This aggregation enhances the visibility and controllability of DERs to system operators and other market participants [5,6]. Specifically, VPP integrates various types of renewable and nonrenewable generators, storages and flexible loads [7,8]. The weakness and strength of different DERs are combined in a complementary way, thus enabling cost-efficient, secure and sustainable energy supply systems [9,10]. Moreover, depending on the net power between production and demand, VPP can be a producer or a consumer (i.e., a prosumer). The profit variability in electricity markets becomes the major concern and VPP should optimize the scheduling of DERs and bid/offer more strategically [11]. Given the market risks in relation to uncertainties, sophisticated risk management strategies are also needed [9].

In the literature, reference [10] has proposed an optimal operation model for VPP participating in both day-ahead (DA) and balancing markets. The model is formulated as a mixed-integer linear programming (MILP) problem, while the conditional value at risk (CVaR) is used as a measure to control the risk of low profit scenarios. In [12], the optimal offering of VPP in the DA and balancing markets is also modeled as an MILP problem, aiming to maximize the expected profit. In [9], the participation of VPP in the DA and real-time (RT) markets is modeled by a stochastic Game theory approach. CVaR is used to compute risks due to uncertainties. In [1], the bidding strategy for VPP in energy and ancillary (spinning reserve services) markets is addressed. The optimal bidding problem is formulated as a single level non-equilibrium model based on the deterministic price-based unit commitment (PBUC). In [13], a stochastic bi-level approach is adopted to address the optimal bidding strategy for VPP. The bi-level problem is formulated as an MILP problem based on the Karush–Kuhn–Tucker optimality conditions and the strong duality theory. Nezamabadi et al. [14] proposed an arbitrage strategy for VPP participating in energy and ancillary markets. The proposed model is based on a security-constrained PBUC, aiming to maximize the profit considering arbitrage opportunities. Similarly, Karimyan et al. [15] used the PBUC method to identify the optimal operation of VPP in energy and reserve markets. The problem is formulated as a mixed-integer nonlinear programming (MINLP) problem, which is solved by the particle swarm optimization (PSO) algorithm. Luo et al. [16] has proposed a two-stage operational planning framework for VPP. The first stage optimizes VPP hourly scheduling in the energy market, while the second stage optimizes the RT control actions using a model predictive control-based dispatch model. Vasirani et al. [17] proposed an agent-based approach to optimal operations of VPP. The model is based on linear programming, and the interactions between the grid and electric vehicle (EV) batteries are investigated.

Most of the above-mentioned references use simplified linear models, and hence the network constraints have been neglected and some important system parameters such as reactive power are not taken into account. On the other hand, some references have not well addressed the risks for VPP participating in pool markets or failed to characterize the correlations between uncertainties (e.g., load, price or wind power availability). By contrast, in this paper, we have proposed an optimal scheduling model for a price-taker VPP participating in the DA and RT balancing markets. The complex network flow constraints and the unit commitments of DERs are explicitly taken into account, thus making the proposed model more applicable in practice. In addition, to comprehend the profit variants in markets, a risk factor is introduced into the optimal scheduling formulation based on the mean-variance Markowitz theory [18]. As a result, value-at-risk (VaR) becomes the model output in this paper, rather than a part of the objective function as reported in [19]. Moreover, the benefits of coordination (VPP aggregation) and the corresponding sensitivity analysis are demonstrated in case studies.

The remainder paper is organized as follows: in Section 2, the problem description is given, followed by proposed optimal scheduling model for VPP in Section 3. Section 4 introduces the applied solution algorithm. Section 5 presents numerical simulations to demonstrate the effectiveness of the proposed approach. Conclusions are given in Section 6.

2. Problem Formulation

As seen in Figure 1, DERs such as battery energy storage system (BESS) units, distributed generation (DG) units (thermal), wind power and solar generators, and end consumers (both interruptible and non-interruptible loads) are physically interconnected to the upstream distribution networks. It should be noted that unlike a microgrid a VPP is not necessarily limited to a geographical location [2]. VPP acts as a representative that manages all transactions of DERs with the independent system operator (ISO). In other words, the VPP needs to submit hourly scheduling plans in energy markets. Thus, VPP is exposed to the market risk and needs to optimize its scheduling based on the entire portfolio’s operational constraints and the system-state estimation [10]. Moreover, end consumers of VPP are supplied with a given retail energy rate, and they are not exposed to the price volatility in spot energy markets [7]. VPP signs a contract with interruptible consumers, in which the upper limit and cost of load curtailment are specified. DG and BESS may be used for trading energy but corresponding costs should be paid. The operational cost of wind power is assumed to be negligible.

The market structure studied in this work is a joint DA and RT electricity markets, which is common in European electricity pools such as Nord Pool [16]. In the DA (spot) market, participants propose their scheduling for each hour of the coming day before the gate closure. After the spot price has been settled and commitments have been made, participants are responsible for deviations due to unpredictable fluctuations in power production or consumption. In the RT market, any deviation from the commitment made in the DA market will be settled by a regulation price. If the actual consumption is more than (or the production is less than) the commitment, the power shortage is purchased at an up-regulation price, which is usually higher than the spot price. If the actual consumption is less than (or the production is more than) the commitment, the power surplus is sold at a down-regulation price, which is usually lower than the spot price [13]. As a result, this market structure encourages participants to reduce their forecast errors and be more strategic in operations [20].

In the above-mentioned market framework, three assumptions are made for the operation of VPP: (i) DERs in VPP are centrally controlled [1]. This means the VPP central controller is responsible for bi-directional communication, monitoring, and control of each VPP component; (ii) VPP can obtain the required data for managing its portfolio (e.g., price and load data). The relevant time-series forecast techniques are available; (iii) VPP is considered as a price-taker, which means its operational strategy dose not influence the market price.

3. Optimal Scheduling of VPP

3.1. Key Models

3.1.1. Price Model

In the DA market, the VPP schedules the hourly power $P_t^{DS}$ (i.e., quantity of power consumption or production) for the next trading day. If $P_t^{DS}>0$, VPP is a producer; if $P_t^{DS}<0$, the VPP is a consumer. Because of the variations of power production and load, there will be a deviation between the scheduling and the forecasted power. The deviation at time $t$ is:

$$\Delta P_t^{\mathrm{Im}b}=P_t^{Fct}-P_t^{DS}$$

(1)

where $P_t^{Fct}$ denotes the forecasted power to be traded.

Thus, the imbalance cost $C_t^{\mathrm{Im}b}$ is calculated as:

$$C_t^{\mathit{Im}b}=\{\begin{array}{l}{\rho }_t^{R+}\cdot \Delta P_t^{\mathit{Im}b},\Delta P_t^{\mathit{Im}b}<0\\ {\rho }_t^{R-}\cdot \Delta P_t^{\mathit{Im}b},\Delta P_t^{Imb}\ge 0\end{array}$$

(2)

where ${\rho }_t^{R+}$ and ${\rho }_t^{R-}$ denote regulation prices for purchasing (up-regulation) and selling (down-regulation) electricity in real-time balancing market at time $t$, respectively.

In addition, in the particular markets, we assume that the regulation price can be calculated as proportions of the DA market price ${\rho }_t^{DA}$:

$$\{\begin{array}{l}{\rho }_t^{R+}=\left(1+{\zeta }_t^{+}\right)\cdot {\rho }_t^{DA}\\ {\rho }_t^{R-}=\left(1-{\zeta }_t^{-}\right)\cdot {\rho }_t^{DA}\end{array}$$

(3)

where ${\zeta }_t^{+}$ and ${\zeta }_t^{-}$ denote the relative differences between the day-ahead price and the up-regulation/down-regulation prices.

Therefore, market price uncertainty can be modeled by a random vector, including random DA price and random differences between DA up and down regulation prices. The distributions of market prices can be obtained by a variety of time-series based forecasting techniques such as autoregressive integrated moving average (ARIMA), which have been well addressed in [21,22].

3.1.2. Interruptible Load (IL) Cost Model

The cost of IL $C_{it}^{IL}$ is assumed to be a function of adjusted load $P_{it}^{IL}$ at time $t$ bus $i$, and is modeled by a quadratic polynomial function as [1]:

$$C_{it}^{IL}=a_{1i}^{IL}\cdot {\left(P_{it}^{IL}\right)}^2+a_{2i}^{IL}\cdot P_{it}^{IL}$$

(4)

where $a_{1i}^{IL}$, $a_{2i}^{IL}$ denote cost coefficients of IL. It should be noted that negative IL (i.e., load increase) also exists for balancing underestimated renewable energy outputs. An example of the upward demand response (DR) could be the controllable electric vehicle (EV) loads through the vehicle-to-grid (V2G) technology.

3.1.3. BESS Cost Model

The operational cost of BESS $C_{it}^{BESS}$ at time $t$ bus $i$ generally refers to maintenance cost [16], which can be modeled by a linear function as:

$$C_{it}^{BESS}={\beta }_i^{BESS}\cdot P_{it}^{BESS}\Delta t+{\beta }_i^{BESS}{E}_{it}^{BESS}{\eta }^L\Delta t$$

(5)

where $P_{it}^{BESS}$ denote the charged or discharged BESS power at time $t$ bus $i$; $E_{it}^{BESS}$ denotes energy stored in BESS; $\Delta t$ denotes a factor that converts power (kW) to energy (kWh), i.e., time duration; ${\eta }^L$ denote leakage loss factor of BESS; and ${\beta }_i^{BESS}$ is the cost coefficient of the BESS lifetime depression, which is calculated as [23]:

$${\beta }_{it}^{BESS}=\frac{I{C}^{BESS}}{E_R^{BESS}\cdot \left(LCN\right)}$$

(6)

where $I{C}^{BESS}$, $E_R^{BESS}$, and $LCN$ denote investment cost, rated energy capacity and total life cycle number of BESS.

3.1.4. DG Cost Model

DG can refer to any kind of DERs though, in this paper DG is only termed as distributed thermal units (e.g., diesel or gas fired power generators), in order to distinguish it from wind power units. The cost of DG $C_{it}^{DG}$ can be modeled by a quadratic function of its real power output [24]:

$$C_{it}^{DG}=a_{1i}^{DG}\cdot P_{DGit}^2+a_{2i}^{DG}\cdot P_{DGit}+a_{3i}^{DG}$$

(7)

where $a_{1i}^{DG}$, $a_{2i}^{DG}$, and $a_{3i}^{DG}$ denote cost coefficients of DG at bus $i$; and $P_{DGit}$ denotes DG output at bus $i$ time $t$.

3.1.5. Model of Other Uncertainties

In this paper, load, market price, and wind speed are considered as uncertainties. The historical data in the DA market is used as correlated scenarios, and hence the correlated probability distributions can be estimated based on the statistical correlations among these uncertainties. More detailed information can be found in [21], which comprehensively discusses a variety of time-series-based methods to generate correlated scenarios. In this paper, the ARIMA method is used to generate correlated time-series wind speed, load and price data. To model the forecast errors produced by the ARIMA model, Monte Carlo (MC) simulations are used to sample scenarios. In general, load forecast errors and market price forecast errors can be modeled by the Gaussian distribution [3]. Wind speed can be modeled by the Weibull distribution [16,25].

3.2. Day-Ahead Scheduling

In the DA market, VPP optimizes its hourly scheduling and unit commitments to maximize the total profit over the scheduling horizon. Meanwhile, the decision made in this stage may affect its operational strategies in the RT market. Thus, it is necessary to consider the decision variables pertaining to the RT market while optimizing the scheduling in the DA market. However, note that the RT market decision variables are not actually implemented in the DA market.

The optimal scheduling model of VPP in the DA market is formulated as:

$$\max f_1=\left\{\begin{array}{l}{\displaystyle \sum _{t=1}^T{\rho }_t^{DA}}{P}_t^{DS}\Delta t+{\displaystyle \sum _{t=1}^T{\rho }_t^{\mathrm{Re}tail}}\left(P_{Dt}-{\displaystyle \sum _{i\in {\mathsf{\Omega }}_D}{P}_{it}^{IL}}\right)\Delta t\\ -{\displaystyle \sum _{t=1}^T{C}_t^{\mathrm{Im}b}}-{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in {\mathsf{\Omega }}_D}{C}_{it}^{IL}}}-{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{BESS}}{C}_{it}^{BESS}}}\\ -{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{DG}}\left(C_{it}^{DG}{\chi }_{it}^{DG}+SU{C}_{it}^{DG}{\chi }_{it}^{SU}+SD{C}_{it}^{DG}{\chi }_{it}^{SD}\right)}}\end{array}\right\}$$

(8)

where ${\rho }_t^{Re tail}$ denotes the retail price agreed between VPP and the customers; $P_{Dt}$ denotes the forecasted load at time $t$; $SU{C}_{it}^{DG}$ and $SD{C}_{it}^{DG}$ denote start-up (usually including cold and hot start-up costs) and shut-down costs of DG, respectively; $T$ denote the total DA scheduling horizon (e.g., 24 h); ${\chi }_{it}^{DG}$, ${\chi }_{it}^{SU}$, and ${\chi }_{it}^{SD}$ are binary variables denoting commitment status, start-up and shut-down decisions for DG at time $t$ bus $i$, respectively; and ${\mathsf{\Omega }}_{BESS}$, ${\mathsf{\Omega }}_{DG}$, and ${\mathsf{\Omega }}_D$ denote sets of BESS, DG and demand, respectively. In Equation (8), the first term represents the cost (revenue) of VPP by purchasing (selling) electricity from (to) the distribution system. The second term represents the revenue by selling electricity to the customers. The remaining terms in Equation (8) represent costs of power imbalance, IL, BESS and DG. The retail price ${\rho }_t^{Re tail}$ in this paper is set according to the residential energy time-of-use (TOU) price in New South Wales, Australia. Specifically, the peak price is 48 cent/kWh, the shoulder price is 19.50 cent/kWh and the off-peak price is 12 cent/kWh.

Given the uncertainties in electricity markets, a decision-making can be risky. Moreover, Equation (8) should be formulated as a probabilistic version. To enhance the computational efficiency, the initial set of scenarios is reduced to a number of representative scenarios using an appropriate scenario reduction technique such as adaptive importance sampling presented in [26]. Moreover, the risk associated with the profit variability is explicitly captured into the model through the mean-variance Markowitz theory [18]. Thus, Equation (8) is rewritten as:

$$\max E\left[O_1\right]+\varpi \cdot {\sigma }_{O_1}$$

(9)

where $E\left[O_1\right]$ denotes the expected profit; ${\sigma }_{O_1}$ denote the standard deviation; $\varpi \in \left[0,+\infty \right)$ is a weighting factor for the incorporation of risk in the objective function [27]. Higher values of $\varpi$ indicate greater risk aversion. Particularly, $\varpi =0$ means an absolute risk-neutral strategy. The calculation of Equation (9) is given in Equations (10) and (11).

$$E\left[O_1\right]={\displaystyle \sum _{k\in {\mathsf{\Omega }}_K}{Pr}_k\cdot f_{1k}}$$

(10)

$${\sigma }_{O_1}^2=E\left[O_1^2\right]-E^2\left[O_1\right]={\displaystyle \sum _{k\in {\mathsf{\Omega }}_K}{Pr}_k\cdot }{f}_{1k}^2-{\left({\displaystyle \sum _{k\in {\mathsf{\Omega }}_K}{Pr}_k\cdot f_{1k}}\right)}^2$$

(11)

where ${Pr}_k$ denotes the probability of scenario $k$; and ${\mathsf{\Omega }}_K$ denotes the set of total scenarios. The probability of scenario ${Pr}_k$ is obtained as follows. We firstly use the ARIMA model to generate correlated uncertainty data such as load, price, and wind speed. The parameters of the ARIMA are obtained based on the historical data in Nord Pool in 2015. Afterwards, the Monte Carlo (MC) simulations are used to generate uncertainty forecast errors produced by the ARIMA model. The probability of each scenario (${Pr}_k$) is obtained based on the histogram of the sampled correlated scenarios in MC simulations. The total number of scenarios is the bin size. The values of scenarios are the mean of samples in the corresponding bin. Lastly, the scenarios are used as model inputs to calculate the objective in Equation (8), which is $f_{1k}$. The larger is the bin size, the more accurate is the result.

For each scenario, the complete constraints of Equation (9) are given below.

3.2.1. Supply-Demand Balance Constraint

$${\displaystyle \sum _{i\in {\mathsf{\Omega }}_{GW}}{P}_{GWkit}}+{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{BESS}}{P}_{kit}^{BESS}+}{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{DG}}{P}_{DGkit}}=P_{Dkt}-{\displaystyle \sum _{i\in {\mathsf{\Omega }}_D}{P}_{kit}^{IL}}+P_{kt}^{DS}+P_{kt}^{Loss},\forall t\in 1:T$$

(12)

where $P_{kt}^{Loss}$ denote power loss of VPP at time $t$ scenario $k$; $P_{GWkit}$, $P_{kit}^{BESS}$, and $P_{DGkit}$ denote outputs of wind power, BESS and DG at time $t$ bus $i$ scenario $k$; $P_{Dkt}$ denotes power demand; $P_{kit}^{IL}$ denotes interrupted load (i.e., incentive-based demand response); and ${\mathsf{\Omega }}_{GW}$ denotes set of wind power.

3.2.2. Wind Power Constraint

$$0\le P_{GWkit}\le {\overline{P}}_{GWi};\forall i\in {\mathsf{\Omega }}_{GW},\forall t\in 1:T$$

(13)

$$\frac{P_{GWkit}}{\sqrt{P_{GWkit}^2+Q_{GWkit}^2}}=\mathrm{constant}$$

(14)

where $\overline{(•)}$ denotes the upper limit; and $Q_{GWkit}$ denotes reactive power output of wind at time $t$ bus $i$ scenario $k$.

3.2.3. DG Constraint

$${\underset{\_}{P}}_{DGi}{\chi }_{kit}^{DG}\le P_{DGkit}\le {\overline{P}}_{DGi}{\chi }_{kit}^{DG},\forall t\in 1:T,\forall i\in {\mathsf{\Omega }}_{DG}$$

(15)

$$\{\begin{array}{l}{P}_{DGkit}-P_{DGkit-1}\le Ram{p}_{DGi}^{Up},if{P}_{DGkit}\ge P_{DGkit-1}\\ P_{DGkit-1}-P_{DGkit}\le Ram{p}_{DGi}^{Down},if{P}_{DGkit-1}\ge P_{DGkit}\end{array}$$

(16)

$$\{\begin{array}{l}\left[T_{kit-1}^{On}-MU{T}_i\right]\cdot \left[{\chi }_{kit-1}^{DG}-{\chi }_{kit}^{DG}\right]\ge 0\\ \left[T_{kit-1}^{Off}-MD{T}_i\right]\cdot \left[{\chi }_{kit-1}^{DG}-{\chi }_{kit}^{DG}\right]\ge 0\end{array}$$

(17)

$$\{\begin{array}{l}{\chi }_{kit}^{DG}-{\chi }_{kit-1}^{DG}\le {\chi }_{kit}^{SU}\\ {\chi }_{kit-1}^{DG}-{\chi }_{kit}^{DG}\le {\chi }_{kit}^{SD}\\ {\chi }_{kit}^{DG}-{\chi }_{kit-1}^{DG}\le {\chi }_{kit}^{SU}-{\chi }_{kit}^{SD}\end{array}$$

(18)

where $Ram{p}_{DGi}^{Up}$ and $Ram{p}_{DGi}^{Down}$ denote ramping up and ramping down limits of DG at bus $i$; $\underset{\_}{(•)}$ and $\overline{(•)}$ denote the lower and upper limits; $T_{kit-1}^{On}$ and $T_{kit-1}^{Off}$ denote number of hours for which DG unit has been on or off at time $t$ bus $i$ scenario $k$; and $MU{T}_i$ and $MD{T}_i$ denote minimum up and down time limits of DG in hours at bus $i$. Equation (17) ensures the minima up and down time limits of DG are satisfied by checking the previous and current operating states of DG. Similarly, in Equation (18), the start-up and shut-down decisions are checked by calculating the state differences of DG in two continuous time intervals.

3.2.4. IL Constraint

$$\frac{\left|P_{kit}^{IL}\right|}{P_{Dkit}}\le \overline{\varsigma },\forall t\in 1:T,\forall i\in {\mathsf{\Omega }}_D$$

(19)

Equation (19) states the maximum ratio of IL, i.e., $\overline{\varsigma }$.

3.2.5. BESS Constraint

$$E_{kit+1}^{BESS}=E_{kit}^{BESS}-P_{kit}^{BESS}\Delta t-\left|P_{kit}^{BESS}\right|{\eta }^C\Delta t-E_{kit}^{BESS}{\eta }^L\Delta t,\forall t\in 1:T,\forall i\in {\mathsf{\Omega }}_{BESS}$$

(20)

$$SO{C}_{kit}=E_{kit}^{BESS}/E_R^{BESS}$$

(21)

$${\underset{\_}{SOC}}_{it}\le SO{C}_{kit}\le {\overline{SOC}}_{it}$$

(22)

$$\{\begin{array}{l}{P}_{kit}^{BESS}\le {\chi }_{kit}^{BESS,Dis}{\overline{P}}_{it}^{BESS,Dis},if{P}_{kit}^{BESS}\ge 0\\ P_{kit}^{BESS}\ge {\chi }_{kit}^{BESS,Chr}{\overline{P}}_{it}^{BESS,Chr},else\end{array}$$

(23)

$${\chi }_{kit}^{BESS,Dis}+{\chi }_{kit}^{BESS,Chr}\le 1$$

(24)

$$\{\begin{array}{l}{E}_{ki0}^{BESS}=E_{iInitial}^{BESS}\\ E_{kiT}^{BESS}\ge E_{iEnd}^{BESS}\end{array}$$

(25)

where ${\eta }^C$ denotes charging/discharging loss factor; ${\eta }^L$ denotes self-discharge loss factor; $E_{kit}^{BESS}$ denotes stored energy in BESS at time $t$ bus $i$ scenario $k$; $E_R^{BESS}$ denotes the rated energy of BESS; and ${\chi }_{it}^{BESS,Dis}$ and ${\chi }_{it}^{BESS,Chr}$ are binary variables denoting discharging and charging decisions of BESS. Equation (20) states energy balance of BESS, including net energy injection, energy losses during charging or discharging, and leakage loss. Equations (21)–(23) represent state-of-charge (SOC) constraints, and maximum charging or discharging power capacity [23]. Constraint (24) means BESS cannot operate in charging and discharging simultaneously. Equation (25) defines the initial and final energy stored in BESS.

3.2.6. Network Constraint

To accurately model the power losses and reactive power, the alternating current (AC) power flow model is used as follow.

$$P_{kit}({\theta }_{kt},V_{kt})-P_{Gkit}+P_{Dkit}=0,\forall t\in 1:T,\forall i\in {\mathsf{\Omega }}_N$$

(26)

$$Q_{kit}({\theta }_{kt},V_{kt})-Q_{Gkit}+Q_{Dkit}=0,\forall t\in 1:T,\forall i\in {\mathsf{\Omega }}_N$$

(27)

$$S_{kijt}({\theta }_{kt},V_{kt})\le {\overline{S}}_{ij},\forall t\in 1:T,\forall i,j\in {\mathsf{\Omega }}_N$$

(28)

$${\underset{\_}{V}}_i\le V_{kit}\le {\overline{V}}_i,\forall t\in 1:T,\forall i\in {\mathsf{\Omega }}_N$$

(29)

where $P_{kit}({\theta }_{kt},V_{kt})$ and $Q_{kit}({\theta }_{kt},V_{kt})$ denote active and reactive power injection time $t$ bus $i$ scenario $k$; $P_{Gkit}$ and $Q_{Gkit}$ denote active and reactive power generation, while $P_{Dkit}$ and $Q_{Dkit}$ denote active and reactive power demand; $S_{kijt}({\theta }_{kt},V_{kt})$ denotes complex power flow between $i-j$ at time $t$, with voltage angle ${\theta }_{kt}$ and amplitude $V_{kt}$; and ${\mathsf{\Omega }}_N$ denotes set of all buses.

3.2.7. Interconnection Exchange Constraint

$${\underset{\_}{P}}_t^{Exch}\le P_{kt}^{DS}\le {\overline{P}}_t^{Exch},\forall t\in 1:T$$

(30)

Equation (30) states that $P_t^{DS}$ is within the upper and lower limits of interconnection exchange between VPP and the upstream grid.

3.2.8. Steady Security Reserve Constraint

$$\begin{array}{l}{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{DG}}{\overline{P}}_{DGi}{\chi }_{it}^{DG}-}{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{DG}}{P}_{DGkit}{\chi }_{it}^{DG}+}{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{GW}}{P}_{GWkit}}\\ +{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{BESS}}{P}_{kit}^{BESS}+{\displaystyle \sum _{i\in {\mathsf{\Omega }}_D}\left({\overline{P}}_{it}^{IL}-P_{kit}^{IL}\right)}}\\ +\left({\overline{P}}_t^{Exch}-P_{kt}^{DS}\right)\ge P_{Dkt}-{\displaystyle \sum _{i\in {\mathsf{\Omega }}_D}{P}_{kit}^{IL}+Rsv(t)}\end{array},\forall t\in 1:T$$

(31)

Equation (31) states the system reserve requirement for static security and adequacy [1]. The system reserve requirement $Rsv(t)$ is assumed to be 1.5 times the total system peak load.

3.3. Real-Time Balancing

The real-time dispatch is to maintain the power balance. Because of the forecast errors of uncertainties, any deviation from the commitment made in the DA market should be balanced. The smaller is the RT dispatch interval, the smaller are the forecast errors. Thus, the frequency of RT clearing is very high (e.g., 5-min). Time-series forecast techniques such as ARIMA are appropriate to capture the persistent behavior of uncertainties [22]. Moreover, the last slight imbalances in each dispatch interval can be handled by automatic generation control and/or emergency DR (i.e., instantaneous control).

The objective of VPP in the RT market is to minimize the imbalance cost and wind power curtailment. In this particular market, we assume that the curtailment of wind power will be published by a monetary penalty coefficient, as given below.

$$\min f_2={\displaystyle \sum _{t=1}^{NT}{\widehat{C}}_t^{\mathrm{Im}b}+}\lambda {\displaystyle \sum _{t=1}^{NT}{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{GW}}\left(P_{GWit}^{Available}-{\widehat{P}}_{GWit}\right)}}$$

(32)

where $NT$ denotes total number of dispatch intervals in the RT market; ${\widehat{C}}_t^{Imb}$ denotes imbalance cost due to production and consumption variations in the RT market; $\lambda$ is the penalty price of wind power curtailment; and $P_{GWit}^{Available}$ is the available wind power at time $t$ bus $i$.

Note that we simply use $\left(\widehat{•}\right)$ to distinguish the variables in the RT market from the variables in the DA market. The detailed calculations of ${\widehat{C}}_t^{Imb}$ are given in Equations (33)–(35).

$$\Delta {\widehat{P}}_t^{Imb}={\widehat{P}}_t^{RT}-P_t^{DS}$$

(33)

$${\widehat{P}}_t^{RT}={\displaystyle \sum _{i\in {\mathsf{\Omega }}_{GW}}{\widehat{P}}_{GWit}}+{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{BESS}}{\widehat{P}}_{it}^{BESS}+}{\displaystyle \sum _{i\in {\mathsf{\Omega }}_{DG}}{\widehat{P}}_{DGit}}-\left({\widehat{P}}_{Dt}+{\widehat{P}}_t^{Loss}-{\displaystyle \sum _{i\in {\mathsf{\Omega }}_D}{\widehat{P}}_{it}^{IL}}\right)$$

(34)

$${\widehat{C}}_t^{Imb}=\{\begin{array}{l}{\rho }_t^{R+}\cdot \Delta {\widehat{P}}_t^{Imb},\Delta P_t^{Imb}<0\\ {\rho }_t^{R-}\cdot \Delta {\widehat{P}}_t^{Imb},\Delta P_t^{Imb}\ge 0\end{array}$$

(35)

The constraints of Equation (32) include Equation (12)–(15), (19)–(23), and (25)–(30).

4. Solution Algorithm

4.1. Solution to the DA Scheduling Problem

The hourly DA scheduling of VPP is a mixed-integer nonlinear programing problem. The model is a here-and-now decision-making process, meaning that decisions are made prior to the knowledge of uncertainty data such as wind speed, load and price data [28]. The decision variables in Equation (9) include $P_t^{DS}$, $P_{it}^{IL}$, $P_{it}^{BESS}$, $P_{DGit}$, ${\chi }_{it}^{DG}$, ${\chi }_{it}^{SU}$, ${\chi }_{it}^{SD}$, ${\chi }_{it}^{BESS,Dis}$, ${\chi }_{it}^{BESS,Chr}$. The mathematical programming methods are not suitable for solving this problem, since the formulated VPP scheduling model is a nonlinear and non-convex optimization problem. By contrast, the heuristic-based optimization methods are free from the problem complexity [1,29]. In this paper, the differential evolution (DE) is employed to solve the upper level optimization model. The computational efficiency of this algorithm has been demonstrated in [30]. The applied solution algorithm is briefly introduced below.

(1) Input data

Input data such as network, wind speed, load, price, etc.

(2) Initialization

Start the DE algorithm, generate an initial population $x^1,\dots x^N$ randomly, and set the objective function value $F{V}^i=f\left(x^i\right)$. Note that, here, $N$ denotes the population size, which is the number of generated candidate solutions (called individuals) in each iteration. The more individuals there are in the population, the more differential vectors are available. As a result, the more directions can be explored. However, the computational complexity per generation increases with the size of the population.

(3) Check constraints and calculate the fitness value

For the population size $i=1,\dots ,N$, perform the following steps:

Section 3.1: Check the operational requirements by running power flow functions. The ON/OFF commitment constraint handling algorithm in [31] is applied to adjust the binary states of DG and BESS to satisfy Constraints (15)–(18) and (20)–(25). Based on the given values of DG output $P_{DGkit}$, BESS output $P_{kit}^{BESS}$ and IL amount $P_{kit}^{IL}$, the power flow function [32] is run to calculate the VPP power losses and the value of exchanged power with the upstream power grid. The exchanged power is equal to the scheduled power $P_{kt}^{DS}$. In addition, other constraints such as the network security and adequacy constraints (Equations (26)–(29)), the interconnection constraint (Equation (30)) and the reserve constraint (Equation (31)) are checked based on running the power flow function. Furthermore, after obtaining the power flow results for all scenarios, the expected profit $E\left[O_1\right]$ and its standard deviation ${\sigma }_{O_1}$ are calculated using (9)–(11).

Section 3.2: If there is any violation of the above-mentioned constraints, a variable penalty is assigned. The variable penalty is defined as a function of the distance from the feasible area [1]. The summation of absolute distances of violated constraints is scaled by a large penalty factor and then is combined with the objective function to constitute the fitness. By doing so, infeasible solutions are not discarded. Instead, their information can be used to guide the algorithm search [1].

(4) Update the best individual and selection

Find and update the best individual of the current generation. This step involves comparing individuals with each other. The fitness sharing technique proposed in [30] is used to select the best one whose fitness value is the highest.

(5) Offspring generation

Produce the offspring generation according to the evolution rule in DE. Typical evolution operations include mutation, recombination and crossover. Detailed mathematical formulations of offspring generation can be found in [30,33].

(6) Termination

If any stopping criterion is met, e.g., the fitness value doses not improve in a number of successive iterations, stop and export the best individual. Otherwise, go back to Step 3.

4.2. Solution to the RT Balancing Problem

After solving the optimal scheduling model in the DA, the ON/OFF states of DERs are determined. VPP needs to reschedule DERs in real-time (e.g., on a five-minute basis) to compensate the energy deviations incurred. The lower level optimal operation model is a nonlinear programming problem with continuous variables. Variables such as ${\widehat{P}}_{it}^{IL}$, ${\widehat{P}}_{it}^{BESS}$, ${\widehat{P}}_{DGit}$ are wait-and-see decisions, meaning that decisions can be made when the information about the future is already known [28]. The flow chart of the applied solution algorithm in two stages is shown in Figure 2. We use the solver AMPL/IPOPT 3.8.0 (900 Sierra Place SE, Albuquerque, NM, USA) to solve the RT balancing. This solver is called Interior Point Optimizer, which is based on a mathematical programming algorithm called the interior point method. This method adopts barrier functions to achieve optimization by going through the middle of the solid defined by the problem rather than around its surface.

5. Case Studies

5.1. Experimental Setting

The proposed VPP scheduling approach is tested on a real 12-bus distribution system of the CSIRO energy center in Newcastle, Australia. The VPP scheduling models are coded in MATLAB^® (R2015b, Chatswood, Australia). As seen in Figure 3, the system is comprised of 4 wind turbines, 3 BESS, 3 DG units, and 11 load buses. The total generation capacity is 6500 kW, including 3000 kW of wind power and 3500 kW of DG. We assume that IL is located at every load bus, and up to 10% of the system load can be adjusted if necessary. The corresponding IL costs are composed of three bands, i.e., early morning (first 6 h), day time (middle 12 h) and evening (last 6 h). The capacity of power interconnection between VPP and the main grid is set as 4000 kW. The power curve of a Vestas V27 wind turbine is used [16]. The rated, cut-in, and cut-out wind speeds are 14.0, 3.5 and 25.0 m/s, respectively. The start-up and shut-down costs of DG are assumed to be $70 and $20, respectively. In the base case, the weighting factor $\varpi$ is set to be 0.1. The correlated scenarios of market price, wind speed, and demand are modeled based on the Nord Pool historical data in 2015, which are publically available on [34]. Note that the data have been scaled down and monetary units are transferred to US$. The historical wind power production data in Denmark are used. The wind power output is calculated using the power-speed curve as [35]:

$$P_{GW}=\{\begin{array}{ll}0,& \upsilon \le {\upsilon }_{In},\upsilon \ge {\upsilon }_{Out}\\ \frac{\upsilon -{\upsilon }_{In}}{{\upsilon }_R-{\upsilon }_{In}}{\overline{P}}_{GW}& {\upsilon }_{In}\le \upsilon \le {\upsilon }_R\\ {\overline{P}}_{GW}& {\upsilon }_R\le \upsilon \le {\upsilon }_{Out}\end{array}$$

(36)

where $P_{GW}$ and $\upsilon$ denote wind power output and wind speed, respectively; ${\upsilon }_{In}$, ${\upsilon }_{Out}$, and ${\upsilon }_R$ denote cut-in, cut-out and rated wind speeds, respectively; and ${\overline{P}}_{GW}$ denotes rated wind power.

5.2. Results and Discussion

To verify the effectiveness of the proposed model, two cases are compared. Case 1: An uncoordinated operation strategy. Individual VPP element operates strategically as a price-taker in the DA market, aiming to maximize their profits. The total profit is the summation of individual’s profit virtually obtained [12]. The imbalance cost is calculated after the realization of scenarios. Case 2: The proposed integrated VPP scheduling approach.

The profit distributions in the DA market for Cases 1 and 2 are illustrated in Figure 4. Compared to Case 1, the expected profit for Case 2 is higher ($4595.50 for Case 1 and $5823.30 for Case 2); while the standard deviation (Std.) for Case 2 is relatively lower. In addition, we use the Value-at-Risk (VaR) at a confidence level of 95% to evaluate the risk of different operational strategies. In other words, the VaR-95% denotes the expected value of the 5% scenarios with lowest profit. As seen in Figure 4, the VaR is $2506.20 and $3485.60 for Cases 1 and 2, respectively. This implies that Case 2 incurs less risk in the DA market.

The iteration convergences of DE for Cases 1 and 2 are illustrated in Figure 5. We can see that the DE solution algorithm requires 381 and 342 iterations for Cases 1 and 2, respectively, corresponding to computational time periods of 1982 and 1671 seconds on a PC with Intel Core i7-6600 CPU @ 2.80 GHZ with 8.00 GB RAM.

Moreover, the energy stored in BESS and the dis/charged power in the DA and RT markets are shown in Figure 6. In both cases, BESS is charged in the early morning (e.g., 2–4 a.m.) or late afternoon (2–4 p.m.). BESS is discharged to satisfy the morning and evening peaks, when electricity prices are relatively higher. In other words, the profits made by BESS are based on the price differences between off-peak and on-peak periods. It is worth mentioning that power deviations in the two markets for Case 1 are smaller than Case 2. For the uncoordinated operation in Case 1, BESS only faces load uncertainty (i.e., load forecast errors). By contrast, for the aggregated VPP operation, BESS is required to cover the risks faced by other types of DERs and hence more power fluctuations are observed in Case 2. For instance, wind availability poses a significant risk to WP and WP is free from the market price volatility (generation cost for VPP is zero). On the contrary, wind availability is not an issue for thermal DG units, but DG operations are greatly influenced by market prices. Meanwhile, load uncertainty has more direct impacts on BESS and IL. Therefore, the coordination of VPP aims to share the risks of individual component, produce less power deviations, and make more profits in the markets.

Furthermore, Figure 7 shows the summation of DERs’ power outputs (WP, BESS, IL and DG), demand, and the energy traded in the DA market. We can see that VPP produces more energy when price and demand are high; while VPP consumes energy in off-peak periods. For instance, the summation of DERs at 4 a.m. is nearly zero for both cases and VPP relies on the imported power from the grid to meet the base demand. Energy generated by WP and DG is mainly used to charge BESS. By contrast, VPP uses DERs (e.g., discharging BESS and using IL) and sells electricity to the grid at 7–8 a.m. and 6–8 p.m. More importantly, the coordination of VPP in Case 2 sells more energy in total, thus making more profits. This indicates that the aggregated production profile in Case 2 can better reduce the uncertainties of DERs and more energy surplus is produced.

The expected profits and imbalance costs in the DA and RT markets for Cases 1 and 2 are given in

Table 1. Profits and costs for VPP components in Case 1.

Cost/Profit	DG	BESS	WP	IL	Total
Expected profit ($)	369.40	2089.50	1679.80	456.80	4595.50
Imbalance cost ($)	32.80	8.70	1426.80	125.20	1593.50
Net profit ($)	336.60	2080.80	253.00	331.60	3002.00

and

Table 2. Profits and costs for VPP components in Case 2

Cost/Profit	DG	BESS	WP	IL	Total
Expected profit ($)	1201.20	2178.90	1856.90	586.30	5823.30
Imbalance cost ($)	31.50	55.00	815.40	184.10	1086.00
Net profit ($)	1169.70	2123.90	1041.50	402.20	4737.30

. It can be seen that more profits can be made for all DERs in Case 2. For DG and WP imbalance costs are lower in Case 2, but BESS and IL incur higher imbalance costs. This is because for the aggregated VPP operation BESS and IL are more frequently used to balance the fluctuations of WP. The imbalance cost of WP is markedly reduced from $1426.80 to $815.40. Overall, the net profits of DERs are much higher in Case 2 than those in Case 1. The total net profit is also higher in Case 2 ($3002.00 in Case 1 and $4737.30 in Case 2).

To demonstrate the impacts of different risk aversion levels, two more cases are added. Cases 1 and 2 are set with a higher weighting factor ($\varpi =0.4$), i.e., Case 1.2 and Case 2.2. Eight indicators are compared in the radar chart in Figure 8. Given the different scales, these indicators are presented in per unit (p.u.) based on the highest values in the four cases, i.e., $5823.30 = 1 p.u. expected profit; $1254.70 = 1 p.u. Std. of profit; $1593.50 = 1 p.u. imbalance cost; $4737.30 = 1 p.u. net profit; 1586.66 kWh = 1 p.u. wind power curtailment (WPC); $3982.30 = 1 p.u. VaR; 50.1463 MWh = 1 p.u. energy traded in DA; 1563.52 kWh = 1 p.u. energy traded in RT. As seen in Figure 8, for a higher weighting factor (more risk-averse), VPP becomes more conservative, thus more energy is traded in the DA market. In the meantime, more WP is curtailed, and VPP is more willing to buy electricity from the grid. Besides, a higher risk-aversion level leads to lower profit, lower Std. of profit, lower imbalance cost, and higher VaR. This reveals the trade-off between profit and risk, and the VPP should make trading decisions based on its risk preference. Moreover, for the uncoordinated strategies (Cases 1 and 1.2), DERs earn less profit but are exposed to higher risks. Therefore, the proposed model can effectively reduce the risks in relation to uncertainties.

The four cases are expanded to another four cases (Cases 1.3, 1.4, 2.3 and 2.4). This is done by increasing the uncertainties by 10% (i.e., the forecast errors of market price, load and wind power are increased by 10%). As seen in

Table 3. Sensitivity analysis for uncoordinated operation.

#	Case 1	Case 1.3	Case 1.2	Case 1.4
Expected profit ($)	4595.50	4265.20	3985.60	3968.50
Std. ($)	1254.70	1365.20	786.20	1058.90
Imbalance cost ($)	1593.50	1896.20	1486.30	1898.20
Net profit ($)	3002.00	2369.00	2499.30	2070.30
WPC (kWh)	1132.65	1685.35	1586.66	1785.85
VaR ($)	2506.20	2235.50	2863.30	2248.50
Energy traded DA (MWh)	43.98	44.26	50.14	51.26
Energy traded RT (kWh)	1563.50	1653.38	1354.61	1458.69

and

Table 4. Sensitivity analysis for integrated VPP operation.

#	Case 2	Case 2.3	Case 2.2	Case 2.4
Expected profit ($)	5823.30	4985.30	5068.00	5012.50
Std. ($)	1209.70	1458.60	893.50	1126.60
Imbalance cost ($)	1086.00	1365.20	689.30	742.50
Net profit ($)	4737.30	3620.10	4378.70	4270.10
WPC (kWh)	956.36	1456.39	1263.28	1268.92
VaR ($)	3485.60	2866.30	3982.30	3569.80
Energy traded DA (MWh)	48.25	48.66	50.12	50.98
Energy traded RT (kWh)	1012.35	1456.36	725.36	996.48

, the increased uncertainties result in higher imbalance costs and lower VaR for both operation strategies, which means higher risks. This implies that the VPP should use more sophisticated forecast techniques to hedge against the profit variation risks if increased uncertainties are expected (e.g., increased WP penetration in the future or increased market volatility). Moreover, in a more uncertain environment, WP is less efficiently used (e.g., WPC is increased from 956.36 kWh in Case 2 to 1456.39 kWh in Case 2.3) and more energy is traded in the RT market. It is worth mentioning that for a more conservative strategy (cases 2.2 and 2.4), the reduction in the net profit caused by the increased uncertainties is the least (the net profits are $4378.70 in Case 2.2 and $4270.10 in Case 2.4). Besides, the net profit in Case 2.4 ($4270.10) is higher than that in Case 2.3 ($3620.10). This proves that the profit reduction due to risk-aversion is minimized by adopting a more conservative integrated VPP operational strategy. In other words, the cost in relation to risk-aversion is effectively reduced.

To investigate the correlation between risk aversion and uncertainty parameters, four more cases are added by combing the low and high values of risk aversion and uncertainty parameters. Specifically, low and high risk aversion values are set at $\varpi =0.05$ and $\varpi =0.7$, respectively, while low and high uncertainty values are set by decreasing and increasing the uncertainties by 10%. Case 2.5 is defined as a low risk aversion and low uncertainty case; Case 2.6 is defined as a low risk aversion and high uncertainty case; Case 2.7 is defined as a high risk aversion and low uncertainty case; and Case 2.8 is defined as a high risk aversion and high uncertainty case. The detailed results of cases 2.5–2.8 are given in

Table 5. Sensitivity analysis for risk aversion and uncertainty parameters.

#	Case 2.5	Case 2.6	Case 2.7	Case 2.8
Expected profit ($)	6525.60	6023.30	4612.70	4472.10
Std. ($)	1337.40	1886.30	986.10	1209.70
Imbalance cost ($)	1093.30	1532.10	832.30	1224.50
Net profit ($)	5432.30	4491.20	3780.40	3247.60
WPC (kWh)	655.18	1308.62	652.32	756.48
VaR ($)	4002.40	3026.20	4523.90	3369.20
Energy traded DA (MWh)	55.26	51.25	42.36	41.23
Energy traded RT (kWh)	701.92	1512.35	422.18	1217.36

. Based on the sensitivity analysis of risk aversion and uncertainty parameters, three findings can be identified in Table 5. First, a high risk aversion value results in marked reductions in expected profit, net profit, and energy trades in both DA and RT markets. This is because that if the VPP is more conservative, the VPP is more likely to use the DERs such as DG and BESS to meet the local demands, rather than trading energy in the markets. Second, a higher uncertainty value results in relatively smaller reductions in expected profit, net profit, and energy trades in the DA market. However, a higher uncertainty value results in marked increases in imbalance cost, wind power curtailment and energy trades in the RT market. This implies that the higher uncertainty value is likely to increase the inaccuracy of DA scheduling plans of the VPP, and hence more energy imbalance is traded. Third, the risk aversion and uncertainty parameters can influence each other. A positive correlation is observed between them, since both higher risk aversion and uncertainty parameters can result in reductions in the objective value (i.e., the net profit).

Furthermore, the proposed VPP scheduling approach is compared to another approach proposed in [16]. As seen in

Table 6. Result comparison to an approach in the literature.

#	Proposed Approach	Approach in [16]
Expected profit ($)	5823.30	5172.70
Std. ($)	1209.70	1463.30
Imbalance cost ($)	1086.00	1073.60
Net profit ($)	4737.30	4099.10
WPC (kWh)	956.36	1726.39
VaR ($)	3485.60	2466.30
Energy traded DA (MWh)	48.25	45.18
Energy traded RT (kWh)	1012.35	1004.88

, the net profits of the VPP are $4737.30 and $4099.10 for the two approaches respectively. This proves that the proposed approach is superior in terms of optimality. More importantly, since the approach in [16] does not include any risk measure in the model, the risk of the VPP is marked higher (VaR for the approach in [16] is $2466.30 ). Note that a lower VaR means a higher risk of profit variability. Using the approach in [16], the VPP cannot understand the trade-off between risk and profit, and identifies a scheduling solution that entails a higher market risk. In addition, the renewable energy is not efficiently utilized, as the wind power curtailment is 1726.39 kWh for the approach in [16].

6. Conclusions

This paper presents an optimal scheduling model for a price-taker VPP participating in the DA and RT markets. In the DA market, the VPP maximizes the expected profit by determining the unit commitments and hourly scheduling of the DERs in its portfolio. To capture the risk of low profit scenarios, a risk factor based on the mean-variance Markowitz theory is incorporated into the objective function. In the RT market, the VPP minimizes the imbalance cost and wind power curtailment by adjusting the outputs of DERs on a five-minute basis. According to the simulation results in case studies, the VPP can help cost-efficient integration of DERs, and the integration is an effective risk-hedging mechanism. Specifically, the aggregated VPP scheduling approach can make more profit surplus, mitigate the impact of uncertainties, and reduce the cost of risk-aversion.