A Probabilistic Reserve Decision-making Method Based on Cumulative Probability Approximation for High-Penetration Renewable Energy Power System

1. Introduction

The demand for operating reserve capacity in power systems depends on the net load reserve capacity requirement and the contingency reserve capacity requirement [1]. The contingency reserve is generally determined based on the maximum unit generation capacity or the maximum DC import capacity, which is a fixed value independent of load power and renewable energy generation [2]. The net load reserve capacity requirement, on the other hand, depends on load forecast errors and renewable energy output forecast errors. The primary focus of operating reserve decisions lies in net load reserve decisions [3]. Historically, China has set operating reserve capacity based on deterministic criteria, typically ranging from 2% to 5% of the system’s peak load [4]. Deterministic criteria are based on load forecasts and renewable energy forecasts to formulate generation plans, with reserve decisions being manually made. This leads to arbitrariness and blindness in the selection of reserve capacity. When the reserve capacity is too low, it becomes difficult to match the significant fluctuations in renewable energy output and load power, thus jeopardizing system security. Conversely, excessive reserve capacity results in the wastage of reserve resources, affecting the economic efficiency of system operation [5].

With the continuous growth of electricity demand and the ongoing progress in renewable energy technologies, the integration of large-scale wind and solar power generation into the power system has led to a significant increase in the penetration of renewable energy in the grid [6]. The large-scale integration of wind and solar power generation increases the uncertainty of generation output, while the existing deterministic load reserve criteria lack enough theoretical support, resulting in overly conservative scenarios or risky decisions in practical applications [7]. Moreover, the combined uncertainties of renewable energy and load significantly increase the volatility of net load, making it difficult for traditional, experience-based reserve decisions to respond quickly and accurately. There is an urgent need for more precise and flexible reserve decision-making methods [8].

At present, probabilistic modeling of renewable energy output forecast errors and load forecast errors serves as an important basis for analyzing the impact of forecast errors and controlling the consequences of these errors. The reserve decision for high-penetration renewable energy systems heavily depends on the accuracy of the net load forecast error probability model [9,10]. The normal distribution is a commonly used model for fitting renewable energy output forecast errors and load forecast errors [11]. For example, literature [12] adopts a normal distribution to fit the probabilistic model of wind power output forecast errors, thereby facilitating reserve decision-making and comparing it with traditional deterministic criteria. Literature [13] fits the forecast errors of load, wind power, and solar power by using normal distribution models to simplify the aggregation of the respective probability density functions. However, normal distribution is typically used to ensure consistency with the sample mean and standard deviation in practical applications and it lacks adequate precision in estimating the small probability intervals of both tails [14]. Furthermore, statistical results show that the forecast errors of renewable energy output do not strictly follow a normal distribution [15,16]. To improve the overall fitting accuracy of the probability model, some studies adopt other parametric distributions for probabilistic modeling of renewable energy output and load forecast errors, such as the beta distribution [17], Laplace distribution [18], Cauchy distribution [19], Weibull distribution [20] and Gaussian mixture models (GMM) [21]. However, the above methods may increase the complexity of the probabilistic decision problem and negatively impacts the efficiency of solving these problems, while the accuracy of fitting the small-probability regions remains difficult to guarantee [22].

Currently, commonly used probabilistic modeling methods focus on improving the overall fitting accuracy of the probability model. However, for reserve decisions in power systems with a high penetration of renewable energy, no reserve is required when the net load does not exceed the forecast. Hence, the reserve decision model is only related to the probability distribution of net load exceeding the forecasted value [23]. In addition, existing probabilistic reserve decisions are typically based on a certain confidence level to determine the reserve capacity. However, the selection of the confidence level lacks theoretical justification. For example, literature [24,25,26] all propose opportunity constraint approaches based on a certain confidence level, but none provide a theoretical basis for determining the confidence level.

Therefore, this paper proposes a probabilistic reserve decision-making method based on cumulative probability approximation for high-penetration renewable energy power systems. Compared with existing technologies, the advantages of the proposed technical solution are mainly reflected in the following aspects:

1)

The method comprehensively considers three factors: reserve cost, expected generation revenue from reserves, and expected loss due to load shedding. It provides optimal conditions for probabilistic reserve decisions and an iterative solving algorithm, achieving the optimal expected overall benefit in probabilistic reserve decision-making.

