Power distribution network voltage out-of-limit probability prediction method considering load fluctuation and photovoltaic output

Abstract

The invention provides a power distribution network voltage out-of-limit probability prediction method considering load fluctuation and photovoltaic output, and relates to the technical field of power distribution network safety. The method comprises the following steps: constructing a time-varying Gumbel-Copula joint distribution function according to a Spearman rank correlation coefficient, obtained by calculation, of a load and photovoltaic output; based on a time-varying Gumbel-Copula joint distribution function, probability distribution of the load and probability distribution of the photovoltaic output are reconstructed through a Hermite polynomial, and parallel probabilistic power flow solving is carried out; and calculating the voltage out-of-limit probability of each node according to the result of the parallelization probabilistic power flow solution. According to the method, dynamic coupling characteristics of different time periods are considered, a time-varying Gumbel-Copula joint distribution function is introduced, probability distribution of loads and probability distribution of photovoltaic output are reconstructed through a Hermite polynomial, the influence of intermittent superposition of load fluctuation and photovoltaic output can be overcome, and real-time prediction of the voltage out-of-limit probability can be accurately completed in extreme weather.

Description

Technical Field

[0001] This invention relates to the field of distribution network safety technology, and in particular to a method for predicting the probability of voltage overruns in distribution networks that takes into account load fluctuations and photovoltaic output. Background Technology

[0002] With the large-scale integration of distributed photovoltaic (PV) power into the distribution network, the intermittency of its output, coupled with load fluctuations, leads to frequent voltage exceedances at nodes. In distribution network operation, there is a challenge in predicting the probability of voltage exceedances due to the dynamic coupling between load fluctuations and the uncertainty of PV output.

[0003] Probabilistic power flow calculation is one of the core technologies in distribution network analysis. Its core lies in establishing the mathematical relationship between random variables such as load and photovoltaic (PV) output and the power flow equations of the power grid, predicting the variation range of key parameters such as voltage through probabilistic statistical methods. Monte Carlo simulation, which obtains statistical results through large-scale random sampling, offers high accuracy but its computational efficiency is insufficient for online applications. Point estimation methods, which approximate the probability distribution using a small number of feature sampling points, significantly improve computational efficiency but require addressing the distortion problem in non-Gaussian distribution modeling. Another important foundation is the random variable correlation modeling technique. Copula functions, as connection functions, can describe the statistical dependence between load and PV output, while time series models are used to capture the temporal fluctuation patterns of random variables. These two well-known techniques together constitute the theoretical basis for the probabilistic analysis of PV integration into the distribution network.

[0004] The distribution network probabilistic power flow calculation method based on static Copula functions and point estimation is one of the technical solutions that is relatively close to this invention. The specific implementation process of this solution includes: first, collecting historical load and photovoltaic data, fitting the variation law of load fluctuations with a normal distribution, and modeling the weather-related characteristics of photovoltaic output with a Beta distribution; then, establishing the joint probability distribution of load and photovoltaic output through Copula functions, using a fixed correlation coefficient to describe the static correlation between the two; finally, applying the classic 2m+1 point estimation method for sampling, inputting the sampled combination into the power flow calculation module, and outputting the mathematical expectation and standard deviation of node voltage as risk assessment indicators.

[0005] The core module of this scheme comprises a three-tiered cascade structure: random variable modeling, correlation coupling based on a static Copula function, and probability calculation based on point estimation. However, the static Copula function relies on only a single fixed coefficient to describe the correlation, failing to reflect the time-varying characteristics of load during midday's strong sunlight period and the decoupling between photovoltaic high synchronization and nighttime conditions. The point estimation process uses only two sampling points, which is insufficiently adaptable to the multi-peak distribution caused by frequent abrupt changes in photovoltaic output under cloudy weather conditions, and the probability prediction bias significantly amplifies under extreme weather conditions. Furthermore, the timeliness problem of the Monte Carlo alternative remains unresolved; the prediction time for a 100-node network is still in the hundreds of seconds, failing to overcome the technical bottleneck of minute-level real-time early warning. While this scheme establishes the basic framework for correlation and probability calculation, the dynamic coupling mechanism and real-time bottleneck constitute key obstacles to the existing technology.

[0006] Specifically, existing technologies use static Copula functions with fixed coefficients to model the correlation between load and photovoltaic output. Because they ignore the dynamic coupling characteristics at different times (such as high synchronization during midday high-irradiance periods and decoupling during cloudy / rainy weather), voltage exceedance events under extreme conditions are systematically underestimated, resulting in a large deviation between predicted values ​​and actual risk levels, and a high rate of missed detections during critical midday periods. This technical deficiency stems from the lack of time-based adaptability of static models, which cannot capture the dynamic correlations during operation mode switching.

[0007] Traditional point estimation methods rely solely on the mean and variance to reconstruct the probability distribution. When photovoltaic power output exhibits a bimodal or multimodal distribution under cloudy weather conditions, low-order statistical moments cannot describe the asymmetric characteristics of sudden power drops, leading to increased calculation errors in the probability of voltage exceeding limits under extreme operating conditions. Essentially, this is a distortion of the mathematical tool's representation of the actual non-Gaussian distribution, directly affecting the reliability of the tail of the risk probability density function.

