Multi-stage robust clearing model considering renewable energy output uncertainty and unit effective reserve calculation in electricity market

With the advancement of China’s carbon policy, the proportion of renewable energy, mainly wind and solar, has increased, which brings greater challenges to the reservation of power systems. Due to the complexity of China’s power systems, adopting the method of the zonal reserve to ensure stable system operation is very difficult. Among the existing clearing results for calculating reservation, the unit reservation will be restricted by the network security constraints, leading to system operation risks. To obtain the unit reservations efficiently without breaking the security constraints, this paper proposes an effective reserve calculation method for engineering implementation. A box robust optimization is further utilized for Security Constrained Unit Commitment in the method to ensure the consumption of renewable energy, and the stability and clearing efficiency of the power system. Furthermore, a data-driven robust random optimization is employed for Security Constrained Economic Dispatch optimizing the economy of power systems. This multi-stage robust optimization model, which has good extensibility, is in line with China’s clearing process. The simulation analysis based on the actual operation data of a province in China and the IEEE 300-bus system verifies the correctness and feasibility of the proposed model and theory.

1 Introduction

With the promotion of China’s policy and the task of carbon peaking and carbon neutralization and the continuous promotion of power system reform (Lu et al., 2022), the development of clean and renewable energy is imperative (Farhani, 2015), and its access proportion in the power system will gradually rise. Although renewable energy such as wind and solar energy is environmentally friendly, the volatility and randomness of output lead to the phenomenon of abandoning wind power and photovoltaic power to ensure the safe operation of the system when the electricity market is cleared. To cope with the challenges brought by the large-scale access of renewable energy to the safe and stable operation of China’s power system, more system operation reserves need to be reserved in the power system dispatching optimization. Therefore, it is of great significance to put forward a set of reasonable reserve calculation methods based on absorbing renewable energy as much as possible.

The positive and reserve capacity of the unit is usually determined by the unit capacity and the current output of the unit. However, in actual operation, if the positive and reserve capacity of the unit is called, it may cause an overload of the relevant transmission lines, and this part of the reserve is usually considered an invalid reserve. To solve this problem, the common method is to determine the reserve demand by using the method of partition reserve (Ma et al., 2007; Ahmadi-Khatir et al., 2013; Zhou, 2017) but the method of partition reserve is only limited to the case where the power supply and load distribution are relatively concentrated. The network structure of China’s power system is complex, and the resources are scattered. It is difficult to ensure the safe and stable operation of the power system through zoning (Wenhuan, 2022) analyzes the existing engineering reserve deduction methods and puts forward a reserve evaluation algorithm, but it still fails to solve the problem of limited reserve in essence.

On the other hand, there are two main methods to deal with the uncertainty of renewable energy in power systems, robust optimization (Ben-Tal and Nemirovski, 1998; Bertsimas and Sim, 2004; Bertsimas et al., 2013) and stochastic optimization (Papavasiliou et al., 2011; Sahin et al., 2013; Ramezanzade et al., 2018; Zhao and Guan, 2018; Zhao, 2021). Stochastic optimization depends on the probability information of renewable energy output. Robust optimization only needs to describe the uncertainty interval of renewable energy output. Although the mathematical model of stochastic optimization is simple, it needs many calculations and is difficult to accurately obtain the probability information of renewable energy output. Although robust optimization has better robustness, it is difficult to convert the model. Different uncertainty sets have different robust equations (Ben-Tal, Ghaoui, Nemirovski), and it is difficult to solve the model. Some scholars proposed to combine robust optimization with stochastic optimization (Su, 2018; Xu et al., 2018; Su, 2019; Zhu et al., 2019), but the transformation theory of such models is often difficult to apply to SCUC problems. (Bertsimas et al., 2013). proposed a two-stage robust optimization problem for SCUC called Benders Decomposition to solve the model, but the convergence and solution speed of the model is difficult to guarantee.

Based on the clearing process of China’s electricity market, this paper adopts different optimization models at different stages to ensure the economy and security of the power system. The highlights of the article are as follows:

1. In the SCUC and SCED models, the operation reserve variables of the unit are introduced, and the operating reserve of the unit can be accurately obtained.

2. Generation shift distribution factor (GSDF) is used to construct network security constraints. In this way, the binary variables generated when the unit operation reserve variable is coupled with the network security constraint can be avoided. The modeling methods in extreme and conventional scenarios are given.

3. According to the clearing process of China’s electricity market, a multi-stage robust optimization model for promoting energy consumption is proposed, which effectively avoids the disadvantage that the decomposition algorithm cannot converge, and takes into account the efficiency, economy, and security of the model.

The rest of this article is organized as follows: In section 2, the model of effective reserve calculation in the electricity market is proposed. In section 3, the multi-stage robust optimization model is proposed. In section 4, the parity of the multi-stage robust optimization model is given. Section 5 is numerical simulation. The conclusion is given in section 6.

2 SCUC model for modeling effective reserve of unit

2.1 SCUC model

The objective function of day-ahead electricity market clearing SCUC is to minimize the system operation cost:

$\min {\displaystyle \sum _{i=1}^G{\displaystyle \sum _{t=1}^T\left[C_{i,t}\left(P_{i,t}\right)+C_{i,t}^U\right]}}+{\displaystyle \sum _{s=1}^{SE}{\displaystyle \sum _{t=1}^{T}M\left[S{L}_s^{+}+S{L}_s^{-}\right]}}(1)$

$\mathrm{s}.\mathrm{t}.{\displaystyle \sum _{i=1}^G{P}_{i,t}}=D_t(2)$

${\displaystyle \sum _{i=1}^G{\alpha }_{i,t}{P}_{i,t}^{\max }}\ge D_t+R_t^U(3)$

${\displaystyle \sum _{i=1}^G{\alpha }_{i,t}{P}_{i,t}^{\min }}\le D_t-R_t^D(4)$

${\alpha }_{i,t}{P}_{i,t}^{\min }\le P_{i,t}\le {\alpha }_{i,t}{P}_{i,t}^{\max },i\in G(5)$

$P_{i,t}-P_{i,t-1}\le \mathrm{\Delta }{P}_i^U{\alpha }_{i,t-1}+P_{i,t}^{\min }({\alpha }_{i,t}-{\alpha }_{i,t-1})+P_{i,t}^{\max }(1-{\alpha }_{i,t}),i\in G(6)$

$P_{i,t-1}-P_{i,t}\le \mathrm{\Delta }{P}_i^D{\alpha }_{i,t}-P_{i,t}^{\min }({\alpha }_{i,t}-{\alpha }_{i,t-1})+P_{i,t}^{\max }(1-{\alpha }_{i,t-1}),i\in G(7)$

$P_{i,t}={\displaystyle \sum _{n=1}^N{P}_{i,t,n}},i\in G(8)$

$C_{i,t}^U={\eta }_{i,t}{C}_i^U,i\in G(9)$

$P_s^{\min }\le {\displaystyle \sum _{i=1}^G{G}_{s-i}{P}_{i,t}}+{\displaystyle \sum _{j=1}^{NT}{G}_{s-j}{T}_{j,t}}-{\displaystyle \sum _{k=1}^K{G}_{s-k}{L}_{k,t}}-S{L}_s^{+}+S{L}_s^{-}\le P_s^{\max },s\in SE(10)$

The objective function 1) represents the minimum generation cost of the power system. Where $C_{i,t}\left(P_{i,t}\right)$、 $C_{i,t}^U$ respectively represent unit output cost and startup cost, $S{L}_s^{+}$ and $S{L}_s^{-}$ represent the slack variables of the branch/section respectively, and $M$ represents the penalty factor, which is usually a large positive integer.Where (2) represents system balance constraint, formula 3, 4 represent the positive and negative reserve constraints of the system, the constraints (5) represent the upper and lower limit constraints of unit output, the constraints (6) and (7) represent the ramp constraints of the unit, the constraints (8) represent unit segment output constraints, the constraints (9) represents the start-up cost of the unit, the constraints (10) represent network security constraint. Where ${\displaystyle \sum P_{i,t}}$ represents unit output and $D$ represents system load, $P^{\max }$ and $P^{\min }$ represent the upper and lower limits of unit output respectively, $R^U$ and $R^D$ respectively represent the positive and negative reserve capacity demand of the power system, $\alpha$ are binary variables, one refers to the startup unit, and 0 refers to the shutdown unit, $\mathrm{\Delta }{P}_i^U$ and $\mathrm{\Delta }{P}_i^D$ respectively represent the ramp speed of the unit, $P_{i,t,m}$ represents the segment output of the unit, and $N$ refers to the number of segmentation. $\eta$ are 0–1 variables, one represents the unit startup, and the rest is 0. $P_s^{\min }$ and $P_s^{\max }$ respectively represents the upper and lower limits of the transmission line $s$, $T$ refers to the power of tie line. The inflow is positive and the outflow is negative, $G_{s-i}$ represents the generator output power transfer distribution factor of line $s$ to node $i$.

2.2 Modeling considering effective reserve variables of units

At present, the calculation formula for the positive reserve of the power system is (3), but the section of the unit may be out of limit due to the provision of the reserve, that is:

$P_s+{\displaystyle \sum _{i=1}^{PRG}{G}_{s-i}(P_{i,t}^{\max }-\stackrel{\sim }{P_{i,t}}})\notin [P_s^{\min },P_s^{\max }],s\in SE(11)$

Where $P_s$ represents the current transmission power of the line, PRG represents the number of units to be provided with reserve and $\stackrel{\sim }{P_{i,t}}$ represents the current unit output. When the above conditions occur, it is considered that the units in the PRG set provide invalid reserve. To solve this problem, reserved variable $P{R}_{i,t}$ for each unit can provide reserve. $P{R}_{i,t}$ represents the maximum available reserve that the unit $i$ can provide at time $t$. Therefore, the positive reserve capacity formula is changed (3) to the following form:

${\displaystyle \sum _{i=1}^{N}P{R}_{i,t}}\ge R_t^U(12)$

When the unit is shut down, the reserve cannot be provided, so $P{R}_{i,t}$ should be coupled with the operating state of the unit,

$0\le P{R}_{i,t}\le {\alpha }_{i,t}{P}_{i,t}^{\max }(13)$

In addition, the upper limit constraint of $P{R}_{i,t}$ shall be less than the maximum output and actual output of the unit, and the following constraints shall be met,

$P{R}_{i,t}\le P_{i,t}^{\max }-P_{i,t}(14)$

Considering the constraints (12) (13), there is

$0\le P{R}_{i,t}\le {\alpha }_{i,t}{P}_{i,t}^{\max }-P_{i,t}(15)$

If considering the startup and shutdown curve of the unit, the unit cannot provide reserve during startup and shutdown, the formula should be modified as follows:

