Probabilistic Power and Energy Balance Risk Scheduling Method Based on Distributed Robust Optimization

0. Introduction

The advancement of new energy technologies is a crucial pathway to realizing the objectives of "carbon peak and carbon neutrality", addressing environmental pollution, and resolving energy crisis issues. The rapid progress in new energy technology has led to a significant increase in the integration of new energy into the power grid, marking a distinctive characteristic of the modern power system. By the end of 2023, China's installed capacity of wind power and photovoltaics reached 440 million kW and 610 million kW respectively. The proportion of new energy installed capacity has reached 36%, making it the second largest power source in China [1].

When analyzing the traditional power system dominated by thermal power generation, the larger controllable range of thermal power generation often leads to an oversight of the temporal coupling relationship of operating constraints, and a segmented description of unit outage limits. Furthermore, with relatively small variance in load forecast errors, deterministic reserve rules effectively mitigate power imbalance risk to a low level [2]. The structural and operational characteristics of high-proportion renewable energy power systems exhibit substantial differences from traditional systems. Factors such as generation, demand, and storage collectively contribute to achieving power balance, and the temporal coupling characteristics result in interdependencies among viable solutions for power equilibrium at various time intervals [3,4]. The significant random fluctuations in wind and solar power have emerged as the primary source of electricity within the system, characterized by intricate and variable operating conditions that defy representation through a limited set of typical scenarios. There is an urgent imperative for probabilistic analysis of power quantity balance. Conventional time-series simulation-based methods rely on specific operational points to analyze power quantity balance, yet these operational points are uncertain under ambiguous conditions, leading to analytical challenges. Mere augmentation of simulated scenarios or introduction of segmented probability balance analysis does not effectively resolve the aforementioned issue. Consequently, traditional methods for analyzing power quantity balance are increasingly inadequate in addressing the need for temporal coupling and probability balance analysis in contexts where wind and solar energy penetration rates are escalating, necessitating a reevaluation of the matter concerning power quantity balance in power systems.

The analysis of power quantity balance is a fundamental cornerstone for designing power distribution and formulating operational strategies, serving as the central and primary concern in establishing new power systems. Presently, research on deterministic power quantity balance has achieved a significant level of maturity, with key areas of emphasis including optimized allocation of regulating resources, balanced analysis across multiple time scales within the system, and rapid solution algorithms for production simulation [5,6]. Document [7] established the groundwork for power and energy balance analysis calculations, introducing a table-based analysis method based on typical days. However, this approach is only suitable for early power systems dominated by thermal power plants. Additionally, Document [8] addressed the power and energy balance of systems with complex power supply structures, considering internal section limitations, external AC/DC power receiving operating conditions, and different resource attributes in each division to achieve inter-regional power and energy balance. With the increasing integration of fluctuating renewable energy sources, scholars are prompted to consider the evolving power and energy balance characteristics of systems with high proportions of renewable energy. Document [9] examines the variability of wind and photovoltaic energy generation across multiple time scales, as well as the complementary characteristics of new energy generation among different provinces, based on actual data from regional and provincial grids. Additionally, it assesses the absorption and supply risks faced by the power system in a high proportion of new energy sources scenario using future grid operation data. Document [10] proposes a power balance mechanism suitable for China's Northwest region that accounts for the uncertainty of new energy sources. To capture the temporal characteristics of new energy generation, time-series production simulation technology is extensively employed to analyze supply-demand balance in high-proportion new energy systems. Document [11] introduces a generator scheduling optimization algorithm that considers fluctuations in wind power generation output by pre-arranging starting modes of generators based on deterministic wind power generation scenarios in the main problem, and adjusting rescheduling plans of thermal power generation units in sub-problems to accommodate the intermittent nature of wind power generation. The final result is achieved through a hierarchical iterative solution method. Document [12] presents a fully parallelized stochastic generator scheduling method that decomposes the stochastic optimization problem into modules for generator configuration, optimal power flow calculation, and connection module; it calculates operation modes of each generator within the system under multiple wind power generation scenarios simultaneously, effectively enhancing simulation speed for large-scale transmission, distribution, transformation networks as well as direct current/alternating current hybrid transmission and distribution networks. There are several commonly used production simulation calculation software available both domestically and internationally, including the power system planning and decision-making software developed by Tsinghua University (referred to as grid optimization planning tools, GOPT) [13], and the new energy evaluation and optimization system (NEOS) [14] developed by State Grid Energy Research Institute Co., Ltd. Nevertheless, existing methods and software for analyzing power generation balance do not fully account for the significant uncertainty associated with high levels of new energy output. The decision-making units' operational modes fail to encompass the output characteristics during extreme weather scenarios, resulting in a heightened risk of compromised power supply reliability and curtailment of renewable energy within the power system.

Stochastic optimization and robust optimization are effective methods for dealing with the high uncertainty of renewable energy output, which can help achieve multi-source coordinated optimization. Among them, stochastic optimization methods aim to minimize the expected value of system operating costs by constructing models that satisfy specific probability level constraints. Document [15] considered the randomness of wind power output and optimized the unit combination based on the reliability index of wind power; Document [16] proposed a two-layer stochastic optimization model that takes into account the potential risks brought by the uncertainty of new energy output and determines the optimal standby capacity for participating in reserve auxiliary services; Document [17] applied stochastic optimization methods to the power market and achieved economic dispatch of traditional fossil-fuel power plants. The robust optimization method, when addressing the uncertainty of renewable energy output, does not necessitate explicit probability distribution of uncertain variables. Instead, it converts the uncertainty optimization problem into a deterministic optimization problem under extreme scenarios. This approach is currently extensively utilized in integrated energy systems [18], unit combinations [19], and optimal power flow [20] optimization models. To tackle the challenges associated with obtaining probability distribution information for wind and solar power output errors and excessively conservative results from robust optimization, experts have proposed the distributionally robust optimization (DRO) method. DRO establishes fuzzy sets for all potential distributions of uncertain variables based on historical data, effectively managing the uncertainty of renewable energy output. Document [21] considers the intra-provincial/inter-provincial coupled scheduling operation mode and establishes an uncertainty set of source-load power matrix based on distributionally robust optimization method, achieving the cross-province and cross-region integration of clean energy. Document [22] fully exploits the flexible regulation capability of the power system, establishing a two-stage distributionally robust optimization model to achieve economic dispatch of flexible ramping reserves. However, these models do not quantitatively evaluate the power balance risk in new-type power systems or provide probabilistic risk assessment indicators that adapt to high proportions of new energy access while considering strong uncertainty characteristics. Document [23], based on a time-varying probability model for clustered wind turbine virtual units, proposes a quantitative method for risk indicators considering the scheduling characteristics of clustered wind turbine virtual units; however, it does not consider the adjustment capability of fossil fuel-based flexible ramping reserves in response to risks. In summary, existing studies on uncertainty balance often characterize balanced risks based on deterministic probability distributions or intervals without providing probabilistic indicators with clear physical meanings for power balance risks and without considering flexibility resources' ability to respond to risks within the system, leading to biased risk assessment results.