[0008] Monte Carlo simulations and unoptimized point estimation methods both require extensive serial computations, with a single prediction for a distribution network with hundreds of nodes taking over 100 seconds, which cannot match the minute-level online decision-making window of the power grid dispatching system. This timeliness bottleneck stems from the lack of parallel acceleration of core computing modules and optimization of sampling strategies, resulting in lagging voltage runaway risk prevention and control in high-penetration photovoltaic access scenarios. Summary of the Invention

[0009] To address the technical problems existing in the prior art, this invention provides a method for predicting the probability of voltage exceedance in distribution networks, taking into account load fluctuations and photovoltaic output, including:

[0010] By dynamically tracking the time-varying coupling characteristics of load and photovoltaic output through a time-varying Copula function parameter sliding window mechanism, the problem of time-related omissions caused by static modeling is solved, significantly reducing the prediction omission rate during the midday peak period. Based on the improved point estimation method of the third-order Hermite polynomial, a skewness and kurtosis calculation layer is introduced to reconstruct the multi-peak non-Gaussian distribution, effectively suppressing prediction errors under extreme conditions and achieving accurate quantification of tail risks. Finally, a joint probabilistic power flow analysis and acceleration framework is designed. Through matrix operation parallelization technology, the prediction time of hundreds of nodes is significantly reduced, opening up a technical channel for minute-level online updates, providing a high-precision and time-efficient voltage safety protection system for active distribution networks.

[0011] Specifically, the technical solution provided by this invention includes the following steps:

[0012] Step S1: Construct the time-varying Gumbel-Copula joint distribution function based on the Spearman rank correlation coefficient between load and photovoltaic output;

[0013] Step S2: Based on the time-varying Gumbel-Copula joint distribution function, the probability distributions of load and photovoltaic output are reconstructed using Hermite polynomials, and parallel probabilistic power flow solutions are performed to obtain the voltage correction amount of each node in the distribution network under each scenario.

[0014] Step S3: For each node:

[0015] In each scenario, the reference voltage is added to the corresponding voltage correction to obtain the actual voltage value. The proportion of the over-limit portion in all actual voltage values ​​is calculated to obtain the voltage over-limit probability of the corresponding node in the distribution network.

[0016] Preferably, step S1 specifically includes:

[0017] Step S11: Calculate the Spearman rank correlation coefficient between load and photovoltaic output using a sliding time window, as shown in the following formula:

[0018]

[0019] In the formula, τ(t) is the Spearman rank correlation coefficient at time t, i is the index of the sample pair within the sliding time window, and d i N represents the rank difference between sample pairs within the sliding time window. w This represents the number of sample pairs within the sliding time window; where each sample pair includes a pair of load data acquisitions and photovoltaic output data acquisitions at the same time.

[0020] Step S12: Update and calculate the Spearman rank correlation coefficient τ(t) according to a preset period, and convert it into the linear correlation coefficient ρ(t), as shown in the following formula:

[0021]

[0022] Step S13: Construct the time-varying Gumbel-Copula joint distribution function C(u,v; θ(t)) based on the linear correlation coefficient ρ(t), as shown in the following formula:

[0023]

[0024] In the formula, u is the CDF function value of the load, v is the PDF function value of the photovoltaic output, and θ(t) is a time-varying parameter at time t that can be dynamically adjusted for tail correlation.

[0025] Furthermore, it also includes the method for calculating the CDF function value;

[0026] Divide the entire day into K time periods. For each time period:

[0027] Construct the state transition probability matrix P of the load (k) The formula is as follows:

[0028] P (k) =[p ij (k) ], k = 1, 2, ..., K;

[0029] In the formula, i and j are the load state numbers, and p ij (k) Let k be the probability of transitioning from load state i to load state j during time period k.

[0030] For the state transition probability matrix P (k) The eigenvalue decomposition formula is as follows:

[0031] P (k) υ=λυ;

[0032] In the formula, υ is the eigenvector and λ is the eigenvalue;

[0033] Solve the following formula to obtain the left eigenvector:

[0034] ((P (k) ) T -I1)π T =0;

[0035] In the formula, T represents the transpose, I1 is the identity matrix of the corresponding dimension, and π is the steady-state distribution matrix, which corresponds to the eigenvector υ with eigenvalue λ = 1.

[0036] Normalize the left eigenvector to obtain the steady-state probability distribution;

[0037] Constructing CDF function F L(x,k), calculate the CDF function value u, as follows:

[0038]

[0039] In the formula, x is the normalized value of the load power, k is the time period, L represents the load, and s j p is the normalized discrete load state value. k (s j ) represents the load state s during time period k. j The probability value.

[0040] Preferably, the entire day is divided into K time periods, specifically including:

[0041] The load curve corresponding to the whole day is divided into K time periods by using the DTW algorithm, and each load data collection is assigned to the nearest cluster according to the DTW distance;

