Quantification of Supply and Demand and Distribution Robust Optimal Scheduling Considering New Energy Power System Flexibility

1. Introduction

The rapid development of renewable energy power generation equipment and grid-connected technology has led to the gradual formation of a new type of power system dominated by new energy sources [1,2]. The stochastic volatility of a high proportion of new energy leads to a significant increase in the difficulty of power system operation and scheduling, and it is difficult for traditional flexibility resources, such as conventional units, to effectively meet the dramatically increased flexibility demand of the system. At the same time, the current situation that the power grid is seriously insufficient in local time, makes the demand-side flexibility gradually become the core of the research on power grid scheduling and operation. The demand-side flexibility has gradually become the core of the research on grid scheduling and operation.

Flexibility, as delineated by the International Energy Agency (IEA), represents the capacity of grid resources to address flexibility demands within specified timeframes [3,4]. Understanding and quantifying this demand is pivotal for analyzing the supply-demand equilibrium of flexibility in modern power systems. While several studies have ventured into quantifying flexibility demand [5,6,7], they have largely overlooked the critical role of uncertainty. To address this stochastic nature of net load, some research efforts have employed joint probability distribution functions [8,9,10]. However, these probability-based methodologies are inherently limited by their reliance on predefined probability distributions, which often fail to capture the precise flexibility demand and regulation capacity required in each time period. Additionally, attempts to quantify uncertainty in flexibility demand using an error coefficient based on forecast values, as demonstrated in [11], although a step forward, remain overly simplistic and insufficient to fully grasp the complexity and dynamics of uncertainty in power systems. Thus, there is a pressing need for more sophisticated approaches to accurately quantify and manage flexibility demand, taking into account the multifaceted nature of uncertainty.

In recent years, robust optimization (RO) and stochastic optimization (SO) have been applied to the study of flexibility supply and demand equilibrium under the influence of uncertainties. Literature [12,13] constructs the operational domain of flexibility supply and demand equilibrium through robust and conditional value-at-risk. Literature [14,15] establishes a comprehensive stochastic optimization model by solving the flexibility margin through the multi-scenario method and determining the benchmark value according to the Bernstein polynomial theory. All the above methods have the problem of being too conservative or risky in the portrayal of flexibility supply and demand. The distributionally robust optimization (DRO) model combines the strengths of both models by using uncertain probability densities to describe uncertainty and thus find the probability distribution of the worst-case scenario. Traditionally, DRO optimizes fuzzy sets of probability distributions based on moment information or distance measures such as KL (Kullback-Leibler ) dispersion [16]. The former is usually a non-deterministic polynomial (NP) hard problem that is difficult to solve, requires the type of distribution of uncertainty to be specified in advance, and is too conservative when the sample data are large; the latter has the risk of losing the significance of the probability fuzzy set when the assumption error is large. Based on this, a data-driven DRO has been emphasized [17]. The model characterizes uncertainty by extracting a large number of historical data samples of power system operation and establishing a comprehensive set of paradigm uncertainty probabilities. Decision making is guided by historical data rather than moment information that characterizes the overall performance, which makes full use of the sample data and has the advantages of simple solution and good economic flexibility [18]. At present, data-driven DRO has been initially applied in the fields of site selection and capacity determination [19] and unit combination [20], but there is no research report on flexibility supply and demand optimization scheduling.

Demand-side resources have emerged as viable alternatives to generation-side resources in providing flexibility balancing, offering operational adequacy and economic benefits [21]. Unfortunately, prior research has overlooked the substantial flexibility supply potential of demand-side resources. The majority of studies examining their involvement in power system optimal scheduling have been confined to the narrow focus of “power balance,” treating it as a static process of supply and demand. This narrow perspective fails to account for the dynamic fluctuations in supply and demand across various time scales, thereby underutilizing the regulatory potential of demand-side flexibility resources. As the operational paradigm of power systems evolves from the traditional one-way model of power supply matching load to a more sophisticated two-way model encompassing source-network-load-storage coordination and interaction, grid scheduling must adapt. It must harness the full spectrum of demand-side resources to achieve “flexibility balance,” ensuring that the system can meet the flexibility demands of each period effectively [22,23]. This evolution underscores the urgency for research that fully integrates and leverages the flexibility capabilities of demand-side resources.

In addition, the existing supply and demand flexibility studies rarely consider the correlation of wind and light outputs, and the traditional correlation analysis method is difficult to effectively reflect the complementary characteristics of new energy outputs, resulting in the inaccurate description of the set of outputs [24]. Based on the above problems, this paper establishes a data-driven distribution robust optimization scheduling model considering flexibility supply and demand balance, and the main work is as follows:

1) The Copula theory is applied to construct the scenariory spatio-temporal correlation output set, and the flexibility demand is quantified by combining the scenario method and the interval method;

2) Considering the flexibility supply capacity of demand-side resources such as electric vehicles, the effectiveness of the model is verified by comparing multiple flexibility regulation methods;

3) The proposed data-driven DRO model can effectively control the robustness of grid dispatch flexibility and economy, and the superiority of the model is verified by comparing with the existing uncertainty algorithms.

2. Power System Supply and Demand Flexibility and Quantitative Model

As a supplement to the abundance, the flexibility of the power system studies the fluctuation of power under a certain time scale. The dual fluctuation of source and load requires the use of flexible resources to match. If the system is not flexible enough, there will be risks of wind curtailment and load shedding.

To quantify the flexibility of the power system, it is necessary to determine the flexibility demand and the flexibility supply capacity of the flexibility resources respectively. Considering the time scale characteristics, directionality and state dependence of flexibility, the explanation is given from the perspectives of flexibility demand, flexibility supply and flexibility balance.

2.1. Quantification of Flexibility Requirements Considering Spatiotemporal Correlation and Uncertainty

2.1.1. Wind and Photovoltaic Power Generation based on Copula Theory

Copula theory states that the joint probability distribution function can be decomposed into the product of the marginal distribution function and the Copula function. Compared with the traditional simple linear correlation analysis, Copula theory can describe the complex and nonlinear correlation between different random variables. The selection of Copula function is related to the accuracy of the wind and photovoltaic power output set. Considering the complementary characteristics between wind and photovoltaic power output, the Frank-Copula function is selected for description [26]. The specific steps are as follows:

1) The wind power output probability density f_{W, t}(x_{W, t}), distribution function F_{W, t}(x_{W, t}), and photovoltaic output probability density f_PV,t(x_{PV, t}), distribution function F_{PV, t}(x_{PV, t}) were calculated based on the historical sample data.

2) The Copula function is constructed to obtain the joint probability distribution H (x_PV, x_W) of the wind-wind output in each period, and generate a set of random values of the edge distribution function of wind-wind variables with spatiotemporal correlation {u_{PV, t}, u_{W, t}}.

3) The set of random values is inverted by the inverse function, that is, $F_{\mathrm{PV},t}^{-1}$ (u_{PV, t}) and $F_{\mathrm{W},t}^{-1}$ (u_{W, t}), and the wind-photovoltaic output set with temporal and spatial dependence is obtained.

2.1.2. Quantifying Flexibility Requirements

