Locational Marginal Pricing Mechanism for Uncertainty Management Based on Improved Multi-ellipsoidal Uncertainty Set

Zongzheng Zhao，Yixin Liu，Li Guo，Linquan Bai，Chengshan Wang; Zongzheng Zhao; Yixin Liu; Li Guo; Linquan Bai; Chengshan Wang

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Locational Marginal Pricing Mechanism for Uncertainty Management Based on Improved Multi-ellipsoidal Uncertainty Set PDF

- ORCID：
Zongzheng Zhao
✉
- ORCID：
Yixin Liu
✉
- ORCID：
Li Guo
✉
- ORCID：
Linquan Bai
✉
- ORCID：
Chengshan Wang
✉

Key Laboratory of Smart Grid of Ministry of Education, Tianjin University, Tianjin 300072, China； Department of Systems Engineering and Engineering Management, University of North Carolina at Charlotte, Charlotte, NC 28223, USA

Updated：2021-08-02

DOI：10.35833/MPCE.2020.000824

Abstract

The large-scale integration of renewable energy sources (RESs) brings huge challenges to the power system. A cost-effective reserve deployment and uncertainty pricing mechanism are critical to deal with the uncertainty and variability of RES. To this end, this paper proposes a novel locational marginal pricing mechanism in day-ahead market for managing uncertainties from RES. Firstly, an improved multi-ellipsoidal uncertainty set (IMEUS) considering the temporal correlation and conditional correlation of wind power forecasting is formulated to better capture the uncertainty of wind power. The dimension of each ellipsoidal subset is optimized based on a comprehensive evaluation index to reduce the invalid region without large loss of modeling accuracy, so as to reduce the conservatism. Then, an IMEUS-based robust unit commitment (RUC) model and a robust economic dispatch (RED) model are established for the day-ahead market clearing. Both the reserve cost and ramping constraints are considered in the overall dispatch process. Furthermore, based on the Langrangian function of the RED model, a new locational marginal pricing mechanism is developed. The uncertainty locational marginal price (ULMP) is introduced to charge the RES for its uncertainties and reward the generators who provide reserve to mitigate uncertainties. The new pricing mechanism can provide effective price signals to incentivize the uncertainty management in the day-ahead market. Finally, the effectiveness of the proposed mechanism is verified via numerous simulations on the PJM 5-bus system and IEEE 118-bus system.

Keywords

Day-ahead market; ellipsoidal uncertainty set; locational marginal pricing; reserve; robust unit commitment

I. Introduction

ALTHOUGH the large-scale integration of renewable energy sources (RESs) alleviates the environmental pressure, they bring great challenges to the power system [

1], [2]. The uncertainty and variability of RES need to be effectively managed in the market operation, calling for new pricing mechanisms for uncertainty management in the electricity markets.

The security-constrained economic dispatch (SCED) model is used in the current market-clearing to optimally dispatch the generators and provide price signals, i.e., the locational marginal price (LMP), to market participants [

3]. The LMP is composed of three components: energy price, congestion price, and loss price. The traditional LMP does not account for uncertainties in the system. However, in today’s electricity market with high penetrations of RESs, it is essential to formulate a pricing mechanism considering uncertainty to guarantee the security and economy of system operation. Recently, some researchers have utilized SCED to design and derive the LMP mechanism for pricing system uncertainty. The SCED model in [4] is formulated as a two-stage stochastic programming problem, where the first stage clears the market and the second stage models the system operation under wind power uncertainty. A marginal pricing mechanism including pool energy prices and balancing energy prices is established. Reference [5] designs a chance-constrained stochastic market, which is capable of providing effective price signals that internalize the uncertainty of renewable generation resources and risk tolerance of the market operators. Reference [6] proposes an energy and reserve pricing mechanism that takes into account the RES stochasticity with the intention to produce more accurate signals to market participants. However,

[4] - [6]

consider a single-period market clearing, ignoring the impact of time-coupling constraints such as ramp rate on market prices. Reference [7] proposes a novel market framework to credit the generation and reserve and to charge the load and uncertainty by the SCED model in the day-ahead market. The thermal generators offer zero prices for their reserve products, and consequently, the marginal prices for clearing reserve in some periods are zero. In addition, this paper considers the ramping constraints in both the dispatch and redispatch processes. Reference [8] puts forward a stochastic market-clearing SCED model with an energy-only pricing scheme to yield LMPs. But it is unable to provide the pricing information for both the reserve and uncertainty. Reference [9] derives an uncertainty-contained LMP (U-LMP), which includes two new uncertainty components, i.e., transmission line overload price and generation violation price. Unlike traditional LMP, the U-LMP differs between different generators and loads even at the same bus due to their uncertainty levels, which exacerbates the difficulty and complexity of market clearing.

On the other hand, the security-constrained unit commitment (SCUC) model is essential to determine the startup and shutdown statuses of thermal generators before SCED calculation. Both robust SCUC (RSCUC) and robust SCED (RSCED) have been widely used because they do not require the detailed probability distribution of uncertain variables or the complex calculation with a large number of scenarios. One of the key obstacles of the robust optimization (RO) methods is the uncertainty set modeling, which directly affects the economy and robustness of the decision-making. Common uncertainty sets include box set, polyhedral set, and ellipsoidal set [

10]. The box set utilizes upper and lower bounds to describe the uncertainty, but ignores the correlation widely existing in RES generation and load demand [11]. As a result, the box set may contain numerous scenarios with extremely low occurrence probability, aggravating the conservativeness of RO model. The polyhedral set is able to capture the correlation, but still contains large invalid space, which cannot effectively reduce the conservativeness of RO model. Compared with the above two uncertainty sets, the ellipsoidal set can better preserve the correlation between variables, and has been proved to be effective in improving modeling accuracy and reducing conservativeness Reference [12] proposes a novel RSCUC model, where ellipsoidal uncertainty set (EUS) is adopted to well fit the correlated wind power, and a novel criterion for budget value selection of the EUS is presented to reduce the conservatism. Reference [13] presents a risk-averse two-stage optimization model. The first-stage problem minimizes the dispatch cost using the wind power forecasting. The second-stage problem minimizes the expected redispatch cost in the worst-case, where the minimum volume enclosing ellipsoid (MVEE) algorithm is used to construct the uncertainty set of correlated wind power. Reference [14] proposes an adjustable RSCED model with wind power uncertainty. The MVEE is employed as a convex hull to solve the scenario-based RO model. However, [13] and [14] ignore the temporal correlation between wind power, which may lead to an extreme climbing phenomenon of thermal generators [12], [15]. Furthermore, the temporal correlation of uncertain variables, e.g., wind power as well as its forecasting error, are normally inversely proportional to their time interval distance. As a result, the high-dimensional ellipsoidal set, e.g., a 24-dimensional ellipsoidal set used in the conventional day-ahead SCUC and SCED problems, may contain a large number of uncertain variables with weak correlation, increasing the volume and hindering the reduction of conservatism of the ellipsoidal set.

In view of the above shortcomings, this paper proposes an RO framework for day-ahead market-clearing based on an improved multi-ellipsoidal uncertainty set (IMEUS). A novel locational marginal pricing mechanism for pricing and managing uncertainties in the electricity market is built. In this work, the IMEUS is proposed to better characterize the uncertainties of wind power to reduce the conservatism of RO. Based on the proposed IMEUS, an RSCUC model and an RSCED model are established to optimize the dispatch scheme for thermal generators and generate price signals for energy and reserve, respectively. The main contributions of this paper are summarized as follows.

1) An IMEUS modeling method that comprehensively considers the temporal correlation of the forecasting errors and the conditional correlation between the forecasting errors and the forecasting values is proposed. Two indexes named integrity index and efficiency index are established to evaluate the performance and determine the optimal dimension of each ellipsoid set. The integrity index is defined as the coverage ratio of the IMEUS to the wind power data, which is used to ensure that more actual data is involved in the IMEUS; while the efficiency index reflects the volume of IMEUS, which is adopted to reduce the invalid region of IMEUS, so as to reduce the conservativeness of the uncertainty set.

2) A day-ahead IMEUS-based RSCUC model and an RSCED model are established. The generation cost for reserve provision of thermal generators and the ramping constraints in the overall dispatch process are considered to guarantee the economic and safe operation of the generators and power system. The optimized power output and reserve capacity of thermal generators can meet the energy and uncertainty demands, and avoid waste of flexible resources.

3) A novel locational marginal pricing mechanism is developed. The LMP is used for pricing energy. The uncertainty locational marginal price (ULMP) is introduced to charge the RES and load for their uncertainties and reward the generators who provide reserve to mitigate uncertainties. This pricing mechanism can provide effective price signals to guarantee the cost recovery of thermal generators and incentivize the uncertainty management.

The remainder of this paper is organized as follows. Section II proposes the IMEUS model. Section III introduces the IMEUS-based RSCUC model. Section IV presents the new locational marginal pricing mechanism. Section V presents the simulation results. Section VI concludes this paper with major findings.

II. IMEUS Model

The uncertainty set modeling of day-ahead wind power includes two stages. The samples denoting the possible realization of wind power in the next day are generated at the first stage, and the Gaussian Copula approach is utilized to generate the samples to consider the temporal correlation of the forecasting error as well as the conditional correlation between the forecasting error and the forecasting value. Then, the generated samples are adopted to construct the IMEUS at the second stage.

A. Gaussian Copula-based Wind Power Sampling

Suppose $T$ is the number of scheduling periods that is usually considered as 24 in day-ahead decision-making with time resolution of 1 hour. Define $x = [x_{1}, x_{2}, \dots, x_{T}]$ , $y = [y_{1}, y_{2}, \dots, y_{T}]$ , and $e = [e_{1}, e_{2}, \dots, e_{T}]$ as the historical data sets of actual value, forecasting value, and forecasting error of wind power, respectively, and $x = y + e$ . The joint distribution of forecasting error and forecasting value can be written as (1) according to Gaussian Copula theory [

16].

\begin{array}{l} f_{e, y} (e, y) = ϕ_{R} (z_{x, 1}, z_{x, 2}, \dots, z_{x, T}, z_{y, 1}, z_{y, 2}, \dots, z_{y, T} | R) \cdot \\ \prod_{t = 1}^{T} f_{x_{t}} (y_{t} + e_{t}) f_{y_{t}} (y_{t}) \end{array}

(1)

where $z_{x, t} = Φ_{0}^{- 1} (F_{x, t} (y_{t} + e_{t}))$ , $z_{y, t} = Φ_{0}^{- 1} (F_{y, t} (y_{t}))$ , $t = 1,2, \dots, T$ ; $f (\cdot)$ is the probability density function (PDF); $F (\cdot)$ is the cumulative distribution function (CDF); $Φ_{0}^{- 1}$ denotes the inverse of the CDF of the standard Gaussian distribution $Φ_{0} (\cdot)$ ; $F_{x, t} (y_{t} + e_{t})$ and $F_{y, t} (y_{t})$ can be calculated from historical data; and $ϕ_{R} (\cdot)$ is the PDF of standard Gaussian distribution with covariance matrix R.

Let $z_{x} = [z_{x, 1}, z_{x, 2}, \dots, z_{x, T}]^{T}$ , $z_{y} = [z_{y, 1}, z_{y, 2}, \dots, z_{y, T}]^{T}$ , and $z = [z_{x}^{T}, z_{y}^{T}]^{T}$ obeys standard Gaussian distribution, i.e., $z ~ N (0, R)$ . R can be calculated as:

R = [\begin{matrix} R_{x x} & R_{x y} \\ R_{y x} & R_{y y} \end{matrix}]

(2)

R_{x x} = [\begin{matrix} 1 & ρ (z_{x, 1}, z_{x, 2}) & \dots & ρ (z_{x, 1}, z_{x, T}) \\ ρ (z_{x, 2}, z_{x, 1}) & ρ (z_{x, 2}, z_{x, 2}) & \dots & ρ (z_{x, 2}, z_{x, T}) \\ ⋮ & ⋮ & ⋮ \\ ρ (z_{x, T}, z_{x, 1}) & ρ (z_{x, T}, z_{x, 2}) & \dots & 1 \end{matrix}]

(3)

R_{x y} = R_{y x}^{T} = [\begin{matrix} ρ (z_{x, 1}, z_{y, 1}) & ρ (z_{x, 1}, z_{y, 2}) & \dots & ρ (z_{x, 1}, z_{y, T}) \\ ρ (z_{x, 2}, z_{y, 1}) & ρ (z_{x, 2}, z_{y, 2}) & \dots & ρ (z_{x, 2}, z_{y, T}) \\ ⋮ & ⋮ & ⋮ \\ ρ (z_{x, T}, z_{y, 1}) & ρ (z_{x, T}, z_{y, 2}) & \dots & ρ (z_{x, T}, z_{y, T}) \end{matrix}]

(4)

R_{y y} = [\begin{matrix} 1 & ρ (z_{y, 1}, z_{y, 2}) & \dots & ρ (z_{y, 1}, z_{y, T}) \\ ρ (z_{y, 2}, z_{y, 1}) & ρ (z_{y, 2}, z_{y, 2}) & \dots & ρ (z_{y, 2}, z_{y, T}) \\ ⋮ & ⋮ & ⋮ \\ ρ (z_{y, T}, z_{y, 1}) & ρ (z_{y, T}, z_{y, 2}) & \dots & 1 \end{matrix}]

(5)

where $ρ (z_{x, i}, z_{y, j}) = 2 s i n (ρ_{r} (z_{x, i}, z_{y, j}) π / 6), i, j = 1,2, \dots, T$ ; $ρ (\cdot)$ and $ρ_{r} (\cdot)$ are the linear correlation coefficient and Spearman correlation coefficient, respectively [

17]. Similarly,

ρ (z_{x, i}, z_{x, j})

and

ρ (z_{y, i}, z_{y, j})

can be obtained in the same way.

The conditional PDF of forecasting error can be formulated as (6) based on (1).

\begin{array}{l} f_{e | y} (e | y) = \frac{f_{e, y} (e, y)}{f_{y} (y)} = \frac{ϕ_{R} (z_{x}^{T}, z_{y}^{T})}{ϕ_{R_{y y}} (z_{y}^{T})} \prod_{t = 1}^{T} f_{x_{t}} (y_{t} + e_{t}) = \\ f_{z_{x} | z_{y}} (z_{x}^{T} | z_{y}^{T}) \prod_{t = 1}^{T} f_{x_{t}} (y_{t} + e_{t}) \end{array}

(6)