[0042] The SBA algorithm is used to solve for the median sequence of each time period to update the cluster center, and the iteration continues until convergence.

[0043] Furthermore, it also includes methods for calculating PDF function values;

[0044] The Beta-ARIMA hybrid photovoltaic power output model is constructed by including a base output layer and a fluctuation correction layer, including:

[0045] The basic output layer is constructed using the following formula:

[0046]

[0047]

[0048] In the formula, y represents the normalized photovoltaic output. Let represent the probability distribution of photovoltaic output under the Beta model, where α and β are the shape parameters of the Beta distribution, B(·) represents the Beta function, μ is the per-unit value of irradiance, and σ is the standard deviation of historical irradiance.

[0049] Construct the fluctuation correction layer using the following formula:

[0050] Δy k =φ1Δy k-1 +…+θ p Δy k-p +ε k +θ1ε k-1 +…+θ q ε k-q ;

[0051] In the formula, Δy k Let p be the difference in photovoltaic output during time period k, and p be the number of time periods before time period k, from φ1 to φ pThe autoregressive coefficients for each time period are θ1 to θ2. q ε represents the moving average coefficient for each time period. k To ε k-q This represents the cloud coverage rate for each time period.

[0052] Construct the probability density function of photovoltaic power output and calculate the PDF function value v, as shown in the following formula:

[0053]

[0054] In the formula, f PV (y,k) is the probability density function of photovoltaic power output during time period k.

[0055] Preferably, step S2 specifically includes:

[0056] Step S21: Inverse sample the Gumbel-Copula joint distribution function to obtain the load or photovoltaic output, as shown in the following formula:

[0057] X = F -1 (C(u,v;θ(t)));

[0058] In the formula, X represents the load or photovoltaic output, and F... -1 (·) indicates inverse sampling;

[0059] Standardize X, the formula is as follows:

[0060]

[0061] In the formula, ξ represents the standardized variable, μ X σ represents the mean of variable X. X The standard deviation of variable X;

[0062] Step S22: Reconstruct the probability distribution Z of variable X using Hermite polynomials, as shown in the following formula:

[0063]

[0064] In the formula, H0 to H3 are the coefficients of the Hermite polynomial, and γ X To measure the skewness of the distribution asymmetry of variable X, κ X Kurtosis is used to measure the sharpness of the distribution of variable X;

[0065] Step S23: Parallelize the probabilistic power flow solution;

[0066] The probability distribution Z is mapped back to the original variable space and converted into power flow input parameters, as shown in the following formula:

[0067]

[0068] In the formula, P Load μ is the power flow input parameter for the load. L Let σ be the mean load. L Z represents the standard deviation of the load. L Let P be the probability distribution of the load. PV For the power flow input parameters of photovoltaic output, μ PV σ is the average photovoltaic output. PV Z represents the standard deviation of photovoltaic power output. PV The probability distribution of photovoltaic power output;

[0069] The power flow equations are matrixed using the Newton-Raphson method through iteration, as shown in the following formula:

[0070] JΔV=ΔS=f(P Load ,P PV ,Z);

[0071] In the formula, J is the Jacobian matrix, ΔV is the voltage correction, ΔS is the power deviation, and f(·) denotes matrixing;

[0072] Top-level parallel solution M group Scene grouping:

[0073]

[0074]

[0075] In the formula, m is the sequence number of the scene group, m = 1, 2, ..., M group ; i is the node index, N is the number of nodes in each scene group, i = 1, 2, ..., N; ΔV (m) For the voltage correction group of the m-th scenario, ΔS (m) The power deviation group for the m-th scene. Let i be the voltage correction value for node i in a group of m scene elements. Let i be the power deviation of node i in a group of m scene groups.

[0076] Furthermore, the top-level parallel solution of M group Before grouping each scenario, the process also includes: performing LU decomposition on the Jacobian matrix J and pre-storing the decomposition results.

[0077] Furthermore, it also includes a dynamic correction method for coefficient H3;

[0078] The formula for calculating kurtosis correction is as follows:

[0079]

[0080] In the formula, Δκ XE[·] represents the kurtosis correction for the sharpness of the distribution of variable X, and E[·] represents the expectation.

[0081] The dynamic correction factor H3 is used to obtain the correction factor. The formula is as follows:

[0082]

[0083] With correction factor Replacement coefficient H3.

[0084] Preferably, step S3 specifically includes:

[0085] Step S31: Statistically calculate the actual voltage values ​​of each node in the M scenarios, and calculate the upper and lower limits of the node's voltage exceedance probability, using the following formula:

[0086]

[0087] M = M group ×N;

[0088]

[0089] In the formula, P up,i Let P be the probability of node i exceeding its upper bound. down,i Let I be the lower bound probability of node i; I(·) is an indicator function, which is 1 when the condition is true and 0 otherwise. V represents the actual voltage value of node i in the m-th scene. base The reference voltage;

[0090] Step S32: Calculate the voltage exceedance probability of the node, using the following formula:

[0091] P VL,i =P up,i +P down,i ;

[0092] In the formula, P VL,i Let be the voltage over-limit probability of node i in the distribution network.