$\begin{array}{l}0\le P{R}_{i,t}\le {\alpha }_{i,t}{P}_{i,t}^{\max }-P_{i,t}-{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{tt=1}^{UD}(P_{i,\max }-P_{i,U}\left(\mathrm{t}\mathrm{t}\right)){\beta }_{i,t-tt+1}}}-\\ {\displaystyle \sum _{i=1}^N{\displaystyle \sum _{tt=1}^{DD}(P_{i,\max }-P_{i,D}\left(\mathrm{t}\mathrm{t}\right)){\gamma }_{i,t+DD-tt+1}}}\end{array}(16)$

$UD$ is the duration of the startup process, calculated to the minimum output of the unit; $DD$ is the duration of the shutdown process, calculated from the minimum output of the unit; $\gamma$ are binary variables, one represents unit shutdown, and the rest is 0.

2.3 Coupling constraints of unit effective reserve and network security constraints

To ensure that the reserve of the unit associated with the section out of limit is 0, it is usually necessary to introduce the if-then constraint. However, for the linearization of logical constraints, binary variables of $SE\ast T$ quantity need to be introduced, which will increase the complexity of the SCUC model by $2^{SE\ast T}$. Here, approximate modeling is given, which does not introduce too many binary variables, making the complexity of the model rise relatively low.

When the positive/negative section exceeds the limit, the corresponding units with positive/negative sensitivity cannot provide the corresponding reserve at this time. Therefore, it has the following forms:

$\begin{array}{l}\text{if}S{L}_s^{+}>0\text{\&\&}{G}_{s-i}>\mathrm{Threshold}\\ P{R}_i=0\\ \text{or}\\ \text{if}S{L}_s^{-}>0\text{\&\&}{G}_{s-i}<-\mathrm{Threshold}\\ P{R}_i=0\end{array}(17)$

(17) needs to be linearized as follows:

$\begin{array}{l}S{L}_s^{+}+S{L}_s^{-}>-M\left(1-{\tau }_s\right)\\ P{R}_i\le M\left(1-{\tau }_s\right),\left|G_{s-i}\right|>\mathrm{Threshold}\end{array}(18)$

${\tau }_s$ is a binary variable. When ${\tau }_s$ is 1, it means that section $s$ is out of limit, else ${\tau }_s$ is 0. Since the current single machine capacity in the market does not exceed 3,000, M here can be set to 3,000. The above method will introduce more binary variables. For each section, binary variables will be introduced at each time, which will greatly increase the difficulty of solving the model. A method without introducing binary variables is proposed as follows:

When considering the reserve network security constraints, the units with positive sensitivity and the units with negative sensitivity are added respectively. Use the following formula:

${\displaystyle \sum _{i=1}^G{G}_{s-i}{P}_{i,t}+{\displaystyle \sum _{i=1}^G\stackrel{+}{G_{s-i}}P{R}_{i,t}}}+{\displaystyle \sum _{j=1}^{NT}{G}_{s-j}{T}_{j,t}}-{\displaystyle \sum _{k=1}^K{G}_{s-k}{L}_{k,t}}-S{L}_s^{+}\le P_s^{\max },s\in SE(19)$

$P_S^{\min }\le {\displaystyle \sum _{i=1}^G{G}_{s-i}{P}_{i,t}+{\displaystyle \sum _{i=1}^G\stackrel{-}{G_{s-i}}P{R}_{i,t}}}+{\displaystyle \sum _{j=1}^{NT}{G}_{s-j}{T}_{j,t}}-{\displaystyle \sum _{k=1}^K{G}_{s-k}{L}_{k,t}}+S{L}_s^{-},s\in SE(20)$

$G_{s-i}^{+}$ and $G_{s-i}^{-}$ represent positive and negative GSDF above the threshold, respectively. The above formula can be abbreviated as follows:

$P_{s,t}^{current}+{\displaystyle \sum _{i=1}^G\stackrel{+}{G_{s-i}}P{R}_{i,t}}-S{L}_s^{+}\le P_s^{\max },s\in SE(21)$

$P_S^{\min }\le P_{s,t}^{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}}+{\displaystyle \sum _{i=1}^G\stackrel{-}{G_{s-i}}P{R}_{i,t}}+S{L}_s^{-},s\in SE(22)$

Where $P_{s,t}^{\text{current}}$ represents the current power flow of the section. For (21) ${\displaystyle \sum _{i=1}^G\stackrel{+}{G_{s-i}}P{R}_{i,t}}\ge 0$, if $P_{s,t}^{\text{current}}>P_s^{\max }$, then $S{L}_s^{+}>0$.In the objective function, $S{L}_s^{+}$ has a penalty term ${\displaystyle \sum _{s=1}^{SE}{\displaystyle \sum _{t=1}^{T}M\left[S{L}_s^{+}+S{L}_s^{-}\right]}}$, and the objective function is minimized, so $S{L}_s^{+}\to 0$. If $P_{s,t}^{\text{current}}>P_s^{\max }$, then $S{L}_s^{+}>0$. It can be inferred that when $P_{s,t}^{\text{current}}>P_s^{\max }$ is established, ${\displaystyle \sum _{i=1}^G\stackrel{+}{G_{s-i}}P{R}_{i,t}}=0$ is established, and because of $\stackrel{+}{G_{s-i}}>0$, $P{R}_{i,t}=0$ is established. For (22), the same principle can be proved. When $P_{s,t}^{\text{current}}<P_s^{\max }$, $P{R}_{i,t}=0$.It is worth noting that the above proof is limited to the case that the model is optimized and the system reserve is sufficient.

When considering the reserve network security constraints of the above formula, the reserve variables associated with the positive and negative reach bounds of the section will only aggravate the forward or negative power flow of the current section. Therefore, there is no need to add the original network security constraints (10) in this case.

In the above formula 19 (20), the influence components of standby variables above the forward and reverse thresholds on section power flow are added respectively, which will reduce the feasible region of the model to a certain extent. If the influence of correlation between reserve variables is considered, the following conditions may occur when the units higher than the positive sensitivity threshold and the units with negative sensitivity threshold are adjusted at the same time:

$P_s+{\displaystyle \sum _{i=1}^{PRG1}{G}_{s-i}^{+}(P_{i,t}^{\max }-\stackrel{\sim }{P_{i,t}}})+{\displaystyle \sum _{i=1}^{PRG2}{G}_{s-i}^{-}(P_{i,t}^{\max }-\stackrel{\sim }{P_{i,t}}})\le P_s,s\in SE(23)$

At this time, the network security constraint of the system needs to be rewritten into the following form:

$\begin{array}{l}{P}_s^{\min }\le {\displaystyle \sum _{i=1}^N{G}_{l-i}{P}_{i,t}+{\displaystyle \sum _{i=1}^N\stackrel{+}{G_{l-i}}P{R}_{i,t}}}+{\displaystyle \sum _{i=1}^N\stackrel{-}{G_{l-i}}P{R}_{i,t}}+{\displaystyle \sum _{j=1}^{NT}{G}_{l-j}{T}_{j,t}}\\ -{\displaystyle \sum _{k=1}^K{G}_{l-k}{L}_{k,t}}-S{L}_s^{+}+S{L}_s^{-}\le P_s^{\max },s\in SE\end{array}(24)$

The influence component of reserve on section power flow can be any value, which can make the section out of limit no longer without considering the reserve variables. In this situation, the original network security constraint (10) needs to be added.

3 Modeling of effective reserve considering renewable energy uncertainty

3.1 SCUC model of box uncertainty set

Because the model of SCUC is complex and contains many binary variables, it is difficult to solve the model. In the SCUC stage, the box robust optimization is adopted to consider the uncertainty of renewable energy and ensure the solution efficiency of the model. 5) is rewritten as follows:

${\alpha }_{i,t}{P}_{i,t}^{\min }\le P_{i,t}+{\upsilon }_i{\displaystyle \sum _{j=1}^{\mathrm{NE}}\left(\mathrm{\Delta }{P}_{j,t}^w\right)}{\zeta }_j\le {\alpha }_{i,t}{P}_{i,t}^{\max },i\in G(25)$

where NE represents the renewable energy unit set, $\mathrm{\Delta }{P}_{j,t}^w$ represents the maximum deviation between the predicted output of renewable energy and the actual output, $\zeta$ represents uncertainty $\zeta \in \left[\mathrm{0,1}\right]$, $\upsilon$ represents the participation factor of conventional units, which is not a variable, and ${\displaystyle \sum {\upsilon }_i=1}$.

The uncertainty of renewable energy is limited to the box uncertainty set as follows:

${\mathrm{ℛ}}^W\left({\text{Γ}}^W,{\mathrm{P}}_{\mathrm{i}}^{\mathrm{w}},\mathrm{\Delta }{\mathrm{P}}_{\mathrm{i}}^{\mathrm{w}}\right):=\{{\tilde{P}}_i^w:\mathit{\exists }{\zeta }_i\in \mathrm{ℝ},{\Vert {\zeta }_i\Vert }_{\mathit{\infty }}\le {\text{Γ}}^W,\mathit{\forall }i\in NE\}(26)$

where ${\text{Γ}}^W\in \left[\mathrm{0,1}\right]$ is the uncertainty level of the total renewable energy of the system.

3.2 SCED model of the robust stochastic optimization model

In the SCED stage, the model is a linear programming (LP) problem, and the efficiency of solving the model is no longer concerned. Therefore, more refined robust stochastic optimization (RSO) can be considered in this stage.

After considering the stochastic, the model becomes as follows:

$\begin{array}{l}\underset{x}{\min }\underset{P\in \mathcal{P}}{\mathrm{s}\mathrm{u}\mathrm{p}}{\mathbb{E}}_{\mathbb{P}}\left\{f\left(x,\zeta \right)\right\}\\ s.t.g\left(x,\zeta \right)\ge 0\\ f\left(x,\zeta \right)\in \mathcal{A}\left(K,\mathcal{I}\right)\end{array}(27)$

where $\parallel$ represents the ambiguity set of the model. When $\parallel$ can accurately describe the probability information, (25) is stochastic optimization. When $\parallel$ only contains uncertainty sets, (25) is robust optimization. In the SCED stage, the ambiguity sets of the RSO model are as follows:

$\mathcal{P}=\left\{\mathbb{P}\in {\mathcal{P}}_0\left({\mathbb{R}}^D+{\mathbb{R}}^L\right):\begin{array}{l}\left(\left(\mathbf{\xi },\mathbf{u}\right),\stackrel{\sim }{w}\right)\mathbb{\sim }\mathbb{P}\\ \begin{array}{l}{\mathbb{E}}_{\mathbb{P}}\left[\mathbf{u}|\stackrel{\sim }{w}\in {\left[S\right]}_k\right]\le \vartheta \\ {\mathbb{E}}_{\mathbb{P}}\left[\mathbf{\xi }\right]=\mu \end{array}\\ \mathbb{P}\left[\left(\mathbf{\xi },\mathbf{u}\right)\in C|\stackrel{\sim }{w}=w\right]=1,\mathbb{P}\left[\stackrel{\sim }{w}=w\right]=p_w\end{array}\right\}(28)$

where $\xi$ is the vector representation of $\zeta$, $u$ is the vector representation of $\mathrm{\Delta }{P}^w$, Where $w$ represents the real scene, $\vartheta$ is a distance scale, adjustable parameter. The specific form of confidence set $C$ is as follows:

$C=\{(\xi ,\mathbf{u}):\mu +\mathbf{u}\ge \mathbf{\xi },\mu +\mathbf{u}\le \xi \}(29)$

In the existing scheduling systems, it is difficult to obtain accurate probability distribution, but it is relatively easy to obtain typical scenarios. Therefore, the data-driven robust stochastic optimization based on data scenarios has better application scenarios.

Here, the method of extending typical scenes to related scenes is used to construct scene sets, as follows:

$\mathcal{A}(K,\mathit{ℐ})=\{f:\begin{array}{c}f\left(\mathbf{x},\mathbf{\xi }\right)=f_0+{\displaystyle \sum _{i\in \mathit{ℐ}}{f}_i(\mathbf{x},\mathbf{\xi })}\\ forsome{f}_0,f_i\left(x\right)\in \mathcal{A}\left(K\right)\end{array}\}(30)$

where $K$ represents the main scene and $\mathit{ℐ}$ represents the sub-scenarios extended from the main node. $\left(K\right)$ is an event-wise static adaptation. The formula is as follows:

$A(K)=\left\{f^{\prime }:\begin{array}{c}x\left(\xi \right)=x^{\tau }\\ \begin{array}{l}\tau =\mathit{ℱ}\left(\xi \right)\\ forsome{x}^{\tau }\in ℛ\end{array}\end{array}\right\}(31)$

$\left(K\right)$ scenes here are independent of each other and have no relevance.

4 Model reconstruction

In the stage of SCUC and SCED, the model cannot be solved directly because the uncertainty variables are considered in the model, so it is necessary to make a reformulation to the model.

4.1 Model reformulation of SCUC

Formulas 25, 26 respectively describe the unit output and box uncertainty set under the renewable energy output uncertainty, To convert it to a solvable form, the display filter is used for conversion, as follows (Bai et al., 2016):

$\max {\gamma }_i{\displaystyle \sum _{i\in \mathcal{W}}\mathrm{\Delta }{p}_i^w{\zeta }_i}={\Vert {\gamma }_i{\displaystyle \sum _{i\in \mathcal{W}}\mathrm{\Delta }{p}_i^w}\Vert }_{P^{\ast }}(32)$

where ${\Vert •\Vert }_{p^{\ast }}$ is the dual norm, $1/p+1/p\ast =1$. Therefore, (25) is transformed into the following form:

${\alpha }_{i,t}{P}_{i,t}^{\min }\le P_{i,t}+{\upsilon }_i{\Vert {\gamma }_i{\displaystyle \sum _{i\in \mathcal{W}}\text{Δ}{p}_i^w}\Vert }_{P^{\ast }}\le {\alpha }_{i,t}{P}_{i,t}^{\max },i\in G(33)$

p* is one norm. If the influence of $PR$ is considered, (33) comes as follows:

${\alpha }_{i,t}{P}_{i,t}^{\min }\le P_{i,t}+{\upsilon }_i{\Vert {\gamma }_i{\displaystyle \sum _{i\in \mathcal{W}}\mathrm{\Delta }{p}_i^w}\Vert }_{P^{\ast }}+P{R}_{i,t}\le {\alpha }_{i,t}{P}_{i,t}^{\max },i\in G(34)$

15) and 16) comes as follows:

$0\le P{R}_{i,t}\le {\alpha }_{i,t}{P}_{i,t}^{\max }-P_{i,t}-{\upsilon }_i{\Vert {\gamma }_i{\displaystyle \sum _{i\in \mathcal{W}}\mathrm{\Delta }{p}_i^w}\Vert }_{P^{\ast }}(35)$

$\begin{array}{l}0\le P{R}_{i,t}\le {\alpha }_{i,t}{P}_{i,t}^{\max }-P_{i,t}-{\upsilon }_i{\Vert {\gamma }_i{\displaystyle \sum _{i\in \mathcal{W}}\mathrm{\Delta }{p}_i^w}\Vert }_{P^{\ast }}-{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{tt=1}^{UD}(P_{i,\max }-P_{i,U}\left(\mathrm{t}\mathrm{t}\right)){\beta }_{i,t-tt+1}}}-\\ {\displaystyle \sum _{i=1}^N{\displaystyle \sum _{tt=1}^{DD}(P_{i,\max }-P_{i,D}\left(\mathrm{tt}\right)){\gamma }_{i,t+DD-tt+1}}}\end{array}(36)$

Now, the robust model of the SCUC stage is transformed into a mixed-integer linear programming (MILP) model. So far, this paper proposes a robust optimization SCUC model ((1)–(2), (4), (6)–(9) (12) (34) (35)/(36)) in multiple scenarios. When binary variables are introduced 18) need to be added; Without introducing binary variables (19)–(20) need to be added in extreme scenarios, and (10) 24) in non-extreme scenarios.

4.2 Model reformulation of SCED

The biggest difference between the SCED model and the SCUC model is that there is an expectation in the objective function of the SCED model, and the expectation function needs to be linearized. Here, the dual theory is used for linearization.

The following transformation theory is given (Su et al., 2019; Chen et al., 2022):