This paper presents a method for probabilistic power balance risk analysis based on distributionally robust optimization. Initially, it integrates the affine factor and the flexible ramping reserve capacity of thermal power to define probabilistic metrics that characterize power balance risk. By leveraging conditional value-at-risk theory, it introduces risk metrics for shedding load and curtailed energy at a specified confidence level. Subsequently, employing the distributionally robust approach to measure Wasserstein distance between different probability distributions, fuzzy sets are constructed for the probability distributions of uncertain variables. Finally, in order to minimize thermal generation operating costs, curtailment risk costs, and shedding load risk costs, this study establishes a distributionally robust optimization scheduling model based on flexible ramping reserves of thermal power.

1. Assessment of Probabilistic Power and Energy Balance Characteristics and Risk Indicators for Supply-Demand Equilibrium

1.1. Analysis of Probability Balance Characteristics

In conventional power systems, achieving power and energy equilibrium relies on 'load balancing source,' where regulating load fluctuations becomes pivotal for maintaining system stability. However, within high-penetration renewable energy frameworks, wind and photovoltaic generation display significant stochastic traits that limit their regulatory capacity. Instead, numerous adaptable loads contribute to supply-demand equilibrium as an emerging form of regulatory asset—ushering in an innovative paradigm characterized by 'interconnected equilibrium across all facets of the power network.' Consequently, analyzing power-energy equilibrium has evolved into assessing adjustable resources alongside regulatory needs to better reflect operational realities. Table 1 illustrates key adjustable resources and regulatory requisites within high-penetration systems.

The probabilistic balance characteristics of power and energy arise from the stochastic variations in non-adjustable loads, wind/solar generation, and other sources. In comparison, the uncertainty associated with wind/solar generation surpasses that of load uncertainty, leading to significant disparities in prediction error levels across different time periods. Fixed reserve capacities struggle to ensure reliable balance across all time periods due to these differences, making the analysis of probabilistic balance characteristics a focal point of concern.

In a significant portion of renewable energy power system planning for a specific region, the determination of system standby capacity is based on the 10% peak load. The probability of standby deficiency for each time slot was calculated using the error distribution of net load forecast for a particular day. The error distribution was color-coded to represent the range of standby deficiency probability, as illustrated in Figure 1. It is evident from the figure that due to substantial variations in net load error distribution across different time slots, there are also considerable fluctuations in standby deficiency probability, with values reaching up to 0.3. Therefore, it is imperative to incorporate probabilistic characteristic analysis into power and energy balance assessments to effectively capture the impact of uncertainty on the system and address reliability and economic challenges stemming from standby deficiencies.

It is evident that the structure and operational characteristics of high-penetration renewable energy power systems differ significantly from traditional systems. Various adjustable resources in the sources, loads, and storage sectors collectively participate in achieving power and energy balance. The temporal coupling characteristics lead to mutual influences on the existence of feasible solutions for power and energy balance across different time periods. Wind/solar power sources with strong stochastic fluctuations have become the primary source of electricity generation for these systems. Their complex and variable operating conditions are challenging to represent using a limited number of typical scenarios, emphasizing the urgent need for probabilistic characteristic analysis in power and energy balance assessments.

1.2. Indicators of Supply and Demand Equilibrium Risk

According to the analysis of probability balance characteristics, it can be inferred that the high level of uncertainty in predicting errors for new energy is the primary factor contributing to power and energy balance risks in a significant portion of new energy power systems. The combined error in wind and solar predictions is represented as follows,

$${\tilde{\psi }}_t={\tilde{\omega }}_t^{\mathrm{wind}}+{\tilde{\omega }}_t^{\mathrm{pv}}$$

(1)

In the equation, ${\tilde{\omega }}_t^{\mathrm{wind}}$/${\tilde{\omega }}_t^{\mathrm{pv}}$ denotes the forecast errors of the wind farm and photovoltaic power station at time t.

The affine strategy is utilized as a method of adjustment for thermal power units to mitigate the forecast errors associated with wind and solar power. The actual output of the thermal power units is represented by an affine function, incorporating the forecast errors of wind and solar power as variables. The specific value of this affine function can only be computed once the actual outputs of wind and solar energy are determined.

$$\begin{array}{l}{\tilde{P}}_{i,t}^{\mathrm{G}}=P_{i,t}^{\mathrm{G}}+{\alpha }_{i,t}^{\mathrm{G}}{\tilde{\psi }}_t\\ =P_{i,t}^{\mathrm{G}}+{\overline{\alpha }}_{i,t}^{\mathrm{G}}\max \{{\tilde{\psi }}_t,0\}+{\underset{\_}{\alpha }}_{i,t}^{\mathrm{G}}\min \{{\tilde{\psi }}_t,0\}\end{array}$$

(2)

$$P_{i,t}^{\mathrm{G}}-R_{i,t}^{\mathrm{down}}\le {\tilde{P}}_{i,t}^{\mathrm{G}}\le P_{i,t}^{\mathrm{G}}+R_{i,t}^{\mathrm{up}}$$

(3)

In the equation, ${\alpha }_{i,t}^{\mathrm{G}}$ represents the participation factor of the thermal power unit; ${\overline{\alpha }}_{i,t}^{\mathrm{G}}$/${\underset{\_}{\alpha }}_{i,t}^{\mathrm{G}}$ represent the upward/downward adjustment factors of thermal power unit i at time t; $R_{i,t}^{\mathrm{up}}$、$R_{i,t}^{\mathrm{down}}$ represent the upward and downward spinning reserve capacities of thermal power unit i at time t.