[0093] Furthermore, step S3 is followed by:

[0094] The adaptive control command for reactive power regulation of the photovoltaic inverter is generated using the following formula:

[0095] Q set =K p ·P VL,i ;

[0096] In the formula, Q set K is the reactive power setpoint for the photovoltaic inverter. p P is the proportional gain coefficient. VL,iLet be the voltage over-limit probability of node i in the distribution network.

[0097] Furthermore, step S3 is followed by:

[0098] The adaptive control command for tap position adjustment of on-load tap-changing transformers is given by the following formula:

[0099]

[0100] In the formula, Δtap is the tap adjustment amount of the on-load tap-changing transformer, ΔV is the voltage correction amount, and V step This is for adjusting the gear position step size.

[0101] Compared with existing technologies, the technical solution provided by this invention considers the dynamic coupling characteristics at different times, introduces a time-varying Gumbel-Copula joint distribution function, and reconstructs the probability distributions of load and photovoltaic output through Hermite polynomials. This can overcome the influence of the intermittent superposition of load fluctuations and photovoltaic output, and accurately predict the probability of voltage over-limit in real time under extreme weather conditions. Attached Figure Description

[0102] Figure 1 This is a flowchart illustrating a distribution network voltage over-limit probability prediction method in one embodiment of the present invention. Detailed Implementation

[0103] The technical solutions provided by the present invention will be further described in detail below through embodiments and accompanying drawings.

[0104] Example 1

[0105] The flowchart of the distribution network voltage over-limit probability prediction method provided in this embodiment is as follows: Figure 1 As shown, data is collected through the distribution network measurement system, and then preprocessed through data cleaning and alignment to achieve the fusion of multi-source heterogeneous data, laying the foundation for subsequent modeling. Intelligent measurement devices, such as PMUs (Phasor Measurement Units), deployed at various nodes of the distribution network collect time-series data such as load and power (including industrial, commercial, and residential load types) and key parameters of photovoltaic output (irradiance, module temperature, cloud cover). Simultaneously, local meteorological data for the next 72 hours is provided by the weather forecasting system. The data cleaning module uses the sliding quartile method to remove outliers and uses cubic spline interpolation to fill in missing data points. The preprocessing process ensures the integrity and continuity of the input data, forming a time-aligned dual-track data stream of load and photovoltaic output (15-minute time resolution), providing high-quality data support for time-segmentation modeling.

[0106] Specifically, the distribution network voltage over-limit probability prediction method provided in this embodiment includes the following four parts.

[0107] I. Principles of Dynamic Spatiotemporal Coupling Modeling

[0108] A piecewise Markov chain is used to describe the intraday load variation. The day is divided into K time periods, and a state transition probability matrix P is defined for each time period. (k) The formula is as follows:

[0109] P (k) =[p ij (k) ], k = 1, 2, ..., K;

[0110] In the formula, p ij (k) Let be the probability of transitioning from load state i to state j within time period k. The state transition matrix is ​​decomposed using eigenvalues, as shown in the following formula:

[0111] P (k) υ=λυ;

[0112] In the formula, υ is the eigenvector and λ is the eigenvalue.

[0113] The steady-state distribution matrix π corresponds to the eigenvector υ with eigenvalue λ = 1. The left eigenvector is calculated by solving the following equation, as shown in the formula below:

[0114] ((P (k) ) T -I1)π T =0;

[0115] In the formula, T represents transpose, and I1 is the identity matrix of the corresponding dimension.