2)

The method determines the parameters of the net load forecast error’s normal distribution model based on the key factors influencing the optimality of reserve decisions. By employing cumulative probability approximation at key points, the accuracy of the probabilistic reserve decision results is enhanced.

3)

By selecting key points for cumulative probability approximation, the proposed method effectively improves the fitting accuracy of the cumulative probability in the small probability interval of the tail.

The rest the paper is arranged as follows: Section 2 presents the procedures for the proposed method. Section 3 discusses the required data and the preprocessing operations. Section 4 plots the cumulative probability distribution curve of the samples. Section 5 performs probabilistic reserve decision modeling. Section 6 estimates the parameters of the normal distribution model for net load forecasting errors. Section 7 iteratively solves the reserve results. Section 8 provides a case analysis. Section 9 concludes the paper.

2. Procedures for the Proposed Method

The flowchart of probabilistic reserve decision method based on cumulative probability approximation for high-penetration renewable energy power systems proposed in this paper is shown in Figure 1, which includes the following steps:

1)

Data loading and preprocessing;

2)

Calculation of the cumulative probability distribution of historical net load forecast error sample data;

3)

Probabilistic reserve decision modeling based on the normal distribution probability model of net load forecast errors;

4)

Estimation of the normal distribution model parameters based on cumulative probability approximation at key points;

5)

Probabilistic reserve decision-making based on the normal distribution model of net load forecast errors.

3. Data Loading and Preprocessing

The data required for the reserve decision-making of high-penetration renewable energy power systems include:

1)

historical actual and forecasted data of renewable energy output;

2)

historical actual and forecasted data of load;

3)

unit cost of reserve capacity;

4)

unit revenue from increased power generation during reserve activation;

5)

unit cost of load shedding.

Data preprocessing includes:

1)

subtracting the actual renewable energy output from the historical actual load data at the same time to obtain the actual net load values;

2)

subtracting the forecasted renewable energy output from the historical forecasted load data at the same time to obtain the forecasted net load values;

3)

subtracting the historical forecasted values from the historical actual net load values at the same time to obtain the net load forecast errors.

4. Cumulative Probability Distribution Curve of the Net Load Forecast Errors

The net load forecast errors are sorted in ascending order. Assuming the sorted net load forecast errors are $x_1,x_2\cdots x_N$, where N represents the number of samples. Each sample corresponds to a probability of ${\displaystyle \frac{1}{N+1}}$. Therefore, the cumulative probabilities for $x_1,x_2\cdots x_N$ are ${\displaystyle \frac{1}{N+1}}$,${\displaystyle \frac{2}{N+1}},$$\cdots$,${\displaystyle \frac{N}{N+1}}$, respectively. Using the net load forecast errors as the abscissa (x-axis) and the sample cumulative probability as the ordinate (y-axis), the sample points are sequentially connected to form the sample cumulative probability curve. Let γ(x) denote the sample cumulative probability function.

Figure 2. The cumulative probability distribution curve of net load forecast errors.

Figure 2. The cumulative probability distribution curve of net load forecast errors.

5. Probabilistic Reserve Decision Model

The objective of probabilistic reserve decision-making is to minimize the total cost, which is the sum of the reserve configuration cost and the expected cost of power shortage per unit time, after subtracting the revenue from reserve activation:

$$\min J\left(R\right)=C_{\mathrm{R}}R+C_{\mathrm{EDNS}}{P}_{\mathrm{EDNS}}\left(R\right)-C_{\mathrm{INC}}{P}_{\mathrm{INC}}\left(R\right)$$

(1)

Where, C_R represents the unit cost of reserve capacity, R represents the reserve decision capacity. C_EDNS represents the unit cost of power shortage expectation, P_EDNS(R) represents the power shortage expectation, which is a function of reserve capacity. C_INC represents the unit benefit of additional generation when reserve capacity is dispatched, and P_INC(R) represents the expected power generation of reserve capacity.

The normal distribution is commonly used to describe the distribution of forecasting errors in renewable energy output and load prediction. Since the sum or difference of normal distributed random variables also follows a normal distribution, the normal distribution can be used to model the net load forecast error. The probability density function of the normal distribution can be expressed as:

$$\rho \left(x|\mu ,\sigma \right)={\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}{\mathrm{e}}^{\left(-{\displaystyle \frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}}\right)}，-\infty <x<\infty$$