where $z_{y} ~ N (0, R_{y y})$ ; and $z_{x} | z_{y}$ is a conditional distribution of the multivariate Gaussian distribution, which obeys Gaussian distribution, i.e., $z_{x} | z_{y} ~ N (μ_{z_{x} | z_{y}}, R_{z_{x} | z_{y}})$ . The expected value vector $μ_{z_{x} | z_{y}}$ and covariance matrix $R_{z_{x} | z_{y}}$ can be obtained as [

18]:

\{\begin{array}{l} μ_{z_{x} | z_{y}} = R_{x y} R_{y y}^{- 1} z_{y} \\ R_{z_{x} | z_{y}} = R_{x x} - R_{x y} R_{y y}^{- 1} R_{y x} \end{array}

(7)

$z_{x} | z_{y}$ characterizes the temporal correlation and conditional correlation between the actual value and forecasting value of wind power. Accordingly, the wind power samples for the next day can be generated based on the latest day-ahead forecasting value. The specific sampling steps are as follows.

1) Obtain CDFs of the actual value and forecasting value of wind power in each time period based on the historical data. Calculate the covariance matrix $R$ via (2)-(5).

2) Calculate $z_{y, t} = Φ_{0}^{- 1} (F_{y, t} (y_{t}))$ based on the latest day-ahead forecasting value of wind power.

3) Calculate the expected value $μ_{z_{x} | z_{y}}$ and covariance matrix $R_{z_{x} | z_{y}}$ of $z_{x} | z_{y}$ through (7).

4) Generate samples of $z_{x}$ by sampling $z_{x} | z_{y}$ and then get the actual value samples $x_{s}$ of wind power through $x_{s, t} = F_{x, t}^{- 1} (Φ_{0} (z_{x, t}))$ .

The above sampling method generates day-ahead actual value samples $x_{s} = [x_{s, 1}, x_{s, 2}, \dots, x_{s, T}]^{T}$ according to the conditional distribution obtained from the historical data and the latest forecasting value. It not only reveals the temporal correlation and conditional correlation of wind power, but also updates the day-ahead samples based on the latest forecasting, showing good adaptability to wind power variation and better modeling accuracy.

B. IMEUS Model

An ellipsoidal uncertainty set (EUS) covering the samples $x_{s}$ with confidence degree $α_{c}$ is formulated as:

(x_{s} - μ_{s})^{T} R_{s}^{- 1} (x_{s} - μ_{s}) \leq C_{α}

(8)

where $μ_{s}$ and $R_{s}$ are the expectation vector and covariance matrix of samples $x_{s}$ , respectively; and $C_{α}$ is a constant corresponding to $α_{c}$ .

$μ_{s}$ and $R_{s}$ can be easily calculated from samples $x_{s}$ . Take the samples into (9) to obtain the CDF of $C$ . Let $P (C \leq C_{α}) = α_{c}$ , and then we can get $C_{α}$ corresponding to $α_{c}$ .

C = (x_{s} - μ_{s})^{T} R_{s}^{- 1} (x_{s} - μ_{s})

(9)

$μ_{s}$ , $R_{s}$ , and $C_{α}$ determine the center, shape, and volume of the EUS, respectively. When $C_{α}$ is large, the EUS covers more samples with a big volume, which may increase the conservativeness of the uncertainty set.

The traditional method considers a T-dimensional EUS to model the temporal correlation of wind power. However, although the high-dimensional EUS is more likely to cover the actual realization of wind power, the weak correlation among distant time periods makes EUS too conservative due to its large volume. On the contrary, a low-dimensional EUS has smaller volume but may not include enough wind power data. As a result, a trade-off between wind power coverage capacity and conservatism should be made in an optimal manner. To this end, the integrity index and efficiency index are put forward to optimize the dimension of EUS. The rolling modeling process is shown in Fig. 1. Firstly, the T-dimensional ellipsoidal set is decomposed into several low-dimensional ellipsoidal subsets chronologically. Define $N_{s}$ as the number of subsets and T_D as the optimal dimension of EUS, $N_{s} = T - T_{D} + 1$ . Define s_n as the index for the serial number of subsets. The $s_{n}^{t h}$ subset $Ω_{s_{n}}$ expressed as (10) considers the temporal correlation among time interval $[s_{n}, s_{n} + T_{D} - 1]$ , $1 \leq s_{n} \leq N_{s}$ , and is used to involve the possible realization of wind power among this time interval. Then, the intersection of these subsets is taken as the IMEUS $Ω$ , which is formulated as (11).

Ω_{s_{n}} = {x_{s_{n}} | (x_{s_{n}} - μ_{s, s_{n}})^{T} R_{s, s_{n}}^{- 1} (x_{s_{n}} - μ_{s, s_{n}}) \leq C_{α, s_{n}}}

(10)

Ω = \cap_{s_{n} = 1}^{N_{s}} Ω_{s_{n}}

(11)

Fig. 1 Modeling process of IMEUS.

The integrity index $ζ$ is defined as the average coverage ratio of the IMEUS $Ω$ to the actual data of wind power for D days, which is used to ensure that more actual data is contained in the IMEUS. If the maximum number of time periods when the IMEUS of the $d^{t h}$ day $(1 \leq d \leq D)$ can cover the actual wind power is $N_{d}^{m a x}$ , $ζ$ can be expressed as:

ζ = \frac{1}{D} \sum_{d = 1}^{D} \frac{N_{d}^{m a x}}{T}

(12)

The efficiency index $η$ evaluates the volume of IMEUS $Ω$ by referring to the box set $Ω_{b o x}$ shown in (13). This index limits the volume of IMEUS to reduce the conservativeness.

Ω_{b o x} = {x_{t} | x_{t} \in [x_{t}^{m i n}, x_{t}^{m a x}]}

(13)

To estimate the volume of IMEUS, $N_{b o x, d}$ samples are generated in box set $Ω_{b o x, d}$ of the $d^{t h}$ day, and the number of these samples in $Ω_{d}$ is recorded as $N_{e l l, d}$ . The efficiency index $η$ is defined as:

η =

1 - \frac{l g (\frac{1}{D} \sum_{d = 1}^{D} N_{e l l, d})}{l g (\frac{1}{D} \sum_{d = 1}^{D} N_{b o x, d})}

(14)

Taking into account the two indexes comprehensively, the comprehensive index $C I$ is defined as:

C I = k ζ + (1 - k) η k \in [0,1]

(15)

Since there are few potential solutions, the exhaustive method is utilized to optimize $T_{D}$ . The pseudocode for the process of optimizing $T_{D}$ is demonstrated in Algorithm 1.

Algorithm 1 : process of optimizing $T_{D}$

Input: wind power samples $x_{s}$ obtained by the method presented in Section II-A, the confidence degree $α_{c}$ which is used to calculate $C_{α, s_{n}}$ in (10), and the coefficient k in (15)

Output: the optimal value of the dimension $T_{D}$ of each low-dimensional ellipsoidal set

1: for $T_{D} = 2 : T$ do

2: for $d = 1 : D$ do

3: Establish all low-dimensional ellipsoidal sets $Ω_{s_{n}}$ shown in (10)

4: Take the intersection of these low-dimensional ellipsoidal sets to form the IMEUS $Ω$ formulated as (11)

5: end for

6: Calculate the integrity index $ζ$ and the efficiency index $η$ which are shown in (12) and (14), respectively

7: Calculate the comprehensive index $C I$ which is shown in (15)

8. end for

9: Select $T_{D}$ with the maximum value of CI as the solution

Through Cholesky decomposition $R_{s, s_{n}}^{- 1} = L_{s, s_{n}}^{T} L_{s, s_{n}}$ , (10) can be transformed into the following form:

Ω_{s_{n}} = \{x_{s_{n}} |(x_{s_{n}} - μ_{s, s_{n}})^{T} L_{s, s_{n}}^{T} L_{s, s_{n}} (x_{s_{n}} - μ_{s, s_{n}}) \leq C_{α, s_{n}}\}

(16)

This formula can be further converted into:

Ω_{s_{n}} = \{x_{s_{n}}| | | C_{α, s_{n}}^{- 1 / 2} L_{s, s_{n}} (x_{s_{n}} - μ_{s, s_{n}}) | |_{2} \leq 1\}

(17)

Replace $x_{s_{n}} = [x_{s_{n}}, x_{s_{n} + 1}, \dots, x_{s_{n} + T_{D} - 1}]^{T}$ with the day-ahead wind power $P_{s_{n}}^{w, I M E U S} = [P_{s_{n}}^{w, I M E U S}, P_{s_{n} + 1}^{w, I M E U S}, \dots, P_{s_{n} + T_{D} - 1}^{w, I M E U S}]^{T}$ , and (17) can be converted into (18) to realize the uncertainty modeling of day-ahead wind power.

Ω_{s_{n}} = \{P_{s_{n}}^{w, I M E U S}| | | C_{α, s_{n}}^{- 1 / 2} L_{s, s_{n}} (P_{s_{n}}^{w, I M E U S} - μ_{s, s_{n}}) | |_{2} \leq 1\}

(18)

In order to realize the adjustability of the robustness, the IMEUS is finally formulated as:

\begin{array}{l} P^{w, I M E U S} : = \{\begin{array}{l} P_{t}^{w, I M E U S}| | | C_{α, s_{n}}^{- 1 / 2} L_{s, s_{n}} (P_{s_{n}}^{w, I M E U S} - μ_{s, s_{n}}) | |_{2} \leq 1 \\ P_{t}^{w, I M E U S} \geq B_{w, t} P_{f, t}^{w} \\ \sum_{t = 1}^{T} B_{w, t} \geq T - Γ^{w} \end{array} \\ t = 1,2, \dots, T, s_{n} = 1,2, \dots, N_{s} \end{array}

(19)

where $P_{f, t}^{w}$ is the day-ahead forecasting value of wind power at time t; and $B_{w, t}$ is a binary variable related to the realization of wind power. In the worst-case scenario of RO, wind power takes the lower bound of IMEUS at time t when $B_{w, t} = 0$ , otherwise, takes the forecasting value $P_{f, t}^{w}$ at time t when $B_{w, t} = 1$ . $Γ^{w}$ is the uncertainty budget denoting the maximum number of periods when wind power is taken at the lower bound of the uncertainty set, which is an integer value between 0 and T. Therefore, the solution is more conservative when $Γ^{w}$ is bigger.

Considering the fact that the forecasting accuracy of load demand is higher than that of wind power, the uncertainty of load is addressed only by box sets in this paper, which can be expressed as:

P^{d, b o x} : = \{\begin{array}{l} P_{t}^{d, b o x} |P_{t}^{d, b o x} = P_{f, t}^{d} + B_{d, t} Δ P_{t}^{d, b o x} \\ \sum_{t = 1}^{T} B_{d, t} \leq Γ^{d} \end{array}

(20)

where $P_{t}^{d, b o x}$ , $P_{f, t}^{d}$ , and $Δ P_{t}^{d, b o x}$ are the uncertain variables related to load demand, forecasting value of load, and the deviation between the load demand and the forecasting value at time $t$ , respectively, $P^{d, b o x} = [P_{1}^{d, b o x}, P_{2}^{d, b o x}, \dots, P_{T}^{d, b o x}]^{T}$ ; $B_{d, t}$ is a binary variable related to the realization of load; and $Γ^{d}$ is the uncertainty budget, denoting the maximum number of periods when load is taken at the upper bound of the uncertainty set.

III. RSCUC

A. RSCUC Model

The RSCUC model optimizes the dispatch scheme of energy and reserve with the lowest cost in the worst-case scenario, and provides unit commitment status and the worst-case uncertainty realization for RSCED calculation. It contains basic dispatch and redispatch processes in this paper. In basic dispatch process, thermal generators provide energy according to the forecasting values of load and wind power. In the redispatch process, reserve is optimized to cope with the deviation between the forecasting value and the worst-case realization of load and wind power. The RSCUC model is formulated as:

\underset{u, v, I, P^{g}, Δ P^{g}}{m i n} \sum_{t = 1}^{T} \sum_{i = 1}^{N_{G}} [C_{i}^{g} (P_{i, t}^{g} + Δ P_{i, t}^{g}) + C_{i}^{S U} u_{i, t} + C_{i}^{S D} v_{i, t}]

(21)

s.t.

\sum_{i = 1}^{N_{G}} P_{i, t}^{g} + \sum_{j = 1}^{N_{W}} P_{j, f, t}^{w} = \sum_{m = 1}^{N_{M}} P_{m, f, t}^{d} \forall t (λ_{t}^{b})

(22)

- P_{i, t}^{g} \geq - I_{i, t} P_{i, m a x}^{g} \forall i, t ({\bar{β}}_{i, t}^{b})

(23)

P_{i, t}^{g} \geq I_{i, t} P_{i, m i n}^{g} \forall i, t ({\underset{̲}{β}}_{i, t}^{b})

(24)

P_{i, t - 1}^{g} - P_{i, t}^{g} \geq - r_{i}^{r u} (1 - u_{i, t}) - r_{i}^{s u} u_{i, t} \forall i, t ({\bar{α}}_{i, t}^{b})

(25)

P_{i, t}^{g} - P_{i, t - 1}^{g} \geq - r_{i}^{r d} (1 - v_{i, t}) - r_{i}^{s d} v_{i, t} \forall i, t ({\underset{̲}{α}}_{i, t}^{b})

(26)

- \sum_{m = 1}^{N_{M}} G S F_{l, m} \cdot P_{m, t}^{i n j} \geq - F_{l} \forall l, t ({\bar{η}}_{l, t}^{b})

(27)

\sum_{m = 1}^{N_{M}} G S F_{l, m} \cdot P_{m, t}^{i n j} \geq - F_{l} \forall l, t ({\underset{̲}{η}}_{l, t}^{b})

(28)

\sum_{i = 1}^{N_{G}} Δ P_{i, t}^{g} = \sum_{m = 1}^{N_{M}} Δ P_{m, t}^{d} - \sum_{j = 1}^{N_{W}} Δ P_{j, t}^{w} (λ_{t}^{r})

(29)

- P_{i, t}^{g} - Δ P_{i, t}^{g} \geq - I_{i, t} P_{i, m a x}^{g} \forall i, t ({\bar{β}}_{i, t}^{r})

(30)

P_{i, t}^{g} + Δ P_{i, t}^{g} \geq I_{i, t} P_{i, m i n}^{g} \forall i, t ({\underset{̲}{β}}_{i, t}^{r})

(31)

- Δ P_{i, t}^{g} \geq - r_{i}^{r u} (1 - u_{i, t}) - r_{i}^{s u} u_{i, t} \forall i, t ({\bar{α}}_{i, t}^{r 1})

(32)

Δ P_{i, t}^{g} \geq - r_{i}^{r d} (1 - v_{i, t}) - r_{i}^{s d} v_{i, t} ({\underset{̲}{α}}_{i, t}^{r 1})

(33)

P_{i, t - 1}^{g} + Δ P_{i, t - 1}^{g} - P_{i, t}^{g} - Δ P_{i, t}^{g} \geq - r_{i}^{r u} (1 - u_{i, t}) - r_{i}^{s u} u_{i, t} \forall i, t ({\bar{α}}_{i, t}^{r 2})

(34)

P_{i, t}^{g} + Δ P_{i, t}^{g} - P_{i, t - 1}^{g} - Δ P_{i, t - 1}^{g} \geq - r_{i}^{r d} (1 - v_{i, t}) - r_{i}^{s d} v_{i, t} \forall i, t ({\underset{̲}{α}}_{i, t}^{r 2})

(35)

- \sum_{m = 1}^{N_{M}} G S F_{l, m} \cdot (P_{m, t}^{i n j} + Δ P_{m, t}^{i n j}) \geq - F_{l} \forall l, t ({\bar{η}}_{l, t}^{r})

(36)

\sum_{m = 1}^{N_{M}} G S F_{l, m} \cdot (P_{m, t}^{i n j} + Δ P_{m, t}^{i n j}) \geq - F_{l} \forall l, t ({\underset{̲}{η}}_{l, t}^{r})

(37)

- \sum_{q = t - U T_{i}^{g} + 1}^{t} u_{i, q}^{t} \geq - I_{i, t} \forall i, t \in [U T_{i}^{g}, T]

(38)

- \sum_{q = t - D T_{i}^{g} + 1}^{t} v_{i, q}^{t} \geq I_{i, t} - 1 \forall i, t \in [D T_{i}^{g}, T]

(39)

u_{i, t} - v_{i, t} = I_{i, t} - I_{i, t - 1} u_{i, t}, v_{i, t}, I_{i, t} \in {0,1}, \forall i, \forall t

(40)

- u_{i, t} - v_{i, t} \geq - 1 \forall i, \forall t

(41)

where i, j, m, and l are the indices for thermal generators (except RES), wind farms, buses, and transmission lines, respectively; $N_{G}$ , $N_{W}$ , and $N_{M}$ are the numbers of thermal generators, wind farms, and buses, respectively; $P_{i, t}^{g}$ is the energy output of thermal generator i at time t in the basic dispatch process; $Δ P_{i, t}^{g}$ is the reserve in the redispatch process; $C_{i}^{g} (\cdot)$ is the operation cost of thermal generator i; $C_{i}^{S U}$ and $C_{i}^{S D}$ are the startup and shutdown costs of thermal generator i, respectively; $u_{i, t}$ is the startup variable (1 for startup, 0 otherwise); $v_{i, t}$ is the shutdown variable (1 for shutdown, 0 otherwise); $I_{i, t}$ is the unit commitment binary variable (1 for online, 0 otherwise); $P_{j, f, t}^{w}$ and $P_{m, f, t}^{d}$ are day-ahead forecasting power of wind farm j and load at bus m (corresponding to $P_{f, t}^{w}$ in (19) and $P_{f, t}^{d}$ in (20), respectively); $G S F_{l, m}$ is the generation shift factor of bus m to line l; $F_{l}$ is the maximum transmission flow of line l; $P_{i, m i n}^{g}$ and $P_{i, m a x}^{g}$ are the lower and upper output limits of thermal generator i, respectively; the variables in brackets of (22)-(37) are the dual variables of corresponding constraints; $U T_{i}^{g}$ and $D T_{i}^{g}$ are the minimum up time and minimum down time, respectively; $r_{i}^{r u}$ and $r_{i}^{r d}$ are the maximum ramp up/down rates of thermal generator i, respectively; $r_{i}^{s u}$ and $r_{i}^{s d}$ are the startup and shutdown ramp rates, respectively; $Δ P_{m, t}^{d}$ and $Δ P_{j, t}^{w}$ are the difference between forecasting value and the worst-case realization of load ( $P_{m, t}^{d}$ ) and wind power ( $P_{j, t}^{w}$ ), respectively; and $P_{m, t}^{i n j}$ and $Δ P_{m, t}^{i n j}$ are the power injection in basic dispatch process and incremental power injection in redispatch process, respectively. $Δ P_{m, t}^{d}$ , $Δ P_{j, t}^{w}$ , $P_{m, t}^{i n j}$ , and $Δ P_{m, t}^{i n j}$ can be expressed as:

Δ P_{m, t}^{d} = P_{m, t}^{d} - P_{m, f, t}^{d}

(42)

Δ P_{j, t}^{w} = P_{j, t}^{w} - P_{j, f, t}^{w}

(43)

P_{m, t}^{i n j} : = \sum_{i \in G (m)} P_{i, t}^{g} + \sum_{j \in W (m)} P_{j, f, t}^{w} - P_{m, f, t}^{d}

(44)

Δ P_{m, t}^{i n j} : = \sum_{i \in G (m)} Δ P_{i, t}^{g} + \sum_{j \in W (m)} Δ P_{j, t}^{w} - Δ P_{m, t}^{d}

(45)

where $G (m)$ and $W (m)$ are the sets of wind farms and loads at bus m, respectively.

The objective function (21) minimizes the total operation costs over the entire dispatch horizon. Constraints (22)-(24) relate to the basic dispatch, redispatch, and unit on/off status, respectively; (22) and (29) ensure the power balance at each bus; (23), (24), (30), and (31) enforce the generation limits of thermal generators; (25), (26) and (32)-(35) denote the ramping constraints; (27), (28), (36), and (37) enforce the line capacity; (38) and (39) impose the minimum up/down time constraints; (40) describes the relationship between the unit commitment, startup and shutdown variables; and (41) ensures that the startup and shutdown of one thermal generator cannot occur at the same time.

In order to provide reserve capacity for system uncertainty, the generators will limit their energy output $P_{i, t}^{g}$ , which may reduce their revenue from energy provision. Therefore, we define the opportunity cost for providing reserve capacity as the reserve cost $C_{i}^{g} (Δ P_{i, t}^{g})$ , which is the product of per-unit generation cost and reserve capacity of generators. This term guarantees the economy of the reserve deployment and the overall RSCUC model. Furthermore, the price signals derived from the reserve cost-contained RSCED model (described in Section IV) ensure that the reserve revenue covers the reserve cost of thermal generators. The rest in the objective function denote the energy cost $C_{i}^{g} (P_{i, t}^{g})$ , startup cost $C_{i}^{S U} u_{i, t}$ , and shutdown cost $C_{i}^{S D} v_{i, t}$ .

In addition, the proposed RSCUC model includes three parts of ramping constraints. The first two parts, namely (25), (26) and (32), (33), present the requirements for the energy output and reserve capacity of thermal generators, respectively. These constraints are also applied in the RSCUC model presented in [

7]. Different from the previous works, this paper adds constraints (34) and (35) as the third part to ensure that the aggregate of energy and reserve of thermal generators does not violate the ramping limits. In practical situation, even if one thermal generator meets the ramping constraints in terms of both energy output and reserve capacity, the sum of energy and reserve may exceed the ramping limits in some extreme scenarios when the reserve capacity is fully deployed by the independent system operator (ISO). The energy output and reserve capacity of thermal generators are mutually restricted. Constraints (34) and (35) guarantee the safe operation of the thermal generators.

The above RSCUC model can be modeled as a two-stage “min-max-min” RO model as:

\{\begin{array}{l} \underset{x}{m i n} c^{T} x + \underset{\begin{array}{l} P_{m}^{d} \in P^{d, b o x} \\ P_{j}^{w} \in P^{w, I M E U S} \end{array}}{m a x} \underset{y \in Ω (x, P_{m}^{d}, P_{j}^{w})}{m i n} b^{T} y \\ s . t . A x + B y \geq d (γ_{1}) \\ C x + D y = g (γ_{2}) \\ E y + F P_{m}^{d} - F P_{j}^{w} \geq h (γ_{3}) \\ G y + P_{m}^{d} - P_{j}^{w} = p (γ_{4}) \end{array}

(46)

where $c$ and $b$ are the cost coefficient of startup/shutdown cost and operation cost of thermal generators corresponding to objective function (21), respectively; $A$ , $B$ , $C$ , $D$ , $E$ , $F$ , and $G$ are the coefficient matrices of the corresponding constraints; $d$ , $g$ , $h$ , and $p$ are the constant column vectors; and $γ_{1}$ , $γ_{2}$ , $γ_{3}$ , and $γ_{4}$ are the corresponding dual variables of these constraints.

The outer “minimization” model is the first-stage problem with the decision variable $x$ , and the inner “max-min” bi-level model is the second-stage problem with the decision variable $y$ and $(P_{m}^{d}, P_{j}^{w})$ , where $x$ , $y$ , and $(P_{m}^{d}, P_{j}^{w})$ are summarized as (47). $Ω (x, P_{m}^{d}, P_{j}^{w})$ denotes the feasible region of $y$ for a given set of $(x, P_{m}^{d}, P_{j}^{w})$ .

\{\begin{array}{l} x = [I_{i}, u_{i}, v_{i}]^{T} \\ y = [P_{i}^{g}, Δ P_{i}^{g}]^{T} \\ P_{m}^{d} = [P_{m, 1}^{d}, P_{m, 2}^{d}, \dots, P_{m, T}^{d}]^{T} \\ P_{j}^{w} = [P_{j, 1}^{w}, P_{j, 2}^{w}, \dots, P_{j, T}^{w}]^{T} \end{array}

(47)

where $I_{i}$ , $u_{i}$ , and $v_{i}$ are the variable vectors related to the decisions of unit commitment, startup, and shutdown of thermal generator i during [0, T], respectively; $P_{i}^{g}$ and $Δ P_{i}^{g}$ are the vectors of energy and reserve of thermal generator i during [0, T], respectively; and $P_{m}^{d}$ and $P_{j}^{w}$ are the variable vectors related to the worst-case realization of load at bus m and wind power j during [0, T], respectively.

Specifically, after bringing (42)-(45) into (22)-(37), the first line of the constraints in (46) includes (23)-(28), (30)-(35), (38), (39), and (41); the second line includes (22) and (40); the third line includes (36) and (37); and the fourth line includes (29).

B. Solution Method

The column-and-constraint generation (CCG) algorithm [

19] is utilized to solve the RSCUC model. The optimal solution is obtained by decomposing the original problem (46) into the master problem (MP) and subproblem (SP) and solving them alternately. The MP can be expressed as:

\{\begin{array}{l} \underset{x}{m i n} c^{T} x + α \\ s . t . α \geq b^{T} y_{k} \\ A x + B y_{k} \geq d \\ C x + D y_{k} = g \\ E y_{k} + F P_{m, k}^{d *} - F P_{j, k}^{w *} \geq h \\ G y_{k} + P_{m, k}^{d *} - P_{j, k}^{w *} = p \end{array}

(48)

where k is the iteration number; $y_{k}$ is the SP solution at the k^th iteration; and $P_{m, k}^{d *}$ and $P_{j, k}^{w *}$ are the worst-case scenario realizations of wind power and load at the k^th iteration, respectively.

The SP can be expressed as:

\{\begin{array}{l} \underset{\begin{array}{l} P_{m}^{d} \in P^{d, b o x} \\ P_{j}^{w} \in P^{w, I M E U S} \end{array}}{m a x} \underset{y \in Ω (x, P_{m}^{d}, P_{j}^{w})}{m i n} b^{T} y \\ s . t . A x + B y \geq d (γ_{1}) \\ C x + D y = g (γ_{2}) \\ E y + F P_{m}^{d} - F P_{j}^{w} \geq h (γ_{3}) \\ G y + P_{m}^{d} - P_{j}^{w} = p (γ_{4}) \end{array}

(49)

For a given $(x, P_{m}^{d}, P_{j}^{w})$ , the inner “min” problem is a linear optimization problem, which is transformed into a “max” form according to the strong duality theory and combined with the outer “max” problem. The transformed single-level “max” model can be expressed as:

\{\begin{array}{l} \underset{P_{m}^{d}, P_{j}^{w}, γ_{1}, γ_{2}, γ_{3}, γ_{4}}{m a x} ({(d - A x)}^{T} γ_{1} + {(g - C x)}^{T} γ_{2} + h^{T} γ_{3} + p^{T} γ_{4} - \\ \underset{}{} (F P_{m}^{d})^{T} γ_{3} + (F P_{j}^{w})^{T} γ_{3} - (P_{m}^{d})^{T} γ_{4} + (P_{j}^{w})^{T} γ_{4}) \\ s . t . B^{T} γ_{1} + D^{T} γ_{2} + E^{T} γ_{3} + G^{T} γ_{4} \leq b \\ γ_{1}, γ_{3} \geq 0 \end{array}

(50)

There are bilinear terms, i.e., $(F P_{m}^{d})^{T} γ_{3}$ , $(F P_{j}^{w})^{T} γ_{3}$ , $(P_{m}^{d})^{T} γ_{4}$ , and $(P_{j}^{w})^{T} γ_{4}$ , in the objective function. These four nonlinear terms can be linearized by the binary expansion method and the big-M approach [

19]. As a result, the SP can be reformulated as a mixed-integer second-order cone programming (MISOCP) problem, which can be expressed as:

\begin{array}{l} \underset{P_{m}^{d}, P_{j}^{w}, γ_{1}, γ_{2}, γ_{3}, γ_{4}}{m a x} \{{(d - A x)}^{T} γ_{1} + {(g - C x)}^{T} γ_{2} + h^{T} γ_{3} + p^{T} γ_{4} + (P_{j}^{w})^{T} γ_{4, m i n} - \sum_{t = 1}^{T} [P_{m, f, t}^{d} (F_{t}^{P o s} γ_{3, t}^{P o s} \overset{}{\underset{}{+}} F_{t}^{N e g} γ_{3, t}^{N e g}) + Δ P_{m, t}^{d, b o x} (F_{t}^{P o s} ε_{i, t}^{d, P o s} + F_{t}^{N e g} ε_{i, t}^{d, N e g}) - \\ (F_{t}^{P o s} P_{j, t}^{w} γ_{3, t, m i n}^{P o s} + F_{t}^{N e g} P_{j, t}^{w} γ_{3, t, m i n}^{N e g} + F_{t}^{P o s} Δ γ_{3}^{P o s} \sum_{i = 0}^{n} 2^{i} ε_{i, t}^{w, P o s} + F_{t}^{N e g} Δ γ_{3}^{N e g} \sum_{i = 0}^{n} 2^{i} ε_{i, t}^{w, N e g}) + (P_{m, f, t}^{d} γ_{4, t} + Δ P_{m, t}^{d, b o x} ε_{i, t}^{d}) - Δ γ_{4} \sum_{i = 0}^{n} \sum_{t = 1}^{T} 2^{i} ε_{i, t}^{w}]\} \end{array}

(51)

s.t.

B^{T} γ_{1} + D^{T} γ_{2} + E^{T} γ_{3} + G^{T} γ_{4} \leq b γ_{1} \geq 0, γ_{3} \geq 0

(52)

\{\begin{array}{l} 0 \leq ε_{i, t}^{d, P o s} \leq B_{d, t} M \\ γ_{3, t}^{P o s} - (1 - B_{d, t}) M \leq ε_{i, t}^{d, P o s} \leq γ_{3, t}^{P o s} \end{array}

(53)

\{\begin{array}{l} 0 \leq ε_{i, t}^{d, N e g} \leq B_{d, t} M \\ γ_{3, t}^{N e g} - (1 - B_{d, t}) M \leq ε_{i, t}^{d, N e g} \leq γ_{3, t}^{N e g} \end{array}

(54)

\{\begin{array}{l} 0 \leq ε_{i, t}^{w, P o s} \leq χ_{i, t}^{P o s} M \\ P_{j, t}^{w} - (1 - χ_{i, t}^{P o s}) M \leq ε_{i, t}^{w, P o s} \leq P_{j, t}^{w} \end{array}

(55)

\{\begin{array}{l} 0 \leq ε_{i, t}^{w, N e g} \leq χ_{i, t}^{N e g} M \\ P_{j, t}^{w} - (1 - χ_{i, t}^{N e g}) M \leq ε_{i, t}^{w, N e g} \leq P_{j, t}^{w} \end{array}

(56)

\{\begin{array}{l} 0 \leq ε_{i, t}^{d} \leq B_{d, t} M \\ γ_{4, t} - (1 - B_{d, t}) M \leq ε_{i, t}^{d} \leq γ_{4, t} \end{array}

(57)

\{\begin{array}{l} 0 \leq ε_{i, t}^{w} \leq χ_{i, t} \\ P_{j, t}^{w} - (1 - χ_{i, t}) M \leq ε_{i, t}^{w} \leq P_{j, t}^{w} \end{array}

(58)

{‖C_{α, s n}^{- 1 / 2} L_{s, s n} (P_{s n}^{w, I M E U S} - μ_{s, s n})‖}_{2} \leq 1 s_{n} = 1, 2, \dots, N_{s}

(59)

\{\begin{array}{l} P_{t}^{w, I M E U S} \geq B_{w, t} P_{f, t}^{w} \\ \sum_{t = 1}^{T} B_{w, t} \geq T - Γ^{w} \end{array}

(60)

\{\begin{array}{l} P_{t}^{d, b o x} = P_{f, t}^{d} + B_{d, t} Δ P_{t}^{d, b o x} \\ \sum_{t = 1}^{T} B_{d, t} \leq Γ^{d} \end{array}

(61)

The detailed linearizing process and the meaning of each notation are illustrated in Appendix A.

After the above transformation, the CCG algorithm-based solution steps are as follows.

1) Give a set of $(P_{m}^{d}, P_{j}^{w})$ as the initial worst-case scenario, set the lower bound of the objective function $L B \to - \infty$ , and the upper bound $U B \to + \infty$ , $k = 1$ .

2) Solve the MP (28), obtain the solution $(x_{k}^{*}, α_{k}^{*}, y_{k}^{*})$ , and update the lower bound as $L B = c^{T} x_{k}^{*} + α_{k}^{*}$ .

3) Fix $x_{k}^{*}$ and solve the SP (51)-(61), obtain the objective function $f_{k}^{*} (x_{k}^{*})$ and the worst-case scenario $(P_{m, k + 1}^{d *}, P_{j, k + 1}^{w *})$ , and update the upper bound as $U B = m i n {U B, c^{T} x_{k}^{*} + f_{k}^{*} (x_{k}^{*})}$ .

4) Set tolerance $ε$ . If $U B - L B \leq ε$ , the iteration is terminated and return the optimal solution $x_{k}^{*}$ , $y_{k}^{*}$ , and $(P_{m, k + 1}^{d *}, P_{j, k + 1}^{w *})$ ; otherwise, create $y_{k + 1}$ and add the following constraints to the MP:

\{\begin{array}{l} α \geq b^{T} y_{k + 1} \\ A x + B y_{k + 1} \geq d \\ C x + D y_{k + 1} = g \\ E y_{k + 1} + F P_{m, k + 1}^{d *} - F P_{j, k + 1}^{w *} \geq h \\ G y_{k + 1} + P_{m, k + 1}^{d *} - P_{j, k + 1}^{w *} = p \end{array}

(62)

5) Update $k = k + 1$ and go to 2).

IV. Locational Marginal Pricing Mechanism

The ISO should deploy reserve appropriately and form a pricing mechanism in a cost-effective manner to guide power entities to improve the efficiency of system operation while ensuring their own profits. This section first expounds the RSCED model, and then derives a novel locational marginal pricing mechanism.

A. RSCED Model

Fixing the unit commitment status and robust uncertainty scenario realization (i.e., $I_{i, t}$ , $u_{i, t}$ , $v_{i, t}$ , $P_{m, t}^{d}$ , and $P_{j, t}^{w}$ ) obtained in the RSCUC model, the RSCED model can be formulated as:

\{\begin{array}{l} \underset{P^{g}, Δ P^{g}}{m i n} \sum_{t = 1}^{T} \sum_{i = 1}^{N_{G}} C_{i}^{g} (P_{i, t}^{g} + Δ P_{i, t}^{g}) \\ s . t . (22) - (37) \end{array}

(63)

This RSCED model is a linear programming model, which can be easily solved by optimization software.

B. Pricing Mechanism and Market Clearing Mechanism

The Lagrangian function $L (P, Δ P, λ, β, α, η)$ is a by-product of the RSCED model. According to the definition, the LMP at bus $m$ can be derived as:

\begin{array}{l} L M P_{m, t} = \frac{\partial L (P, Δ P, λ, β, α, η)}{\partial P_{m, f, t}^{d}} = λ_{t}^{b} - \sum_{l} G S F_{l, m} \cdot ({\bar{η}}_{l, t}^{b} - {\underset{̲}{η}}_{l, t}^{b}) - \\ \sum_{l} G S F_{l, m} \cdot ({\bar{η}}_{l, t}^{r} - {\underset{̲}{η}}_{l, t}^{r}) \end{array}

(64)

The ULMP is defined as the marginal cost corresponding to the unit increment of forecasting deviation of net load at bus m, which can be derived as:

U L M P_{m, t} = \frac{\partial L (P, Δ P, λ, β, α, η)}{\partial (Δ P_{m,, t}^{d} - \sum_{j \in W (m)} Δ P_{j, t}^{w})} = λ_{t}^{r} - \sum_{l} G S F_{l, m} \cdot ({\bar{η}}_{l, t}^{r} - {\underset{̲}{η}}_{l, t}^{r})

(65)

For the thermal generator i at bus m, LMP and ULMP can also be derived from KKT condition [

8]:

\begin{array}{l} L M P_{m, t} = \frac{\partial C_{i}^{g} (P_{i, t}^{g} + Δ P_{i, t}^{g})}{\partial P_{i, t}^{g}} + {\bar{β}}_{i, t}^{b} - {\underset{̲}{β}}_{i, t}^{b} + {\bar{α}}_{i, t}^{b} - {\underset{̲}{α}}_{i, t}^{b} - \\ {\bar{α}}_{i, t + 1}^{b} + {\underset{̲}{α}}_{i, t + 1}^{b} + {\bar{β}}_{i, t}^{r} - {\underset{̲}{β}}_{i, t}^{r} + {\bar{α}}_{i, t}^{r 2} - {\underset{̲}{α}}_{i, t}^{r 2} - {\bar{α}}_{i, t + 1}^{r 2} + {\underset{̲}{α}}_{i, t + 1}^{r 2} \end{array}

(66)

\begin{array}{l} U L M P_{m, t} = \frac{\partial C_{i}^{g} (P_{i, t}^{g} + Δ P_{i, t}^{g})}{\partial Δ P_{i, t}^{g}} + {\bar{β}}_{i, t}^{r} - {\underset{̲}{β}}_{i, t}^{r} + {\bar{α}}_{i, t}^{r 1} - \\ {\underset{̲}{α}}_{i, t}^{r 1} + {\bar{α}}_{i, t}^{r 2} - {\underset{̲}{α}}_{i, t}^{r 2} - {\bar{α}}_{i, t + 1}^{r 2} + {\underset{̲}{α}}_{i, t + 1}^{r 2} \end{array}

(67)

The detailed derivation processes of (64), (65) and (66), (67) are given in Appendix B and Appendix C, respectively.

It can be seen from (64)-(67) that the marginal prices are affected by various factors such as line capacity and generator capacity [

20], [21], which belong to non-time-coupling factors. On the other hand, the time-coupling factors, mainly the generator ramping constraints, have been proved to cause high marginal prices possibly [22]-[24]. Different from the previous works, this paper proves that ramping constraints may lead to low prices (lower than the generation cost of thermal generators), which is illustrated below.

The dual variables ${\bar{α}}_{i, t}^{b}$ , ${\underset{̲}{α}}_{i, t}^{b}$ , ${\bar{α}}_{i, t + 1}^{b}$ , ${\underset{̲}{α}}_{i, t + 1}^{b}$ , ${\bar{α}}_{i, t}^{r 1}$ , ${\underset{̲}{α}}_{i, t}^{r 1}$ , ${\bar{α}}_{i, t}^{r 2}$ , ${\underset{̲}{α}}_{i, t}^{r 2}$ , ${\bar{α}}_{i, t + 1}^{r 2}$ , and ${\underset{̲}{α}}_{i, t + 1}^{r 2}$ in (66) and (67) correspond to the ramping constraints (25), (26) and (32), (35). When one of these constraints is binding in the day-ahead economic dispatch, the corresponding dual variable is greater than 0. The signs before ${\underset{̲}{α}}_{i, t}^{b}$ , ${\bar{α}}_{i, t + 1}^{b}$ , ${\underset{̲}{α}}_{i, t}^{r 1}$ , ${\underset{̲}{α}}_{i, t}^{r 2}$ , and ${\bar{α}}_{i, t + 1}^{r 2}$ are negative. Therefore, when these dual variables are greater than 0, the prices may decrease. According to the time period when these constraints are binding, it can be discussed in two cases as follows.

1) When the thermal generator i reaches the ramp up limit from time t to time $t + 1$ , the ramping constraint (25) or (34) is binding at time $t + 1$ . At this time, ${\bar{α}}_{i, t + 1}^{b} > 0$ or ${\bar{α}}_{i, t + 1}^{r 2} > 0$ .

2) When the thermal generator i reaches the ramp down limit from time $t - 1$ to time t, the ramping constraint (26), (33), or (35) is binding at time t. At this time, ${\underset{̲}{α}}_{i, t}^{b} > 0$ , ${\underset{̲}{α}}_{i, t}^{r 1} > 0$ , or ${\underset{̲}{α}}_{i, t}^{r 2} > 0$ .

Both cases will decrease LMP and ULMP at time t. The former case happens at the beginning of ramp up stage (ramp up from time t to time $t + 1$ , time t is the beginning time), and the latter case happens at the end of ramp down stage (ramp down from time $t - 1$ to time t, time t is the ending time).

Formulas (29)-(37) are the constraints for the thermal generators to provide reserve capacity against system uncertainties. The related dual variables denote the shadow prices of system resources, such as the energy, generation capacity of thermal generators as well as the line capacity. The expressions of (65) and (67) reveal the influence of uncertainty on ULMP obviously. On the other hand, LMP in (64) is a partial differential of Lagrangian function with respect to the forecasting value of load ( $P_{m, f, t}^{d}$ ), which is a deterministic value. However, the economic dispatch model co-optimizes the basic dispatch process and redispatch process, resulting in the mutual influence of energy scheduling (deal with deterministic forecasting values) and the reserve scheduling (against uncertainties). Accordingly, (64) and (66) that include the dual variables corresponding to redispatch constraints (29)-(37) reflect the impact of system uncertainty on LMP. In a word, both LMP and ULMP internalize the system uncertainty.

Based on the definition, LMP is used for price energy, and ULMP is used for price uncertainty and reserve. The day-ahead market clearing scheme is as follows.

1)　Market Clearing for Thermal Generators

The energy revenue of thermal generator i at bus m is the product of the energy output and LMP, i.e., $L M P_{m, t} P_{i, t}^{g}$ .

The sources that provide reserve to address system uncertainties will get paid for their generation reserve. The revenue for reserve provision of thermal generator i at bus m is the product of its reserve capacity and ULMP, $i . e .$ , $U L M P_{m, t} \cdot Δ P_{i, t}^{g}$ . The term $\partial C_{i}^{g} (\cdot) / \partial Δ P_{i, t}^{g}$ of ULMP in (67) related to the reserve cost in the RSCED model ensures that thermal generators can receive enough revenue commensurate with their reserve cost.

2)　Market Clearing for Wind Farms

The energy revenue of wind farm j at bus m is the product of its power forecasting value and LMP, i.e, $L M P_{m, t} \cdot P_{j, f, t}^{w}$ .

The uncertainty sources need to pay for uncertainties they bring to the system based on the ULMP. Since $Δ P_{j, t}^{w} \leq 0$ in the worst-case scenario, the reserve payment of wind farm j at bus m is $- U L M P_{m, t} \cdot Δ P_{j, t}^{w}$ , which is determined by the maximum forecasting deviation and the ULMP.

3)　Market Clearing for Loads

The energy payment of the load at bus m is the product of its power forecasting value and LMP, i.e., $L M P_{m, t} \cdot P_{m, f, t}^{d}$ . The load at bus m will pay $U L M P_{m, t} \cdot Δ P_{m, t}^{d}$ to the ISO, denoted as the product of the maximum forecasting deviation and the ULMP.

The uncertainties increase the system cost and the payment of uncertainty sources, which are ultimately embodied in the changes of LMP and ULMP. The proposed pricing mechanism can stimulate uncertainty sources to improve the forecasting accuracy and provide effective price signals to incentivize the uncertainty management in the day-ahead market.

V. Case Study

In this section, the effectiveness of the proposed method is verified by the PJM 5-bus system and IEEE 118-bus system. All the simulations have been implemented in MATLAB 2017a with YALMIP interface and CPLEX 12.9.0 in a computational environment with Intel(R) Core(TM) i7-10700F CPU running at 2.9 GHz with 16 GB RAM. The CCG gap is set to be $ε = 0.1 % \times L B$ for both the PJM 5-bus system and IEEE 118-bus system.

A. PJM 5-bus System

In this section, the confidence degree and uncertainty budget are set to be 90% and 24 for all uncertainty sets, respectively.

1)　Comparison of Different Uncertainty Sets

A group of hourly data in 2019 is taken from a 400 MW wind farm in China for research. The temporal correlation coefficient of forecasting error is shown in Fig. 2(a), where the colors indicate the correlation strength and the abscissa/ordinate represent the time interval. The forecasting error shows obvious temporal correlation which is inversely proportional to the time interval distance. This correlation is related to the wind power fluctuation and affects the ramp rate of thermal generators [

25]. The conditional correlation between forecasting error and forecasting value is presented in Fig. 2(b). The forecasting error increases and the error distribution becomes more dispersed as the forecasting value becomes larger. The formulated Gaussian Copula function establishes the relationship between forecasting error and forecasting value, so that the samples of the possible realization of wind power can be updated based on the latest day-ahead forecast power. The temporal correlation and conditional correlation are important features of wind power, which need to be fully considered in uncertainty set modeling.

Fig. 2 Temporal correlation and conditional correlation of wind power forecasting. (a) Temporal correlation coefficient of forecasting error of wind power. (b) Conditional correlation between forecasting error and forecasting value of wind power.

The relationship between the dimension T_D of each ellipsoidal subset and the evaluation indices is illustrated in Fig. 3 when k is set to be 0.3. The integrity index increases monotonously with the increase of the T_D, whereas the efficiency index decreases gradually after $T_{D} \geq 5$ as more invalid region is contained in the IMEUS. The comprehensive index performs best when $T_{D} = 6$ , which is selected as the optimal solution for $T_{D}$ .

Fig. 3 Relationship between dimension of each ellipsoidal subset and evaluation indices.

To verify the advantage of the IMEUS, three uncertainty set modeling methods are compared. Method 1 utilizes the box set that adds error interval on both sides of the forecasting value in each period, which can actually be regarded as the polyhedral set [

26]. Method 2 uses the traditional high-dimensional ellipsoidal set proposed in [10], [12]. Method 3 applies the proposed IMEUS model. The box set is denoted as (20). The traditional high-dimensional ellipsoidal set is a special case of (10), (11) with T_D valued as 24. The proposed IMEUS is formulated as (10), (11) with the optimal T_D. The data from January to November are used as historical data to build the joint distribution of forecasting value and forecasting error, the forecasting value data in the last month is applied to generate the corresponding IMEUS, and the actual wind power data in the last month is used to test the performance. The performance of modeling methods for different uncertainty sets is summarized in Table I. The coverage ratio denotes the total proportion of the actual data in the uncertainty sets, and the average width is the average difference between the upper and lower bounds of the uncertainty sets.

TABLE I Performance of Modeling Methods for Different Uncertainty Sets

Method	Coverage ratio (%)	Average width (MW)
Method 1	90.11	95.84
Method 2	99.85	130.87
Method 3	90.07	80.53

It can be seen that the traditional ellipsoidal set performs better in terms of actual data coverage ratio. However, it contains too much invalid region, resulting in a high degree of conservativeness. On the premise of ensuring the coverage ratio, the IMEUS has the narrowest interval, which can reduce the conservativeness. Although the box set has a similar coverage ratio with IMEUS, its conservativeness is higher than that of the IMEUS. The reason is that the IMEUS takes into account the temporal correlation and conditional correlation to updates the set via the latest day-ahead forecasting data, whereas the box set has a fixed error interval in different days, which is poor in adaptability.

Figure 4 demonstrates the set realization of three methods on a specified day, where the above conclusions can also be drawn. These results indicate that the IMEUS can effectively reduce the conservativeness on the premise of ensuring the modeling accuracy.

Fig. 4 Wind power uncertainty sets of three modeling methods.

The PJM 5-bus system is shown in Fig. 5. The thermal generators originally located at bus A are combined and thermal generator parameters are shown in Table II [

27]. The forecasting total load is distributed to the loads at buses B, C, and D with a ratio of 3:3:4. It is assumed that there is no forecasting error for the load at bus D, and the maximum forecasting error for the loads at bus B and C are 10% and 5% of their forecasting values, respectively. The wind farm mentioned above is located at bus D.

Fig. 5 PJM 5-bus system.

TABLE II Generator Parameters

Generator name	$r^{r u} / r^{s u} / r^{r d} / r^{s d}$ (MW/hour)	C^su ($)	C^sd ($)	$U T^{g}$ (hour)	$D T^{g}$ (hour)
A ( $A l t a + P a r k$ City)	25	360	40	4	3
C (Solitude)	60	500	80	4	4
D (Sundance)	25	300	50	2	2
E (Brighton)	80	550	90	3	3

The proposed RSCUC model is implemented with the three uncertainty set modeling methods for wind power mentioned in Section V-A, respectively. The load is modeled by box set in all cases. The daily average performances of the 31-day (in November) simulations are listed in Table III.

TABLE III RSCUC Results with Different Uncertainty Modeling Methods

Method	The maximum ramp up of wind power (MW)	The maximum ramp down of wind power (MW)	Reserve capacity (MW)	RSCUC cost ($)	Average calculation time (s)	Average iteration number
Method 1	82.17	87.00	1487.31	254193.84	$< 1$
Method 2	91.00	111.08	1915.49	269945.57	2.05	3.16
Method 3	76.00	85.10	1328.62	249342.03	1.91	3.03

Thanks to the improved conservativeness, the extreme wind power fluctuation scenarios with low probability of occurrence are eliminated in the IMEUS. As a result, the ramp-up and ramp-down rate levels of wind power are the lowest in the RSCUC model based on the IMEUS (Method 3). In addition, the IMEUS-based RSCUC can also reduce the reserve demand, and the RSCUC cost is significantly reduced by 1.91% as compared with Method 1 and 7.63% with Method 2. To sum up, the IMEUS reduces the RO conservativeness and the demand for flexible resources, and improves the economy and reliability. On the other hand, the RSCUC calculation of Method 1 consumes less than 1 s and no iterative calculation is required. In Method 1, load and wind power are both modeled by box sets. The upper and lower bounds of the box sets are taken as the worst-case scenario for the load and wind power, respectively. Given the worst-case scenario, the RSCUC is transformed into a single level “min” model, which does not require iterative calculation and takes very little time. The RSCUC based on Method 2 and Method 3 are similar in calculation time and iteration number, which are around 2 s and 3 iterations, respectively.

Furthermore, the performance of uncertainty sets is verified through an out-of-sample evaluation, which utilizes the actual wind power data in the last month as the out-of-sample values to obtain the real benefits of each uncertainty modeling method after the uncertainty is realized. Assume that the wind farm will be punished by the ISO if the actual output is lower than the lower bound of its uncertainty set, as it will lead to additional balancing costs; while the wind farm will be compensated if the actual output is larger than the lower bound of its uncertainty set. The real benefits of wind farm with different penalty coefficients and compensation coefficients are demonstrated in Table IV, where the penalty coefficients are all set to be 2. The two coefficients denote the ratio of the penalty price and compensation price to the day-ahead LMP, respectively.

TABLE IV OUT-OF-SAMPLE RESULTS

Compensation coefficient	Extra profit ( $p e n a l t y + c o m p e n s a t i o n$ ) ($)			Total profit (extra $p r o f i t + p r o f i t$ in day-ahead market) ($)
Compensation coefficient	Method 1	Method 2	Method 3	Method 1	Method 2	Method 3
0.25	1266.70	5961.36	1103.88	27683.28	26002.80	34681.45
0.50	4717.93	12747.90	4487.69	31134.51	32789.34	38065.26
0.75	8169.16	19534.44	7871.51	34585.74	39575.88	41449.08

It can be seen that the extra profit of the wind farm with Method 3 is the lowest among all cases. However, the total profit of Method 3 is the highest in all cases. The reason is that the profit in day-ahead market of Method 3 is significantly higher than the other two methods as the box set and high-dimensional ellipsoidal set are too conservative, which requires higher reserve payments. This out-of-sample test reveals that the proposed IMEUS with lower conservativeness has better economic performance.

To sum up, the proposed IMEUS can reduce the conservativeness of the uncertainty set and improve the accuracy of uncertainty modeling, thus improving the economy of RO.

2)　Comparison of Different RSCUC/RSCED Models

In Section V-A-1), the RSCUC model based on Method 3 (box set for load, IMEUS for wind power) has obtained the worst-case scenario realization of wind power and load, which is fixed in this section to compare the following three different RSCUC/RSCED models.

1) Model 1: the proposed RSCUC/RSCED models in this paper.

2) Model 2: the RSCUC/RSCED models proposed in [

7]. Different from Model 1, the ramping constraints (34), (35) and the reserve cost (

C_{i}^{g} (Δ P_{i, t}^{g})

in (21)) are not considered in this model. The derived ULMP does not include the cost term

\partial C_{i}^{g} (\cdot) / \partial Δ P_{i, t}^{g}

in (66).

3) Model 3: all ramping constraints, i.e., (25), (26) and (32)-(35) are excluded, while other aspects are the same as Model 1.