[0116] The eigenvectors are normalized, and the steady-state probability distribution p is obtained. k (s j The cumulative distribution function (CDF) F, which represents the discrete probability distribution, is then output. L (x,k). In the piecewise Markov chain model, the load is treated as a discrete-continuous mixed variable, and its CDF is constructed piecewise over time. The CDF function value u is calculated using the following formula:

[0117]

[0118] In the formula, x is the normalized value (or per-unit value) of the load power, determined by the power reference value and the actual measured load power value, k is the time period, L represents the load, and s j The normalized discrete load state values, the steady-state probability distribution pk (s j The load is in state s during time period k, calculated using the Markov chain transition matrix. j The probability value.

[0119] This embodiment uses historical load curves as input, normalizes them to eliminate differences in absolute dimensions between days, and introduces DTW (Dynamic Time Warping) distance metric to improve the K-means clustering algorithm instead of Euclidean distance. The day is divided into six typical time periods based on load level (e.g., 00:00-06:00 trough, 06:00-10:00 morning peak, etc.; in some embodiments, the number of typical time periods can be redefined according to actual needs). Samples are first assigned to the nearest cluster according to DTW distance, and then the SBA (Simple Backoff Algorithm) algorithm is used to solve for the median sequence to update the cluster centers. This process is repeated multiple times until convergence. Time period boundaries are automatically adjusted through cross-validation of cluster centers to avoid subjective bias from manual division.

[0120] A Beta-ARIMA hybrid photovoltaic power output model was constructed by combining meteorological data such as irradiance, ambient temperature, and cloud cover, which includes a base output layer and a fluctuation correction layer.

[0121] The basic output layer is constructed as follows:

[0122]

[0123] In the formula, y represents the normalized photovoltaic output. Let represent the photovoltaic output probability distribution under the Beta model, where α and β are the shape parameters of the Beta distribution, and B(·) represents the Beta function.

[0124] A conversion model between irradiance and output is established using the Beta distribution, as shown in the following formula:

[0125]

[0126] In the formula, μ is the per-unit value of irradiance, determined by the actual irradiance and the standard irradiance, and σ is the standard deviation of historical irradiance, used to reflect the stability of local weather. Both are updated continuously from meteorological data. The conversion model between irradiance and photovoltaic output also includes parameters such as the temperature efficiency coefficient and the solar incidence angle. The base output layer mainly describes the impact of weather trends on photovoltaic output.

[0127] The volatility correction layer describes short-term volatility using an Autoregressive Integrated Moving Average (ARIMA) model, as shown in the following formula:

[0128] Δy k =φ1Δy k-1 +…+θ p Δy k-p +ε k +θ1ε k-1 +…+θ q ε k-q ;

[0129] In the formula, Δy k Let p be the difference in photovoltaic output during time period k, and p be the number of time periods before time period k, from φ1 to φ p The autoregressive coefficients for each time period are θ1 to θ2. q ε represents the moving average coefficient for each time period. k To ε k-q This represents the cloud coverage rate for each time period.

[0130] The final output is the probability density function. The PDF (Probability Density Function) value v is calculated using the following formula:

[0131]

[0132] In the formula, f PV (y,k) is the probability density function of photovoltaic power output during time period k.

[0133] II. Dynamic Coupling Mechanism of Time-Varying Copula Functions

[0134] The Spearman rank correlation coefficient between load and photovoltaic output is calculated using a sliding time window, as shown in the following formula:

[0135]

[0136] In the formula, τ(t) is the Spearman rank correlation coefficient at time t, i is the index of the sample pair within the sliding time window, and d i N represents the rank difference between sample pairs within the sliding time window. w This represents the number of sample pairs within the sliding time window; where each sample pair includes a pair of load acquisitions and photovoltaic output acquisitions at the same time.

[0137] Calculate the rank correlation coefficient τ(t) of the data within the window according to a preset period (usually 15 minutes), and convert it into the linear correlation coefficient ρ(t), as shown in the following formula:

[0138]

[0139] The Spearman rank correlation coefficient τ(t) reflects nonlinear correlation, while the linear correlation coefficient ρ(t) provides the linear correlation parameter for the Copula function.

[0140] The time-varying Gumbel-Copula joint distribution function C(u,v; θ(t)) is constructed based on the linear correlation coefficient ρ(t), as shown in the following formula:

[0141]

[0142]

[0143] In the formula, u is the CDF function value of the load, i.e., u = F L (x,t), where v is the PDF function value of the photovoltaic output, i.e., v = f PV (y,t), θ(t) is a time-varying parameter at time t that can be dynamically adjusted to tail correlation and accurately capture the strong coupling effect during the midday high irradiance period.

[0144] III. Improved Probabilistic Power Flow Calculation Framework

[0145] The load or photovoltaic output is obtained by inverse sampling of the Gumbel-Copula joint distribution function, as shown in the following formula:

[0146] X = F -1 (C(u,v;θ(t)));

[0147] In the formula, X represents the load or photovoltaic output, and F... -1 (·) indicates inverse sampling;

[0148] The random load or photovoltaic output is standardized using the following formula:

[0149]

[0150] In the formula, ξ represents the standardized variable, μ X σ represents the mean of variable X. X This represents the standard deviation of variable X; where variable X is the load or photovoltaic output.

[0151] The formula is derived from the Hermite polynomial reconfiguration of the load and the photovoltaic output probability distribution Z, as follows:

[0152]

[0153] In the formula, H0 to H3 are the coefficients of the Hermite polynomial. The probability distribution Z refers to the marginal distribution of the individual random variables of load and photovoltaic output. The reconstruction of the marginal distribution is a necessary prerequisite for solving the non-Gaussian property.

[0154] The coefficients in the above formula are calculated as follows:

[0155]

[0156] In the formula, H0 to H3 are the coefficients of the Hermite polynomial, and γ X To measure the skewness of the distribution asymmetry of variable X, κ X Kurtosis is used to measure the sharpness of the distribution of variable X.

[0157] By innovatively introducing two higher-order statistical moments, skewness compensation coefficient H2 and kurtosis compensation coefficient H3, and using the kurtosis correction factor Δκ... X Dynamic correction factor H3.

[0158] The formula for calculating kurtosis correction is as follows:

[0159]

[0160] In the formula, Δκ X E[·] represents the kurtosis correction for the sharpness of the distribution of variable X, and E[·] represents the expectation.

[0161] The dynamic correction factor H3 is used to obtain the correction factor. The formula is as follows:

[0162]

[0163] With correction factor The replacement coefficient H3 accurately fits the bimodal distribution of photovoltaic power output under cloudy weather conditions.

[0164] Next, parallel probabilistic power flow solution is performed:

[0165] The probability distribution Z is mapped back to the original variable space and converted into power flow input parameters, as shown in the following formula:

[0166]

[0167] In the formula, P Load μ is the power flow input parameter for the load. L Let σ be the mean load. L Z represents the standard deviation of the load. L Let P be the probability distribution of the load. PV For the power flow input parameters of photovoltaic output, μ PV σ is the average photovoltaic output. PV Z represents the standard deviation of photovoltaic power output. PV The probability distribution of photovoltaic power output;

[0168] The power flow equations are matrixed using the Newton-Raphson method through iteration. The iterative process is described below:

[0169] JΔV=ΔS=f(P Load ,P PV,Z);

[0170] In the formula, J is the Jacobian matrix, ΔV is the voltage correction, ΔS is the power deviation, and f(·) denotes matrixing.

[0171] This embodiment implements the parallel acceleration strategy through the following two steps:

[0172] 1. Pre-store the LU decomposition results of the Jacobian matrix J in GPU memory;

[0173] 2. Top-level view of M group Solve each scenario in groups in parallel:

[0174]

[0175] In the formula, m is the sequence number of the scene group, m = 1, 2, ..., M group ; i is the node index, N is the number of nodes in each scene group, i = 1, 2, ..., N; ΔV (m) For the voltage correction group of the m-th scenario, ΔS (m) The power deviation group for the m-th scene. Let i be the voltage correction value for node i in a group of m scene elements. Let i be the power deviation of node i in a group of m scene groups.

[0176] The middle layer allocates several (N) CPU computing nodes to each scene group. The bottom-level nodes initiate GPU parallel computing power flow equations (single instruction stream, multiple data stream). Parallelization of matrix operations reduces the computational complexity from O(MN). 3 ) decreased to O(N 3 )).

[0177] IV. Voltage Exceedance Probability Aggregation

[0178] The node voltage V at node i (i = 1, 2, ..., N) i The actual voltage values ​​of M scenarios are statistically analyzed, and the upper and lower limits of the node are calculated using the following formula:

[0179]

[0180] M = M group ×N;

[0181]

[0182] In the formula, P up,i Let P be the probability of node i exceeding its upper bound. down,i Let I be the lower bound probability of node i; I(·) is an indicator function, which is 1 when the condition is true and 0 otherwise. V represents the actual voltage value of node i in the m-th scene. base The reference voltage is 0.9 pu and 1.1 pu. In this embodiment, 0.9 pu and 1.1 pu are used as the upper and lower limits of the node voltage, which can be adjusted according to the relevant IEEE standards based on the operating environment.

[0183] The total probability of exceeding the limit is calculated using the following formula:

[0184] P VL,i =P up,i +P down,i ;

[0185] In the formula, P VL,i Let be the total over-limit probability of node i, that is, the voltage over-limit probability of node i in the distribution network, to realize the prediction of the voltage over-limit probability of the distribution network.

[0186] The adaptive control commands for the relevant devices are generated by outputting the total out-of-limit probability.

[0187] 1. Reactive power regulation of photovoltaic inverters, the formula is as follows:

[0188]

[0189] In the formula, Q set K is the reactive power setpoint for the photovoltaic inverter. p This is the proportional gain coefficient.

[0190] 2. The formula for adjusting the tap position of an on-load tap-changing transformer is as follows:

[0191]

[0192] In the formula, Δtap is the tap adjustment amount of the on-load tap-changing transformer, ΔV is the voltage correction amount, and V step This is the step size for gear adjustment.

[0193] As can be seen from the above embodiments and accompanying drawings, in this invention:

[0194] (1) The invention of a time-varying Copula function parameter sliding window mechanism using dynamic coupling technology allows for automatic updates of the load-photovoltaic output correlation coefficient via real-time data streams in each preset period, overcoming the limitation of traditional static models that cannot track changes in all-weather operation. Utilizing the rank correlation coefficient and the time-varying Gumbel-Copula joint distribution function, it accurately captures the strong midday correlation and nighttime decoupling effect. Specifically, the rank correlation coefficient is calculated in real-time through a sliding time window, driving adaptive updates of the Copula function parameters. This significantly reduces the time-specific omission rate and greatly improves the prediction accuracy for midday high-irradiance periods.

[0195] (2) Non-Gaussian Distribution Reconstruction: A third-order statistical moment compensation layer is embedded in the classical point estimation method. By using skewness and kurtosis as dual parameters to collaboratively correct the distribution structure of random variables, the problem of distorted modeling of multi-peak distribution of photovoltaic power output under cloudy weather conditions is overcome. An adaptive fitting method is designed to dynamically adjust the compensation intensity based on the data. A third-order Hermite polynomial compensator is constructed by introducing skewness and kurtosis to reconstruct the asymmetric multi-peak distribution. The prediction error of cloudy weather conditions is significantly reduced, and the reliability of quantifying the tail risk of voltage over-limit is improved.

[0196] (3) A three-level acceleration system for probabilistic power flow computation is constructed using a hierarchical parallel acceleration architecture: the top level allocates scenario groups according to risk sensitivity, the middle level schedules tasks across multiple CPU nodes, and the bottom level implements GPU parallelization based on pre-stored matrix templates. The design of the three-level parallel architecture, combined with the pre-stored reuse of Jacobian matrices, reduces the computation time to the minute level.

[0197] In summary, compared with the prior art, the technical solution provided by this invention takes into account the dynamic coupling characteristics of different time periods, introduces a time-varying Gumbel-Copula joint distribution function, and reconstructs the probability distributions of load and photovoltaic output through Hermite polynomials. This can overcome the influence of the intermittent superposition of load fluctuations and photovoltaic output, and accurately predict the probability of voltage over-limit in real time under extreme weather conditions.

[0198] Preferably, the rank correlation coefficient and the time-varying Gumbel-Copula joint distribution function can accurately capture the strong midday correlation and the decoupling effect at night, which helps to improve the prediction accuracy and precision. The skewness and kurtosis dual parameters work together to correct the distribution structure of random variables, which can overcome the technical problem of distortion in the modeling of multi-peak distribution of photovoltaic power output under cloudy weather, and help to improve the prediction accuracy and precision. The hierarchical parallel architecture helps to improve the prediction efficiency.

Claims

1. A method for probabilistic prediction of voltage excursion in a distribution network accounting for load fluctuations and photovoltaic output, characterized in that, The method comprises the following steps: Step S1: constructing a time-varying Gumbel-Copula joint distribution function according to the Spearman rank correlation coefficient of the load and the photovoltaic output; Step S2: based on the time-varying Gumbel-Copula joint distribution function, reconstructing the probability distribution of the load and the photovoltaic output respectively through Hermite polynomials, and performing parallelized probability power flow calculation to obtain the voltage correction of each node in the distribution network under each scenario; Step S3: for each node: In each scenario, the reference voltage is added to the corresponding voltage correction to obtain the actual voltage value, the proportion of the number of out-of-limit parts in all actual voltage values is calculated to obtain the voltage out-of-limit probability of the corresponding node in the distribution network.

2. The power distribution network voltage excursion probability forecasting method of claim 1, wherein, Step S1 specifically comprises: Step S11: calculating the Spearman rank correlation coefficient of the load and the photovoltaic output through a sliding time window, and the formula is as follows: In the formula, τ(t) is the Spearman rank correlation coefficient at time t, i is the serial number of the sample pair in the sliding time window, d i is the rank difference of the sample pair in the sliding time window, N w is the number of sample pairs in the sliding time window; wherein the sample pair comprises a pair of load collection quantity and photovoltaic output collection quantity at the same time; Step S12: updating and calculating the Spearman rank correlation coefficient τ(t) every period, and converting it into a linear correlation coefficient ρ(t), and the formula is as follows: Step S13: constructing a time-varying Gumbel-Copula joint distribution function C(u,v; θ(t)) according to the linear correlation coefficient ρ(t), and the formula is as follows: In the formula, u is the CDF function value of the load, v is the PDF function value of the photovoltaic output, and θ(t) is the time-varying parameter at time t, which can dynamically adjust the tail correlation.

3. The power distribution network voltage excursion probability forecasting method of claim 1, wherein, Further comprising a calculation method of the CDF function value; Divide the whole day into K time periods, and for each time period: Constructing the state transition probability matrix P of the load (k) The formula is as follows: P (k) = [p ij (k) ], k = 1, 2,..., K; where i and j are the indices of the load states, p ij (k) Pij(k) is the probability of moving from load state i to load state j for time period k; The state transition probability matrix P (k) Eigenvalue decomposition is performed, and the formula is as follows: P (k) υ = ly; In the formula, υ is the eigenvector, and λ is the eigenvalue; Solve the following formula to obtain the left eigenvector: ((P (k) ) T -I1)π T =0; In the formula, T represents transposition, I1 is a unit matrix corresponding to the dimension; π is a steady-state distribution matrix corresponding to the eigenvector υ with the eigenvalue λ = 1; Normalize the left eigenvector to obtain the steady-state probability distribution; Constructing the CDF function F L (x, k), compute the CDF function value u, as follows: where x is a normalized value of the load power, k is a time period, L represents the load, s j is a normalized discrete load state value, p k (s j ) is a probability value of the load being in state s j in the time period k.

4. The power distribution network voltage excursion probability forecasting method of claim 3, wherein, The whole day is divided into K time periods, specifically comprising: The corresponding load curve of the whole day is divided into K time periods through the DTW algorithm, and each load collection quantity is distributed to the nearest cluster according to the DTW distance; The median sequence of each time period is solved through the SBA algorithm to update the cluster center, and iteration is performed until convergence.

5. The power distribution network voltage excursion probability forecasting method of claim 2, wherein, Further comprising a calculation method of the PDF function value; The Beta-ARIMA hybrid photovoltaic output model comprises a basic output layer and a fluctuation correction layer, comprising: The basic output layer is constructed, and the formula is as follows: where y is the normalized photovoltaic output, is the photovoltaic output probability distribution under the Beta model, a and β are the Beta distribution shape parameters, B(·) denotes the Beta function, μ is the unit value of irradiance, and σ is the standard deviation of historical irradiance; The fluctuation correction layer is constructed, and the formula is as follows: Δy k = φ1Δy k-1 + φ2Δy p + … + φnΔy k-p + ε k = θ1ε k-1 + θ2ε q + … + θmε k-q ; where Δy k is the difference value of photovoltaic output at time period k, p is the number of time periods before time period k, φ1 to φ p are the autoregressive coefficients corresponding to each time period, θ1 to θ q are the moving average coefficients corresponding to each time period, and ε k to ε k-q are the cloud cover rates corresponding to each time period. The probability density function of the photovoltaic output is constructed, and the PDF function value v is calculated, and the formula is as follows: where f PV (y, k) is the probability density function of the photovoltaic power output at time period k.

6. The method for probabilistic prediction of power distribution network voltage excursion according to claim 2, wherein, Step S2 specifically comprises: Step S21: inverse sampling of the Gumbel-Copula joint distribution function to obtain the load or photovoltaic output, and the formula is as follows: X = F -1 (C(u, v; θ(t))); where X is the load or photovoltaic power output, F -1 (·) denotes inverse sampling; Standardize X, and the formula is as follows: where ξ denotes a standardized variable, μ X denotes the mean of the variable X, σ X denotes the standard deviation of the variable X; Step S22: reconstruct the probability distribution Z of the variable X through Hermite polynomials, and the formula is as follows: where H0to H3are coefficients of Hermite polynomials, γ X is skewness, which measures the asymmetry of the distribution of the variable X X is kurtosis, which measures the sharpness of the distribution of the variable X Step S23: parallelized probability power flow calculation; Map the probability distribution Z back to the original variable space and convert it into a power flow input parameter, and the formula is as follows: where P Load is the load's tidal input parameter, μ L is the load's mean, σ L is the load's standard deviation, Z L is the load's probability distribution, P PV is the photovoltaic output's tidal input parameter, μ PV is the photovoltaic output's mean, σ PV is the photovoltaic output's standard deviation, Z PV is the photovoltaic output's probability distribution; The Newton-Raphson method is used to matrix the power flow equation through iteration, and the formula is as follows: JAV = AS = f(P Load ,P PV ,Z). In the formula, J is the Jacobian matrix, ΔV is the voltage correction, ΔS is the power deviation, and f(·) represents matrixing. Top-level parallel solution of M group a scene group: wherein m is the sequence number of the scene group, m = 1, 2,..., M group ; i is the node sequence number, N is the number of nodes in each scene group, i = 1, 2,..., N; AV (m) is the voltage correction amount group of the mth scene group, AS (m) is the power deviation amount group of the mth scene group, AV i (m) is the voltage correction amount of node i in the m scene groups, AS i (m) is the power deviation amount of node i in the m scene groups.

7. The power distribution network voltage excursion probability forecasting method of claim 6, wherein, The top layer solves M group Before solving the scene group, the method further comprises: performing LU decomposition on the Jacobian matrix J and pre-storing the decomposition result.

8. The power distribution network voltage excursion probability forecasting method of claim 6, wherein, The dynamic correction method of the coefficient H3 is also included; The kurtosis correction amount is calculated, and a formula is as follows: where Δκ X is the kurtosis correction for the sharpness of the distribution of the variable X, and E[·] denotes the expectation. The dynamic correction coefficient H3 is obtained by correcting the coefficient H2 The formula is as follows: with a correction factor replacement coefficient H3.

9. The power distribution network voltage excursion probability forecasting method of claim 6, wherein, The step S3 specifically includes: In the step S31, the actual voltage values of the nodes in the M scenes are counted, and the over-limit probability and the under-limit probability of the nodes are calculated, and formulas are as follows: where P up,i is the over-limit probability of node i, P down,i is the under-limit probability of node i; I(·) is an indicator function, which is 1 if the condition is true, otherwise 0; V i (m) is the actual voltage value of node i in the mth scenario, V base is the reference voltage; In the step S32, the voltage over-limit probability of the nodes is calculated, and a formula is as follows: P VL,i = P up,i + P down,i ; In the formula, P VL,i is the probability of voltage excursion of node i in the distribution network.

10. The method for probabilistic prediction of power distribution network voltage excursion according to claim 1, wherein, After the step S3, the following steps are further included: The adaptive control instruction of the reactive power regulation of the photovoltaic inverter is generated, and a formula is as follows: Q set = K p · P VL,i ; In the formula, Q set is the reactive power set value of the photovoltaic inverter, K p is the proportional gain coefficient, P VL,i is the voltage out-of-limit probability of node i in the power distribution network.

11. The method for probabilistic prediction of power distribution network voltage excursion according to claim 1, wherein, After the step S3, the following steps are further included: The adaptive control instruction of the tap adjustment of the on-load voltage regulating transformer is generated, and a formula is as follows: In the formula, Δtap is the tap adjustment amount of the on-load tap changer, ΔV is the voltage correction amount, V step is the tap adjustment step.

Images