(2)

Here, x represents the random variable, which refers to the net load forecast error. μ denotes the mean of the normal distribution, and σ represents the standard deviation of the normal distribution, with $\sigma >0$.

If the system’s net load forecast error exceeds the reserve capacity of the system, the portion of the load beyond the reserve capacity will be lost, and the system’s power shortage expectation is:

$$\begin{array}{l}{P}_{\mathrm{EDNS}}\left(R\right)={\displaystyle {\int }_R^{\infty }\left(x-R\right)}d\left(x|\mu ,\sigma \right)dx\\ ={\displaystyle {\int }_R^{\infty }\left(x-R\right)}{\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}{\mathrm{e}}^{\left(-{\displaystyle \frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}}\right)}dx\\ ={\displaystyle {\int }_R^{\infty }\frac{x-\mu }{\sqrt{2\pi }\sigma }{\mathrm{e}}^{\left(-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right)}dx}+\left(\mu -R\right){\displaystyle {\int }_R^{\infty }\frac{1}{\sqrt{2\pi }\sigma }{\mathrm{e}}^{\left(-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right)}dx}\\ ={\displaystyle {\int }_{\frac{{\left(R-\mu \right)}^2}{2{\sigma }^2}}^{\infty }\frac{\sigma }{\sqrt{2\pi }}{\mathrm{e}}^{\left(-t\right)}dt}+\left(\mu -R\right){\displaystyle {\int }_{\frac{R-\mu }{\sigma }}^{\infty }\frac{1}{\sqrt{2\pi }}{\mathrm{e}}^{\left(-\frac{z^2}{2}\right)}dz}\\ ={\displaystyle \frac{\sigma }{\sqrt{2\pi }}}{\mathrm{e}}^{\left(-{\displaystyle \frac{{\left(R-\mu \right)}^2}{2{\sigma }^2}}\right)}+\left(\mu -R\right){\displaystyle {\int }_{-\infty }^{\frac{\mu -R}{\sigma }}\frac{1}{\sqrt{2\pi }}{\mathrm{e}}^{\left(-\frac{z^2}{2}\right)}dz}\\ =\sigma \varphi \left({\displaystyle \frac{\mu -R}{\sigma }}\right)+\left(\mu -R\right)\Phi \left({\displaystyle \frac{\mu -R}{\sigma }}\right)\end{array}$$

(3)

Where, R represents the decision reserve capacity. P_EDNS(R) denotes the power shortage expectation, which is a function of the reserve capacity decision.

$$t={\displaystyle \frac{{(x-\mu )}^2}{{\sigma }^2}}$$

(4)

$$z={\displaystyle \frac{x-\mu }{\sigma }}$$

(5)

Where ϕ(⋅) represents the probability density function of the standard normal distribution, and Φ(⋅) represents the cumulative distribution function of the standard normal distribution.

The standard normal distribution refers to a normal distribution with a mean of zero and a standard deviation of one, and its probability density function is:

$$\varphi \left(x\right)={\displaystyle \frac{1}{\sqrt{2\pi }}}{\mathrm{e}}^{\left(-{\displaystyle \frac{x^2}{2}}\right)}，-\infty <x<\infty$$

(6)

The expected generation power of the decision reserve capacity consists of two parts: when the net load forecast error x is greater than zero but less than the decision reserve capacity R, the generation power of the reserve capacity is x; when the net load forecast error x exceeds the decision reserve capacity R, the generation power of the reserve capacity is R. Therefore, the expected energy generation of the decision reserve capacity per unit time is:

$$\begin{array}{l}{P}_{\mathrm{INC}}\left(R\right)={\displaystyle {\int }_0^{R}x}\rho \left(x|\mu ,\sigma \right)dx+R{\displaystyle {\int }_R^{\infty }\rho \left(x|\mu ,\sigma \right)dx}\\ ={\displaystyle {\int }_0^{R}x}{\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}{\mathrm{e}}^{\left(-{\displaystyle \frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}}\right)}dx+R{\displaystyle {\int }_0^R\frac{1}{\sqrt{2\pi }\sigma }{\mathrm{e}}^{\left(-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right)}dx}\\ ={\displaystyle {\int }_0^R\frac{x-\mu }{\sqrt{2\pi }\sigma }{\mathrm{e}}^{\left(-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right)}dx}+(\mu +R){\displaystyle {\int }_0^R\frac{1}{\sqrt{2\pi }\sigma }{\mathrm{e}}^{\left(-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right)}dx}\\ ={\displaystyle {\int }_{\frac{{\mu }^2}{2{\sigma }^2}}^{\frac{{\left(R-\mu \right)}^2}{2{\sigma }^2}}\frac{\sigma }{\sqrt{2\pi }}{\mathrm{e}}^{\left(-t\right)}dt}+(\mu +R){\displaystyle {\int }_{\frac{-\mu }{\sigma }}^{\frac{R-\mu }{\sigma }}\frac{1}{\sqrt{2\pi }}{\mathrm{e}}^{\left(-\frac{z^2}{2}\right)}dz}\\ ={\displaystyle \frac{\sigma }{\sqrt{2\pi }}}\left[{\mathrm{e}}^{\left(-{\displaystyle \frac{{\mu }^2}{2{\sigma }^2}}\right)}-{\mathrm{e}}^{\left(-{\displaystyle \frac{{\left(R-\mu \right)}^2}{2{\sigma }^2}}\right)}\right]+(\mu +R){\displaystyle {\int }_{\frac{-\mu }{\sigma }}^{\frac{R-\mu }{\sigma }}\frac{1}{\sqrt{2\pi }}{\mathrm{e}}^{\left(-\frac{z^2}{2}\right)}dz}\\ =\sigma \left[\varphi \left({\displaystyle \frac{\mu }{\sigma }}\right)-\varphi \left({\displaystyle \frac{\mu -R}{\sigma }}\right)\right]+(\mu +R)\left[\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)-\Phi \left(-{\displaystyle \frac{\mu }{\sigma }}\right)\right]\end{array}$$

(7)

Where P_INC(R) represents the expected generation power per unit time of the decision reserve capacity, which is a function of the decision reserve capacity.

The necessary condition for the extremum of Equation (1) is:

$${\displaystyle \frac{dJ\left(R\right)}{dR}}=C_{\mathrm{R}}+C_{\mathrm{EDNS}}{\displaystyle \frac{d{P}_{\mathrm{EDNS}}\left(R\right)}{dR}}-C_{\mathrm{INC}}{\displaystyle \frac{d{P}_{\mathrm{INC}}\left(R\right)}{dR}}=0$$

(8)

Due to:

$$\begin{array}{l}{\displaystyle \frac{d{P}_{\mathrm{EDNS}}\left(R\right)}{dR}}\\ ={\displaystyle \frac{\sigma }{\sqrt{2\pi }}}\exp \left(-{\displaystyle \frac{{\left(R-\mu \right)}^2}{2{\sigma }^2}}\right)\left({\displaystyle \frac{\mu -R}{{\sigma }^2}}\right)-\Phi \left({\displaystyle \frac{\mu -R}{\sigma }}\right)-\left(\mu -R\right){\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}\exp \left(-{\displaystyle \frac{{\left(R-\mu \right)}^2}{2{\sigma }^2}}\right)\\ =-\Phi \left({\displaystyle \frac{\mu -R}{\sigma }}\right)=\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)-1\end{array}$$

(9)

$$\begin{array}{l}{\displaystyle \frac{d{P}_{\mathrm{INC}}\left(R\right)}{R}}={\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}\exp \left(-{\displaystyle \frac{{\left(R-\mu \right)}^2}{2{\sigma }^2}}\right)\left(R-\mu \right)\\ +\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)-\Phi \left(-{\displaystyle \frac{\mu }{\sigma }}\right)+(\mu +R){\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}\exp \left(-{\displaystyle \frac{{\left(R-\mu \right)}^2}{2{\sigma }^2}}\right)\\ ={\displaystyle \frac{2R}{\sqrt{2\pi }\sigma }}\exp \left(-{\displaystyle \frac{{\left(R-\mu \right)}^2}{2{\sigma }^2}}\right)+\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)-\Phi \left(-{\displaystyle \frac{\mu }{\sigma }}\right)\\ ={\displaystyle \frac{2R}{\sigma }}\varphi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)+\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)-\Phi \left(-{\displaystyle \frac{\mu }{\sigma }}\right)\end{array}$$

(10)

By substituting Equations (9) and (10) into Equation (8), we obtain:

$${\displaystyle \frac{dJ\left(R\right)}{dR}}=C_{\mathrm{R}}-C_{\mathrm{EDNS}}+C_{\mathrm{EDNS}}\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)-C_{\mathrm{INC}}\left[{\displaystyle \frac{2R}{\sigma }}\varphi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)+\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)-\Phi \left(-{\displaystyle \frac{\mu }{\sigma }}\right)\right]=0$$

(11)

By rearranging Equation (9), we obtain:

$$\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)={\displaystyle \frac{C_{\mathrm{EDNS}}-C_{\mathrm{R}}+C_{\mathrm{INC}}\left[\frac{2R}{\sigma }\varphi \left(\frac{R-\mu }{\sigma }\right)-\Phi \left(-\frac{\mu }{\sigma }\right)\right]}{\left(C_{\mathrm{EDNS}}-C_{\mathrm{INC}}\right)}}$$

(12)

6. Parameter Estimation of the Normal Distribution Model by Using Key Points

From Equation (12), we can find that the factors influencing the probabilistic reserve decision include the cumulative probability at a net load forecast error of 0, as well as the cumulative probability and probability density at a net load forecast error of R.

Since the normal distribution has only two parameters, the mean and standard deviation, it has only two degrees of freedom. Therefore, the parameters of the normal distribution can be determined by the two points on the cumulative probability distribution curve.

For actual power systems, $\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)$ is close to 1 and $\varphi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)$ is a smaller positive number since C_EDNS>>C_INC and C_EDNS>>C_R. The key factors influencing the probabilistic reserve decision are the cumulative probability at a net load forecast error of 0 and at a net load prediction error of R.

To improve the accuracy of the reserve decision as much as possible, the cumulative probabilities of the normal distribution probability model for net load forecast errors at 0 and R are set to be equal to the corresponding sample cumulative probabilities at those positions, i.e.:

$$\Phi \left(-{\displaystyle \frac{\mu }{\sigma }}\right)=\gamma \left(0\right)$$

(13)

Where $\gamma (\ast )$ represents the cumulative probability function of net load forecast error.

$$\Phi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)=\gamma \left(R\right)$$

(14)

From Equations (12), (13) and (14) we have:

$$\gamma \left(R\right)={\displaystyle \frac{C_{\mathrm{EDNS}}-C_{\mathrm{R}}+C_{\mathrm{INC}}\left[\frac{2R}{\sigma }\varphi \left(\frac{R-\mu }{\sigma }\right)-\gamma \left(0\right)\right]}{\left(C_{\mathrm{EDNS}}-C_{\mathrm{INC}}\right)}}$$

(15)

Since R is a decision variable, it needs to be determined after the reserve decision and remains an unknown quantity at the time of parameter estimation for the normal distribution model. Considering that $\varphi \left({\displaystyle \frac{R-\mu }{\sigma }}\right)$ is a smaller positive number, the initial value of the decision variable R can be calculated according to Equation (16):

$$\gamma \left(R^{\left(0\right)}\right)={\displaystyle \frac{C_{\mathrm{EDNS}}-C_{\mathrm{R}}-\gamma \left(0\right)C_{\mathrm{INC}}}{\left(C_{\mathrm{EDNS}}-C_{\mathrm{INC}}\right)}}$$

(16)

Where the superscript 0 indicates the initial value.

The right end term of Equation (16) can be obtained from the original parameters and the sample data. Thus, R⁽⁰⁾ can be obtained by linear interpolation of the sample cumulative probability curves.

Comparison of Equation (15) and Equation (16) shows that:

$$\gamma \left(R\right)\ge \gamma \left(R^{\left(0\right)}\right)$$

(17)

Since the sample cumulative probability curve is monotonically increasing:

$$R\ge R^{\left(0\right)}$$

(18)

Equation (18) shows that R⁽⁰⁾ is a lower bound for the probabilistic reserve decision.

By applying R⁽⁰⁾ from the alternative expression (14) of R, we obtain:

$$\Phi \left({\displaystyle \frac{R^{\left(0\right)}-\mu }{\sigma }}\right)=\gamma \left(R^{\left(0\right)}\right)$$

(19)

It can be obtained from Equation (13) that:

$$-{\displaystyle \frac{\mu }{\sigma }}={\Phi }^{-1}\left(\gamma \left(0\right)\right)$$