Model 2 is used to explore the influence of the ramping constraints (34), (35) and the reserve cost $C_{i}^{g} (Δ P_{i, t}^{g})$ on the RSCUC results. Model 3 is designed to explore the influence of all ramping constraints on the RSCUC results and to prove that ramping constraints are one of the reasons leading to low prices.

Figure 6 and Table V indicate the dispatching scheme. Each period includes three histograms, corresponding to the results of three models, respectively. In Fig. 6, “Gen” represents the generator. It can be seen that the load is mainly supported by Gen A and Gen E due to their low generation costs. In Model 1, the residual capacity of Gen A with the cheapest offer is limited after providing energy. Therefore, the reserve is mainly provided by Gen E with large residual capacity and low offer price. In Model 2, the reserve cost is not considered; thus, Gen C and Gen D with higher costs supply a great quantity of reserve capacity.

Fig. 6 Dispatching results of three RSCUC models.

TABLE V ENERGY OUTPUT AND RESERVE RESULTS

Model	Energy output (MWh)				Reserve (MW)
Model	Gen A	Gen C	Gen D	Gen E	Gen A	Gen C	Gen D	Gen E
Model 1	4859.01	554.52	3.00	6345.45	99.94	64.11	12.90	1442.65
Model 2	4921.00	282.15	14.27	6544.57	25.00	1310.89	231.52	52.19
Model 3	5040.00	136.99	0.00	6585.01	0.00	327.93	0.00	1291.68