Combined with the affine factor and the total forecast error of wind and solar power, ${\alpha }_{i,t}^{\mathrm{G}}{\tilde{\psi }}_t$ and ${\underset{\_}{\alpha }}_{i,t}^{\mathrm{G}}{\tilde{\psi }}_t$ respectively represent the system's upward and downward regulation requirements. When the thermal power's upward ramping reserve regulation capacity is insufficient to cover the fluctuation adjustment demand of forecasted wind and solar outputs, there exists a risk of load shedding in the system. Similarly, when the thermal power's downward ramping reserve regulation capacity is inadequate to cover the fluctuation adjustment demand of forecasted wind and solar outputs, there exists a risk of curtailment as illustrated in Figure 2.

If the forecast error of wind and solar power is a stochastic variable, then the expression for the system's probabilistic electricity balance risk is as follows.

$$\left\{\begin{array}{l}\underset{P_t^{\psi }\in {\Re }_t^{\psi }}{\inf }{P}_t^{\psi }[{\overline{\alpha }}_{i,t}^{\mathrm{G}}{\tilde{\psi }}_t\le R_{i,t}^{\mathrm{up}}]\ge 1-{\gamma }_1={\beta }_1\\ \underset{P_t^{\psi }\in {\Re }_t^{\psi }}{\inf }{P}_t^{\psi }[-R_{i,t}^{\mathrm{down}}\le {\underset{\_}{\alpha }}_{i,t}^{\mathrm{G}}{\tilde{\psi }}_t]\ge 1-{\gamma }_2={\beta }_2\end{array}\right.$$

(4)

In equation (4), ${\gamma }_1$/${\gamma }_2$ represents the allowable loss-of-load/curtailment probability of the system. Equation (4) describes that, even when the wind and solar prediction error distribution $P_t^{\psi }$ is in the worst distribution of ${\Re }_t^{\psi }$, there is a probability less than ${\gamma }_1$/${\gamma }_2$ that the constraint cannot be satisfied, as shown in Figure 2.

According to the expression for probabilistic electricity balance risk in equation (4), random variables for load shedding and curtailed power can be defined as follows,

$${\tilde{P}}_t^{\mathrm{cut}}={\displaystyle \sum _{i=1}^G{({\overline{\alpha }}_{it}^{\mathrm{G}}{\tilde{\psi }}_t-R_{it}^{\mathrm{up}})}^{+}}$$

(5)

$${\tilde{S}}_t^{\mathrm{re}}={\displaystyle \sum _{i=1}^G{(-R_{it}^{\mathrm{down}}-{\underset{\_}{\alpha }}_{it}^{\mathrm{G}}{\tilde{\psi }}_t)}^{+}}$$

(6)

In the equation, ${\tilde{P}}_t^{\mathrm{cut}}$ represents the random variable for load shedding, and ${\tilde{S}}_t^{\mathrm{re}}$ represents the random variable for curtailed power.

In accordance with the Value at Risk (VaR) definition, the threshold expressions for the random variables of load shedding and curtailed power can be delineated at a specific confidence level.

$${\theta }_t^{\mathrm{cut}}=\min [b_1:\mathbb{P}\{{\tilde{P}}_t^{\mathrm{cut}}\le b_1\}\ge 1-{\gamma }_1={\beta }_1]$$

(7)

$${\theta }_t^{\mathrm{S}}=\min [b_2:\mathbb{P}\{{\tilde{S}}_t^{\mathrm{re}}\le b_2\}\ge 1-{\gamma }_2={\beta }_2]$$

(8)

The variables ${\theta }_t^{\mathrm{cut}}$ and ${\theta }_t^{\mathrm{S}}$ denote the threshold values of the random variables corresponding to load shedding and curtailed power at time t, respectively.

The physical interpretation of ${\theta }_t^{\mathrm{cut}}$ is the minimum threshold for load shedding power at a confidence level of β1, indicating the potential imbalance between supply and demand when the dispatch instruction for thermal power reserve exceeds ${\displaystyle \sum _{i=1}^G{R}_{it}^{\mathrm{up}}}$ with a probability greater than β₁ at time t. The physical interpretation of ${\theta }_t^{\mathrm{cut}}$ aligns closely with this concept.

In the optimization strategy for balancing probabilistic power supply and demand in the power system, the objective is to identify the optimal balance point for dispatching thermal power standby capacity. This takes into consideration both safety and economic efficiency, aiming to minimize the likelihood of significant load shedding or curtailment during actual operation while ensuring effective implementation of dispatch instructions. The conditional value-at-risk (CVaR) theory is utilized to quantify the risk associated with large load shedding or curtailment caused by prediction errors in the probability distribution of new energy sources on both upper and lower tails. This study quantifies risk indicators for supply and demand balance in a high proportion new energy power system based on CVaR, expressing risk indicators for load shedding power and curtailment power at confidence levels β1 and β2 as average VaR values,

$${\phi }_t^{\mathrm{cut}}={\displaystyle \frac{1}{1-{\beta }_1}}{\displaystyle {\int }_{{\tilde{P}}_t^{\mathrm{cut}}\ge {\theta }_t^{\mathrm{cut}}}{\tilde{P}}_t^{\mathrm{cut}}p({\tilde{\psi }}_t)d{\tilde{\psi }}_t}$$

(9)

$${\phi }_t^{\mathrm{S}}={\displaystyle \frac{1}{1-{\beta }_2}}{\displaystyle {\int }_{{\tilde{S}}_t^{\mathrm{re}}\ge {\theta }_t^S}{\tilde{S}}_t^{\mathrm{re}}p({\tilde{\psi }}_t)d{\tilde{\psi }}_t}$$

(10)

The variables ${\phi }_t^{\mathrm{cut}}$ and ${\phi }_t^{\mathrm{S}}$ represent the risk indicators for load shedding power and curtailed wind power, respectively, at time t. $p({\tilde{\psi }}_t)$ denotes the probability distribution function of wind and solar forecasting errors.

${\phi }_t^{\mathrm{cut}}$ can be interpreted physically as the conditional expectation value of the load shedding power exceeding its threshold ${\theta }_t^{\mathrm{cut}}$, and a similar physical interpretation applies to ${\phi }_t^{\mathrm{S}}$.

In fact, by minimizing ${\phi }_t^{\mathrm{cut}}$/${\phi }_t^{\mathrm{S}}$, a small ${\theta }_t^{\mathrm{cut}}$/${\theta }_t^{\mathrm{S}}$ is also obtained due to the constraint that ${\theta }_t^{\mathrm{cut}}$/${\theta }_t^{\mathrm{S}}$ does not exceed ${\phi }_t^{\mathrm{cut}}$/${\phi }_t^{\mathrm{S}}$. An auxiliary function $F_{\beta ,t}^{\mathrm{cut}}$/$F_{\beta ,t}^{\mathrm{S}}$ is introduced to represent ${\phi }_t^{\mathrm{cut}}$/${\phi }_t^{\mathrm{S}}$, which is defined as follows,

$$F_{\beta ,t}^{\mathrm{cut}}=b_1+{\displaystyle \frac{1}{1-{\beta }_1}}{\displaystyle {\int }_{{\tilde{\psi }}_t\in \Re }{[{\tilde{P}}_t^{\mathrm{cut}}-b_1]}^{+}p({\tilde{\psi }}_t)d{\tilde{\psi }}_t}$$

(11)

$$F_{\beta ,t}^{\mathrm{S}}=b_2+{\displaystyle \frac{1}{1-{\beta }_2}}{\displaystyle {\int }_{{\tilde{\psi }}_t\in \Re }{[{\tilde{S}}_t^{\mathrm{re}}-b_2]}^{+}p({\tilde{\psi }}_t)d{\tilde{\psi }}_t}$$

(12)

In this formula, ${[\cdot ]}^{+}$ represents $\max \{\cdot ,0\}$. $F_{\beta ,t}^{\mathrm{cut}}$/$F_{\beta ,t}^{\mathrm{S}}$ can be estimated using sample data of random variables ${\tilde{\psi }}_t^j$ (j=1,2,...M) as follows,

$${\tilde{F}}_{\beta ,t}^{\mathrm{cut}}=b_1+{\displaystyle \frac{1}{M(1-{\beta }_1)}}{\displaystyle \sum _{j=1}^M{[{\tilde{P}}_t^{\mathrm{cut}}({\tilde{\psi }}_t^j)-b_1]}^{+}}$$

(13)

$${\tilde{F}}_{\beta ,t}^{\mathrm{S}}=b_2+{\displaystyle \frac{1}{M(1-{\beta }_2)}}{\displaystyle \sum _{j=1}^M{[{\tilde{S}}_t^{\mathrm{re}}({\tilde{\psi }}_t^j)-b_2]}^{+}}$$

(14)

Furthermore, let,

$$H_{\beta ,t}^{\mathrm{cut}}=\max \{b_1,{\displaystyle \frac{1}{1-{\beta }_1}}{\tilde{P}}_t^{\mathrm{cut}}-{\displaystyle \frac{{\beta }_1{b}_1}{1-{\beta }_1}}\}$$

(15)

$$H_{\beta ,t}^{\mathrm{S}}=\max \{b_2,{\displaystyle \frac{1}{1-{\beta }_2}}{\tilde{S}}_t^{\mathrm{re}}-{\displaystyle \frac{{\beta }_2{b}_2}{1-{\beta }_2}}\}$$

(16)

Moreover, we can acquire,

$$F_{\beta ,t}^{\mathrm{cut}}=b_1+{\displaystyle \frac{1}{1-{\beta }_1}}{E}_{{\widehat{P}}_M}({[{\tilde{P}}_t^{\mathrm{cut}}-b_1]}^{+})=E_{{\widehat{P}}_M}(H_{\beta ,t}^{\mathrm{cut}})$$

(17)

$$F_{\beta ,t}^{\mathrm{S}}=b_2+{\displaystyle \frac{1}{1-{\beta }_2}}{E}_{{\widehat{P}}_M}({[{\tilde{S}}_t^{\mathrm{re}}-b_2]}^{+})=E_{{\widehat{P}}_M}(H_{\beta ,t}^{\mathrm{S}})$$

(18)

2. Probabilistic Balancing Analysis Model for New Energy Power Systems Based on Distributed Robust Optimization

2.1. A New Energy Power System Balance Analysis Model Based on Risk Assessment Indicators

(1)

The objective function

The objective function aims to minimize the cost of thermal power generation, standby capacity expenses, and the potential loss costs under extreme distribution conditions for system load shedding and wind curtailment.

$$\min ({\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i=1}^G(F_{i,t}^{\mathrm{G}}+F_{i,t}^{\mathrm{R}})}}+\underset{\mathbb{P}\in \mathbb{B}\epsilon ({\widehat{\mathbb{P}}}_M)}{\sup }{E}_{\mathbb{P}}(CL+CE))$$

(19)

Among them,

$$\left\{\begin{array}{l}{F}_{i,t}^{\mathrm{G}}=a_i{(P_{i,t}^g)}^2+b_i{P}_{i,t}^g+c_i\\ F_{i,t}^{\mathrm{R}}={\rho }_i^{\mathrm{G}}{R}_{i,t}^{\mathrm{up}}+{\rho }_i^{\mathrm{G}}{R}_{i,t}^{\mathrm{down}}\end{array}\right.$$

(20)

$$\begin{array}{l}CL=c^{\mathrm{l}}{F}_{\beta ,t}^{\mathrm{cut}}=c^{\mathrm{l}}(b_1+{\displaystyle \frac{1}{1-{\beta }_1}}{E}_{{\widehat{\mathbb{P}}}_M}({[{\tilde{P}}_t^{\mathrm{cut}}-b_1]}^{+}))\\ =c^{\mathrm{l}}{E}_{{\widehat{\mathbb{P}}}_M}(H_{\beta ,t}^{\mathrm{cut}})\end{array}$$

(21)

$$\begin{array}{l}CE=c^{\mathrm{e}}{F}_{\beta ,t}^{\mathrm{S}}=b_2+{\displaystyle \frac{1}{1-{\beta }_2}}{E}_{{\widehat{\mathbb{P}}}_M}({[{\tilde{S}}_t^{\mathrm{re}}-b_2]}^{+})\\ =E_{{\widehat{\mathbb{P}}}_M}(H_{\beta ,t}^{\mathrm{S}})\end{array}$$

(22)

In the equation, ${\rho }_i^{\mathrm{G}}$ denotes the cost of spinning reserve capacity for thermal power units; ai, bi, ci denote the coefficients of a polynomial; $P_{i,t}^{\mathrm{G}}$ denotes the power output of unit i at time t; CL represents load shedding costs; CE represents curtailment costs; cl is the penalty cost for unit capacity load shedding and ce is the penalty cost for unit capacity curtailment.

(2)

Equations governing power flow

$$\begin{array}{l}{P}_{i,t}={\displaystyle \sum _{j\in B}{G}_{i,j}{V}_{i,t}{V}_{j,t}\cos {\theta }_{i,j,t}}+\\ {\displaystyle \sum _{j\in B}{B}_{i,j}{V}_{i,t}{V}_{j,t}\sin {\theta }_{i,j,t}}\begin{array}{cc},& i\in B\end{array}\end{array}$$

(23)

$$P_{i,t}=P_{i,t}^g+P_{i,t}^{\mathrm{wind}}+P_{i,t}^{\mathrm{pv}}-P_{i,t}^{\mathrm{d}}$$

(24)

$$\begin{array}{l}{Q}_{i,t}=-{\displaystyle \sum _{j\in B}{B}_{i,j}{V}_{i,t}{V}_{j,t}\cos {\theta }_{i,j,t}}+\\ {\displaystyle \sum _{j\in B}{G}_{i,j}{V}_{i,t}{V}_{j,t}\sin {\theta }_{i,j,t}}\begin{array}{cc},& i\in B\end{array}\end{array}$$

(25)

$$Q_{i,t}=Q_{i,t}^{\mathrm{r}}-Q_{i,t}^{\mathrm{d}}$$

(26)

$$P_{i,j,t}^{\mathrm{l}}=-V_{i,t}^2{G}_{i,j}+V_{i,t}{V}_{j,t}(G_{i,j}\cos {\theta }_{i,j,t}+B_{i,j}\sin {\theta }_{i,j,t})$$

(27)

$$-{\overline{P}}^{\mathrm{l}}\le P_{i,j,t}^{\mathrm{l}}\le {\overline{P}}^{\mathrm{l}}\begin{array}{cc},& (i,j)\end{array}\in L$$

(28)

$${\underset{¯}{V}}_i\le V_{i,t}\le {\overline{V}}_i$$

(29)

$${\underset{¯}{Q}}_i^{\mathrm{r}}\le Q_{i,t}^{\mathrm{r}}\le {\overline{Q}}_i^{\mathrm{r}}$$

(30)

In the equations, B and L represent the sets of buses and lines respectively; $P_{i,t}^{\mathrm{wind}}$ and $P_{i,t}^{\mathrm{pv}}$ denote the forecasted output of wind farms and photovoltaic fields at time period t; $P_{i,t}^g$ represents the thermal power output at node i; $P_{i,t}^{\mathrm{d}}$ denotes the conventional load in the system; $Q_{i,t}^{\mathrm{d}}$ represents reactive power demand in the system; $Q_{i,t}^{\mathrm{r}}$,${\underset{¯}{Q}}_i^{\mathrm{r}}$,${\overline{Q}}_i^{\mathrm{r}}$ respectively represent reactive power generation, lower bound of reactive power output, upper bound of reactive power output; $P_{i,j,t}^{\mathrm{l}}$ represents line power flow; $V_{i,t}$,${\overline{V}}_i$,${\underset{¯}{V}}_i$respectively represent node voltage, upper bound of node voltage, lower bound of node voltage.

(3)

Limitations on the Commencement and Termination of Operations for Thermal Power Units.

$$\left\{\begin{array}{l}[t_{i(t-1)}^{\mathrm{on}}-T_{i\min }^{\mathrm{on}}][X_{i(t-1)}-X_{it}]\ge 0\\ [t_{i(t-1)}^{\mathrm{off}}-T_{i\min }^{\mathrm{off}}][X_{it}-X_{i(t-1)}]\ge 0\end{array}\right.$$

(31)

In the equation, $T_{i\min }^{\mathrm{on}}$/$T_{i\min }^{\mathrm{off}}$ represents the minimum start-up/shut-down time for unit i; $t_{i(t-1)}^{\mathrm{on}}$/$t_{i(t-1)}^{\mathrm{off}}$ represents the duration of operation/shut-down from unit i to time period t-1; $X_{it}$ is a binary integer variable representing the start-up/shut-down status of unit i in time period t, with a value of 1 indicating start-up and a value of 0 indicating shut-down.

(4)

Restrictions on the Output Limitation of Thermal Power Units

$$X_{it}{P}_i^{\mathrm{G},\min }\le P_{i,t}^{\mathrm{G}}\le X_{it}{P}_i^{\mathrm{G},\max }$$

(32)

In the equation, $P_i^{\mathrm{G},\max }$ and $P_i^{\mathrm{G},\min }$represent the upper and lower bounds of unit i's output, respectively.

(5)

Limitations on the Spinning Reserve Capacity of Thermal Power Units

$$\left\{\begin{array}{l}{P}_{i,t}^{\mathrm{G}}-R_{i,t}^{\mathrm{down}}\ge X_{it}{P}_i^{\mathrm{G},\min }\begin{array}{cc},& R_{i,t}^{\mathrm{down}}\ge 0\end{array}\\ P_{i,t}^{\mathrm{G}}+R_{i,t}^{\mathrm{up}}\le X_{it}{P}_i^{\mathrm{G},\max }\begin{array}{cc},& R_{i,t}^{\mathrm{up}}\ge 0\end{array}\end{array}\right.$$

(33)

In the equation, $R_{i,t}^{\mathrm{up}}$ and $R_{i,t}^{\mathrm{down}}$ represent the upper and lower spinning reserve capacity of thermal power unit i in time period t, respectively.

(6)

Limitations on the Ramp Rate of Thermal Power Units

$$\left\{\begin{array}{l}{P}_{i,t}^{\mathrm{G}}+R_{i,t}^{\mathrm{up}}-(P_{i,t-1}^{\mathrm{G}}-R_{i,t-1}^{\mathrm{down}})\le r_i^{\mathrm{up}}\\ P_{i,t-1}^{\mathrm{G}}+R_{i,t-1}^{\mathrm{up}}-(P_{i,t}^{\mathrm{G}}-R_{i,t}^{\mathrm{down}})\le r_i^{\mathrm{down}}\end{array}\right.$$

(34)

In the equation, $r_i^{\mathrm{up}}$ and $r_i^{\mathrm{down}}$ represent the ramp-up rate and ramp-down rate of thermal power unit i, respectively.

2.2. Wasserstein Fuzzy Sets for Characterizing Uncertainty in Wind and Solar Resources

In the objective function of the model presented in Equation (19), the expected value of system contingency load shedding and curtailment condition risk loss is approximated based on the empirical distribution ${\widehat{\mathbb{P}}}_M$ of wind and solar prediction errors. However, this approach only considers historical sample data and cannot accurately represent the true distribution $\mathbb{P}$ of wind and solar prediction errors. The fundamental challenge in robust optimization with distribution lies in constructing a fuzzy set that can authentically reflect the uncertain prediction error distribution of new energy sources. In this paper, we propose a fuzzy set based on Wasserstein distance, which effectively measures differences between probability distributions. Specifically, given two distributions $\mathbb{Q}$ and $\mathbb{P}$ within a compact support spaceΞ, the Wasserstein distance is defined as,

$$W(\mathbb{Q},\mathbb{P})=\underset{\pi \in m(E\times E)}{\inf }\{{\displaystyle {\int }_{E\times E}\mathrm{d}({\xi }^q,{\xi }^p)\Pi }(\mathrm{d}{\xi }^q,\mathrm{d}{\xi }^p)\}$$

(35)

$$(\mathrm{d}{\xi }^q,\mathrm{d}{\xi }^p)=‖{\xi }^q-{\xi }^p‖$$

(36)

In the equation, ${\xi }^q$ and ${\xi }^p$ represent stochastic variables; $\Pi$ denotes the joint distribution of $\mathbb{Q}$ and $\mathbb{P}$; Equation (36) defines the distance between ${\xi }^q$ and ${\xi }^p$; $‖·‖$ represents a norm operation, with this paper utilizing the $\infty$-norm to enhance model solvability. The fuzzy set based on Wasserstein distance is defined as,

$${\mathbb{B}}_{\epsilon }(\mathbb{Q})=\{\mathbb{P}|W(\mathbb{P},\mathbb{Q})\le \epsilon \}$$

(37)

The variable $\epsilon$ represents a constant derived from historical data, while ${\mathbb{B}}_{\epsilon }(\mathbb{Q})$ denotes a Wasserstein ball with a radius of $\epsilon$ and centered at $\mathbb{Q}$.

The fuzzy set based on the Wasserstein distance should be centered around a specific distribution $\mathbb{Q}$, as indicated by Equation (37). However, due to the unknown nature of the error distribution in wind and solar power generation, this center should be determined based on relevant historical data to ensure alignment with the actual distribution of uncertain output errors. Assuming we have a finite historical sample dataset $\left\{{\psi }^1,{\psi }^2,\cdots ,{\psi }^M\right\}$ for an uncertain variable $\psi$ following a certain distribution $\mathbb{P}$, we can construct an empirical distribution ${\widehat{\mathbb{P}}}_M$ to approximate the true distribution. Given that the true distribution typically does not deviate significantly from the empirical one, we select the empirical distribution as the center of our Wasserstein ball for estimating the unknown true distribution. The expression for this empirical distribution is as follows,

$${\widehat{\mathbb{P}}}_M={\displaystyle \frac{1}{M}}{\displaystyle \sum _{j=1}^M{\delta }_{{\psi }^j}}$$

(38)

In the equation, M represents the sample size, ${\delta }_{{\psi }^j}$ is the Dirac measure (if ${\psi }^j\in A$, then ${\delta }_{{\psi }^j}(A)=1$; otherwise ${\delta }_{{\psi }^j}(A)=0$, A denotes any set belonging to $\sigma (\Xi )$, and $\Xi$ signifies the support set of $\psi$. As M approaches infinity, ${\widehat{\mathbb{P}}}_M$ should tend toward the true distribution $\mathbb{P}$. In other words, as more data becomes available, the "distance" between ${\widehat{\mathbb{P}}}_M$ and $\mathbb{P}$ will become smaller. The Wasserstein fuzzy set centered around an empirical distribution ${\widehat{\mathbb{P}}}_M$ can be expressed as,

$${\mathbb{B}}_{\epsilon }({\widehat{\mathbb{P}}}_M)=\{\mathbb{P}|W(\mathbb{P},{\widehat{\mathbb{P}}}_M)\le \epsilon \}$$

(39)

In this equation, the parameter A can regulate the level of conservatism in the model, which significantly influences its performance based on the Wasserstein distance. A potential choice is suggested in literature [24],

$$\epsilon =D\sqrt{{\displaystyle \frac{2}{M}}\log ({\displaystyle \frac{1}{1-\eta }})}$$

(40)

In the equation, $\eta$ represents the confidence level, and D denotes the diameter of the support points of the random variable. Nevertheless, literature [25] suggests that $\epsilon$ derived from Equation (40) is excessively conservative and proposes a new alternative as follows,

$$P\{W(\mathbb{P},{\widehat{\mathbb{P}}}_M)\le \epsilon \}\ge 1-\exp (-M{\displaystyle \frac{{\epsilon }^2}{c^2}})$$

(41)

The calculation method for $\epsilon$ can be derived by setting the right-hand side of Equation (41) equal to the confidence level $\eta$ of the given fuzzy set.

$$\epsilon =C\sqrt{{\displaystyle \frac{1}{M}}\ln ({\displaystyle \frac{1}{1-\eta }})}$$

(42)

In the equation, C represents a constant. If $C=\sqrt{2}D$ is true, then Equation (42) transforms into Equation (40). The value of C can be estimated using historical data, and the calculation formula is as follows,

$$\begin{array}{l}C\le 2{\min }_{\rho >0}({\displaystyle \frac{1}{2\rho }}(1+\ln {\mathbb{E}}_{\mathbb{P}}[e^{\rho {‖\lambda -\widehat{\mu }‖}_1^2}]))\\ \approx {\min }_{\rho >0}2\sqrt{{\displaystyle \frac{1}{2\rho }}(1+\ln ({\displaystyle \frac{1}{M}}{\displaystyle \sum _{j=1}^M{e}^{\rho {‖\lambda -\widehat{\mu }‖}_1^2}}))}\end{array}$$

(43)

In this equation, $\widehat{\mu }$ represents the average of the historical sample, specifically defined as $\widehat{\mu }={\displaystyle \frac{1}{M}}{\displaystyle \sum _{j=1}^M{\psi }^j}$.

3. Transformation of the Model

3.1. Linearization of Trends

In the aforementioned model, the power flow equations for AC systems exhibit non-convex and nonlinear characteristics, rendering direct solutions impractical. Consequently, linearizing the power flow equations for AC systems becomes imperative. The expansion of Equations (23) and (25) results in,

$$\begin{array}{l}{P}_{i,t}=g_{i,i}{V}_{i,t}^2+{\displaystyle \sum _{j\in B,j\ne i}[g_{i,j}{V}_{i,t}(V_{i,t}-V_{j,t}\cos {\theta }_{i,j,t})}-\\ b_{i,j}{V}_{i,t}{V}_{j,t}\sin {\theta }_{i,j,t}]\begin{array}{cc},& i\in B\end{array}\end{array}$$

(44)

$$\begin{array}{l}{Q}_{i,t}=-b_{i,i}{V}_{i,t}^2-{\displaystyle \sum _{j\in B,j\ne i}[b_{i,j}{V}_{i,t}(V_{i,t}-V_{j,t}\cos {\theta }_{i,j,t})}+\\ g_{i,j}{V}_{i,t}{V}_{j,t}\sin {\theta }_{i,j,t}]\begin{array}{cc},& i\in B\end{array}\end{array}$$

(45)

Approximate the non-linear terms in Equations (44) and (45) as outlined below,

$$\left\{\begin{array}{l}{V}_{i,t}^2\approx V_{i,t}\\ V_{i,t}{V}_{j,t}\sin {\theta }_{i,j,t}\approx {\theta }_{i,t}-{\theta }_{j,t}\\ V_{i,t}(V_{i,t}-V_{j,t}\cos {\theta }_{i,j,t})\approx V_{i,t}-V_{j,t}\end{array}\right.$$

(46)

The ultimate outcome produces the subsequent decoupled linear power flow.

$$\left\{\begin{array}{l}{P}_{i,t}={\displaystyle \sum _{j\in B}{G}_{i,j}{V}_{j,t}}-{\displaystyle \sum _{j\in B}{B}_{i,j}^{\text{'}}{\theta }_{j,t}}\\ Q_{i,t}=-{\displaystyle \sum _{j\in B}{B}_{i,j}{V}_{j,t}}-{\displaystyle \sum _{j\in B}{G}_{i,j}^{\text{'}}{\theta }_{j,t}}\end{array}\right.\begin{array}{cc},& i\in B\end{array}$$

(47)

Likewise, the formula for line power is as follows.

$$P_{i,j,t}^{\mathrm{l}}=G_{i,j}(V_{i,t}-V_{j,t})-B_{i,j}({\theta }_{i,t}-{\theta }_{j,t})\begin{array}{cc},& (i,j)\end{array}\in L$$

(48)

3.2. The Dual Transformation of a Distribution Robust Model

Distributionally robust optimization combines the principles of stochastic optimization and robust optimization, typically framed as a minimax problem. The model is outlined as follows,

$$\underset{x}{\inf }\underset{P\in {\Omega }_P(\tilde{\omega })}{\sup }{E}_P\{f(x,\tilde{\omega })\}$$

(49)

In this formula, ${\Omega }_P(\tilde{\omega })$ represents the fuzzy set of uncertain quantity $\tilde{\omega }$; $\underset{P\in {\Omega }_P(\tilde{\omega })}{\sup }$ represents the worst-case distribution of uncertain quantity ${\Omega }_P(\tilde{\omega })$ within fuzzy set $\underset{P\in {\Omega }_P(\tilde{\omega })}{\sup }$; and $\underset{x}{\inf }$ represents the lower bound of the expected value of $f(x,\tilde{\omega })$ under the worst-case distribution of $\tilde{\omega }$.

Distributionally robust optimization utilizes fuzzy set ${\Omega }_P(\tilde{\omega })$ to represent the higher-order uncertainty of uncertain quantity $\tilde{\omega }$, encompassing both the intrinsic uncertainty of the quantity and the ambiguity in its probability distribution information.

The extreme distribution problem of model (49) can be reformulated in accordance with the strong duality theory [27] as follows,

$$\begin{array}{l}\underset{x}{\inf }\underset{P\in {\Omega }_P(\tilde{\omega })}{\sup }{E}_P\{f(x,\tilde{\omega })\}=\\ \underset{\kappa \ge 0}{\inf }\left\{\kappa \cdot \epsilon +{\displaystyle \frac{1}{M}}{\displaystyle \sum _{m=1}^M\underset{\tilde{\omega }\in \Xi }{\sup }}\left(f(x,\tilde{\omega })-\kappa ‖\tilde{\omega }-{\widehat{\omega }}_m‖\right)\right\}\end{array}$$

(50)

The model (50) can be further reformulated as follows,

$$\left\{\begin{array}{l}\underset{\kappa \ge 0}{\inf }\kappa \cdot \epsilon +{\displaystyle \frac{1}{M}}{\displaystyle \sum _{m=1}^M{\tau }_m}\\ \mathrm{s}.\mathrm{t}.f(x,\overline{\omega })-\kappa \cdot \left(\overline{\omega }-{\widehat{\omega }}_m\right)\le {\tau }_m\\ f(x,\underset{\_}{\omega })+\kappa \cdot \left(\underset{\_}{\omega }-{\widehat{\omega }}_m\right)\le {\tau }_m\\ f\left(x,{\widehat{\omega }}_m\right)\le {\tau }_m\end{array}\right.$$

(51)

In this equation, κ represents a dual variable; $\overline{\omega }$ and $\underset{\_}{\omega }$ denote the vectors formed by the boundaries of ${\tilde{\omega }}^i$.

4. Simulation Study

The proposed method was validated by testing it on the IEEE118 node system with an enhanced structure as a case study. The system structure utilized in the case study is depicted in Figure A1 of Appendix A. It comprises 11 thermal power units (G1 - G11) with specific parameters detailed in Table A1 of Appendix A. The installed capacity for wind power is 900 MW, connected to node 25, while the installed capacity for photovoltaic power is 1500 MW, connected to node 65. New energy accounts for a substantial 47% of the total installed capacity, signifying a significant proportion within the system. Additionally, each node's maximum load in the system has been limited to 75% of that in the original IEEE118 node system.

Utilizing the actual and predicted values of historical wind and solar power output, this study modeled the empirical distribution of errors in wind and solar power output. Subsequently, 1000 samples of errors in wind and solar power output were randomly generated at each time point based on the historical empirical distribution.

4.1. Systematic Probabilistic Analysis of Supply and Demand Balance

The paper conducts a probabilistic risk analysis of power generation balance in a high-proportion new energy power system, utilizing proposed load shedding and curtailed power conditions. For instance, the system dispatch optimization results at a 95% confidence level are illustrated in Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12.

Figure 3 depicts the voltage levels at system nodes at a specific time, demonstrating that the system's voltage levels are within permissible limits. This validates the effectiveness of the power flow model proposed in this paper.

As depicted in Figure 4, the output of new energy sources reaches its minimum value in the early morning. Due to the lower limit on the output of new energy sources, downward fluctuations are constrained by uncertainty. Consequently, at this time, the total upward ramping reserve is at its daily nadir.

The system encounters insufficient climbing reserve capacity around the 19th hour, leading to a potential risk of load shedding. This means that relying solely on the regulation capacity of thermal power plants at this time is no longer adequate to manage the uncertainty of new energy output. Therefore, it is advisable to utilize a risk value for the probabilistic supply-demand balance cut-load condition as a guide for configuring storage, load demand response, and other flexible resources.

The risk assessment results for curtailed power under different time periods, as illustrated in Figure 7, indicate that insufficient downward ramping reserve capacity is experienced by the system during the 14th hour, leading to a significant risk of curtailed power.

Moreover, by integrating the operational status and output outcomes of each thermal power unit depicted in Figure 8 and Figure 9, it is evident that during the time intervals of 0-10 and 17-24 hours, thermal power units G1-G4 are in operation. Conversely, during the period of 10-17 hours, units G5, G6, G10, and G11 are active. The operational state of these units is predominantly influenced by their operating cost coefficients. Consequently, during the time slots of 0-10 and 17-24 hours, reliance primarily rests on thermal power units G1-G4 to furnish reserves for managing uncertainties in new energy output; whereas during the period of 10-17 hours reliance mainly shifts to G5, G6, G10, and G11 for providing upward ramping reserves.

For instance, during the 10-17 time period, the active power generated by thermal power units G5 and G6 is notably below their installed capacity. Consequently, they are able to provide a higher upward ramping reserve capacity. The output and reserve capacity configuration of unit G6 is depicted in Figure 11.

The output and reserve configuration of thermal power unit G1 are depicted in Figure 10. During the peak load period from 17 to 24 hours, in the decision-making process of power system dispatch, G1 operates at full capacity to optimize economic benefits, primarily providing downward ramping reserve capacity.

Figure 12 illustrates the correlation between the total output of the power plant and the ramp reserve. During the 1-8 hour period in the early morning, there is minimal fluctuation in the output of new energy sources, leading to a small allocation of reserve capacity by the system. However, during the 9-24 hour period, significant fluctuations in new energy source output necessitate a larger allocation of reserve capacity to manage uncertainty. At 14 hours, when the system's thermal power output reaches its lowest point, it provides ample up reserve capacity but insufficient down reserve capacity.

4.2. Analyzing the Influence of Varying Confidence Levels on Scheduling Outcomes

In order to validate the proposed distributed robust coordinated optimization scheduling model under various levels of confidence in load shedding and curtailed wind power risk, this study sets six distinct confidence levels. The corresponding optimized scheduling results are obtained and presented in Figure 13.

Based on Figure 13, it is evident that as the confidence level increases while maintaining the same thermal power regulation capacity, the associated risks also escalate. The expenses related to curtailed wind power and load shedding gradually mount. At a confidence level of 95%, the costs for curtailed wind power and load shedding amount to 36,500 yuan and 41,000 yuan respectively.

The higher the confidence level, the greater the need for adjustable resources that the system must allocate in order to manage fluctuations in new energy power and mitigate uncertainty at its source. Simultaneously, system reserve costs also increase as confidence levels rise, reaching 86,900 yuan when the confidence level is 95%.

To assess the scheduling performance of the proposed distributed robust coordinated optimization scheduling model under various levels of fuzzy set confidence, this study establishes six distinct fuzzy set confidence levels. The corresponding optimized scheduling results are computed and presented in Figure 14.

Based on the findings illustrated in Figure 14, it is evident that a higher fuzzy set confidence level corresponds to a larger Wasserstein ball radius and increased uncertainty in new energy output. Consequently, as the fuzzy set confidence level rises, so does the expected cost of tail risk for the system. At a 95% confidence level, the costs of curtailment and load shedding are RMB 36,500 and RMB 41,000 respectively.

4.3. Analyzing Scheduling Outcomes Employing Various Optimization Techniques

The distributed robust optimization method proposed in this paper is further compared with the robust optimization method and stochastic optimization. In the robust optimization method, the new energy output error at each time is used as the boundary of the sample for new energy output errors at each time; in distributed robust optimization, a conditional risk confidence level of β=80% and a fuzzy set confidence level of η=95% are set; for the stochastic optimization model, the confidence parameter value is set to $\beta (—)$=80%. Subsequently, based on historical empirical distribution of wind and solar power output errors, 2000 samples are randomly generated to simulate external random scenarios for experimentation. The curtailment probability and load shedding probability for each optimization method are then statistically analyzed. The results of these three methods are presented in Table 2, while adjustments made for reserve allocation are illustrated in Figure 15.

According to Table 2, it is evident that the distributed robust optimization method based on conditional value-at-risk proposed in this paper has led to an approximate increase of 81.4% and 0.7% in total ramping requirements and total scheduling costs respectively, compared to the stochastic optimization method. However, it has significantly reduced the probabilities of load shedding and curtailment. In contrast, while the robust optimization method can further reduce supply-demand balance risk, it also results in an approximately 7.2% increase in total ramping requirements compared to our proposed approach. The comparative results suggest that our proposed method effectively balances safety and economic demands by enhancing optimization robustness without incurring excessive additional costs due to overly conservative measures.

5. Conclusions

In order to effectively evaluate the potential risk of power and energy balance disruption resulting from the variability of wind and solar power generation, as well as mitigate the likelihood of curtailment and supply shortages in modern power systems, this study introduces a probabilistic index for assessing power and energy balance risk based on conditional value at risk (CVaR) theory. Additionally, it develops a robust optimization dispatch model that incorporates flexible ramping reserves from coal-fired power plants. The key findings are as follows,

• Tailored safety and economic risk indicators for modern power systems' power and energy balance have been proposed, which can be used in scheduling decisions to achieve operational equilibrium between safety and economic risks.
• By fully considering the flexible ramping reserves provided by coal-fired power plants, the results of probabilistic risk assessment can guide the rational planning of flexible resources such as energy storage and demand-side response within the system.