(20)

It can be obtained from Equation (19) that:

$${\displaystyle \frac{R^{\left(0\right)}-\mu }{\sigma }}={\Phi }^{-1}\left(\gamma \left(R^{\left(0\right)}\right)\right)$$

(21)

By solving Equations (20) and (21) simultaneously, we obtain:

$$\sigma ={\displaystyle \frac{R^{\left(0\right)}}{{\Phi }^{-1}\left(\gamma \left(R^{\left(0\right)}\right)\right)-{\Phi }^{-1}\left(\gamma \left(0\right)\right)}}$$

(22)

$$\mu ={\displaystyle \frac{R^{\left(0\right)}{\Phi }^{-1}\left(\gamma \left(0\right)\right)}{{\Phi }^{-1}\left(\gamma \left(0\right)\right)-{\Phi }^{-1}\left(\gamma \left(R^{\left(0\right)}\right)\right)}}$$

(23)

The solution of Equations (22) and (23) requires the use of the cumulative probability curve of the standard normal distribution to obtain the quantiles from the cumulative probability.

The cumulative probability table of the standard normal distribution can be obtained through numerical integration. The cumulative probability of the standard normal distribution possesses the following properties:

$$\Phi \left(0\right)=0.5$$

(24)

$$\Phi \left(-z\right)=1-\Phi \left(z\right)$$

(25)

Due to $\Phi (5)\approx 0.9999997$, it is enough in practice to construct a cumulative probability table of the standard normal distribution for the range [0, 5]. For values exceeding 5, the cumulative probability can be approximated as 1, which meets the accuracy requirements for engineering applications. When constructing the cumulative probability table of the standard normal distribution via numerical integration, it is recommended to use a small integration step size (e.g., 0.0001) to ensure numerical accuracy. When using the table lookup method, linear interpolation can be applied between the data intervals to obtain the cumulative probability corresponding to a specified value. The integral formula for the cumulative probability within the range [0, 5] is as follows:

$$\Phi \left(z\right)=0.5+{\displaystyle {\int }_0^z\frac{1}{\sqrt{2\pi }}{\mathrm{e}}^{\left(-\frac{s^2}{2}\right)}}ds$$

(26)

The numerical integration formula for calculating the cumulative probability of the standard normal distribution is as follows:

$$\Phi \left((k+1)\Delta z\right)=\Phi \left(k\Delta z\right)+{\displaystyle \frac{\varphi \left(k\Delta z\right)+\varphi \left((k+1)\Delta z\right)}{2}}\Delta z$$

(27)

Where $\Delta z$ represents the integration step size and k denotes the number of integration steps.

7. Probabilistic Reserve Decision-Making

After determining the mean $\mu$ and standard deviation $\sigma$ of the forecast error normal distribution model, the reserve decision variable R can be obtained based on the optimality conditions of Equation (12).

Since the right-hand side of Equation (12) contains the decision variable R, which poses challenges for direct numerical solution, an iterative method based on Equation (28) is constructed to solve this.

$$\Phi \left({\displaystyle \frac{R^{(i+1)}-\mu }{\sigma }}\right)={\displaystyle \frac{C_{\mathrm{EDNS}}-C_{\mathrm{R}}+C_{\mathrm{INC}}\left[\frac{2R^{(i)}}{\sigma }\varphi \left(\frac{R^{(i)}-\mu }{\sigma }\right)-\Phi \left(-\frac{\mu }{\sigma }\right)\right]}{\left(C_{\mathrm{EDNS}}-C_{\mathrm{INC}}\right)}}$$

(28)

In the equation, the superscript i represents the iteration number.

The solution of Equation (28) also requires the use of the standard normal distribution cumulative probability table to obtain the quantile from the cumulative probability.

If the deviation between two consecutive iteration results is smaller than the predefined allowable error, the iteration converges, and R⁽ⁱ⁺¹⁾ is taken as the probabilistic reserve decision result. The convergence criterion is set as follows:

$$\left|R^{(i+1)}-R^{(i)}\right|\le \epsilon$$

(29)

Where $\epsilon$ represents the allowable error.

8. Case Analysis

8.1. Case Test of Belgian Transmission Network

Firstly, a case study analysis is conducted by using the relevant data from Belgian transmission network, published by the Elia Group [27]. The frequency distribution histogram of the net load forecast errors for a specific year of this transmission network is shown in Figure 3. The cumulative probability distribution curve of the net load forecast errors is shown in Figure 4.

It can be observed that situations where the net load significantly exceeds expectations (with net load forecast errors greater than 1000 MW) are concentrated in the tail cumulative probability of less than 2%.

Let the unit cost of load shedding, C_EDNS, be 1000 $/(MW·h) [28], the unit cost of reserve capacity, C_R, be 20 $/(MW·h) [29], and the unit revenue from additional generation during reserve activation, C_INC, be 5 $/(MW·h) [30]. According to Equation (16), it can be obtained:

R⁽⁰⁾=1141.14MW

Further, based on Equations (22) and (23), the parameters of the normal distribution probability model for net load forecast errors can be obtained:

μ=228.19

σ=428.53

The iterative solution is carried out according to Equation (28), setting the permissible deviation of the results from two adjacent iterations to be 0.01. The programming is carried out on the MATLAB R2022b software platform, taking 0.3201 seconds. After 4 iterations, the result converges. The final reserve decision-making value was:

R=1152.40MW

It can be found that the final value R obtained for the reserve decision-making is larger than the initial value R⁽⁰⁾, which is consistent with the conclusion in the theoretical derivation.

We use two key points on the sample cumulative probability distribution curve to obtain the two key parameters of the net load forecast error normal distribution model: the mean (μ) and the standard deviation (σ). If the sample mean ($\overline{x}$) and standard deviation (s) are directly taken as the parameters of the normal distribution model:

$$\overline{x}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N{x}_i}$$

$$s={\displaystyle \frac{1}{N-1}}\sqrt{{\displaystyle \sum _{i=1}^N{(x_i-\overline{x})}^2}}$$

The two obtained normal distribution models and the frequency distribution histogram of the net load forecast errors from Belgian transmission network, are plotted together, and the result is shown in Figure 5:

The normal distribution obtained by using the method proposed in this paper is denoted as Model 1. The normal distribution obtained by using the sample mean and standard deviation is denoted as Model 2. From the Figure 5, it can be observed that Model 1 fits the small probability interval of the tail much more accurately than Model 2. Furthermore, when the normal distribution parameters of Model 2 are directly substituted into the reserve decision-making iteration model, the obtained reserve decision result is 997.87 MW, indicating that the reserve result obtained by Model 2 is riskier.

8.2. Case Test of Guangdong Power Grid in China

The second case study analysis is conducted by using the relevant data from Guangdong Province, China. All the parameter values are the same as those in Section 8.1.

Model 1 (the normal distribution model obtained by using the method proposed in this paper, μ=1210.27, σ=2753.80) takes 4 iterations and 0.4967 seconds, yielding a reserve result of 7135.77 MW. Model 2 (the normal distribution obtained by using the sample mean and standard deviation) takes 5 iterations and 0.3422 seconds, yielding a reserve result of 7758.21 MW. The two obtained normal distribution models and the frequency distribution histogram of net load forecasting errors from Guangdong Province, China, are plotted together, as shown in Figure 6.

From the Figure 6, it can be observed that Model 1 fits the small probability interval of the tail much more accurately than Model 2. Additionally, the reserve decision result obtained by Model 2 is more conservative.

8.3. Result Analysis

Based on the above, the probabilistic reserve decision-making method proposed in this paper has the following mathematical advantages:

1)

The normal distribution model obtained through fitting the key points provides a better fit for the small probability interval of the tail.

2)

Since the normal distribution is still used to fit the net load forecast error variable, compared to other parameter distributions, the method has a lower computational complexity and is easier to compute.

9. Conclusion

This paper proposes a probabilistic reserve decision-making method for high-penetration renewable energy power systems based on cumulative probability approximation. It comprehensively considers three factors: reserve cost, expected generation revenue from reserves, and expected loss due to load shedding. While using a normal distribution model to fit the net load forecast error variable, the method improves the fitting accuracy in the small probability interval of the tail, demonstrating good practical value.

It is worth noting that this paper does not consider the limitations of transmission line capacity in the power system. Moreover, with the increasing frequency of extreme weather in recent years, the power grid is prone to N-k line failures during such events. Therefore, considering transmission line capacity and N-k line failures in reserve decision making will be a direction for future research.