The above results show that the proportion of energy provided by each generator in Model 1 and Model 2 is close. However, the two models differ greatly in reserve provision. In Model 1, Gen E with lower cost supplies 89.07% reserve. In Model 2, Gen C with higher cost provides 80.94% reserve. Table VI illustrates the day-ahead cost for energy and reserve of the three models. It can be seen that the reserve cost of Model 2 is significantly higher than that of Model 1, as the reserve is provided by thermal generators with higher cost. In addition, the total cost of Model 1 is higher than that of Model 3. The reason is that Model 1 considers the ramping constraints which limit the output of low-cost generators.

TABLE VI RSCUC Cost Results

Model	Energy cost ($)	Reserve cost ($)	Total cost (objective value) ($)
Model 1	216550.14	32791.89	249342.03
Model 2	213742.05	50006.87	263748.92
Model 3	211409.90	35671.73	247081.63

The proposed RSCUC model tightens the ramping constraints. Taking hours 17-18 as an example, the energy output of Gen E increases from 307 MW to 376.87 MW, the reserve capacity increases from 54.64 MW to 64.78 MW, and the total output will increase from 361.64 MW to 441.64 MW if the reserve is fully dispatched in Model 1. Both the energy output and reserve do not violate the ramping constraints, yet their sum just triggers the limit of 80 MW. In Model 2, however, the total output of the energy and reserve increases from 305.00 MW to 393.56 MW, exceeding the ramping limit. A more serious situation occurs in Model 3, where both the energy output (from 313.00 MW to 423.00 MW) and the overall output (from 367.64 MW to 516.56 MW) violate the constraints. Accordingly, the proposed RSCUC/RSCED model ensures the safe operation of generators by considering ramping constraints comprehensively.

3)　Pricing Mechanism and Day-ahead Market Clearing

The LMP and ULMP derived from the three RSCED models mentioned in Section V-A-2) are shown in Fig. 7 and Fig. 8, respectively. In Model 1, the prices are high during heavy load periods in hours 18-22. The prices are different at different buses in each period of hours 20-23 due to the network congestion on line D-E. In addition, there are some differences between LMP and ULMP. For example, LMP is 15 $/MWh higher than ULMP in hour 2 (40 $/MWh vs 25 $/MWh). It can be explained by (66) and (67), and the marginal prices at bus A is taken as an example firstly. $\partial C_{i}^{g} (\cdot) / \partial P_{i, t}^{g}$ and $\partial C_{i}^{g} (\cdot) / \partial Δ P_{i, t}^{g}$ of Gen A (at bus A) are both equal to its generation cost, i.e., 15 $/MWh. Besides, Gen A reaches the output limit in hour 2 and ramp down limit from hour 2 to hour 3, and consequently, the results of ${\underset{̲}{α}}_{A, 3}^{b}$ and ${\bar{β}}_{A, 2}^{r}$ are equal to 15 $/MWh and 10 $/MWh, respectively. Then the marginal prices at bus A are obtained by introducing these nonzero terms into (66) and (67). The prices of other buses can be derived in the same way.

Fig. 7 LMP derived from three RSCED models. (a) Model 1. (b) Model 2. (c) Model 3.

Fig. 8 ULMP derived from three RSCED models. (a) Model 1. (b) Model 2. (c) Model 3.