The flexibility demand of power grid comes from the random fluctuation of net load obtained from wind-wind load synthesis. In the Data Driven DRO, the need for power system flexibility is quantified using a combination of scenario and interval methods, thereby reducing conservatism while ensuring flexible operation in the day-ahead phase. In this paper, the probabilistic distance based scenario reduction technique is used to obtain typical scenarios from the sample data .

Scenario reduction method based on probability distance is a kind of rapid previous scenario reduction technology, which can find a given number of typical scenarios in a large number of historical sample scenarios, and obtain the occurrence probability of each scenario. Compared with the typical scenario obtained by the clustering method, which is the cluster center, the scenarios obtained by the cutting method are all from the real scenario of the sample, which makes the typical scenario more authentic. The reduction steps are as follows:

Step 1: Determine the cut scenario $s^{(k^{\ast })},k^{\ast }\in (1,\cdots ,\mathrm{K})$. Calculate the geometric distance d(s⁽ⁿ⁾,s^(m)). between any two scenarios in the history sample S. Considering the probability distance D⁽ⁿ⁾ between scenario probability p⁽ⁿ⁾ and Euclidean distance d, the scenario s^(k’) with the smallest sum of probability distance from the remaining scenario s^(k*) is selected.

$$s^{(k^{\prime })}=\underset{n\ne m,m\in [1,\cdots ,K]}{\min }\left\{\arg \sum _{m=1}^K{p}^{(n)}d\left(s^{(n)},s^{(m)}\right)\right\}.$$

(1)

Step 2: Change the probability p^{(k ')} of the scenario closest to the probability s^{(k ')} of the deleted scenario

$$p^{(-k^{\prime })}=p^{(-k^{\prime })}+p^{(-k^{\ast })}$$

(2)

Step 3: Check whether the number of remaining scenarios meets requirements. Otherwise, return to Step 1.

Taking the load as an example, the load model considering the maximum fluctuation error for the kth scenario can be expressed as:

$$\left\{\begin{array}{c}{P}_{\mathrm{L},t,k}^{\max }=(1+{\epsilon }_{\mathrm{L}})P_{\mathrm{L},t,k}^0\\ P_{\mathrm{L},t,k}^{\min }=(1-{\epsilon }_{\mathrm{L}})P_{\mathrm{L},t,k}^0\end{array}\right.$$

(3)

where: $P_{\mathrm{L},t,k}^0$ is the predicted value of load at time t in the kth scenario; $P_{\mathrm{L},t,k}^{\max }$ and $P_{\mathrm{L},t,k}^{\min }$ are the upper and lower limits of load fluctuation at time t in the kth scenario respectively. ε_L is the maximum predictive error coefficient of the load.

Under the given time scale τ, the flexibility demand for scenario k can be divided into three cases: when $P_{\mathrm{L},t,k}^0$≤$P_{\mathrm{L},t+\tau ,k}^{\min }$, there is only upward flexibility demand; When $P_{\mathrm{L},t+\tau ,k}^{\max }$≤$P_{\mathrm{L},t,k}^0$, there is only downward flexibility demand; When $P_{\mathrm{L},t+\tau ,k}^{\min }$<$P_{\mathrm{L},t,k}^0$<$P_{\mathrm{L},t+\tau ,k}^{\max }$, there is bidirectional flexibility demand, that is, there is both upward flexibility demand and downward flexibility demand. Considering the coupling relationship of time, load fluctuation can be represented by the first-order difference of adjacent periods, as shown in equation (4), and the quantification principle of flexibility demand is shown in Figure 2.

$$\left\{\begin{array}{c}{F}_{\mathrm{L},t,k}^{\mathrm{de},\mathrm{up}}=\max (0,P_{\mathrm{L},t+\tau ,k}^{\max }-P_{\mathrm{L},t,k}^0)\\ F_{\mathrm{L},t,k}^{\mathrm{de},\mathrm{dn}}=\max (0,P_{\mathrm{L},t,k}^0-P_{\mathrm{L},t+\tau ,k}^{\min })\end{array}\right.$$

(4)

Where: $F_{\mathrm{L},t,k}^{\mathrm{de},\mathrm{up}}$and$F_{\mathrm{L},t,k}^{\mathrm{de},dn}$are respectively the upward flexibility demand and downward flexibility demand generated by load at time t in the kth scenario.

Figure 1. Quantification diagram of flexibility demands .

Figure 1. Quantification diagram of flexibility demands .

The process of quantifying the flexibility demand generated by wind power and photovoltaic is similar to that of load, but it should be noted that as a power source its downward fluctuation actually produces an upward flexibility demand, as shown in equations (5) and (6).

$$\left\{\begin{array}{l}{P}_{\mathrm{PV},t,k}^{\max }=(1+{\epsilon }_{\mathrm{PV}})P_{\mathrm{PV},t,k}^0\\ P_{\mathrm{PV},t,k}^{\min }=(1-{\epsilon }_{\mathrm{PV}})P_{\mathrm{PV},t,k}^0\\ F_{\mathrm{PV},t,k}^{\mathrm{de},\mathrm{up}}=\max (0,P_{\mathrm{PV},t,k}^0-P_{\mathrm{PV},t+\tau ,k}^{\min })\\ F_{\mathrm{PV},t,k}^{\mathrm{de},\mathrm{dn}}=\max (0,P_{\mathrm{PV},t+\tau ,k}^{\max }-P_{\mathrm{PV},t,k}^0)\end{array}\right.$$

(5)

$$\left\{\begin{array}{l}{P}_{\mathrm{W},t,k}^{\max }=(1+{\epsilon }_{\mathrm{W}})P_{\mathrm{W},t,k}^0\\ P_{\mathrm{W},t,k}^{\min }=(1-{\epsilon }_{\mathrm{W}})P_{\mathrm{W},t,k}^0\\ F_{\mathrm{W},t,k}^{\mathrm{de},\mathrm{up}}=\max (0,P_{\mathrm{W},t,k}^0-P_{\mathrm{W},t+\tau ,k}^{\min })\\ F_{\mathrm{W},t,k}^{\mathrm{de},\mathrm{dn}}=\max (0,P_{\mathrm{W},t+\tau ,k}^{\max }-P_{\mathrm{W},t,k}^0)\end{array}\right.$$

(6)

where: $P_{\mathrm{PV},t,k}^{\max }$ and $P_{\mathrm{PV},t,k}^{\min }$ is the upper and lower limits of PV at time t in the kth scenario respectively;$P_{W,t,k}^{\max }$ and $P_{W,t,k}^{\min }$ are respectively the upper and lower limits of wind power output at time t in the kth scenario.$P_{\mathrm{PV},t,k}^0$and $P_{\mathrm{W},t,k}^0$are the predicted values of PV and wind power output at time t in the kth scenario, respectively.$F_{\mathrm{PV},t,k}^{\mathrm{de},\mathrm{up}}$ and $F_{\mathrm{PV},t,k}^{\mathrm{de},\mathrm{dn}}$ are respectively the upward and downward flexibility requirements generated by PV at time t in the kth scenario.$F_{\mathrm{W},t,k}^{\mathrm{de},\mathrm{up}}$ and $F_{\mathrm{W},t,k}^{\mathrm{de},\mathrm{up}}$ are the upward and downward flexibility demands of wind power at time t in the kth scenario, respectively. ε_PV is the maximum prediction error coefficient of PV; ε_W is the maximum prediction error coefficient of wind power.

$$\left\{\begin{array}{l}{F}_{t,k}^{\mathrm{de},\mathrm{up}}=F_{\mathrm{L},t,k}^{\mathrm{de},\mathrm{up}}+F_{\mathrm{PV},t,k}^{\mathrm{de},\mathrm{up}}+F_{\mathrm{W},t,k}^{\mathrm{de},\mathrm{up}}\\ F_{t,k}^{\mathrm{de},\mathrm{dn}}=F_{\mathrm{L},t,k}^{\mathrm{de},\mathrm{dn}}+F_{\mathrm{PV},t,k}^{\mathrm{de},\mathrm{dn}}+F_{\mathrm{W},t,k}^{\mathrm{de},\mathrm{dn}}\end{array}\right.$$

(7)

Where: $F_{t,k}^{\mathrm{de},\mathrm{up}}$and$F_{t,k}^{\mathrm{de},\mathrm{dn}}$ are the upward and downward flexibility requirements generated by the net load at time t in scenario k, respectively.

2.2. Analysis of Flexibility Resources and Flexibility Supply Capabilities

Existing literature has focused on the power supply side's participation in the flexibility supply and demand balance, while ignoring the flexibility supply capacity on the demand side. This section considers thermal power units and pumped storage power stations as traditional flexibility resources on the power supply side, and conventional fixed loads can provide some flexibility through load transfer. As a time-coupled resource, electric vehicles are limited in flexibility by the adjustable power and capacity in adjacent time periods, and need to be modeled and analyzed differently from ordinary loads.

2.2.1. Conventional Units

Conventional units in the power system are mainly thermal power units. Although the total installed capacity of thermal power units is large, their ramp rate is low and they are limited by dispatching instructions, so their adjustment ability is poor and they can only provide limited flexibility. Flexibility supply can be expressed as:

$$F_{\mathrm{G},t}^{\mathrm{su},\mathrm{up}}=\min (R_{\mathrm{G}}^{\mathrm{up}}\tau ,P_{\mathrm{G}}^{\max }-P_{\mathrm{G},t})$$

(8)

$$F_{\mathrm{G},t}^{\mathrm{su},\mathrm{dn}}=\min (R_{\mathrm{G}}^{\mathrm{dn}}\tau ,P_{\mathrm{G},t}-P_{\mathrm{G}}^{\min })$$

(9)

where:$F_{\mathrm{G},t}^{\mathrm{su},\mathrm{up}}$and $F_{\mathrm{G},t}^{\mathrm{su},\mathrm{dn}}$are the upward and downward flexibility supplies of the thermal power unit at time t, respectively; $R_{\mathrm{G}}^{\mathrm{up}}$ and $R_{\mathrm{G}}^{\mathrm{dn}}$ are the upward and downward ramp rates of the thermal power unit, respectively; $P_{\mathrm{G}}^{\max }$ and $P_{\mathrm{G}}^{\min }$ are the maximum and minimum values of the technical output P_G,t of the thermal power unit, respectively.

2.2.2. Pumped Storage Power Station

Pumped storage power stations have flexible operation modes and good two-way flexibility supply capabilities at multiple time scales, which can share the peak demand for conventional units. Its flexibility supply is affected by the time scale and reservoir energy storage, which can be expressed as:

$$F_{\mathrm{P},t}^{\mathrm{su},\mathrm{up}}=\min \left(P_{\mathrm{P}}^{\max }-P_{\mathrm{P},t},{\displaystyle \frac{W_{\mathrm{P}}^{\max }-W_{\mathrm{P},t}}{\tau }}\right)$$

(10)

$$F_{\mathrm{P},t}^{\mathrm{su},\mathrm{dn}}=\min \left(P_{\mathrm{P},t}-P_{\mathrm{P}}^{\min },{\displaystyle \frac{W_{\mathrm{P},t}-W_{\mathrm{P}}^{\min }}{\tau }}\right)$$

(11)

where: $F_{\mathrm{P},t}^{\mathrm{su},\mathrm{up}}$and $F_{\mathrm{P},t}^{\mathrm{su},\mathrm{dn}}$ are the upward and downward flexibility supply of the pumped-storage power station at time t, respectively; P_P,t is the pumped power of the pumped-storage power station at time t; $P_{\mathrm{P}}^{\max }$ and $P_{\mathrm{P}}^{\min }$ are the maximum and minimum pumped power of the pumped-storage power station, respectively; W_P,t is the energy storage of the reservoir at time t; $W_{\mathrm{P}}^{\max }$and $W_{\mathrm{P}}^{\min }$ are the upper and lower limits of the reservoir energy storage, respectively.

2.2.3. Conventional Transferable Load

The transferable load has a certain flexibility in the grid dispatching process. Residents can appropriately change their electricity consumption plans according to incentives, providing two-way flexibility supply[7], which can be expressed as:

$$F_{\mathrm{TL},t}^{\mathrm{su},\mathrm{up}}=P_{\mathrm{TL},t}-P_{\mathrm{TL}}^{\min }$$

(12)

$$F_{ \mathrm{TL},t}^{ \mathrm{su},\mathrm{dn}}=P_{\mathrm{TL}}^{\max }-P_{\mathrm{TL},t}$$

(13)

where: $F_{\mathrm{TL},t}^{\mathrm{su},\mathrm{up}}$and $F_{\mathrm{TL},t}^{\mathrm{su},\mathrm{dn}}$ are the upward and downward flexibility supplies of the transferable load at time t, respectively; P_TL,t is the power of the transferable load at time t; $P_{\mathrm{TL}}^{\max }$ and $P_{\mathrm{TL}}^{\min }$ are the upper and lower limits of the transferable load power, respectively.

2.2.4. Cluster Electric Vehicles

Electric vehicles have the characteristics of small capacity and large scale. Only through the participation of electric vehicle aggregators as a whole in optimized scheduling can they exert flexibility. Cluster electric vehicles have a fast response time and can provide a large amount of two-way flexible supply. Its supply capacity is mainly constrained by the power of charging piles and the travel habits of car owners, which can be expressed as:

$$F_{\mathrm{EV},t}^{\mathrm{su}.\mathrm{up}}=\min \left(P_{\mathrm{EV},t}-P_{\mathrm{EV},t}^{\min },{\displaystyle \frac{E_{\mathrm{EV},t}-E_{\mathrm{EV},t}^{\min }}{\tau }}\right)$$

(14)

$$F_{\mathrm{EV},t}^{\mathrm{su},\mathrm{dn}}=\min \left(P_{\mathrm{EV},t}^{\max }-P_{\mathrm{EV},t},{\displaystyle \frac{E_{\mathrm{EV},t}^{\max }-E_{\mathrm{EV},t}}{\tau }}\right)$$

(15)

where: $F_{\mathrm{EV},t}^{\mathrm{su},\mathrm{up}}$ and $F_{\mathrm{EV},t}^{\mathrm{su},\mathrm{dn}}$ are the upward and downward flexibility supply of cluster electric vehicles at time t, respectively; P_EV,t and E_EV,t are the charging and discharging power and electric energy of cluster electric vehicles at time t, respectively; $P_{\mathrm{EV},t}^{\max }$ and $P_{\mathrm{EV},t}^{\min }$ are the upper and lower limits of the charging and discharging power of cluster electric vehicles at time t, respectively; $E_{\mathrm{EV},t}^{\max }$ and $E_{\mathrm{EV},t}^{\min }$ are the upper and lower limits of the electric energy of cluster electric vehicles at time t, respectively.

The feasible domain of charging and discharging power of EVs is not only constrained by the upper and lower limits of its own power, but also affected by electric energy. In this paper, we utilize the extreme scenario method, integrate the parameters of EVs and charging piles, as well as the travel habits of vehicle owners, and establish a general model for single EVs with a feasible domain of half-planar polyhedral sets.

Figure 2. Modeling principle of EV.

Figure 2. Modeling principle of EV.

2.3. Flexibility Balance

The flexibility balance of the power system refers to the ability of the total flexibility supply of flexibility resources to meet the flexibility demand of net load. By integrating the flexibility supply capacity of various resources, we can get:

$$\left\{\begin{array}{l}{F}_t^{\mathrm{su},\mathrm{up}}=F_{\mathrm{G},t}^{\mathrm{su},\mathrm{up}}+F_{\mathrm{P},t}^{\mathrm{su},\mathrm{up}}+F_{\mathrm{T}1.,t}^{\mathrm{su},\mathrm{up}}+F_{\mathrm{EV},t}^{\mathrm{su},\mathrm{up}}\\ F_t^{\mathrm{su},\mathrm{dn}}=F_{\mathrm{G},t}^{\mathrm{su},\mathrm{dn}}+F_{\mathrm{P},t}^{\mathrm{su},\mathrm{dn}}+F_{\mathrm{T}1.,t}^{\mathrm{su},\mathrm{dn}}+F_{\mathrm{EV},t}^{\mathrm{su},\mathrm{dn}}\end{array}\right.$$

(16)

where: $F_t^{\mathrm{su},\mathrm{up}}$ and $F_t^{\mathrm{su},\mathrm{dn}}$ are the upward and downward flexibility supply of the power grid at time t.

The flexibility margin of the power grid is:

$$\left\{\begin{array}{c}\mathsf{\Delta }{F}_t^{\mathrm{up}}=F_t^{\mathrm{su},\mathrm{up}}-F_t^{\mathrm{de},\mathrm{up}}\\ \mathsf{\Delta }{F}_t^{\mathrm{dn}}=F_t^{\mathrm{su},\mathrm{dn}}-F_t^{\mathrm{de},\mathrm{dn}}\end{array}\right.$$

(17)

where: $F_t^{\mathrm{de},\mathrm{up}}$ and$F_t^{\mathrm{de},\mathrm{dn}}$ are the upward and downward flexibility demands of the power grid at time t, respectively; Δ$F_t^{\mathrm{up}}$ and Δ$F_t^{\mathrm{dn}}$ are the upward flexibility margin and downward flexibility margin of the power grid, respectively.

When the margin is less than 0, there will be a flexibility deficit, that is:

$$F_t^{\mathrm{defc},\mathrm{up}}=\left\{\begin{array}{c}\mid \mathsf{\Delta }{F}_t^{\mathrm{up}}\mid \hspace{1em}\mathsf{\Delta }{F}_t^{\mathrm{up}}⩽0\\ 0\hspace{1em}\hspace{1em}\hspace{1em} \hspace{1em}\mathsf{\Delta }{F}_t^{\mathrm{up}}>0\end{array}\right.$$

(18)

$$F_t^{\mathrm{defc},\mathrm{dn}}=\left\{\begin{array}{c}\mid \mathsf{\Delta }{F}_t^{\mathrm{dn}}\mid \hspace{1em}\mathsf{\Delta }{F}_t^{\mathrm{dn}}⩽0\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathsf{\Delta }{F}_t^{\mathrm{dn}}>0\end{array}\right.$$

(19)

where: $F_t^{\mathrm{defc},\mathrm{up}}$ and$F_t^{\mathrm{defc},\mathrm{dn}}$ are the upward and downward flexibility deficits of the grid at time t, respectively.

After calculating the flexibility supply of various resources in the power system, the flexible adjustment factor can be defined as shown in formula (20).

$$\left\{\begin{array}{c}{\alpha }_{m,t}^{\mathrm{up}}={\displaystyle \frac{F_{m,t}^{\mathrm{su},\mathrm{up}}}{F_t^{\mathrm{de},\mathrm{up}}}}\\ {\alpha }_{m,t}^{\mathrm{dn}}={\displaystyle \frac{F_{m,t}^{\mathrm{su},\mathrm{dn}}}{F_t^{\mathrm{de},\mathrm{up}}}}\end{array}\right.\hspace{1em}m\in \{\mathrm{G},\mathrm{P},\mathrm{TL},\mathrm{EV}\}$$

(20)

where: ${\alpha }_{m,t}^{\mathrm{up}}$ and ${\alpha }_{m,t}^{\mathrm{dn}}$ are the upward and downward flexible adjustment factors of various flexibility resources, respectively, indicating the proportion of flexibility resources in meeting the grid flexibility demand at time t, reflecting the contribution of different flexibility resources to ensuring the flexibility balance of the power system.

When ${\alpha }_{\mathrm{G},t}^{\mathrm{up}}+{\alpha }_{\mathrm{P},t}^{\mathrm{up}}+{\alpha }_{\mathrm{TL},t}^{\mathrm{up}}+{\alpha }_{\mathrm{EV},t}^{\mathrm{up}}=1$, it means that the upward flexibility supply capacity of the flexibility resources can fully track the upward flexibility demand of the power system; when ${\alpha }_{\mathrm{G},t}^{\mathrm{up}}+{\alpha }_{\mathrm{P},t}^{\mathrm{up}}+{\alpha }_{\mathrm{TL},t}^{\mathrm{up}}+{\alpha }_{\mathrm{EV},t}^{\mathrm{up}}=1$, it means that the downward flexibility supply capacity of the flexibility resources can fully track the downward flexibility demand of the power system. At this time, the power system flexibility has no shortage or redundancy.

It should be noted that the flexibility supply of flexible resources includes ramp capacity and spare capacity, so the capacity spare that flexible resources can provide belongs to flexible capacity.

2.4. Model Transformation

The power balance constraints must be met at each moment:

$$P_{\mathrm{L},t}^0+P_{\mathrm{TL},t}+P_{\mathrm{EV},t}=P_{\mathrm{W},t}^0+P_{\mathrm{PV},t}^0+P_{\mathrm{P},t}+P_{\mathrm{G},t}$$

(21)

where: $P_{\mathrm{L},t}^0$ , $P_{\mathrm{W},t}^0$and $P_{\mathrm{PV},t}^0$ are the predicted values of load, wind power and photovoltaic power of the power system at time t, respectively.

2.4.1. Conventional Unit Constraints

Since this paper studies the day-ahead time scale, it is assumed that the unit is in a normally open state and the start and stop of the unit are ignored. The constraints of the thermal power unit include ramp and power constraints, which are expressed as:

$$-R_G^{\mathrm{dn}}\tau ⩽P_{G,t+\tau }-P_{G,t}⩽R_G^{\mathrm{up}}\tau$$

(22)

$$P_{\mathrm{G}}^{\min }⩽P_{\mathrm{G},t}⩽P_{\mathrm{G}}^{\max }$$

(23)

where: P_G,t is the power generation capacity of the thermal power unit at time t; $P_{\mathrm{G}}^{\max }$and$P_{\mathrm{G}}^{\min }$ are the upper and lower limits of the thermal power unit output.

2.4.2. Pumped Storage Power Station Constraints

Pumped storage hydropower stations are similar to energy storage unit models and can be divided into two states: power generation and pumping. They cannot be operated simultaneously during operation. For conventional pumped storage power stations, the power generation can be continuously changed, while the pumping power cannot be adjusted to a constant value, and the following constraints are met:

$$0⩽P_{\mathrm{P}.t}^{\mathrm{G}}⩽U_{\mathrm{P}.t}{P}_{\mathrm{P}}^{\max }$$

(24)

$$P_{\mathrm{P}.t}^{\mathrm{P}}=-(1-U_{\mathrm{P}.t})P_{\mathrm{P}}^{\min }$$

(25)

$$P_{\mathrm{P}.t}=P_{\mathrm{P}.t}^{\mathrm{G}}-P_{\mathrm{P}.t}^{\mathrm{P}}$$

(26)

where: $P_{\mathrm{P}}^{\max }$ and $P_{\mathrm{P}}^{\min }$ are the upper and lower limits of the pumped storage power station’s pumped power; U_P.t is an operating state variable, which is 1 for power generation and 0 for pumping; $P_{\mathrm{P}.t}^{\mathrm{G}}$ and $P_{\mathrm{P}.t}^{\mathrm{P}}$ are the power generation and pumping power of the pumped storage power station; P_P.t is the pumped power of the power station at time t, with power generation being a positive value and pumping being a negative value.

At the same time, the reservoir of the pumped storage power station stores water at the beginning and end of each regulation cycle to maintain the balance of the reservoir storage volume and ensure its cyclical participation in power grid regulation. The pumped power of the pumped storage power station and the reservoir energy storage satisfy the following relationship:

$${\eta }_{\mathrm{P}}\sum _{t=1}^T{P}_{\mathrm{P},t}^{\mathrm{G}}\tau -{\displaystyle \frac{1}{{\eta }_{\mathrm{P}}}}\sum _{t=1}^T{P}_{\mathrm{P},t}^{\mathrm{P}}\tau =0.$$

(27)

$$W_{\mathrm{P}}^{\min }⩽W_{\mathrm{P}}^0-{\eta }_{\mathrm{P}}\sum _{t=1}^T{P}_{\mathrm{P}A}^{\mathrm{G}}\tau +{\displaystyle \frac{1}{{\eta }_{\mathrm{P}}}}\sum _{t=1}^T{P}_{\mathrm{P},t}^{\mathrm{P}}\tau ⩽W_{\mathrm{P}}^{\max }.$$

(28)

where: ηP is the pumping efficiency of the pumped storage power station; T is the total dispatching period; W is the initial energy storage of the reservoir; W and W are the upper and lower limits of the energy storage of the power station reservoir during the regulation process.

2.4.3. Conventional Transferable Load Constraints

When conventional transferable loads perform demand response, the adjustment of the power consumption plan is affected by the power demand, and the power consumption characteristics that are satisfied can be expressed as:

$$\sum _{t=1}^T{P}_{\mathrm{TL},t}=D_{\mathrm{TL}}.$$

(29)

$$P_{\mathrm{TL}}^{\min }⩽P_{\mathrm{TL},t}⩽P_{\mathrm{TL}}^{\max }$$

(30)

where: D_TL is the total power demand of the transferable load during the period of grid dispatch; P_TL,t is the power of the transferable load at time t; $P_{\mathrm{TL}}^{\max }$ and $P_{\mathrm{TL}}^{\min }$ are the upper and lower limits of the power.

3. Based on a Data-Driven DRO Model

3.1. Deterministic Model

This paper introduces flexibility deficiency penalty to quantify flexibility and constructs a flexible dispatch model with the lowest total grid operation cost as the objective function. The objective function without considering new energy and load uncertainty is:

$$F_{\mathrm{total}}=\min (F_1+F_2)$$

(31)

where: F₁ is the operating cost; F₂ is the penalty cost for the lack of flexibility. According to the policy of giving priority to new energy access to the grid, this paper assumes that all wind and photovoltaic power are absorbed and do not participate in the optimization and dispatching of the power grid.

The grid operation costs include the coal consumption costs of thermal power units, the peak regulation costs of pumped storage, the transferable load compensation costs, and the discharge loss compensation costs of cluster electric vehicles.

$$\left\{\begin{array}{l}{F}_1=\sum _{t=1}^T(C_{\mathrm{G},t}+C_{\mathrm{P},t}+C_{\mathrm{TL},t}+C_{\mathrm{EV},t})\\ C_{\mathrm{G},t}=a{P}_{\mathrm{G},t}^2+b{P}_{\mathrm{G},t}+c\\ C_{\mathrm{P},t}=K_{\mathrm{P}}(P_{\mathrm{P},t}^{\mathrm{G}}+P_{\mathrm{P},t}^{\mathrm{P}})\\ C_{\mathrm{TL},t}=K_{\mathrm{TL}}|P_{\mathrm{TL},t}-P_{\mathrm{TL},t}^{\ast }|\\ C_{\mathrm{EV},t}=K_{\mathrm{EV}}{P}_{\mathrm{EV},t}^{\mathrm{dis}}\end{array}\right.$$

(32)

where: C_G,t C_p,t,C_TL,t and C_EV,t are the coal consumption cost of thermal power units, peak-shaving cost of pumped storage, compensation cost of transferable loads and discharge loss cost of cluster electric vehicles in the power grid at time t, respectively; T is the dispatch period; a, b and c are the cost coefficients of thermal power units; $P_{\mathrm{P},t}^{\mathrm{G}}$ and $P_{\mathrm{P},t}^{\mathrm{P}}$ are the power generation and pumping power of the pumped storage power station at time t, respectively; K_P is the unit peak-shaving cost of pumped storage; K_TL is the unit compensation cost of transferable loads participating in demand response; K_EV is the unit discharge compensation cost of cluster electric vehicles; $P_{\mathrm{EV},t}^{\mathrm{dis}}$ is the discharge power of cluster electric vehicles at time t; $P_{\mathrm{TL},t}^{*}$ is the expected electricity demand of residents at time t. By introducing auxiliary variables P_TL1,t and P_TL2,t and adding constraints as shown in Equations (33) and (34), the transferable load compensation cost is equivalent to a linear form as shown in Equation (35).

$$P_{\mathrm{TL},t}-P_{\mathrm{TL},t}^{\ast }+P_{\mathrm{TL},1,t}-P_{\mathrm{TL},2,t}=0$$

(33)

$$\left\{\begin{array}{c}{P}_{\mathrm{TL}1,t}⩾0\\ P_{\mathrm{TL}2,t}⩾0\end{array}\right.$$

(34)

$$C_{\mathrm{TL},t}=K_{\mathrm{TL}}(P_{\mathrm{TL},1,t}+P_{\mathrm{TL},2,t})$$

(35)

The coal consumption cost function of the thermal power unit is piecewise linearized to improve the model calculation speed, and the coal consumption cost of the thermal power unit becomes:

$$C_{\mathrm{G},t}=\sum _{s=1}^s{K}_s{P}_{\mathrm{G},t,s}+C_{\mathrm{G},0}$$

(36)

where: K_s is the slope of each segment after the original quadratic function is piecewise linearized; S is the number of segments; C_G,0 is the cost of coal consumption generated when the thermal power unit operates at the minimum output $P_{\mathrm{G}}^{\min }$; P_G,t,s is the segmented output of the unit, satisfying formula (38).

$$C_{\mathrm{G},0}=a{\left(P_{\mathrm{G}}^{\min }\right)}^2+b{P}_{\mathrm{G}}^{\min }+c$$

(37)

$$\left\{\begin{array}{l}0⩽P_{\mathrm{G},t,s}⩽{\displaystyle \frac{P_{\mathrm{G}}^{\max }-P_{\mathrm{G}}^{\min }}{m}}\\ P_{\mathrm{G},t}=\sum _{s=1}^s{P}_{\mathrm{G},t,s}+P_{\mathrm{G}}^{\min }\end{array}\right.$$

(38)

When the power system flexibility margin is insufficient, it will be penalized. The flexibility deficiency penalty cost is expressed as the sum of the upward flexibility deficiency penalty cost and the downward flexibility deficiency penalty cost, as shown in formula (39).

$$\left\{\begin{array}{l}{F}_2=\sum _{t=1}^T{C}_{\mathrm{defc},t}\\ C_{\mathrm{defc},t}=K_{\mathrm{defc}}(F_t^{\mathrm{defc},\mathrm{up}}+F_t^{\mathrm{defc},\mathrm{dn}})\end{array}\right.$$

(39)

where: C_defc,t is the flexibility deficit penalty cost of the power grid at time t; K_defc is the flexibility deficit penalty factor.

3.2. Two-Stage DRO Model

In order to find the optimal flexible scheduling solution under the probability distribution of the worst scenario, this section establishes a data-driven DRO model based on the historical sample data of load and new energy. The DRO involved in this article is a two-stage min-max-min three-layer optimization problem. The first stage min problem uses the output of each flexibility unit on the power supply side as the decision variable to achieve the optimal economic operation goal. The second stage max-min problem finds the probability distribution of the worst scenario that maximizes the minimum value of the demand-side resource call cost and the grid flexibility shortage penalty cost. The expression is as follows:

$$\underset{x}{\min }\left(f_1(x)+\underset{p_k\in \mathsf{\Omega }}{\max }\sum _{k=1}^K{p}_{k}r(x,{\xi }_k)\right)$$

(40)

$$\mathrm{s}.\mathrm{t}.\hspace{1em}h(x)⩽0$$

(41)

$$\left\{\begin{array}{l}r(x,{\xi }_k)=\underset{y_k\in U(x,{\xi }_k)}{\min }{f}_2(x,{\xi }_k,y_k)\\ \mathrm{s}.\mathrm{t}.\hspace{1em}g(x,{\xi }_k,y_k)⩽0\end{array}\right.$$

(42)

$$f_1(x)=\sum _{t=1}^T(C_{\mathrm{G},t}+C_{\mathrm{P},t})$$

(43)

$$f_2(y)=\sum _{t=1}^T(C_{\mathrm{TL},t}+C_{\mathrm{EV},t}+C_{\mathrm{defc},t})$$

(44)

where: x is the variable in the first stage; y is the variable in the second stage; K is the number of discrete scenarios; p_k is the probability of the occurrence of the kth scenario; r(x, ξ_k) is the objective function value of the inner minimum problem; ξ_k is the kth scenario; y_k is the variable in the second stage under the kth scenario; f₁ is the output cost of each unit on the power supply side; f₂ is the resource call cost on the demand side and the penalty cost of the grid flexibility shortage; Ω is the comprehensive norm fuzzy set; U (x, ξ_k) is the feasible domain of the optimization variable y_k when a set (x, ξ_k) is given; h (x) ≤ 0 is a constraint related only to the first stage, without considering the uncertainty of the probability of new energy and load scenarios, and is a deterministic constraint; g (x, p_k, y_k) ≤ 0 is the constraint in the second stage.

Since the uncertain probability set of load and new energy is difficult to obtain, this paper adopts a data-driven method to obtain K discrete typical scenarios (ξ₁, ξ₂, …, ξ_K) in M samples of historical data through a probability-based reduction method to characterize the possible values of uncertain information, and uses this probability distribution as the initial probability distribution. Considering that the discrete values of each reduction scenario are still uncertain, in order to ensure that the value of the scenario probability fluctuates within a reasonable range, a comprehensive norm constraint centered on the above initial probability distribution is constructed to restrict the probability distribution of the uncertain scenario, making it closer to the real scenario. The comprehensive norm constraint includes 1-norm and ∞-norm, as shown in Equation (45) and Equation (46), respectively. The constructed comprehensive norm fuzzy set is shown in Equation (47).

$${\begin{array}{c}{p}_k-p_k^0\end{array}}_1=\sum _{k=1}^K\left|\begin{array}{c}{p}_k-p_k^0\end{array}\right|⩽{\theta }_1$$

(45)

$${\begin{array}{c}{p}_k-p_k^0\end{array}}_{\infty }=\underset{1⩽k⩽K}{\max }\left|\begin{array}{c}{p}_k-p_k^0\end{array}\right|⩽{\theta }_{\infty }$$

(46)

$$\begin{array}{l}\mathsf{\Omega }=\left\{p_k|p_k⩾0,k=1,2,\cdots ,K;{\displaystyle \sum _{k=1}^K}{p}_k=1;\right.\\ \left.{\displaystyle \sum _{k=1}^K}|p_k-p_k^0|⩽{\theta }_1;\underset{1\le k\le K}{\max }|p_k-p_k^0|⩽{\theta }_{\infty }\right\}\end{array}$$

(47)

where: ‖⋅‖1 is the 1-norm; ‖⋅‖_∞ is the ∞-norm; $p_k^0$ is the probability distribution value of the initial scenario; θ_∞ and θ₁ are the deviation values allowed by the probability distribution under the ∞-norm and 1-norm constraints, respectively. It can be seen from the literature [27] that the probability distribution {p_k} satisfies the confidence constraints of equations (48) and (49).

$$\mathrm{Pr}\left\{\sum _{k=1}^K|p_k-p_k^0|⩽{\theta }_1\right\}⩾1-2K{\mathrm{e}}^{{\displaystyle \frac{-2M{\theta }_1}{K}}}$$

(48)

$$\mathrm{Pr}\left\{\underset{1⩽k⩽K}{\max }|p_k-p_k^0|⩽{\theta }_{\infty }\right\}⩾1-2K{\mathrm{e}}^{-2M{\theta }_{\infty }}$$

(49)

Let the right sides of equations (48) and (49) be the uncertainty probability confidence levels α₁ and α_∞, respectively. Under a confidence level of 95%, equations (48) and (49) guarantee that there is at least a 95% probability that a fuzzy distribution exists in a given set [27]. θ_∞ and θ₁ can be calculated by equation (50).

$$\left\{\begin{array}{c}{\theta }_1={\displaystyle \frac{K}{2M}}\ln {\displaystyle \frac{2K}{1-{\alpha }_1}}\\ {\theta }_{\infty }={\displaystyle \frac{1}{2M}}\ln {\displaystyle \frac{2K}{1-{\alpha }_{\infty }}}\end{array}\right.$$

(50)

The deviation values of θ_∞ and θ₁ indicate the maximum value by which the DRO scenario probability can deviate from the initial scenario probability. The larger the value, the more conservative the robust model is, and vice versa.

4. Two-Stage Model Solving Algorithm

Two-stage distributed robust optimization is usually based on the idea of zero-sum game, and uses the column and constraint generation (C&CG) algorithm or the cutting plane (Benders) algorithm to solve the original problem, decoupling it into a main problem and sub-problems for iteration. Compared with the Benders algorithm, the C&CG algorithm continuously introduces auxiliary variables and constraints related to the sub-problems in the process of solving the main problem, which speeds up the convergence speed.

Decomposing the two-stage model (i.e., (40) to (42)), the master problem (MP) is:

$$\left\{\begin{array}{l}\underset{x}{\min }(f_1(x)+\theta )\\ \mathrm{s}.\mathrm{t}.\hspace{1em}h(x)⩽0\\ \theta ⩾\sum _{k=1}^K{p}_k^{(l)}{f}_2(x,{\xi }_k,y_k^{(l)})\\ g(x,{\xi }_k,y_k^{(l)})⩽0\end{array}\right..$$

(51)

where: l is the number of iterations; θ is the relaxation amount of the subproblem; y is the auxiliary variable related to the subproblem introduced at the l-th iteration; p is the worst scenario probability of wind and photovoltaic load obtained by solving the subproblem at the l-th iteration, which is a known quantity.

The purpose of the sub-problem is to find the worst scenario probability distribution that maximizes the minimum value of the demand-side resource call cost and the grid flexibility deficiency penalty cost. When the variables in the first stage are given, the sub-problem can be expressed as:

$$\left\{\begin{array}{l}\underset{p_k\in \mathsf{\Omega }}{\max }\left(\sum _{k=1}^K{p}_k\underset{y_k\in U(x,{\xi }_k)}{\min }{f}_2(x^{\ast },{\xi }_k,y_k)\right)\\ \mathrm{s}.\mathrm{t}.\hspace{1em}g(x^{\ast },{\xi }_k,y_k)⩽0\end{array}\right..$$

(52)

where: x* is the solution to the main problem, and is a known quantity here.

Although the form of the subproblem is an NP-hard max-min double-layer optimization problem, it is not difficult to find that the inner min problems in each scenario are independent of each other, and the inner min problems can be solved simultaneously using parallel computing methods. Therefore, the subproblem can be rewritten as two single-layer optimization problems that are solved in sequence, as shown in formula (53).

$$\left\{\begin{array}{l}f(x^{\ast },{\xi }_k)=\underset{y_k\in U(x,{\xi }_k)}{\min }{f}_2(x^{\ast },{\xi }_k,y_k)\\ \underset{p_k\in \mathsf{\Omega }}{\max }\sum _{k=1}^K{p}_{k}f(x^{\ast },{\xi }_k)\end{array}\right..$$

(53)

Equation (53) is a mixed integer programming, which can be quickly solved using the Gurobi solver. The result {p_k} obtained is the probability distribution of the worst scenario, which is passed to the main problem for iterative calculation.

The specific solution process is as follows:

• Start
• Collect historical wind-photovoltaic-load data
  Gather historical sample data for wind, photovoltaic, and load.
• Construct copula temporal and spatial correlation output set
  Use the collected data to build the Copula temporal and spatial correlation output set.
• Initialize iteration parameters
  Set the lower bound L_B of the second-stage cost to -∞ and the upper bound U_B to +∞.
  Set the iteration counter l to 1.
• Main problem optimization
  Solve the main optimization problem considering the worst-case scenario probability distribution.
  Update the upper bound U_B and lower bound L_B based on the solution.
• Check iteration termination conditions
  If the termination conditions (e.g., convergence of U_B and L_B, maximum number of iterations reached) are met, proceed to the next step.
  Otherwise, go back to the Main Problem Optimization step.
• Update worst-case scenario probability distribution
  Solve the subproblem to update the probability distribution of the worst-case scenario.
• Obtain optimal solution
  After the iteration terminates, obtain the optimal solution for the wind-photovoltaic-load management problem.
• End

5. Results and Discussion

The algorithms were solved through MATLAB 2018b platform and Gurobi solver. The time scale for flexible dispatch before the day in this paper is taken as 1 h. The wind and light load data samples are randomly generated using normal distribution on the basis of the predicted values, and the number of samples is taken as 500, and the probability distribution of the initial scenario is obtained by probabilistic distance-based scenario cutting technique, and the number of scenarios after cutting is 5, and the typical scenario curves are shown in the literature [19]. The maximum fluctuation coefficients of ε_L, ε_W and ε_PV for the load, wind power and PV in the day before yesterday are 0.05, 0.10 and 0.08, respectively, and the parameters of each unit of the calculation example are shown in Table 1.

5.1. Power System Flexibility Balancing Results

In this section, four flexibility regulation scenarios are set up to analyze the impact of the participation of load-side resources in the flexibility supply and the inclusion of flexibility balancing constraints on the flexibility margin of the power system during the day-ahead scheduling phase.

Scheme 1: Load-side resources do not participate in demand response (transferable loads are considered as fixed loads and Evs are charged in a disorderly manner), and flexibility balancing constraints are not considered.

Scheme 2: Load-side resources participate in demand response (orderly charging and discharging of Evs) without considering flexibility balancing constraints.

Scheme 3: Load-side resources participate in demand response, taking into account flexibility balancing constraints.

Scheme 4: Load-side resources do not participate in demand response, considering flexibility balancing constraints.

Since Scheme 1 and 2 do not consider flexibility constraints, the flexibility demand and flexibility supply of the power system are then calculated at the end of the optimal scheduling. The Scheme 3 grid flexibility margins are shown in Figure 3. The net load curve of Scenario 3 is flatter in the 16:00-22:00 hours, i.e., the uphill and downhill sections of the evening peak, reducing the flexibility demand of the grid, so that there are sufficient upward and downward flexibility margins in all hours except 19:00.

Table 2 shows the operating costs and flexibility balances for each Scheme. Scheme 1 does not consider demand response and does not optimize the flexibility supply and demand, so there is a large flexibility deficit in the grid. Schemes 2 and 3 take into account the transferability of loads, and although they require some call subsidy costs, the total operating cost and flexibility deficit are lower than in Scenario 1. Scenario 3 takes into account the flexibility balancing constraints in comparison to Scheme 2, and requires more demand-side resources to be called upon, so the total operating cost is slightly higher, but the flexibility deficit is significantly reduced, which ensures that the system is flexible with a small reduction in the economic performance. The flexibility of Scheme 4 is only supplied by the power side, although the flexibility balancing constraint is considered on the basis of power balancing, it only works on the power side resources, and the adjustment range Is not large, so there Is still a large amount of flexibility shortfall, and the operating cost Is close to that of Scheme 1.

5.2. Analysis of Flexibility Adjustment Capacity for Flexibility Resources

In this paper, a flexibility adjustment factor is introduced in the participation of various types of flexibility resources in maintaining the flexibility balance of the grid to reflect the degree of contribution of different flexibility resources to ensure the flexibility balance of the grid. The larger the flexibility adjustment factor of a flexibility resource at time t is, the stronger the flexibility adjustment capability of the resource is. When the sum of the flexibility adjustment factors of all types of resources is less than 1, it means that the power system has a flexibility shortage at time t. Conversely, when the sum of the flexibility adjustment factors is greater than 1, it indicates that flexibility is left in reserve.

Figure 4 shows the contribution of different flexibility resources to take up the flexibility needs of the grid in scheme 3, where the red dashed line indicates that the flexibility resources fully track the flexibility needs of the grid.

In Figure 4(a), the sum of the flexibility regulation factors at the 19:00 hour is below the red line, which means that there is a shortage of flexibility at that hour, which is consistent with the analysis results in Section 4.1. From Fig. 4(a) and (b), it can be seen that the flexibility regulation ability of each type of flexibility resource is different in each time period. The flexibility of thermal units is limited by capacity and ramp-up, and can only withstand baseload, so the downward flexibility supply capability is better than the upward flexibility supply capability as a whole. During 10:00-12:00 load trough hours, thermal units operate at lower output levels and can only provide a small amount of downward flexibility supply. During the 13:00-20:00 net load rise period, thermal power units are only supplied with ramping capacity for upward flexibility. Pumped storage power plants pump water and generate electricity at sharp net load fluctuations, with strong flexibility, effectively compensating for the lack of climbing capacity of thermal power units in a short period.

Cluster EVs should take into account the travel needs of their owners when participating in vehicle-network interactions: choosing to charge at maximum power during the 12:00-16:00 net load trough hour provides a lot of upward flexibility; discharging during the 19:00-21:00 evening peak hour provides only a small amount of downward flexibility. It should be noted that EVs are also discharged near the 10:00 hour due to the risk of insufficient downward flexibility of the grid at that hour, and the charging and discharging behavior of cluster EVs is adequately managed due to the flexibility balance constraints considered in this paper. The transferable loads are limited by electricity demand and capacity, and can provide less flexibility by transferring peak loads, whose downward flexibility supply is slightly larger than the upward flexibility supply.

5.3. Conservativeness Analysis of Data-Driven DRO Models

Given uncertainty probability confidence levels α₁ and α_∞ of 0.5 and 0.99, respectively, the effect of different historical sample sizes M and the number of typical scenarios cut K on the degree of conservatism of the DRO model is investigated.

As can be seen in Table 3, as the number of historical samples increases, the cost of operation and flexibility deficit of each unit of the grid decreases. This is due to the fact that the increase in the number of historical samples makes the deviation values θ_∞ and θ₁ allowed by the probability distribution decrease accordingly, and the conservatism of the DRO model is alleviated.

Table 4 shows the effect of the number of typical scenarios on the model conservatism after sample data reduction, and it is easy to find that the conservatism degree increases with the number of scenarios. This is because as the number of typical scenarios becomes larger, the resolution of the cut scenarios becomes larger, and the probability that the extreme data in the original large amount of sample data become a class of typical scenarios alone becomes higher. Typical scenarios are portrayed more and more finely, closer to the actual fluctuation range of wind and photovoltaic loads, and the robustness of the model is improved.

The conservatism of the DRO model in this paper is reflected in the fact that the probability of occurrence of severe scenarios of net load fluctuations becomes higher. Both the static equilibrium of power supply and demand are reflected and the fluctuating process of power supply and demand, i.e., the flexibility equilibrium process, is highlighted. The small change in the cost of thermal units when the model conservatism increases is related to their tracking of base load. The increased total operating costs mainly come from the demand response compensation costs of load-side resources, which results in a smaller increase in the flexibility deficit of the power system, proving the validity of the proposed model.

5.4. Comparative Analysis with Other Uncertainty Optimization Methods

The data-driven DRO model (the model in this paper), which takes into account the correlation of wind and scenery output, is compared with the data-driven DRO model without considering the correlation of wind and scenery output, the traditional two-stage RO model, and the SO model. The uncertainty probability confidence level α₁ and α_∞ of the data-driven DRO model are 0.5 and 0.99, respectively, and the number of samples M and the number of cut scenarios K are taken as 500 and 5, respectively. The results of the four models are shown in Table 5.

As can be seen from Table 5, the four models have a small difference in the operating cost of each unit on the power side and a large difference in the demand response cost. the RO model is the most robust considering the worst peak climbing scenario with uncertain information, so it has the largest cost and flexibility shortfall in all aspects. The data-driven DRO model synthesizes the advantages of the RO model and the SO model, and uses comprehensive paradigm constraints to construct a probability distribution fuzzy set, which has better economy and flexibility, and the economy and robustness in the day-ahead scheduling phase are taken as a compromise. When the temporal and spatial correlation of wind and photovoltaic output is described by Copula theory, the calculated grid flexibility demand is smaller due to the consideration of wind and photovoltaic complementarity, and the model proposed in this paper is less conservative compared to the data-driven DRO model that does not consider the wind and photovoltaic correlation. At the same time, since the data-driven DRO model does not require pairwise processing in the subproblem processing compared to the RO model, the parallel computation method it adopts reduces the computation time significantly.

6. Conclusions

This paper proposes a data-driven distributed robust optimization (DRO) model that considers flexibility balance. By comparing simulation examples, the following conclusions are obtained:

1) By combining the interval method with the scenario-based approach within the DRO framework, we quantify flexibility demand while fully accounting for net load fluctuations across various scenarios. Under the constraint of flexibility supply and demand balance, load-side resources such as electric vehicles and conventional transferable loads participate in demand response, providing significant ramping capacity. This effectively reduces net load fluctuations. By considering the flexibility supply capabilities of demand-side resources, the overall flexibility of the power grid is improved, particularly during morning and evening peak hours with high ramping demands. Flexibility margins are notably enhanced, ensuring power system flexibility while only marginally compromising economic performance.

2) The data-driven DRO model thoroughly considers the uncertainty of renewable energy and load fluctuations based on historical data samples. When the degree of conservatism increases, the power grid needs to offer more economic incentives to load-side resources to encourage their participation in demand response. These incentives motivate load-side resources to provide additional flexibility, addressing the challenges posed by adverse net load ramping scenarios. By offering economic compensation, the grid ensures that load-side resources are willing to contribute to flexibility provision, thereby enhancing the overall resilience of the power system against uncertainty.

3) The DRO model not only achieves a compromise between economy and robustness in the day-ahead scheduling stage but also demonstrates superior resilience against uncertainty fluctuations in real-time operations. By leveraging historical data and considering a distribution of possible outcomes, the DRO model is able to better manage uncertainty and risk. Compared to Stochastic Optimization (SO) models, the DRO model exhibits better overall economic performance. It provides a more robust solution that can withstand uncertainty while maintaining economic efficiency, making it a valuable tool for power system flexibility management in the face of increasing renewable energy penetration and load variability.