$\begin{array}{l}\underset{\mathbb{P}\left(\xi \right)\in P}{\mathrm{s}\mathrm{u}\mathrm{p}}{\mathbb{E}}_{\mathbb{P}}\left[r^T\left(\stackrel{\sim }{w}\right)G_m\left(\stackrel{\sim }{w}\right)\stackrel{\sim }{z}+h_m\left(\stackrel{\sim }{w}\right)\right]\\ \Downarrow \\ \mathrm{i}\mathrm{n}\mathrm{f}\gamma \\ s.t.\gamma \ge \alpha p+{\displaystyle \sum }_{k\in \left[K\right]}{\beta }_k{\mu }_k\forall p\in \mathfrak{P},\frac{{\mu }_k}{{\displaystyle \sum }_{w\in {\mathfrak{E}}_k}{p}_w}\in Q_k,k\in \left[K\right]\\ {\alpha }_w+{\displaystyle \sum }_{k\in {\mathfrak{K}}_s}{\beta }_{k}z\ge r\left(w\right)G_m\left(w\right)z+h_m\left(w\right)\forall z\in Z_w,w\in \left[S\right]\\ \gamma \in \mathbb{R},\alpha \in {\mathbb{R}}^S,{\beta }_k\in {\mathbb{R}}^{I_z}\forall k\in \left[K\right]\end{array}(37)$

where $K$ represents scenarios set.

In the SCED stage, the start and stop of the unit have been determined, so the cost of this part is constant. The objective function is as follows:

$\min {\displaystyle \sum _{i=1}^G{\displaystyle \sum _{t=1}^T\left[C_{i,t}\left(P_{i,t}\right)\right]}}+{\displaystyle \sum _{s=1}^{SE}{\displaystyle \sum _{t=1}^{T}M\left[S{L}_s^{+}+S{L}_s^{-}\right]}}(38)$

In the SCED stage, considering the influence of ambiguity sets, the objective function becomes the following form:

$\underset{x}{\min }\underset{P\in \parallel }{\sup }{\displaystyle \sum _{i=1}^G{\displaystyle \sum _{t=1}^T\left[C_{i,t}\left(P_{i,t}\right)\right]}}(39)$

In the SCED stage, the section is usually adjusted to the state of not exceeding the limit, that is, the SL variables are 0. To solve it, need to transform the expectation into a linearly solvable form.

Therefore, the SCED model can be reformulated as follows:

$\begin{array}{l}\underset{\mathbb{P}\left(\xi \right)\in P}{\max }{\mathbb{E}}_{\mathbb{P}}\left[{\displaystyle \sum }_{i=1}^G{\displaystyle \sum }_{t=1}^T{C}_{i,t}\left(P_{i,t}\right)\right]\\ =\{\begin{array}{l}\mathrm{i}\mathrm{n}\mathrm{f}\gamma \\ s.t.\gamma \ge {\alpha }^{T}p+{\displaystyle \sum }_{k\in \left[K\right]}{\beta }_k^T{\sigma }_k\forall p\in P,\frac{{\sigma }_k}{{\displaystyle \sum }_{w\in {ϵ}_k}{p}_w}\in Q_k,k\in \left[K\right]\\ \alpha +{\displaystyle \sum }_{k\in K_s}{\beta }_k^T\Delta \xi \ge {\displaystyle \sum }_{i=1}^G{\displaystyle \sum }_{t=1}^T{C}_{i,t}\left(P_{i,t}\right)+b^T\Delta \xi \\ \forall z\in Z_w,w\in \left[S\right]\\ \gamma ,\alpha ,\beta \in \mathbb{R}\end{array}\end{array}(40)$

Now, the model reconstruction of SCUC and SCED has been completed. The model can be solved directly by calling the solver.

The clearing process of the electricity market after considering the uncertainty set is shown in Figure 1 and the following process:

Step 1: Solution of the SCUC model considering the box uncertainty set and effective reserve modeling. Determine the startup and shutdown state of the unit.

Step 2: According to the typical scenario of renewable energy output. Constructing uncertainty sets of event expressions.

Step 3: Adjust the boundary data cleared by the power system, the planned output of the unit and tie line, etc.

Step 4: Calculate the RSO-SCED model. Determine whether the calculation result satisfies the standard issued by the scheduling agency, if not, go back to step 2, if so, publish the result.

FIGURE 1

FIGURE 1. Day-ahead electricity market clearing process.

5 Numerical simulation

In this section, the simulation analysis is based on the IEEE 300-bus system and the actual operation data of a province in China. All simulation analysis is based on MATLAB or C++. All simulation analyses were performed on a laptop. Processor model: i7-8565u, 1.8GHz, 16 g memory.

To test the engineering application value of effective reserve modeling, the test is carried out in C++ based on the actual electricity market operation data. Due to the lack of data on typical daily operation scenarios of renewable energy in the actual operation data, the simulation analysis for multi-stage robust optimization is carried out based on MATLAB with the IEEE 300-bus system, and the typical daily operation scenarios of renewable energy are generated by the Monte Carlo method. In the IEEE 300-bus system, Renewable energy is connected at nodes 10, 30, 52, 57, 68, 81, 169, 175, 194, 265, and 271 with each capacity of 80 MW.

5.1 Comparative analysis of effective reserve calculation models

PRSCUC is used to represent the SCUC model after considering the modeling (19) (20) of reserve variables, BSCUC represents a model that calculates the effective reserve by introducing binary variables (18), and sPRSCUC is used to represent the SCUC model after slack the feasible region (22) (10), the SCUC model after the reserve variables can affect each other. SCED uses the same naming logic in Figure 2.

FIGURE 2

FIGURE 2. Comparison value of objective functions of different models under different typical scenarios.

Where * means that the model can find a feasible solution in 1800s. As can be seen from Table 1, The BSCUC model with binary variables will impose a great burden on the calculation efficiency, and often the model cannot be solved normally. In actual operation, this situation is not allowed. PRSCUC and sPRSCUC models will make the model speed longer, such as scenario 1, scenario 4, and scenario 8. However, there are also cases where the model solving speed is significantly improved, such as scenario seven This phenomenon is caused by the introduction of new continuous variables, which increases the relevant constraints of the SCUC model and reduces the feasible region. The solver may spend more time and find the optimal solution faster when using branch and bound to solve MILP problems. However, the order of magnitude of the model solution time has not changed, and it belongs to an acceptable fluctuation range from the perspective of engineering applications.

TABLE 1

TABLE 1. Comparison results of calculation time of different models under different typical scenarios.

Figure 2 shows the changes in objective function values of several models in different scenarios. In the SCUC stage, because the slack penalty factor of network security constraints in the objective function value is large, the objective function values of the three models are quite different. In the SCED stage, the slack penalty factor of network security constraints is small, and the difference between the objective function values of different models is not obvious. But in general, the objective function values of the SPR. model and PR. model are higher than those of the original model. This is because the refinement considers the effective reserve value of each unit and the feasible region of the SPR. model and PR. model is smaller than that of the original model, so the objective function value is on the high side as a whole. In scenarios 1, 4, and 8, the difference between the optimization objective functions of the three models is very small, in these scenarios, there is less power load, the transmission pressure of the power grid is small, and most of the transmission lines are in the light load state.

Table 2 shows the efficiency of the system under different models, and the specific calculation formula is 1. Both PRSCUC and SPRSCUC models have high efficiency. Among them, it is worth noting that in scenario 6, in which the system load is high and the transmission lines are generally overloaded, so it is difficult to provide sufficient effective reserve. After reducing the feasible region, the refined modeling part can only provide 71% of the effective reserve. However, considering the cooperation between units, the PRSCUC model can satisfy almost all scenarios.

TABLE 2

TABLE 2. Comparison of effective reserve of different models.

5.2 Comparative analysis of multi-stage robust optimization and conventional robust optimization

To verify the influence of several robust optimizations on the clearing model, a comparative analysis is carried out based on IEEE300. CM represents the clearing model, and PCM represents the clearing model that considers the uncertainty of renewable energy output to improve the reserve demand. ROCM represents the clearing model under box robust optimization. RSOCM indicates that box type robust optimization is adopted in the SCUC stage and the RSO clearing model is adopted in the SCED stage. The reserve demand is 10% of the system load. Considering the extreme scenario that all renewable energy sources are unable to provide output in PCM, the reserve demand increases the total capacity of renewable energy output based on the original demand. The distance parameter $\vartheta$ of the RSOCM model is 1.

As shown in Figure 3, RSOCM approximates from ROCM results to CM results as the number of samples increases. However, when the sample volume is greater than 1,500, the objective function changes slightly, because the uncertainty of 1,500 renewable energy output under the current system has been described as accurately as possible. Although cm has better optimization results, the system has low-risk resistance and poor robustness of clearing results. PCM has the best robustness because it considers the higher reserve demand, but the system economy is poor.

FIGURE 3

FIGURE 3. Comparison of objective functions of different models under different sample volumes.

Since the differences between the above models are all in the SCED stage, a comparative analysis is made only for the calculation speed of RSOCM and other models in the SCED stage. The details are as follows.

Figure 4 shows the time from modeling to the solution of different models. Since the model itself is small and LP problem, the modeling time is considered. Except for the RSOCM model, the solution time of other models is very close, about 1.25s. The RSOCM model takes longer to calculate, but the overall magnitude does not change.

FIGURE 4

FIGURE 4. Comparison of calculation time of different models under different sample volumes.

5.3 Comparative analysis of different models under different renewable energy penetration rates

With the development of the power system, the proportion of renewable energy access will gradually increase. Here, simulation analysis is made for the power system under different renewable energy penetration rates, as follows:

Figure 5 shows the changes in objective functions of several models with the increasing penetration of clean energy. As can be seen from the figure, since the cost of clean energy is much lower than that of conventional units, the objective function value of the CM model is declining. After considering the impact of the reserve, the objective function value of the PCM model is improved. The optimization result of PCM is similar to that of ROCM, but the difference is small. In addition, with the increase of renewable energy permeability, the value of the objective function of the RSOCM model decreases slowly. It is proved that the RSOCM model performs well when clean energy has a high penetration rate.

FIGURE 5

FIGURE 5. Comparison of objective functions of different models under different renewable energy penetration rates.

6 Conclusion

In this paper, a calculation method of system effective reserve is proposed and embedded into the robust optimization model of SCUC and SCED. Different robust optimization methods are adopted in SCUC and SCED respectively to ensure the security and economy of the power system while the solving efficiency requirement is met. The effectiveness of the proposed model is verified by an actual operation case and the IEEE 300-node case. The simulation analysis, it shows that the proposed model considering reserve calculation can well ensure the safe operation of the power system. The proposed multi-stage robust optimization model matches the existing clearing process of the electricity market, and the model has good extensibility. In the future, the uncertainty set in different stages can be changed based on the needs of the operating entity to ensure the economy and security of the system operation.

Data availability statement

The data presented in the study are deposited in the Matpower repository, accession number https://github.com/MATPOWER/matpower/tree/master/data.

Author contributions

QD is responsible for writing the whole article and model derivation. ZC is responsible for technical guidance. XT, XX and HC are responsible for the layout adjustment and code debugging of the article. XS is responsible for the overall idea of the article.

Conflict of interest

Authors QD, ZC, XT, XX, HC and XS were employed by Electric Power Dispatching Control Center of Guangdong Power Grid Co, Ltd and Beijing Tsintergy Technology Co, Ltd

The author declares that this research was funded by key scientific and technological projects of Guangdong Power Grid Corporation 036000KK52210008 (GDKJXM20200713). The funders participated in the following research: research design, collection, analysis and writing.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ahmadi-Khatir, A., Bozorg, M., and Cherkaoui, R. (2013). Probabilistic spinning reserve provision model in multi-control zone power system. IEEE Trans. Power Syst. 28, 2819–2829. doi:10.1109/tpwrs.2013.2243923

CrossRef Full Text | Google Scholar

Bai, X., Qu, L., and Qiao, W. (2016). Robust AC optimal power flow for power networks with wind power generation. IEEE Trans. Power Syst. 31 (5), 4163–4164. doi:10.1109/tpwrs.2015.2493778

CrossRef Full Text | Google Scholar

Ben-Tal, A., Ghaoui, L. E., and Nemirovski, A. Robust optimization. princeton university press princeton nj.

Ben-Tal, A., and Nemirovski, A. (1998). Robust convex optimization. Math. OR. 23 (4), 769–805. doi:10.1287/moor.23.4.769

CrossRef Full Text | Google Scholar

Bertsimas, D., Litvinov, E., Sun, X. A., Zhao, J., and Zheng, T. (2013). Adaptive robust optimization for the security constrained unit commitment problem. IEEE Trans. Power Syst. 28 (1), 52–63. doi:10.1109/tpwrs.2012.2205021

CrossRef Full Text | Google Scholar

Bertsimas, D., and Sim, M. (2004). The price of robustness. Operations Res. 52, 35–53. doi:10.1287/opre.1030.0065

CrossRef Full Text | Google Scholar

Chen, Z., Melvyn, S., and Xiong, P. (2022). Optimization online. Robust stochastic optimization made easy with RSOME. Available: http://www.optimization-online.org/DB_FILE/2017/06/6055.pdf.

Google Scholar

Farhani, S. (2015). Renewable energy consumption, economic growth and CO2 emissions: Evidence from selected MENA countries. Work. Pap.

Google Scholar

Lu, T., Zhang, W., Wang, Y., Xie, H., and Ding, X. (2022). Medium- and long-term trading strategies for large electricity retailers in China's electricity market. Energies 15, 3342. doi:10.3390/en15093342

CrossRef Full Text | Google Scholar

Ma, X., Sun, D., and Cheung, K. (2007). Energy and reserve dispatch in a multi-zone electricity market. IEEE Trans. Power Syst. 14 (3), 913–919. doi:10.1109/59.780903

CrossRef Full Text | Google Scholar

Papavasiliou, A., Oren, S. S., and O'Neill, R. P. (2011). Reserve requirements for wind power integration: A scenario-based stochastic programming framework. IEEE Trans. Power Syst. 26 (4), 2197–2206. doi:10.1109/tpwrs.2011.2121095

CrossRef Full Text | Google Scholar

Ramezanzade, M., Jafari-Nokandi, M., and Barforoshi, T. (2018). Scheduling of generation and reserve of thermal generation resources considering load and wind uncertainty in presence of energy storage and demand response. Tehran, Iran: Faculty of Electrical & Computer Engineering.2

Google Scholar

Sahin, C., Shahidehpour, M., and Erkmen, I. (2013). Allocation of hourly reserve versus demand response for security-constrained scheduling of stochastic wind energy. IEEE Trans. Sustain. Energy 4 (1), 219–228. doi:10.1109/tste.2012.2213849

CrossRef Full Text | Google Scholar

Su, X. (2018). “A novel data-driven robust method applied to OPF with uncertain wind power,” in 2018 2nd IEEE conference on energy internet and energy system integration (EI2), 1–5.

CrossRef Full Text | Google Scholar

Su, X., Bai, X., Liu, C., Zhu, R., and Wei, C. (2019). Research on robust stochastic dynamic economic dispatch model considering the uncertainty of wind power. IEEE Access 7 ({}), 147453–147461. doi:10.1109/access.2019.2946460

CrossRef Full Text | Google Scholar

Su, X. (2019). “Data-driven robust dispatch of integrated electricity-gas energy systems considering uncertainty of wind power,” in 2019 IEEE power & energy society general meeting (PESGM), 1–5.

CrossRef Full Text | Google Scholar

Wenhuan, H. (2022). Research on deduction algorithm for limited power system reserve capacity network. Electr. Eng. (09), 12–14.

Google Scholar

Xu, H., Liu, Y., and Sun, H. (2018). Distributionally robust optimization with matrix moment constraints: Lagrange duality and cutting plane methods. Math. Program. 169 (2), 489–529. doi:10.1007/s10107-017-1143-6

CrossRef Full Text | Google Scholar

Zhao, C., and Guan, Y. (2018). Data-driven risk-averse stochastic optimization with Wasserstein metric. Operations Res. Lett. 46 (2), 262–267. doi:10.1016/j.orl.2018.01.011

CrossRef Full Text | Google Scholar

Zhao, D., X. (2021). Statistical measure for risk-seeking stochastic wind power offering strategies in electricity markets. J. Mod. Power Syst. Clean Energy, 1–6.

Google Scholar

Zhou, M. (2017). “Multi-area generation-reserve joint dispatch approach considering wind power cross-regional accommodation,” in Power and energy systems (Beijing, China: CSEE Journal of).

CrossRef Full Text | Google Scholar

Zhu, R., Wei, H., and Bai, X. (2019). Wasserstein metric based distributionally robust approximate framework for unit commitment. IEEE Trans. Power Syst. 34, 2991–3001. doi:10.1109/tpwrs.2019.2893296

CrossRef Full Text | Google Scholar