The three RSCED models produce different prices due to the different considerations of reserve cost and ramping constraints. The LMPs of Model 1 and Model 2 deviate in hours 6-7, 14-15, 20, and 22-23. For instance, LMP of Model 1 is 5 $/MWh higher than that of Model 2 in hour 7 (25 $/MWh vs 20 $/MWh). The cheap Gen A and Gen E in Model 1 provide all the reserve in this period, and yet the expensive Gen C and Gen D in Model 2 also provide reserve due to the lack of reserve cost. The full use of the output capacity of Gen E makes the sum of energy and reserve reach the ramping limit (from 214.41 MW in hour 6 to 294.41 MW in hour 7) in Model 1. As a result, ${\bar{α}}_{i, t}^{r 2}$ of Gen E is equal to 5 $/MWh in Model 1 with constraint (34) binding, whereas there is no such term in Model 2. Thus, LMP at bus E in Model 1 and Model 2 are 25 $/MWh and 20 $/MWh, respectively, according to (66) as $\partial C_{i}^{g} (\cdot) / \partial P_{i, t}^{g} = 20$ $/MWh for Gen E. The LMP of other buses is the same as that of bus E as there is no congestion.

The ULMP in Model 1 almost matches the reserve cost of generators. However, the ULMP in Model 2 are all equal to 0 $/MWh. The reason is that the reserve cost term $\partial C_{i}^{g} (\cdot) / \partial Δ P_{i, t}^{g}$ is considered in the ULMP expression in Model 1, but is ignored in Model 2. The generator revenue/profit from reserve provision in Model 1 and Model 2 is $33677.20/$885.30 and $0/ $$ - 50006.87$ , respectively. The proposed model provides enough revenue for reserve provision of generators, whereas Model 2 cannot cover the reserve cost, resulting in profit loss. Model 3 considers the reserve cost but ignores all the ramping constraints. As a result, the price variation related to ramping constraint is not reflected in Model 3. LMP and ULMP are 20 $/MWh in hours 1-18 and 24, as there is no network congestion and Gen E is the marginal generator. In hours 19-23, Gen C is activated as the marginal generator due to the congestion of line D-E, and the price at each bus goes up consequently. Model 3 guarantees the generator profits, but the lack of ramping constraints increases the operation risk of the system. In conclusion, the proposed pricing mechanism can generate effective price signals to reflect the actual operation safety of the system and ensure the cost recovery of generators simultaneously.

Figures 7 and 8 reveal that low electricity prices (less than the offer prices of generators at corresponding buses) occur in hours 3, 6, 16, and 17 when the generators trigger ramping limits in Model 1. There are no low prices when ramping limits are not considered in Model 3. The following uses the proposed RSCUC/RSCED model (Model 1) to analyze the low electricity prices. Low prices lead to profit loss of generators in that period. Hours 15-19 coupled with ramping constraints are selected for analysis. The profits except Gen D (all zero values) are presented in Table VII (at the left of the solidus “/”). Gen A and Gen E suffer losses due to low prices in hours 16 and 17, but their total profits ($6250 and $4800, respectively) are positive in hours 15-19. Gen A and Gen E have an incentive to reduce output to minimize profit loss. Suppose that Gen A reduces 1 MW energy in hour 16 and Gen E reduces 1 MW reserve in hour 17. To meet the system constraints, other generators follow to adjust energy and reserve in hours 15-19. After the adjustment, the profit change of generators is shown in Table VII (at the right of the solidus “/”).

TABLE VII Generator Profits and Profit Changes in Hours 15-19

Time (hour)	Energy profit ($)			Reserve profit ($)
Time (hour)	Gen A	Gen C	Gen E	Gen A	Gen C	Gen E
15	925/-5	0/0	0/0	125/0	0/0	0/0
16	-4000/25	0/-80	-6810/30	125/0	0/0	0/0
17	2775/-15	0/0	3070/-10	0/0	0/-30	-1093/20
18	3150/-15	0/0	3769/91	0/0	0/0	648/-111
19	3150/-167	0/0	4569/91	0/167	0/0	648/-111

Gen A in hour 16 and Gen E in hour 17 reduce the profit loss by cutting down energy and reserve, respectively. However, due to the output adjustment in the rest hours, the profits of Gen A and Gen C decrease by $10 and $110, respectively, and the profit of Gen E keeps unchanged. The system costs before/after adjustment are $56703.06/$56823.15 in hours 15-19. Thus, the generators should not deviate from the optimal plan of system in terms of both the system cost and generator revenue. Another verification is that low prices may appear at the beginning of ramp up stage (hours 6 and 16) and the end of ramp down stage (hours 3 and 17). From the previous discussion, it can be seen that low prices at time t can be caused by a generator reaching the ramp down limit from time $t - 1$ to t or the ramp up limit from time t to $t + 1$ . In order to gain enough revenue by exporting more energy and reserve in high price periods, e.g., time $t - 1$ and $t + 1$ , the generators constrained by the ramping constraints still provide a certain amount of energy and reserve in low price periods, $e . g .$ , time t, even though they may bear some losses. When these periods are considered as a whole, each generator can obtain a considerable positive profit.

Table VIII shows the day-ahead market clearing results in Model 1. As the energy output provided by all the generators is far more than the reserve capacity, the generator profit mainly comes from energy supply, accounting for 97.62% of the total profit. The loads and wind farms should pay tariff for their uncertainties. The reserve payment of wind farm with strong uncertainty accounts for 33.54% of its energy revenue. Therefore, it is urgent to improve forecasting accuracy to strengthen the regulation ability.

TABLE VIII Day-ahead Market Clearing Results

Clearing result	Load payment ($)	Wind revenue ($)	Generator profit ($)
Energy	335618.28	55625.40	36419.64
Reserve	14512.86	-22047.83	885.30
Total	350131.13	33577.56	37304.94

B. IEEE 118-bus System

The performance of the proposed method is also tested in the IEEE 118-bus system, which includes 118 buses, 54 traditional thermal generators, and 186 lines. The detailed parameters are given in [

7]. This system contains 10 uncertain loads at buses 11, 15, 49, 54, 56, 59, 60, 62, 80, and 90. Five wind farms are located at buses 10, 26, 49, 65, and 80. The load and wind power are modeled by box set and the IMEUS model, respectively.

1)　Sensitivity Analysis of Load Deviation

The uncertainty budgets are all set to be 24, and the confidence degree of the IMEUS is set to be 90%. The forecasting deviation of uncertain load is set to be 0.2, 0.25, and 0.3 of the forecasting values, respectively, and the simulation results are demonstrated in Table IX. It can be seen that the energy cost/revenue, reserve cost/revenue, and the total profits of generators increase with the load deviation level. This is because the width of the box set of load becomes larger along with the increase of forecasting deviation, thus more energy and reserve are required to deal with the uncertainty. Meanwhile, the electricity prices as well as the payment/revenue of load and wind farms increase accordingly. In addition, it can be found that the reserve revenues of generators in Table IX are smaller than the total reserve payment of load and wind farm, which is caused by the difference of ULMP in different bus, also known as the congestion residual. Therefore, the serious system uncertainties make electricity prices go up, which influences the revenue and payment of each participant, not only the uncertainty sources themselves. On the other hand, the CCG algorithm for solving the RSCUC problem can terminate in a finite number of iterations (about 6 iterations in these cases) and takes about 20 minutes.

TABLE IX Day-ahead Calculation Results

Load deviation	System cost ($)	Generators					Load		Wind		LMP/ULMP average ($/MWh)	Average calculation time (s)	Average iteration
Load deviation	System cost ($)	Energy cost ($)	Energy revenue ($)	Reserve cost ($)	Reserve revenue ($)	Total profit ($)	Energy payment ($)	Reserve payment ($)	Energy revenue ($)	Reserve payment ($)	LMP/ULMP average ($/MWh)	Average calculation time (s)	Average iteration
0.20	2.03×10⁶	1.82×10⁶	2.48×10⁶	2.10×10⁵	2.47×10⁵	6.95×10⁵	2.70×10⁶	1.59×10⁵	2.10×10⁵	9.22×10⁴	20.18	1110.92	6.19
0.25	2.07×10⁶	1.83×10⁶	2.51×10⁶	2.40×10⁵	2.93×10⁵	7.30×10⁵	2.73×10⁶	2.01×10⁵	2.14×10⁵	9.36×10⁴	20.45	1270.03	6.45
0.25	2.11×10⁶	1.83×10⁶	2.68×10⁶	2.83×10⁵	3.54×10⁵	9.29×10⁵	2.93×10⁶	2.60×10⁵	2.28×10⁵	9.97×10⁴	21.94	1435.66	6.61

2)　Sensitivity Analysis of Confidence Degree and Uncertainty Budget

When the uncertainty budgets are all set to be 24, the results with different confidence degrees are shown in Table X. It can be seen that with the increase of confidence degree, the reserve demand in the worst-case scenario will increase significantly, which will cause the rise of average electricity prices, and accordingly, lead to the increase in the energy cost and reserve cost of market participants. It should be noted that when the confidence degree changes from 90% to 100%, the unit commitment outcomes and price results have a significant increase compared with that from 70% to 90%. The reason is that when the confidence degree is close to 100%, there will be some extreme scenarios with very low probability in the actual situation, which make the solution too conservative.

TABLE X RSCUC Results with Different Confidence Degrees

Confidence degree (%)	Reserve demand (MW)	System cost ($)	Average LMP/ULMP ($/MWh)	Load		Wind
Confidence degree (%)	Reserve demand (MW)	System cost ($)	Average LMP/ULMP ($/MWh)	Energy payment ($)	Reserve payment ($)	Energy revenue ($)	Reserve payment ($)
70	4595.43	1.86×10⁶	18.43	2.46×10⁶	4.78×10⁴	1.94×10⁵	3.73×10⁴
80	5863.60	1.89×10⁶	18.58	2.47×10⁶	5.92×10⁴	1.96×10⁵	5.03×10⁴
90	8135.70	1.93×10⁶	18.88	2.52×10⁶	8.03×10⁴	1.99×10⁵	7.40×10⁴
100	19843.47	2.17×10⁶	22.33	3.05×10⁶	2.42×10⁵	2.36×10⁵	2.10×10⁵

When the confidence degrees of the box set and the IMEUS are all set to be 90%, the simulation results under different uncertainty budgets are shown in Table XI.

TABLE XI RSCUC Results Under Different Uncertainty Budgets

Uncertainty budget	Reserve demand (MW)	System cost ($)	Average LMP ($/MWh)	Average ULMP ($/MWh)	Load		Wind
Uncertainty budget	Reserve demand (MW)	System cost ($)	Average LMP ($/MWh)	Average ULMP ($/MWh)	Energy payment ($)	Reserve payment ($)	Energy revenue ($)	Reserve payment ($)
0	0.00	1.78×10⁶	18.17	0.00	2.42×10⁶	0	1.92×10⁵	0
6	2186.46	1.82×10⁶	18.44	5.18	2.47×10⁶	2.44×10⁴	1.94×10⁵	2.10×10⁴
12	4230.36	1.86×10⁶	18.61	10.09	2.49×10⁶	4.64×10⁴	1.96×10⁵	3.92×10⁴
18	6254.20	1.90×10⁶	18.67	14.75	2.50×10⁶	6.56×10⁴	1.97×10⁵	5.75×10⁴
24	8135.71	1.93×10⁶	18.88	18.88	2.52×10⁶	8.03×10⁴	1.99×10⁵	7.40×10⁴

It can be seen that when the uncertainty budgets are set to be 0, the forecasting uncertainty of the system is ignored. In this situation, the reserve demand, ULMP, the reserve payment of the uncertain load and wind farms are all equal to 0. The system cost and LMP are also the lowest among the results with different uncertainty budgets. With the increase of uncertainty budgets, the stronger uncertainty raises the demand for reserve capacity, which leads to the increase of unit commitment cost and the prices. Consequently, the energy payment/revenue and the reserve payment of the load and wind farm increase.

In conclusion, the confidence degree and the uncertainty budget have similar impact on unit commitment outcomes and prices. High confidence degree and uncertainty budget indicate high uncertainty, which leads to the increase of unit commitment cost, reserve demand, and electricity prices.

3)　Comparison Between IMEUS Method and MVEE Method

In this section, the proposed method is further compared with the MVEE-based method. Two cases are carried out for comparison. Case 1 uses the MVEE method to construct the ellipsoidal uncertainty set of wind power of the five wind farms in each hour, and 24 ellipsoidal sets are established in the day-ahead electricity market. Case 2 applies the proposed IMEUS method to construct a 24-hour uncertainty set for each wind farm, thus 5 IMEUSs are needed. The confidence degrees are all set to be 90%, and the uncertainty budget of each IMEUS is set to be 24. Based on the uncertainty sets of the two cases, RSCUC and day-ahead market clearing results are shown in Table XII.

TABLE XII RSCUC and Day-ahead Market Clearing Results

Case	Reserve demand (MW)	System cost ($)	Average LMP/ULMP ($/MWh)	Energy revenue (wind) ($)	Reserve payment (wind) ($)	Total revenue (wind) ($)
1	9482.88	1.95×10⁶	18.97	2.00×10⁵	9.98×10⁴	1.00×10⁵
2	8135.70	1.93×10⁶	18.88	1.99×10⁵	7.40×10⁴	1.25×10⁵

It can be seen that the results of Case 2 are superior to those in Case 1, except for the energy revenue of the wind farms due to the low prices in Case 2. The linear correlation coefficients between the wind power of the five wind farms are between -0.3 and 0.3, showing a very weak spatial correlation [

28]. Therefore, although the MVEE method has the ability to consider the spatial correlation, the influence of the spatial correlation on uncertainty set modeling can be ignored in this case study. Accordingly, the main difference between the two cases is that IMEUS takes into account the temporal correlation and conditional correlation of wind power forecasting, whereas MVEE does not consider, which may worsen the extreme events of wind power and increase the reserve demand, RSCUC costs, and prices. Consequently, Case 2 performs better in RSCUC and market clearing results than Case 1.

To further verify the impact of temporal correlation on the extreme ramp events [

12], [15], the relevant calculation is carried out in Case 1 and Case 2, which is shown in Fig. 9. The extreme ramp events denote the wind power fluctuating between the minimum and maximum of the wind power uncertainty set. It can be seen that the proposed IMEUS method performs better than the MVEE method due to the consideration of temporal correlation. In addition, 14 generators reach the climbing limit for 48 times in the IMEUS-based RSCUC, whereas 17 generators reach the climbing limit for 51 times in the MVEE-based RSCUC. To sum up, IMEUS can relieve the extreme ramp events of wind power and extreme climbing phenomenon of generators due to the ability to internalize temporal correlation.

Fig. 9 Extreme ramp-up/ramp-down event. (a) Ramp-up. (b) Ramp-down.

VI. Conclusion

This paper proposes a novel locational marginal pricing mechanism in day-ahead market for managing uncertainties. The improved multi-ellipsoidal uncertainty set is proposed to better characterize the uncertainties of wind power to reduce the conservativeness of RO problems. The robust unit commitment and economic dispatch models are presented to optimize the dispatch scheme for thermal generators and generate price signals for energy, reserve and uncertainty, respectively. The following conclusions can be drawn from the simulations.

1) The proposed multi-ellipsoidal uncertainty set can reduce the conservativeness of uncertainty set and the robust unit commitment/economic dispatch models, and improve the economy and reliability of system operation.

2) The proposed robust unit commitment/economic dispatch models strengthen the ramping constraints to realize the safe operation of generators. The proposed uncertainty LMP ensures the cost recovery of generators, thus can motivate the participants to provide reserve services. The system operation cost is reduced by the proposed methods.

3) The generator ramping constraint is one of the reasons for the low marginal prices. Low prices may cause the profit loss of generators in that period, but the generator profits of the whole day can be guaranteed.

The uncertainty-contained electricity market should have two sequential market mechanisms, i.e., day-ahead market and real-time market, which are more tightly linked. In the future, our work will focus on designing a real-time market mechanism to complement the proposed day-ahead market mechanism based on the results already obtained in day-ahead market clearing.

Appendix

Appendix A

Appendix A will present the linearizing process of the subproblem. There are bilinear terms, i.e., $(F P_{m}^{d})^{T} γ_{3}$ , $(F P_{j}^{w})^{T} γ_{3}$ , $(P_{m}^{d})^{T} γ_{4}$ , and $(P_{j}^{w})^{T} γ_{4}$ , in the objective function of (30). These four nonlinear terms can be linearized by the binary expansion method and the big-M approach which is illustrated below.

1) Suppose $F = [(F^{P o s})^{T}, (F^{N e g})^{T}]^{T}$ . $F^{P o s}$ and $F^{N e g}$ are diagonal matrices of T order corresponding to constraints (36) and (37), respectively, which can be expressed as:

\{\begin{array}{l} F^{P o s} = d i a g (F_{1}^{P o s}, F_{2}^{P o s}, \dots, F_{T}^{P o s}) \\ F^{N e g} = d i a g (F_{1}^{N e g}, F_{2}^{N e g}, \dots, F_{T}^{N e g}) \end{array}

(A1)

$γ_{3}$ includes $2 T$ elements, which can be divided into:

\{\begin{array}{l} γ_{3} = [γ_{3,1}, γ_{3,2}, \dots, γ_{3,2 T}]^{T} = [(γ_{3}^{P o s})^{T}, (γ_{3}^{N e g})^{T}]^{T} \\ γ_{3}^{P o s} = [γ_{3,1}, γ_{3,1}, \dots, γ_{3, T}]^{T} = [γ_{3,1}^{P o s}, γ_{3,2}^{P o s}, \dots, γ_{3, T}^{P o s}]^{T} \\ γ_{3}^{N e g} = [γ_{3, T + 1}, γ_{3, T + 2}, \dots, γ_{3,2 T}]^{T} = [γ_{3,1}^{N e g}, γ_{3,2}^{N e g}, \dots, γ_{3, T}^{N e g}]^{T} \end{array}

(A2)

Thus, $(F P_{m}^{d})^{T} γ_{3}$ can be transformed into:

\begin{array}{l} (F P_{m}^{d})^{T} γ_{3} = (F^{P o s} P_{m}^{d})^{T} γ_{3}^{P o s} + (F^{N e g} P_{m}^{d})^{T} γ_{3}^{N e g} = \\ \sum_{t = 1}^{T} [P_{m, f, t}^{d} (F_{t}^{P o s} γ_{3, t}^{P o s} + F_{t}^{N e g} γ_{3, t}^{N e g}) + \\ Δ P_{m, t}^{d, b o x} (F_{t}^{P o s} B_{d, t} γ_{3, t}^{P o s} + F_{t}^{N e g} B_{d, t} γ_{3, t}^{N e g})] \end{array}

(A3)

By introducing the dummy variable $ε_{i, t}^{d, P o s} = B_{d, t} γ_{3, t}^{P o s}$ and $ε_{i, t}^{d, N e g} = B_{d, t} γ_{3, t}^{N e g}$ , the bilinear term can be linearized by the use of big-M approach, i.e.,

\{\begin{array}{l} (F P_{m}^{d})^{T} γ_{3} = \sum_{t = 1}^{T} [P_{m, f, t}^{d} (F_{t}^{P o s} γ_{3, t}^{P o s} + F_{t}^{N e g} γ_{3, t}^{N e g}) + \\ Δ P_{m, t}^{d, b o x} (F_{t}^{P o s} ε_{i, t}^{d, P o s} + F_{t}^{N e g} ε_{i, t}^{d, N e g})] \\ 0 \leq ε_{i, t}^{d, P o s} \leq B_{d, t} M \\ γ_{3, t}^{P o s} - (1 - B_{d, t}) M \leq ε_{i, t}^{d, P o s} \leq γ_{3, t}^{P o s} \\ 0 \leq ε_{i, t}^{d, N e g} \leq B_{d, t} M \\ γ_{3, t}^{N e g} - (1 - B_{d, t}) M \leq ε_{i, t}^{d, N e g} \leq γ_{3, t}^{N e g} \end{array}

(A4)

where $M$ is a large real number.

2) $(F P_{j}^{w})^{T} γ_{3}$ can be transformed into:

\begin{array}{l} (F P_{j}^{w})^{T} γ_{3} = (F^{P o s} P_{j}^{w})^{T} γ_{3}^{P o s} + (F^{N e g} P_{j}^{w})^{T} γ_{3}^{N e g} = \\ \sum_{t = 1}^{T} (F_{t}^{P o s} P_{j, t}^{w} γ_{3, t}^{P o s} + F_{t}^{N e g} P_{j, t}^{w} γ_{3, t}^{N e g}) \end{array}

(A5)

Through the binary expansion method, $γ_{3, t}^{P o s}$ and $γ_{3, t}^{N e g}$ can be expanded as:

γ_{3, t}^{P o s} = γ_{3, t, m i n}^{P o s} + Δ γ_{3}^{P o s} \sum_{i = 0}^{n} 2^{i} χ_{i, t}^{P o s}

(A6)

γ_{3, t}^{N e g} = γ_{3, t, m i n}^{N e g} + Δ γ_{3}^{N e g} \sum_{i = 0}^{n} 2^{i} χ_{i, t}^{N e g}

(A7)

where $γ_{3, t, m i n}^{P o s}$ and $γ_{3, t, m i n}^{N e g}$ are the lower bounds of $γ_{3, t}^{P o s}$ and $γ_{3, t}^{N e g}$ , respectively; $Δ γ_{3}^{P o s}$ and $Δ γ_{3}^{N e g}$ are the interval of the binary expansion; and $χ_{i, t}^{P o s}$ and $χ_{i, t}^{N e g}$ are 0-1 binary variables.

Thus, $(F P_{j}^{w})^{T} γ_{3}$ can be transformed as:

\begin{array}{l} (F P_{j}^{w})^{T} γ_{3} = \sum_{t = 1}^{T} (F_{t}^{P o s} P_{j, t}^{w} γ_{3, t, m i n}^{P o s} + F_{t}^{N e g} P_{j, t}^{w} γ_{3, t, m i n}^{N e g} + \begin{matrix}  \end{matrix} \\ F_{t}^{P o s} Δ γ_{3}^{P o s} \sum_{i = 0}^{n} 2^{i} χ_{i, t}^{P o s} P_{j, t}^{w} + F_{t}^{N e g} Δ γ_{3}^{N e g} \sum_{i = 0}^{n} 2^{i} χ_{i, t}^{N e g} P_{j, t}^{w}) \end{array}

(A8)

By introducing $ε_{i, t}^{w, P o s} = χ_{i, t}^{P o s} P_{j, t}^{w}$ and $ε_{i, t}^{w, N e g} = χ_{i, t}^{N e g} P_{j, t}^{w}$ , the bilinear term can be linearized, i.e.,

\{\begin{array}{l} (F P_{j}^{w})^{T} γ_{3} = \sum_{t = 1}^{T} (F_{t}^{P o s} P_{j, t}^{w} γ_{3, t, m i n}^{P o s} + F_{t}^{N e g} P_{j, t}^{w} γ_{3, t, m i n}^{N e g} + \begin{matrix}  \end{matrix} \\ F_{t}^{P o s} Δ γ_{3}^{P o s} \sum_{i = 0}^{n} 2^{i} ε_{i, t}^{w, P o s} + F_{t}^{N e g} Δ γ_{3}^{N e g} \sum_{i = 0}^{n} 2^{i} ε_{i, t}^{w, N e g}) \\ 0 \leq ε_{i, t}^{w, P o s} \leq χ_{i, t}^{P o s} M \\ P_{j, t}^{w} - (1 - χ_{i, t}^{P o s}) M \leq ε_{i, t}^{w, P o s} \leq P_{j, t}^{w} \\ 0 \leq ε_{i, t}^{w, N e g} \leq χ_{i, t}^{N e g} M \\ P_{j, t}^{w} - (1 - χ_{i, t}^{N e g}) M \leq ε_{i, t}^{w, N e g} \leq P_{j, t}^{w} \end{array}

(A9)

3) $γ_{4}$ has $T \times$ 1 elements $γ_{4, t}$ , $t = 1,2, \dots, T$ , which corresponds to constraint (29). By introducing $ε_{i, t}^{d} = B_{d, t} γ_{4, t}$ , $(P_{m}^{d})^{T} γ_{4}$ can be linearized as:

\{\begin{array}{l} (P_{m}^{d})^{T} γ_{4} = \sum_{t = 1}^{T} (P_{m, f, t}^{d} γ_{4, t} + Δ P_{m, t}^{d, b o x} ε_{i, t}^{d}) \\ 0 \leq ε_{i, t}^{d} \leq B_{d, t} M \\ γ_{4, t} - (1 - B_{d, t}) M \leq ε_{i, t}^{d} \leq γ_{4, t} \end{array}

(A10)

4) $γ_{4}$ can be expanded as:

γ_{4, t} = γ_{4, t, m i n} + Δ γ_{4} \sum_{i = 0}^{n} 2^{i} χ_{i, t}

(A11)

where $γ_{4, t, m i n}$ is the lower bound of $γ_{4, t}$ ; $Δ γ_{4}$ is the interval of the binary expansion; and $χ_{i, t}$ is a 0-1 binary variable.

Similar to the linearization of $(F P_{j}^{w})^{T} γ_{3}$ , $(P_{j}^{w})^{T} γ_{4}$ can be linearized as:

\{\begin{array}{l} (P_{j}^{w})^{T} γ_{4} = (P_{j}^{w})^{T} γ_{4, m i n} + Δ γ_{4} \sum_{i = 0}^{n} \sum_{t = 1}^{T} 2^{i} ε_{i, t}^{w} \\ 0 \leq ε_{i, t}^{w} \leq χ_{i, t} \\ P_{j, t}^{w} - (1 - χ_{i, t}) M \leq ε_{i, t}^{w} \leq P_{j, t}^{w} \end{array}

(A12)

where $γ_{4, m i n} = [γ_{4,1, m i n}, γ_{4,2, m i n}, \dots, γ_{4, T, m i n}]^{T}$ is a constant vector; and the dummy variable $ε_{i, t}^{w} = χ_{i, t} P_{j, t}^{w}$ .

As a result, the SP can be reformulated as a mixed-integer second-order cone programming (MISOCP) problem, which can be expressed in (51)-(61).

Appendix B

Appendix B will present the derivation process of prices based on Lagrangian function. The Lagrangian function is expressed as:

\begin{array}{l} L (P, Δ P, λ, β, α, η) = \sum_{t = 1}^{T} \sum_{i = 1}^{N_{G}} C_{i}^{g} (P_{i, t}^{g} + Δ P_{i, t}^{g}) + \sum_{t = 1}^{T} λ_{t}^{b} (\sum_{m = 1}^{N_{M}} P_{m, f, t}^{d} - \sum_{i = 1}^{N_{G}} P_{i, t}^{g} - \sum_{j = 1}^{N_{W}} P_{j, f, t}^{w}) + \sum_{t = 1}^{T} \sum_{i = 1}^{N_{G}} [{\bar{β}}_{i, t}^{b} (P_{i, t}^{g} - I_{i, t} P_{i, m a x}^{g}) + {\underset{̲}{β}}_{i, t}^{b} (I_{i, t} P_{i, m i n}^{g} - P_{i, t}^{g})] + \\ \sum_{t = 1}^{T} \sum_{i = 1}^{N_{G}} \{{\bar{α}}_{i, t}^{b} [P_{i, t}^{g} - P_{i, t - 1}^{g} - r_{i}^{r u} (1 - u_{i, t}) - r_{i}^{s u} u_{i, t}] + {\underset{̲}{α}}_{i, t}^{b} [P_{i, t - 1}^{g} - P_{i, t}^{g} - r_{i}^{r d} (1 - v_{i, t}) - r_{i}^{s d} v_{i, t}]\} + \\ \sum_{t = 1}^{T} \sum_{l = 1}^{N_{L}} \{{\bar{η}}_{l, t}^{b} [G S F_{l, m} (\sum_{i \in G (m)} P_{i, t}^{g} + \sum_{j \in W (m)} P_{j, f, t}^{w} - P_{m, f, t}^{d}) - F_{l}] - {\underset{̲}{η}}_{l, t}^{b} [G S F_{l, m} (\sum_{i \in G (m)} P_{i, t}^{g} + \sum_{j \in W (m)} P_{j, f, t}^{w} - P_{m, f, t}^{d}) + F_{l}]\} + \\ \sum_{t = 1}^{T} λ_{t}^{r} (\sum_{m = 1}^{N_{M}} Δ P_{m, t}^{d} - \sum_{j = 1}^{N_{W}} Δ P_{j, t}^{w} - \sum_{i = 1}^{N_{G}} Δ P_{i, t}^{g}) + \sum_{t = 1}^{T} \sum_{i = 1}^{N_{G}} [{\bar{β}}_{i, t}^{r} (P_{i, t}^{g} + Δ P_{i, t}^{g} - I_{i, t} P_{i, m a x}^{g}) + {\underset{̲}{β}}_{i, t}^{r} (I_{i, t} P_{i, m i n}^{g} - P_{i, t}^{g} - Δ P_{i, t}^{g})] + \\ \sum_{t = 1}^{T} \sum_{i = 1}^{N_{G}} \{{\bar{α}}_{i, t}^{r 1} [Δ P_{i, t}^{g} - r_{i}^{r u} (1 - u_{i, t}) - r_{i}^{s u} u_{i, t}] + {\underset{̲}{α}}_{i, t}^{r 1} [- Δ P_{i, t}^{g} - r_{i}^{r d} (1 - v_{i, t}) - r_{i}^{s d} v_{i, t}]\} + \\ \sum_{t = 1}^{T} \sum_{i = 1}^{N_{G}} \{{\bar{α}}_{i, t}^{r 2} [P_{i, t}^{g} + Δ P_{i, t}^{g} - P_{i, t - 1}^{g} - Δ P_{i, t - 1}^{g} - r_{i}^{r u} (1 - u_{i, t}) - r_{i}^{s u} u_{i, t}] + {\underset{̲}{α}}_{i, t}^{r 2} [P_{i, t - 1}^{g} + Δ P_{i, t - 1}^{g} - P_{i, t}^{g} - Δ P_{i, t}^{g} - r_{i}^{r d} (1 - v_{i, t}) - r_{i}^{s d} v_{i, t}]\} + \\ \sum_{t = 1}^{T} \sum_{l = 1}^{N_{L}} \{{\bar{η}}_{l, t}^{r} [G S F_{l, m} (\sum_{i \in G (m)} P_{i, t}^{g} + \sum_{j \in W (m)} P_{j, f, t}^{w} - P_{m, f, t}^{d} \sum_{i \in G (m)} Δ P_{i, t}^{g} + \sum_{j \in W (m)} Δ P_{j, t}^{w} - Δ P_{m, t}^{d}) - F_{l}] - \\ {\underset{̲}{η}}_{l, t}^{r} [G S F_{l, m} (\sum_{i \in G (m)} P_{i, t}^{g} + \sum_{j \in W (m)} P_{j, f, t}^{w} - P_{m, f, t}^{d} + \sum_{i \in G (m)} Δ P_{i, t}^{g} + \sum_{j \in W (m)} Δ P_{j, t}^{w} - Δ P_{m, t}^{d}) + F_{l}]\} \end{array}

(B1)

LMP is a partial differential of Lagrangian function with respect to the forecasting value of load, which can be derived as:

\begin{array}{l} L M P_{m, t} = \frac{\partial L (P, Δ P, λ, β, α, η)}{\partial P_{m, f, t}^{d}} = λ_{t}^{b} - \\ \sum_{l} G S F_{l, m} \cdot ({\bar{η}}_{l, t}^{b} - {\underset{̲}{η}}_{l, t}^{b}) - \sum_{l} G S F_{l, m} \cdot ({\bar{η}}_{l, t}^{r} - {\underset{̲}{η}}_{l, t}^{r}) \end{array}

(B2)

The ULMP is defined as the marginal cost corresponding to the unit increment of forecasting deviation of net load at bus m, which can be derived as:

U L M P_{m, t} = \frac{\partial L (P, Δ P, λ, β, α, η)}{\partial (Δ P_{m,, t}^{d} - \sum_{j \in W (m)} Δ P_{j, t}^{w})} = λ_{t}^{r} - \sum_{l} G S F_{l, m} ({\bar{η}}_{l, t}^{r} - {\underset{̲}{η}}_{l, t}^{r})

(B3)

APPENDIX C

Appendix C will present the derivation process of prices based on Karush-Kuhn-Tucker (KKT) condition. According to the KKT condition $\partial L (P, Δ P, λ, β, α, η) / \partial P_{i, t}^{g} = 0$ [

8], it can be obtained as:

\begin{array}{l} \frac{\partial C_{i}^{g} (P_{i, t}^{g} + Δ P_{i, t}^{g})}{\partial P_{i, t}^{g}} + {\bar{β}}_{i, t}^{b} - {\underset{̲}{β}}_{i, t}^{b} + {\bar{α}}_{i, t}^{b} - {\underset{̲}{α}}_{i, t}^{b} - {\bar{α}}_{i, t + 1}^{b} + \\ {\underset{̲}{α}}_{i, t + 1}^{b} + {\bar{β}}_{i, t}^{r} - {\underset{̲}{β}}_{i, t}^{r} + {\bar{α}}_{i, t}^{r 2} - {\underset{̲}{α}}_{i, t}^{r 2} - {\bar{α}}_{i, t + 1}^{r 2} + {\underset{̲}{α}}_{i, t + 1}^{r 2} = \\ λ_{t}^{b} - \sum_{l} G S F_{l, m} \cdot ({\bar{η}}_{l, t}^{b} - {\underset{̲}{η}}_{l, t}^{b}) - \sum_{l} G S F_{l, m} \cdot ({\bar{η}}_{l, t}^{r} - {\underset{̲}{η}}_{l, t}^{r}) \end{array}

(C1)

The right side of the equal sign of this equation is exactly the LMP at bus m and time $t$ , i.e., $L M P_{m, t}$ . Therefore, for the generator i at bus m, LMP can also be derived as (C2).

\begin{array}{l} L M P_{m, t} = \frac{\partial C_{i}^{g} (P_{i, t}^{g} + Δ P_{i, t}^{g})}{\partial P_{i, t}^{g}} + {\bar{β}}_{i, t}^{b} - {\underset{̲}{β}}_{i, t}^{b} + {\bar{α}}_{i, t}^{b} - {\underset{̲}{α}}_{i, t}^{b} - \\ {\bar{α}}_{i, t + 1}^{b} + {\underset{̲}{α}}_{i, t + 1}^{b} + {\bar{β}}_{i, t}^{r} - {\underset{̲}{β}}_{i, t}^{r} + {\bar{α}}_{i, t}^{r 2} - {\underset{̲}{α}}_{i, t}^{r 2} - {\bar{α}}_{i, t + 1}^{r 2} + {\underset{̲}{α}}_{i, t + 1}^{r 2} \end{array}

(C2)

Similarly, (C3) can be obtained according to the KKT condition $\partial L (P, Δ P, λ, β, α, η) / \partial Δ P_{i, t}^{g} = 0$ .

\begin{array}{l} \frac{\partial C_{i}^{g} (P_{i, t}^{g} + Δ P_{i, t}^{g})}{\partial Δ P_{i, t}^{g}} + {\bar{β}}_{i, t}^{r} - {\underset{̲}{β}}_{i, t}^{r} + {\bar{α}}_{i, t}^{r 1} - {\underset{̲}{α}}_{i, t}^{r 1} + {\bar{α}}_{i, t}^{r 2} - \\ {\underset{̲}{α}}_{i, t}^{r 2} - {\bar{α}}_{i, t + 1}^{r 2} + {\underset{̲}{α}}_{i, t + 1}^{r 2} = λ_{t}^{r} - \sum_{l} G S F_{l, m} \cdot ({\bar{η}}_{l, t}^{r} - {\underset{̲}{η}}_{l, t}^{r}) \end{array}

(C3)

The right side of the equal sign of (C3) is exactly the ULMP at bus m and time $t$ , i.e., $U L M P_{m, t}$ . Therefore, for the generator i at bus m, ULMP can also be derived as (C4).

\begin{array}{l} U L M P_{m, t} = \frac{\partial C_{i}^{g} (P_{i, t}^{g} + Δ P_{i, t}^{g})}{\partial Δ P_{i, t}^{g}} + {\bar{β}}_{i, t}^{r} - {\underset{̲}{β}}_{i, t}^{r} + \\ {\bar{α}}_{i, t}^{r 1} - {\underset{̲}{α}}_{i, t}^{r 1} + {\bar{α}}_{i, t}^{r 2} - {\underset{̲}{α}}_{i, t}^{r 2} - {\bar{α}}_{i, t + 1}^{r 2} + {\underset{̲}{α}}_{i, t + 1}^{r 2} \end{array}

(C4)

REFERENCES

Z. Ullah, G. Mokryani, F. Campean et al., “Comprehensive review of VPPs planning, operation and scheduling considering the uncertainties related to renewable energy sources,” IET Energy Systems Integration, vol. 1, no. 3, pp. 147-157, Sept. 2019. [Baidu Scholar]

A. Akrami, M. Doostizadeh, and F. Aminifar, “Power system flexibility: an overview of emergence to evolution,” Journal of Modern Power Systems and Clean Energy, vol. 7, no. 5, pp. 987-1007, Sept. 2019. [Baidu Scholar]

X. Fang and M. Cui, “Analytical Model of Day-ahead and Real-time Price Correlation in Strategic Wind Power Offering,” Journal of Modern Power Systems and Clean Energy, vol. 8, no. 5, pp. 1024-1027, Sept. 2020. [Baidu Scholar]

J. M. Morales, A. J. Conejo, K. Liu et al., “Pricing electricity in pools with wind producers,” IEEE Transactions on Power Systems, vol. 27, no. 3, pp. 1366-1376, Aug. 2012. [Baidu Scholar]

Y. Dvorkin, “A chance-constrained stochastic electricity market,” IEEE Transactions on Power Systems, vol. 35, no. 4, pp. 2993-3003, Jul. 2020. [Baidu Scholar]

R. Mieth, J. Kim, and Y. Dvorkin, “Risk- and variance-aware electricity pricing,” Electric Power Systems Research, vol. 189, pp. 1-8, Dec. 2020. [Baidu Scholar]

H. Ye, Y. Ge, M. Shahidehpour et al., “Uncertainty marginal price, transmission reserve, and day-ahead market clearing with robust unit commitment,” IEEE Transactions on Power Systems, vol. 32, no. 3, pp. 1782-1795, May 2017. [Baidu Scholar]

J. Kazempour, P. Pinson, and B. F. Hobbs, “A stochastic market design with revenue adequacy and cost recovery by scenario: benefits and costs,” IEEE Transactions on Power Systems, vol. 33, no. 4, pp. 3531-3545, Jul. 2018. [Baidu Scholar]

X. Fang, B. Hodge, E. Du et al., “Introducing uncertainty components in locational marginal prices for pricing wind power and load uncertainties,” IEEE Transactions on Power Systems, vol. 34, no. 3, pp. 2013-2024, May 2019. [Baidu Scholar]

Y. Guan and J. Wang, “Uncertainty sets for robust unit commitment,” IEEE Transactions on Power Systems, vol. 29, no. 3, pp. 1439-1440, May 2014. [Baidu Scholar]

A. R. Malekpour and A. Pahwa, “Stochastic networked microgrid energy management with correlated wind generators,” IEEE Transactions on Power Systems, vol. 32, no. 5, pp. 3681-3693, Sept. 2017. [Baidu Scholar]

Y. Tian, W. Wu, K. Wang et al., “Robust transmission constrained unit commitment under wind power uncertainty with adjustable conservatism,” IET Generation, Transmission & Distribution, vol. 14, no. 5, pp. 824-832, Mar. 2020. [Baidu Scholar]

X. Xu, Z. Yan, M. Shahidehpour et al., “Data-driven risk-averse two-stage optimal stochastic scheduling of energy and reserve with correlated wind power,” IEEE Transactions on Sustainable Energy, vol. 11, no. 1, pp. 436-447, Jan. 2020. [Baidu Scholar]

T. Ding, J. Lv, R. Bo et al., “Lift-and-project MVEE based convex hull for robust SCED with wind power integration using historical data-driven modeling approach,” Renewable Energy, vol. 92, pp. 415-427, Jul. 2016. [Baidu Scholar]

J. Xu, B. Wang, Y. Sun et al., “A day-ahead economic dispatch method considering extreme scenarios based on wind power uncertainty,” CSEE Journal of Power and Energy Systems, vol. 5, no. 2, pp. 224-233, Jun. 2019. [Baidu Scholar]

M. T. Bina and D. Ahmadi, “Stochastic modeling for the next day domestic demand response applications,” IEEE Transactions on Power Systems, vol. 30, no. 6, pp. 2880-2893, Nov 2015. [Baidu Scholar]

H. Hult and F. Lindskog, “Multivariate extremes, aggregation and dependence in elliptical distributions,” Advances in Applied Probability, vol. 34, no. 3, pp. 587-608, Sept. 2002. [Baidu Scholar]

A. Jamalizadeh and N. Balakrishnan, “Conditional distributions of multivariate normal mean-variance mixtures,” Statistics & Probability Letters, vol. 145, pp. 312-316, Feb. 2019. [Baidu Scholar]

T. Ding, S. Liu, W. Yuan et al., “A two-stage robust reactive power optimization considering uncertain wind power integration in active distribution networks,” IEEE Transactions on Sustainable Energy, vol. 7, no. 1, pp. 301-311, Nov. 2015. [Baidu Scholar]

E. Ela and D. Edelson, “Participation of wind power in LMP-based energy markets,” IEEE Transactions on Sustainable Energy, vol. 3, no. 4, pp. 777-783, Oct. 2012. [Baidu Scholar]

M. Nicolosi, “Wind power integration and power system flexibility–an empirical analysis of extreme events in Germany under the new negative price regime,” Energy Policy, vol. 38, no. 11, pp. 7257-7268, Nov. 2010. [Baidu Scholar]

R. Khatami and M. Parvania, “Continuous-time locational marginal price of electricity,” IEEE Access, vol. 7, pp. 129480-129493, Sept. 2019. [Baidu Scholar]

K. H. Abdul-Rahman, H. Alarian, M. Rothleder et al., “Enhanced system reliability using flexible ramp constraint in CAISO market,” in Proceedings of 2012 IEEE PES General Meeting, San Diego, USA, Jul. 2012, pp. 1-6. [Baidu Scholar]

M. Milligan, E. Ela, D. Lew et al., “Operational analysis and methods for wind integration studies,” IEEE Transactions on Sustainable Energy, vol. 3, no. 4, pp. 612-619, Oct. 2012. [Baidu Scholar]

B. Zhou, G. Geng, and Q. Jiang, “Hierarchical unit commitment with uncertain wind power generation,” IEEE Transactions on Power Systems, vol. 31, no. 1, pp. 94-104, Jan. 2016. [Baidu Scholar]

B. Vatani, B. Chowdhury, S. Dehghan et al., “A critical review of robust self-scheduling for generation companies under electricity price uncertainty,” International Journal of Electrical Power & Energy Systems, vol. 97, pp. 428-439, Apr. 2018. [Baidu Scholar]

F. Li and R. Bo, “Small test systems for power system economic studies,” in Proceedings of 2010 IEEE PES General Meeting, Minneapolis, USA, Jul. 2010, pp. 1-4. [Baidu Scholar]

S. Xu, H. Sun, and J. Wu, “Cross-correlation analysis in mixed traffic flow time series,” International Journal of Modern Physics B, vol. 25, no. 13, pp. 1823-1832, May 2011. [Baidu Scholar]

Address:No.19 Chengxin Avenue, Jiangning District, Nanjing 211106, China

E-mail: mpce@alljournals.cn

Tel:86-25-81093060

Fax:86-25-81093040

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher