Photovoltaic micro-grid voltage out-of-limit risk detection method and system

Abstract

The application discloses a photovoltaic micro-grid voltage out-of-limit risk detection method and system, relates to the field of voltage detection, and accurately simulates the random volatility and time correlation of photovoltaic output by constructing an Ornstein-Uhlenbeck random process model to generate photovoltaic output sequences and Brown motion increments in line with actual physical characteristics. Based on this, the system voltage response is obtained by combining nonlinear power flow simulation, and the Malliavin analysis theory is introduced to generate a Malliavin weight sample set by constructing a Malliavin derivative matrix and covariance analysis, and then the weight set is used for voltage probability density estimation and out-of-limit risk calculation. By introducing the Malliavin weight factor, the sensitivity transfer relationship from random source disturbance to system voltage response is established, and the weight correction technology is used to realize high-precision approximation of the voltage probability distribution (especially the tail out-of-limit area) under a small sample size.

Description

Technical Field

[0001] This application relates to the field of voltage detection, and more specifically, to a method and system for detecting voltage over-limit risks in photovoltaic microgrids. Background Technology

[0002] Microgrids, as an important means of integrating distributed energy resources and improving power supply reliability, have become a key direction in the development of modern power systems. Photovoltaic power generation, as one of the main forms of distributed energy, exhibits significant intermittency and randomness in its output, posing a severe challenge to the safe and stable operation of microgrids. Voltage stability is one of the core indicators for the safe operation of microgrids; drastic fluctuations in photovoltaic output can easily lead to node voltage exceeding limits, causing power quality problems and even threatening system security. Therefore, constructing a scheme that can accurately describe the random characteristics of photovoltaic power and efficiently detect voltage exceeding risk is of significant practical importance for ensuring the reliable operation of microgrids.

[0003] Currently, existing technologies for assessing and detecting voltage exceedance risks mainly employ Monte Carlo simulation, point estimation, and semi-invariant methods. While Monte Carlo simulation is simple in principle and widely applicable, it suffers from inherent drawbacks such as slow convergence and low computational efficiency. Its estimation error is inversely proportional to the square root of the sample size, often requiring extremely large-scale sampling for low-probability voltage exceedance events, leading to a heavy computational burden and failing to meet the timeliness requirements of practical engineering. Point estimation approximates the system response by calculating at specific probability points, improving computational efficiency, but it suffers from significant errors when dealing with non-Gaussian and strongly nonlinear systems, making it difficult to accurately describe the complex stochastic characteristics of photovoltaic power output and the nonlinear features of the system, especially lacking accuracy in assessing low-probability events at the tail of the distribution. Furthermore, while semi-invariant methods can be used for probabilistic power flow calculations, their approximation effect is poor when dealing with non-Gaussian and multimodal distributions, and their complex derivation process limits their widespread application in engineering. Existing methods struggle to accurately characterize photovoltaic stochastic processes (such as the Ornstein-Uhlenbeck process) and system nonlinear responses while ensuring computational efficiency, making it difficult to balance the accuracy and real-time performance of voltage over-limit risk detection. Summary of the Invention

[0004] To address the problems in the prior art, according to one aspect of this application, a method for detecting voltage over-limit risks in photovoltaic microgrids is provided, comprising: The process involves acquiring microgrid network structure parameters, photovoltaic (PV) configuration parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters. The Ornstein-Uhlenbeck process parameters include PV output regression velocity parameters, long-term mean, and volatility parameters. The numerical calculation parameters include time step, number of PV nodes, total number of samples, and total simulation duration. Stochastic process modeling and sample sequence generation are performed on the microgrid network structure parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters to obtain PV output time series samples and Brownian motion increment sequences. Based on the microgrid network structure parameters and PV configuration parameters, nonlinear power flow simulation is performed on the PV output time series samples to obtain a voltage amplitude sample set. Based on the voltage amplitude sample set and Brownian motion increment sequences, a Malliavin weighted sample set is determined. Finally, voltage probability density estimation and over-limit risk calculation are performed on the Malliavin weighted sample set and voltage amplitude sample set to obtain the voltage over-limit risk probability.

[0005] According to another aspect of this application, a photovoltaic microgrid voltage over-limit risk detection system is provided, comprising: a parameter acquisition module for acquiring microgrid network structure parameters, photovoltaic configuration parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters, wherein the Ornstein-Uhlenbeck process parameters include photovoltaic output regression rate parameters, long-term mean and volatility parameters, and the numerical calculation parameters include time step, number of photovoltaic nodes, total number of samples, and total simulation duration; and a modeling and sample sequence generation module for processing microgrid network structure parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters. The system employs stochastic process modeling and sample sequence generation based on parameters and numerical calculation parameters to obtain photovoltaic power output time series samples and Brownian motion increment sequences. A simulation module performs nonlinear power flow simulation on the photovoltaic power output time series samples based on microgrid network structure parameters and photovoltaic configuration parameters to obtain a voltage amplitude sample set. A weighted sample set generation module determines the Malliavin weighted sample set based on the voltage amplitude sample set and Brownian motion increment sequences. A voltage exceedance risk probability generation module estimates the voltage probability density and calculates the exceedance risk from the Malliavin weighted sample set and voltage amplitude sample set to obtain the voltage exceedance risk probability.

[0006] Compared with existing technologies, this application provides a method and system for detecting voltage over-limit risks in photovoltaic microgrids. By constructing an Ornstein-Uhlenbeck stochastic process model, it accurately simulates the random fluctuations and time correlation of photovoltaic power output, generating a photovoltaic power output sequence and Brownian motion increment that conform to actual physical characteristics. Based on this, it combines nonlinear power flow simulation to obtain the system voltage response and introduces Malliavin analytical theory. By constructing the Malliavin derivative matrix and performing covariance analysis, it generates a Malliavin weighted sample set, which is then used for voltage probability density estimation and over-limit risk calculation. This effectively solves the problems of slow convergence and heavy computational burden in the assessment of low-probability events using existing Monte Carlo simulation methods, while also overcoming the shortcomings of point estimation methods in handling strong nonlinearity and non-Gaussian distributions. By introducing Malliavin weighting factors, a sensitivity transfer relationship from random source disturbances to the system voltage response is established. Using weight correction techniques, high-precision approximation of the voltage probability distribution (especially the tail over-limit region) can be achieved with a small sample size, thereby significantly improving the computational efficiency of photovoltaic microgrid voltage over-limit risk detection while ensuring assessment accuracy. Attached Figure Description

[0007] The above and other objects, features and advantages of this application will become more apparent from the more detailed description of the embodiments of this application in conjunction with the accompanying drawings.

[0008] Figure 1 This is a flowchart of a photovoltaic microgrid voltage over-limit risk detection method according to an embodiment of this application.

[0009] Figure 2 This is a schematic diagram of data flow in the photovoltaic microgrid voltage over-limit risk detection method according to an embodiment of this application.

[0010] Figure 3 This is a flowchart of step S2 in the photovoltaic microgrid voltage over-limit risk detection method according to an embodiment of this application.

[0011] Figure 4 This is a block diagram of a photovoltaic microgrid voltage over-limit risk detection system according to an embodiment of this application. Detailed Implementation

[0012] Embodiments of this disclosure will now be described in more detail with reference to the accompanying drawings. It should be understood that the drawings and embodiments of this disclosure are for illustrative purposes only and are not intended to limit the scope of protection of this disclosure.

[0013] To address the technical challenges mentioned above, this application proposes a method for detecting voltage over-limit risks in photovoltaic microgrids. Figure 1This is a flowchart of a photovoltaic microgrid voltage over-limit risk detection method according to an embodiment of this application. Figure 2 This is a schematic diagram of the data flow in the photovoltaic microgrid voltage over-limit risk detection method according to an embodiment of this application. Figure 1 and Figure 2 As shown, the photovoltaic microgrid voltage over-limit risk detection method according to an embodiment of this application includes: S1, acquiring microgrid network structure parameters, photovoltaic configuration parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters, wherein the Ornstein-Uhlenbeck process parameters include photovoltaic output regression velocity parameters, long-term mean, and volatility parameters, and the numerical calculation parameters include time step, number of photovoltaic nodes, total number of samples, and total simulation duration; S2, performing stochastic process modeling and sample sequence generation on the microgrid network structure parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters to obtain photovoltaic output time series samples and Brownian motion increment sequences; S3, performing nonlinear power flow simulation on the photovoltaic output time series samples based on the microgrid network structure parameters and photovoltaic configuration parameters to obtain a voltage amplitude sample set; S4, determining a Malliavin weight sample set based on the voltage amplitude sample set and Brownian motion increment sequences; S5, performing voltage probability density estimation and over-limit risk calculation on the Malliavin weight sample set and voltage amplitude sample set to obtain the voltage over-limit risk probability.

[0014] In step S1, the microgrid network structure parameters, photovoltaic (PV) configuration parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters are obtained. The Ornstein-Uhlenbeck process parameters include PV output regression rate parameters, long-term mean, and volatility parameters. The numerical calculation parameters include time step, number of PV nodes, total number of samples, and total simulation duration. It should be understood that a microgrid is a complex physical network containing multi-source heterogeneous devices, and its operating state is influenced by the dual coupling of physical topology constraints and stochastic weather conditions. To accurately capture the voltage sensitivity characteristics to PV output using the Malliavin analytical method in subsequent analysis, it is first necessary to establish a mathematical description benchmark highly consistent with the actual physical scenario. A lack of accurate network architecture parameters will cause power flow calculations to deviate from actual physical constraints, while coarse PV stochastic process parameters cannot truly reflect the mean regression characteristics and spatiotemporal correlation of PV output, thus rendering the generated Malliavin weights meaningless. Therefore, this step aims to perform a full-dimensional parameterization and initialization of the microgrid's physical topology, energy configuration, and stochastic dynamic characteristics. This not only determined the boundary conditions and time resolution for the numerical solution of stochastic differential equations, but also provided the necessary data foundation for the subsequent construction of a high-precision probability density estimation model, ensuring the consistency and reliability of the voltage over-limit risk detection results at both the physical and statistical levels.

[0015] In an exemplary operation, the specific process of step S1 is as follows: In the initial stage of implementing the photovoltaic microgrid voltage over-limit risk detection method, step S1 is mainly responsible for constructing the parameter foundation for the entire numerical simulation and risk analysis. The specific implementation process of this step is a mapping process from physical entities to mathematical model parameters, mainly covering four core dimensions: microgrid network structure parameter acquisition, photovoltaic configuration parameter setting, Ornstein-Uhlenbeck (OU) process parameter identification, and numerical calculation parameter configuration. These parameters together constitute a structured initialization dataset, which serves as the direct input for subsequent stochastic process modeling and nonlinear power flow simulation.

[0016] First, regarding the acquisition of microgrid network structure parameters, this process aims to digitally reconstruct the microgrid's topological connections and electrical characteristics. During implementation, it is necessary to read the line impedance data (resistance R and reactance X) and node connection information of the target microgrid, thereby constructing the system node admittance matrix. This matrix is ​​a complex sparse matrix, and its dimension depends on the total number of system nodes. Simultaneously, it is necessary to obtain the baseline active and reactive power distribution data for each load node, which is determined based on historical load curves or typical daily load forecast data. Based on this, a key implementation step is to construct the implicit spatial correlation coefficient matrix between photovoltaic (PV) nodes. Since geographically proximate PV power plants are affected synchronously by meteorological factors such as cloud cover, this spatial correlation must be quantified. Specifically, by collecting historical output time-series data of each PV power plant within the target area, the correlation between each pair of nodes is calculated using the Pearson correlation coefficient formula, thus forming a K×K dimensional symmetric matrix (K being the number of PV nodes). For example, in a microgrid with three photovoltaic grid-connected points, if node A and node B are close to each other, the correlation coefficient calculated from historical data may be 0.85, while the correlation coefficient between nodes A and C, which are far apart, may be only 0.3. The resulting matrix will accurately reflect this spatial coupling characteristic and provide a mathematical basis for the subsequent generation of correlated Brownian motion increments.

[0017] Secondly, obtaining the photovoltaic (PV) configuration parameters mainly involves clarifying the physical connection points of the PV power source in the microgrid and its rated power limit. During implementation, the node numbers where PV power generation devices are installed need to be specified in the system topology, and the rated capacity (in kW or MW) of each connection point needs to be defined. These parameters determine how the normalized values ​​generated by subsequent stochastic processes are mapped to the actual active power values ​​injected into the grid. For example, this microgrid contains three PV grid-connected points. Node 12 is configured to connect a 200kW PV system, node 25 a 150kW PV system, and node 5 a 30kW PV system. This configuration information will be stored as a PV location index vector and a capacity vector.

[0018] Next, the parameters of the Ornstein-Uhlenbeck process were obtained and modeled. The OU process was chosen to simulate the stochastic characteristics of photovoltaic power output because it exhibits mean reversion, effectively describing the physical phenomenon of photovoltaic power output fluctuating around a solar irradiance baseline, rather than diverging like ordinary Brownian motion. This process is described by stochastic differential equations, and its core parameters include the regression velocity parameter. Long-term mean and volatility parameters In practice, these parameters are obtained through statistical analysis of high-resolution historical photovoltaic irradiance or power output data for the region. Maximum likelihood estimation (MLE) is used to fit the historical time series: regression velocity parameters... This reflects how quickly photovoltaic output recovers to its long-term average after being disturbed (such as by cloud cover); the larger the value, the faster the recovery. (Long-term average) It represents the average output level over a specific time period; volatility parameter This quantifies the severity of random fluctuations. For example, to determine the photovoltaic power output characteristics under typical cloudy weather conditions, a regression velocity parameter is determined after fitting historical data. The normalized long-term mean is 5.0. The volatility parameter is 0.6. The value is 0.15. These parameters will directly determine the dynamic trajectory of the subsequently generated stochastic differential equations.

[0019] Finally, the numerical computation parameters are configured to establish the spatiotemporal resolution of the simulation. This includes setting the time step. Total simulation duration Number of photovoltaic nodes and the total number of samples in the Monte Carlo simulation The choice of time step needs to meet the stability requirements of discretized solution of stochastic differential equations, and should be set to a small value, such as 0.01 hours or less, to capture rapid photovoltaic fluctuations. The total simulation duration is set to 24 hours to cover a complete daily operating cycle. The total number of samples determines the balance between the accuracy of probability density estimation and computational cost. Although the Malliavin method has a variance reduction effect, a certain sample size (such as 1000 to 5000) is still required to ensure statistical properties. For example, a time step of 0.01 h, a simulation duration of 24 h (i.e., 2400 time steps), 3 photovoltaic nodes, and a total sample size of 2000 can be set.

[0020] In step S2, stochastic process modeling and sample sequence generation are performed on the microgrid network structure parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters to obtain photovoltaic power output time series samples and Brownian motion increment sequences. Correspondingly, the output power of photovoltaic power generation is neither constant nor completely random fluctuation, but exhibits significant mean regression characteristics and time correlation. That is, after being disturbed by meteorological factors such as cloud cover, photovoltaic output tends to recover to a certain long-term average level. If only static probability models or simple normal distribution random number generation methods are used, the physical process of the dynamic evolution of photovoltaic output over time cannot be truly reproduced, resulting in a lack of continuity constraints in the time dimension in the generated samples. Consequently, the voltage over-limit risk calculated subsequently deviates significantly from the actual operating situation. Therefore, this step aims to establish a stochastic process model that can accurately describe this dynamic characteristic, transforming the static parameters set in the previous steps into a physically meaningful dynamic time series, thereby ensuring that the risk detection results can capture the dynamic impact of photovoltaic fluctuations on the voltage stability of the microgrid.

[0021] In one exemplary operation Figure 3 This is a flowchart of step S2 in the photovoltaic microgrid voltage over-limit risk detection method according to an embodiment of this application. Figure 3 As shown, step S2 involves performing stochastic process modeling and sample sequence generation on the microgrid network structure parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters to obtain photovoltaic power output time series samples and Brownian motion increment sequences. This includes: S21, determining the simulation initialization state set based on the numerical calculation parameters and Ornstein-Uhlenbeck process parameters, wherein the simulation initialization state set includes the initial power output vector and discrete event grid sequence; S22, performing correlation decomposition and correlation transformation on the spatial correlation coefficient matrix between implicit photovoltaic nodes in the microgrid network structure parameters based on the numerical calculation parameters to obtain the Brownian motion increment sequence; and S23, performing explicit Euler iteration and power sequence synthesis on the simulation initialization state set based on the Ornstein-Uhlenbeck process parameters and the Brownian motion increment sequence to obtain photovoltaic power output time series samples.

[0022] In the above exemplary operation, the specific process of step S2 is as follows: Sub-step S21 is mainly responsible for constructing the mathematical framework of the stochastic differential equation and completing the time domain and state domain initialization of the numerical simulation. This process is to transform the Ornstein-Uhlenbeck (OU) process parameters and numerical calculation parameters determined in step S1 into specific mathematical model instances and discretized computational grids. The generated simulation initialization state set will serve as the starting point for subsequent iterative solutions. First, the equation definition based on the Ornstein-Uhlenbeck process parameters is implemented. Based on the photovoltaic power output regression velocity parameters obtained in step S1... Long-term mean and volatility parameters A continuous form of a stochastic differential equation describing the dynamic characteristics of photovoltaic power output is constructed. This equation is expressed as: In this mathematical model, The photovoltaic normalized output state variable at time t is calculated using a discrete iterative formula based on the Ornstein-Uhlenbeck stochastic process model, utilizing the output value, regression velocity, long-term mean, and volatility parameters from the previous time step, combined with the Brownian motion increment. The term represents the drift, describing how photovoltaic output shifts from its mean back to its long-term mean. The deterministic trend of regression, including the regression velocity parameter This determines the rate of response, reflecting the physical inertia of photovoltaic power output recovering to stability after being disturbed by meteorological disturbances; The diffusion term represents the effect of random fluctuations, where... For volatility parameters, The increments of the standard Brownian motion (Wiener process) follow a mean of 0 and a variance of . The normal distribution is established. By establishing this equation, a clear dynamic mechanism is given to the photovoltaic power output sample. For example, using the parameters in step S1 =5.0、 =0.6、 =0.15, this equation determines that the output trajectory will fluctuate around the central value of 0.6, and the intensity of the fluctuation is controlled by 0.15. Next, time-domain discretization is performed. Since computers cannot directly process continuous-time differential equations, it is necessary to rely on the total simulation time in the numerical calculation parameters. With calculation step size Construct a discrete time architecture. This process begins by calculating the total number of time steps. The calculation formula is: Subsequently, a monotonically increasing discrete-time grid sequence is generated. ,in =0, at any time Continuing from the previous embodiment, the total simulation duration is set. =24 hours, time step =0.01 hours, then the total number of steps is calculated. =2400. The generated discrete event grid sequence is a vector containing 2401 time points [0.00, 0.01, 0.02, ..., 24.00]. This sequence defines the time base for subsequent Brownian motion increment generation and Euler iterative solution, ensuring that the output state of all photovoltaic nodes evolves synchronously on the same time axis. Next, state initialization is performed to determine the initial output vector. To start the numerical iteration process, the initial state of all photovoltaic nodes at time t=0 needs to be set. Based on the long-term mean parameters of the Ornstein-Uhlenbeck process... This initial value is set so that the system is in statistical equilibrium at the initial moment. For systems containing... For a microgrid with 3 photovoltaic nodes, an initial output vector of dimension K×1 is constructed. If in step S1 the number of photovoltaic nodes K=3, and the long-term mean of each node is... =0.6, then the generated initial output vector for This vector represents the normalized power output value of all photovoltaic access points at the instant the simulation begins. It will serve as the first known state variable in the iterative formula for subsequent calculations. Finally, the simulation initialization state set is encapsulated. The generated initial power vector is combined with the discrete event grid sequence and packaged to form the simulation initialization state set. This set clarifies the starting point and path scale of the stochastic process simulation.

[0023] Sub-step S22 is crucial, serving as a bridge between the preceding and following steps. Its core task is to transform the spatial correlation parameters determined in step S1 into a numerical sequence that can be used to solve stochastic differential equations. In an exemplary operation, step S22, based on the numerical calculation parameters, performs correlation decomposition and correlation transformation on the spatial correlation coefficient matrix between implicit photovoltaic nodes in the microgrid network structure parameters to obtain a Brownian motion increment sequence, including: S221, generating an independent random number matrix that follows a standard normal distribution, the dimension of which is... ,in, S222: Perform correlation decomposition on the spatial correlation coefficient matrix between the implicit photovoltaic nodes to obtain the lower triangular matrix; S223: Perform linear transformation on the independent random number matrix based on the lower triangular matrix to obtain the Brownian motion increment sequence.

[0024] First, step S221 is implemented. The purpose of this step is to construct the original driving force of the stochastic process, namely, a white noise sequence generated purely by probabilistic statistical laws without any spatial correlation information. During implementation, the number of photovoltaic nodes is read based on the numerical calculation parameters determined in steps S1 and S21. and the total number of discrete time steps In this embodiment, the microgrid contains K=3 photovoltaic nodes, labeled as node 12, node 25, and node 5, respectively. The total simulation duration is 24 hours, and the time step is 0.01 hours, resulting in a total of N=2400 time steps. Based on these two dimensional parameters, a standard normal distribution random number generator (e.g., using the Box-Muller transform algorithm or the Ziggurat algorithm) is invoked to generate a matrix with dimensions [3, 2400]. In this matrix, each element (Where i represents the i-th photovoltaic node and j represents the j-th time step) are all independently sampled random numbers, and follow a standard normal distribution N(0,1) with a mean of 0 and a variance of 1. The matrix at this time... This represents three random noise sources that are completely independent of each other under ideal conditions. The sequence in the first row of the matrix represents only the independent random fluctuations of node 12, and is statistically uncorrelated with the fluctuations of node 25 in the second row (i.e., the expected covariance is 0). This independent random number matrix... It will be used as the basis input for subsequent linear transformations.

[0025] Next, proceed to step S222. The core of this step lies in analyzing the spatial correlation coefficient matrix obtained in step S1. This extracts the mathematical structure that can combine independent random variables into correlated random variables. Spatial correlation coefficient matrix. It is a K×K dimension real symmetric positive definite matrix, whose diagonal elements are all 1 (representing the perfect correlation of the nodes themselves), and the off-diagonal elements are all 1. This represents the Pearson correlation coefficient between node i and node j. Continuing with the example above where K=3, consider the spatial correlation coefficient matrix calculated in step S1 based on historical data. as follows: The matrix shows a strong correlation (0.8) between node 12 and node 25, while the correlation between node 12 and node 5 is weak (0.5). To inject this correlation structure into independent random numbers, the matrix needs to be modified... Perform Cholesky decomposition. The Cholesky decomposition aims to find a unique lower triangular matrix. such that the equation is satisfied The calculation process follows a recursive algorithm, proceeding row by row and column by column. First, the matrix is ​​calculated. The first column element: ; ; Next, calculate the elements in the second column: ; Finally, calculate the elements in the third column: After the above decomposition and calculation, the lower triangular matrix is ​​obtained. for: The matrix This is the required transformation operator, which will be used as the weight coefficient matrix for the next stage of linear transformation.

[0026] Finally, step S223 is implemented. This step is the substantive operational stage for generating correlations, physically mixing spatial correlations into the random sequence through matrix multiplication. During implementation, the 3×3 lower triangular matrix obtained in step S222 is... Left multiply the 3×2400 independent random number matrix generated in step S221 That is, to perform operations This is generated This refers to a Brownian motion increment sequence containing spatial correlation, which is strictly speaking a standardized sequence of correlated random variables. When solving SDEs, it usually needs to be multiplied by... However, in the context of this step, it represents a stochastic driving source with a correlated structure. Specifically, we analyze how this linear transformation introduces correlation. For the resulting matrix... The values ​​of the three nodes at any time step (i.e., a column vector) in the vector. , , The calculation is as follows: For the first node (the first row of the matrix), the calculation formula is... It can be seen that the random increment of the first node completely retains the original independent random number and is not affected by other nodes. For the second node (second row of the matrix), the calculation formula is... The key point here is that the increment of node 25 is not only determined by its own independent random number. The decision was also made by adding 0.8 times the random number of node 12. It is precisely because of and All of them contain This common component leads to a statistical correlation between the two. According to the variance calculation rules, the correlation coefficient recovers to 0.8 × 1 + 0 = 0.8, which is consistent with the original matrix. Consistent. For the third node (third row of the matrix), the calculation formula is: The increment of node 5 is mixed with... , and Three components. Among them... Its existence makes it associated with node 12. and The coexistence of these elements makes them related to node 25. This is achieved through analysis of the matrix. The above transformation is performed on all 2400 columns, ultimately yielding a Brownian motion increment sequence matrix of dimension [3, 2400]. In this matrix, each row represents the random fluctuation sequence of a photovoltaic node over 24 hours, and each column represents the instantaneous random increment of each node at the same moment. This sequence possesses the following statistical properties: the sequence of any single node still follows a standard normal distribution on the time axis (which guarantees the basic properties of Brownian motion), but the values ​​between different nodes at any given moment are no longer independent, but strictly obey the spatial correlation coefficient matrix set in step S1. .

[0027] As the core execution stage of stochastic process modeling, sub-step S23 is mainly responsible for transforming the mathematical model into a specific numerical trajectory. Based on the Ornstein-Uhlenbeck process parameters determined in step S1, the simulation initialization state set constructed in step S21, and the spatially correlated Brownian motion increment sequence generated in step S22, this process uses the explicit Euler-Maruyama numerical iteration method to deduce the dynamic evolution of photovoltaic power output step by step, ultimately synthesizing a complete photovoltaic power output time series sample. During implementation, the discrete-time grid sequence generated in step S21 is first invoked. and the initial moment Photovoltaic output vector At the same time, the dimension generated in step S22 is read. Brownian motion increment sequence Each column corresponds to a random increment vector for a time step. Iterative computation begins at t=0, and for each discrete time step k, from 0 to N-1, it utilizes the state from the previous time step. Calculate the state at the next time step The iterative formula adopts the Euler-Marueille scheme: In the formula, This is the normalized vector of photovoltaic power output at the current moment. The drift term reflects the deterministic trend of mean regression. For the regression velocity parameter, This is the long-term mean; The diffusion term reflects the influence of random fluctuations, where For volatility parameters, This is the scaling factor for the time scale. This is the Brownian motion increment vector corresponding to that moment. A specific explanation will be given using the parameters set in steps S1 and S22 as an example. The number of photovoltaic nodes K=3, and the time step... =0.01 hours, Ornstein-Uhlenbeck parameter set as regression rate =5.0, long-term mean =0.6, volatility =0.15. In the first time step (t=0 to t=0.01), the initial state is known. Read the first Brownian motion increment vector generated in step S22. If its value is Substituting into the iterative formula: Drift term: 5.0 × (0.6 - 0.6) × 0.01 = 0. Since the initial value is exactly equal to the mean, there is currently no regression trend. Diffusion term: Therefore, the state at the next time step t=0.01 After completing the single-step calculation, the new state vector will be obtained. As a known quantity, the above process is repeated for the next time step calculation until all N=2400 time steps have been iterated. Throughout the iteration process, the photovoltaic output is not only pulled towards the mean of 0.6 by the regression force, but also disturbed by spatially correlated random terms, causing the output trajectories of the three nodes to exhibit both random fluctuations and a synchronous changing trend. Finally, a sequence encapsulation operation is performed. The N+1 state vectors generated during the iteration process are arranged in chronological order to construct a photovoltaic output time series matrix of dimension [K, N+1]. Each row of this matrix represents the power change curve of a specific photovoltaic node over the entire simulation period, and each column represents the photovoltaic output state of the entire system at a specific moment. This synthesized photovoltaic output time series sample truly reflects the stochastic dynamic characteristics of distributed power sources in the microgrid.

[0028] In step S3, based on the microgrid network structure parameters and photovoltaic (PV) configuration parameters, a nonlinear power flow simulation is performed on the PV output time series samples to obtain a voltage amplitude sample set. It is understandable that the previous steps successfully simulated the dynamic trajectory of PV output over time using stochastic differential equations, generating a power time series that reflects the influence of meteorological factors such as cloud cover. However, simple active power values ​​cannot directly reflect the voltage safety state of the microgrid because the node voltage response in the power network is not a simple linear relationship with the injected power, but is subject to nonlinear coupling constraints from various physical factors such as network topology, line impedance, and reactive power balance. To accurately assess the specific impact of PV fluctuations on grid voltage stability, it is necessary to map the statistical power samples to the physical electrical state. This step aims to construct a physical bridge connecting the stochastic power source and the grid state. By solving nonlinear power flow equations that strictly follow Kirchhoff's laws, the randomly changing PV injected power is transformed into specific node voltage amplitude responses, thereby providing physically consistent sample data support for subsequent voltage exceedance risk quantification.

[0029] In an exemplary operation, step S3, based on the microgrid network structure parameters and photovoltaic configuration parameters, performs nonlinear power flow simulation on the photovoltaic output time series samples to obtain a voltage amplitude sample set, including: S31, based on the photovoltaic configuration parameters, performing net injected power mapping and equation system construction on the microgrid network structure parameters and photovoltaic output time series samples to obtain instances of nonlinear equation systems to be solved; S32, performing iterative solution using the Newton-Raphson method on the instances of nonlinear equation systems to be solved to obtain the converged system state vector; S33, extracting the voltage response from the converged system state vector to obtain a voltage amplitude sample set.

[0030] In the above exemplary operation, the specific process of step S3 is as follows: As a preprocessing stage for nonlinear power flow simulation, sub-step S31 is mainly responsible for fusing the static physical parameters of the microgrid with the dynamic photovoltaic power output samples to construct a mathematically solvable nonlinear equation system instance. The implementation process first involves data reading and integration. It is necessary to call the microgrid network structure parameters configured in step S1, among which the core data includes system load data and node admittance matrix parameters. The system load data comes from the historical load curves recorded by the smart meters installed in the microgrid or the short-term load forecast results based on typical daily operating scenarios. It defines the basic electricity demand (including active load) of each node at different times. and reactive load For example, in a microgrid scenario in an industrial park, the basic active power load of node 5 during the simulation period is obtained through the energy management platform. The value is a fixed 20kW. Meanwhile, the node admittance matrix is ​​a complex matrix describing the electrical topology characteristics of the microgrid, calculated from line resistance, reactance, and ground susceptance, and its elements... Includes the conductance between nodes i and j and susceptance Information. Subsequently, node power synthesis is performed to calculate the net injected power. For each discrete time step t in the photovoltaic output time series sample, where t ranges from 1 to M, it is necessary to determine the actual net power injected into the grid by each node. For example, in the photovoltaic output time series generated in step S2, the normalized output value of node 5 at time t=1 is 0.6, and the photovoltaic installed capacity of this node is extracted as 30kW from the photovoltaic configuration parameters set in step S1. Then, the actual photovoltaic power generation at that time... The result is 0.6 × 30 = 18 kW. At this point, the photovoltaic power generation is algebraically superimposed with the base load of this node: Using the above values, the net injected active power at node 5 is... =18-20=-2kW, indicating that the node at this moment acts as a load node absorbing 2kW of power from the grid. For a pure load node without photovoltaic access, its net injected power is the negative value of the load power; for a photovoltaic node without load, it is the positive value of the photovoltaic output. This process traverses all nodes of the microgrid, generating the power boundary condition vector of the entire network at the current moment. Next, the balance equations are established. Based on the complex domain expression of Kirchhoff's Current Law (KCL), a system of nonlinear algebraic equations describing the power balance of the nodes is constructed. For each node i in the microgrid, its injected power must be equal to the sum of the power flowing out through all lines connected to that node. This physical constraint is described by the following power equation in polar coordinates: In this formula, This represents the net injected active power of node i calculated in the previous step, i.e., the known quantity; and These represent the voltage magnitudes at node i and node j, respectively (unknown quantities to be solved, unless it is a saturated node or a PV node). This represents the total number of system nodes. and These are the real part (conductance) and imaginary part (susceptance) of the element in the i-th row and j-th column of the nodal admittance matrix, respectively. This represents the voltage phase angle difference between node i and node j, i.e. (phase angle) (Also an unknown quantity to be solved). This equation profoundly reveals the relationship between active power P, voltage amplitude U, and phase angle. The strong nonlinear coupling between nodes means that voltage changes at any node will affect the power distribution of the entire network through sine and cosine functions. Finally, the computational task is instantiated. The 2S nonlinear equations (including active and reactive power) constructed for the S nodes of the entire network are encapsulated into a single mathematical object to be solved. At this point, the left side of the equation set represents the estimated power calculated based on network parameters and unknown voltage variables, while the right side represents the actual net injected power (a specified value) obtained through power mapping. To initiate subsequent numerical iterations, initial guesses for the state variables need to be set. A flat-start strategy is generally adopted, assuming that the initial voltage amplitude of all PQ nodes is 1.0 per unit (pu) and the initial phase angle is 0 degrees. This complete mathematical instance, containing the equation structure, parameter coefficients (G, B), boundary conditions (P, Q), and initial value vectors, is the instance of the nonlinear equation set to be solved.

[0031] Sub-step S32, as the core computational step in nonlinear power flow simulation, primarily functions as a mathematical solver. This step receives the instance of the nonlinear equations to be solved from step S31, utilizes the efficient convergence characteristics of the Newton-Raphson iterative algorithm to find an accurate voltage solution that satisfies the power balance constraints, and finally outputs the converged system state vector. The implementation process begins with the construction and calculation of the Jacobian matrix. Based on the initial voltage state set in step S31 (usually a flat-start state, i.e., voltage amplitude of 1.0 per unit and phase angle of 0 degrees) or the voltage variables updated in the previous iteration, the partial derivatives of the power equations with respect to the state variables are calculated. For a microgrid containing S nodes, the Jacobian matrix... It is a block matrix whose elements consist of four submatrices: These partial derivative values ​​reflect the active power at the current operating point. and reactive power voltage phase angle and amplitude The sensitivity is assessed. For example, for node 5 mentioned in step S31, the partial derivative of its active power equation with respect to the voltage magnitude of adjacent nodes is calculated. This value quantifies the rate of change of active power at node 5 when the voltage of adjacent nodes changes slightly. This matrix linearizes the nonlinear power equation near the current state point. Next, the correction is calculated. The core of this step is to quantify the difference between the current state and the target state. First, the power imbalance vector is calculated. The elements therein are determined by the net injection power of each node. That is, the combined photovoltaic and load power in step S31 minus the estimated power calculated by substituting the current voltage into the power flow equation. To obtain, that is For example, in the first iteration, the specified net injected power of node 5 is -2kW (which could be reduced to a per-unit value of -0.02pu), while the estimated power calculated using the initial flat-start voltage is 0. Therefore, the active power imbalance at this node is -0.02. Subsequently, a system of linear correction equations is constructed. In this formula, The Jacobian matrix calculated above, Let be the correction vector for the voltage phase angle and magnitude to be solved. Let be the power imbalance vector. Solve the linear equations using Gaussian elimination or LU decomposition to obtain the voltage state adjustment required to bring the power imbalance towards zero. Then, perform state updates and convergence checks. Update the current system state using the obtained correction vector; the calculation formula is as follows: as well as After obtaining the new voltage state, it is necessary to determine whether the calculation has reached a stable state. This is done by checking whether the maximum absolute value in the power imbalance vector is less than the preset convergence accuracy, such as... This is implemented by... If the maximum imbalance is 0.002, which is greater than the set threshold, it means that the current solution has not yet met the physical constraint accuracy requirements, and iteration needs to continue. Finally, perform loop iteration until convergence. If the above convergence criterion is not met, the updated voltage state... and As a new baseline, we return to the first step to recalculate the Jacobian matrix and power imbalance. Due to the quadratic convergence property of the Newton-Raphson method, the error can be reduced to an acceptable range after 3 to 5 iterations. The iteration terminates when the convergence condition is met, at which point the voltage magnitude vector U and phase angle vector obtained from the last update are... This is the steady-state solution of the system that satisfies Kirchhoff's laws at that time step. This solution is encapsulated as a converged system state vector, which contains the precise voltage response data of each node in the microgrid under the current photovoltaic output scenario.

[0032] Sub-step S33, as the final stage of nonlinear power flow simulation, is primarily responsible for data cleaning, filtering, and structured reorganization. This process transforms discrete single-step calculation results into a continuous spatiotemporal response matrix, serving as a crucial data interface connecting physical simulation and statistical analysis. The implementation begins with feature separation. From the converged system state vector output in step S32, only the core indicator representing the risk of voltage exceedance—the node voltage magnitude vector—is extracted. The converged system state vector is a composite vector containing 2S elements, where S is the total number of nodes, and the first S elements may represent the phase angle of each node. The last S elements represent the voltage amplitude. For example, in a microgrid with 3 nodes, the state vector might be... After iteratively solving step S32, the specific value of the state vector at a certain time t is obtained as follows: (Phase angles are in radians, and voltages are in per-unit values). At this point, the feature separation operation will use index filtering to remove the first three phase angle data points, retaining only the last three voltage amplitude data points, to generate the voltage amplitude vector for the current moment. This operation removes phase angle information that has no direct impact on limit exceedance judgment, significantly reducing data dimensionality. Subsequently, a time-series reassembly operation is performed. Since nonlinear power flow simulation is performed step-by-step, each time step results in an independent column vector. To analyze the voltage evolution and statistical characteristics over time, these discrete vectors need to be assembled in time order. The total simulation duration is set to... =24 hours, step size =0.01 hours, total N=2400 time steps. In the simulation loop, after each time step t is completed, the calculation and extraction are performed. Then, it is appended as a column vector to the resulting matrix. The t-th column. For example, when t=1, The first column is When t=2, spliced ​​in Next, this process is similar to stitching together individual snapshots into a continuous film, ensuring the continuity and order of the voltage data in the time dimension. Finally, the set construction and output are completed. After traversing all M sample paths and all time steps under each path, a multi-dimensional voltage amplitude sample set will be formed. In a single random process simulation (one sample path), this set is represented as a K×N matrix, where K is the number of nodes and N is the number of time steps, and the element in the i-th row and j-th column is... The voltage amplitude of the i-th node at time j was precisely recorded. If we consider the total number of samples M in the Monte Carlo simulation, such as the 2000 samples set in step S1, the final voltage amplitude sample set will be a three-dimensional tensor or a list containing 2000 of the aforementioned matrices. This structured dataset fully preserves the spatiotemporal response characteristics of the microgrid under photovoltaic random fluctuations.

[0033] In step S4, the Malliavin weighted sample set is determined based on the voltage amplitude sample set and the Brownian motion increment sequence. Specifically, since voltage over-limit risk usually focuses on specific critical moments (such as periods of severe photovoltaic fluctuations or peak load), this step first selects the critical assessment moments from the three-dimensional voltage amplitude sample set (number of samples × number of nodes × time step) output in step S3. Node number to watch The voltage amplitude at that node at that moment is defined as the objective random variable. The subsequent calculations of the Malliavin derivative and weights are all performed with respect to this scalar. It was carried out.

[0034] It is understandable that in the stochastic dynamics analysis of photovoltaic microgrids, while obtaining a large number of Monte Carlo samples of voltage amplitude can depict the general outline of the voltage distribution, traditional frequency statistics methods often face problems of extremely slow convergence and excessive variance for low-probability out-of-limit events located at the tail of the distribution. To accurately reconstruct the continuous probability density function of voltage with a limited sample size, especially to accurately capture the probabilistic characteristics of extreme risk regions, it is necessary to introduce mathematical tools capable of analytically quantifying the differential relationship between the input random source and the output voltage response. The Malliavin analytical method, by establishing the differential operation rules for random variables, can transform the expectation calculation of the non-smooth indicator function in probability density estimation into the expectation calculation of the smooth weights. Therefore, this step aims to construct a Malliavin derivative matrix describing the sensitivity relationship between the Brownian motion increment (input source) and the voltage amplitude samples (output response) generated previously, and further derive the Malliavin weight sample set containing information on the dynamic evolution of the system.

[0035] In an exemplary operation, step S4, determining the Malliavin weight sample set based on the voltage amplitude sample set and the Brownian motion increment sequence, includes: S41, constructing the Malliavin derivative matrix based on the voltage amplitude sample set and the Brownian motion increment sequence; S42, performing covariance analysis on the Malliavin derivative matrix to obtain the Malliavin weight sample set.

[0036] In an exemplary operation, the specific process of step S4 is as follows: In the overall process of determining the Malliavin weight sample set, sub-step S41 is a key stage in constructing the core differential structure, aiming to numerically solve the sensitivity matrix of the microgrid voltage response to the underlying stochastic driving source. In an exemplary operation, step S41, based on the voltage amplitude sample set and the Brownian motion increment sequence, constructs the Malliavin derivative matrix, including: S411, injecting a random source perturbation into the Brownian motion increment sequence to obtain the perturbed Brownian motion increment sequence; S412, performing perturbation response propagation on the perturbed Brownian motion increment sequence to obtain the voltage amplitude sample set under perturbation; S413, performing Malliavin derivative numerical difference calculation on the perturbed voltage amplitude sample set and the voltage amplitude sample set to obtain the Malliavin derivative matrix. This process essentially constructs an analysis engine in the computer capable of quantifying the butterfly effect, that is, calculating how a small meteorological perturbation of a photovoltaic node at a certain moment propagates along the time axis and grid topology, ultimately affecting the voltage state of the entire grid.

[0037] First, S411 is implemented. The purpose of this step is to artificially introduce a controllable deviation with a clear direction and extremely small amplitude into the original random input space, serving as a probe to detect the sensitivity of the system. Here, if a full ergodic perturbation is applied to all historical moments, the computational load would be extremely large (e.g., thousands of simulations would be required for thousands of time steps). Therefore, to balance computational efficiency and accuracy, this step employs a sparse sampling strategy. Since the photovoltaic power output model uses an Ornstein-Uhlenbeck process, which has mean regression characteristics, the impact of historical perturbations on the current state decays exponentially over time. Therefore, it is not necessary to traverse all historical moments; instead, the impact is based on the photovoltaic power output regression velocity parameters. Calculate the effective relevant time window (e.g., the evaluation time). Before A time step, for example, taking (Only within the time interval) The Brownian motion increment sequence is perturbed and injected internally. For Previous historical moments were considered to have an impact on the current voltage. The sensitivity has decayed to zero, so its perturbation calculation is ignored.

[0038] Based on the above strategy, the first step in implementation is to define a deterministic perturbation function. In discrete numerical computation scenarios, to accurately distinguish the impact of different time steps and different photovoltaic nodes on voltage, the unit impulse function is selected as the orthogonal basis function. This means that when calculating the derivative of the random source with respect to the k-th photovoltaic node at the s-th time step, the deterministic perturbation function... It is a Brownian motion increment sequence matrix A matrix of the same dimension (dimension [K,N]) is used, where only the element in the k-th row and s-th column is 1, and all other elements are 0. This setting ensures the locality and orthogonality of the perturbation, allowing the subsequently calculated derivative to accurately point to specific spatiotemporal coordinates. Next, the amplitude of the small perturbation is set. This is a pre-defined, extremely small positive number used to simulate the limiting process in the definition of differential. In numerical computation, The selection of needs to balance truncation error and rounding error, and is set to . or For example, if we choose = This means that the random increment will be shifted by a factor of one million. Then, a perturbation superposition operation and new path generation are performed. The original Brownian motion increment sequence generated in step S2 is read. For example, the current process is the first Monte Carlo sample path. Based on the selected perturbation direction and amplitude , calculating the perturbed Brownian motion increment sequence , the calculation formula is . From the perspective of data structure, this operation is to take the original matrix and make a copy, then add the value at the k-th row and s-th column by . For example, in a certain operation, the sensitivity of the concerned node 12 at t=10.00 hours is focused on, the original has a value of 0.015 at this position, then the perturbed matrix has a value of 0.015001 at this position, while the values of the other 2399 time steps and other nodes in the matrix remain completely unchanged. This set of new matrices represents a random input scenario that experiences a tiny meteorological change at a specific moment (such as an instantaneous slight increase in light intensity). For cases where a complete derivative matrix needs to be calculated, the computer will traverse the time steps s within the valid correlation time window (that is, satisfying ≤s≤ ), and all nodes k, sequentially generate the corresponding perturbation sequences, and obtain the perturbed Brownian motion increment sequence .

[0039] Immediately proceed to S412. This step is a numerical implementation process of the chain rule, which aims to observe how the aforementioned tiny input perturbation passes through the complex nonlinear physical system and is finally manifested as a voltage change at the output end. This process is divided into two physical phases: output reproduction and power flow recalculation. The first phase is output reproduction. The perturbed Brownian motion increment sequence generated in step S411 is substituted back into the discrete iteration formula of the Ornstein-Uhlenbeck (OU) stochastic differential equation determined in step S23. Since the OU process is a memory-enabled autoregressive process (the state at the next moment depends on the current state), the tiny perturbation injected at the moment t=s will not only change the photovoltaic output at this moment , but also through the regression term is transmitted to all subsequent output states at moments t>s. Specifically, using the formula , the photovoltaic output trajectory of the entire time period is re-deduced starting from t=0. In the time period of t<s, due to and being identical, the reproduced output trajectory is completely consistent with the original trajectory; once the time advances to t=s, the output state undergoes a tiny jump; thereafter, in the time period of t>s, although returns to the original value, due to having already deviated, this deviation will remain in the subsequent output sequence in the form of exponential decay afterwards. Finally, a brand new set of photovoltaic output sequences with perturbation traces is obtained The second stage is power flow recalculation. The perturbed photovoltaic output sequence obtained from the above recalculation will be... As a new boundary condition, it is again input into the nonlinear power flow simulation module in step S3. This means that the power balance equations need to be reconstructed and the Newton-Raphson solver invoked for the affected time steps (i.e., all times t ≥ s). Although the perturbation is extremely small, mathematically it results in a power injection vector on the right-hand side of the equations. The solver will search for new voltage equilibrium points based on the microgrid's network structure parameters (node ​​admittance matrix). Since the power flow equations are nonlinear, small increments in power will be nonlinearly mapped to small increments in voltage amplitude. Through this process, a sample set of voltage amplitudes under disturbances is extracted. This set records the voltage response trajectories of each node in the system after a slight bias in the input source. For example, a disturbance at node 12 at 10 o'clock may cause the voltage at the same node to change from 0.980000 pu to 0.980002 pu at 10:05.

[0040] Finally, step S413 is implemented, using the finite difference method to quantify the differences in the physical response into mathematical derivative values, thereby completing the construction of the sensitivity matrix. During implementation, the difference operation is performed first. For each observed time step t, the focus is on the portion t ≥ s and each voltage node i, reading the original voltage amplitude generated in step S3 and the perturbation voltage amplitude generated in step S412. The difference between the two is then calculated. This difference represents the voltage response increment caused purely by that specific disturbance. A ratio calculation is then performed. The aforementioned voltage change... Divide by the preset disturbance amplitude According to the definition of calculus, when As it approaches 0, the ratio converges to the voltage variable. Increment of Brownian motion The directional derivative, i.e., the Malliavin derivative. For example, if the voltage increment is... The disturbance amplitude is If the Malliavin derivative at that point is 2.0, then the voltage at that moment has a linear amplification effect of 2 times on the random source. Finally, matrix construction is performed. The Malliavin derivative matrix is ​​a structured high-dimensional data volume. For a single Monte Carlo simulation, a large matrix of dimension [K×N, K×N] needs to be constructed, or more commonly, a [K, t] matrix is ​​constructed for each observation time t, where rows represent photovoltaic nodes (random source locations) and columns represent historical time steps. The elements in the matrix... The system records the specific impact of random fluctuations at time s on the voltage at time t. Due to the causality of the physical system, when s > t, future disturbances do not affect the past, therefore the corresponding derivative element is 0. By iterating through all selected disturbance directions, i.e., changing... For s and k in the matrix, repeat the process from S411 to S413 above to gradually fill each element of the matrix. The resulting Malliavin derivative matrix accurately depicts the overall picture of the stochastic sensitivity of the photovoltaic microgrid under the current operating state: it not only tells us which factors affect the voltage, but also quantifies the specific degree to which this effect decays or amplifies over time (for example, the weight of the effect of morning cloud cover on the voltage at noon is almost zero, while the weight of the current moment's irradiance change on the voltage is the greatest).

[0041] In the overall process of determining the Malliavin weight sample set, step S42 serves as a bridge connecting sensitivity analysis and probability assessment. Its core task is to compress the high-dimensional Malliavin derivative matrix and transform it into scalar statistical weights. In an exemplary operation, step S42, performing covariance analysis on the Malliavin derivative matrix to obtain the Malliavin weight sample set, includes: S421, determining the Malliavin covariance based on the Malliavin derivative matrix; S422, regularizing and inverting the Malliavin covariance to obtain the regularized inverse covariance; S423, generating the Malliavin weight sample set based on the regularized inverse covariance and the Malliavin derivative matrix. This process transforms physical-level sensitivity information into statistical-level correction coefficients, which is crucial for achieving high-precision assessment with small samples.

[0042] First, implement step S421. The purpose of this step is to quantify the total uncertainty of the target voltage variable caused by random input fluctuations over the entire simulation period. During implementation, read the Malliavin derivative matrix generated in step S41. For a specific photovoltaic node i in a single Monte Carlo simulation at a specific time t, the voltage amplitude... Its Malliavin derivative matrix is ​​a row vector (or time series) that records the Brownian motion increments of all historical times s from 0 to time t with respect to the current voltage. derivative According to the Malliavin definition of calculus, the objective variable... Malliavin covariance It is the squared norm of its derivative process in Hilbert space. In discrete numerical computation, this integral operation is transformed into a simple summation of squares. The specific formula is: In the formula, This represents the total number of time steps, corresponding to the simulation duration. , For time step, for The Malliavin derivative value at time step T. For example, step S41 calculates the derivative of the voltage at node 12 at time T=12 with respect to the previous three times (t=0, 6, 12) as 0.1, 0.5, and 0.8 respectively (0 for the other times for simplicity), and the time step size. =0.01. The Malliavin covariance of this node at that time is calculated as follows: The value 0.009 represents the cumulative effect of the total fluctuation energy of the input source transferred to the output voltage terminal under this sample path. For the case of multidimensional random variables, the Malliavin covariance is a matrix, but in this embodiment, for the analysis of a single voltage amplitude, the result is a non-negative scalar.

[0043] Next, step S422 is performed. This step aims to generate scaling factors for subsequent weight calculations while addressing potential numerical instability issues. Mathematically, calculating the Malliavin weights requires the inverse of the covariance. However, in actual numerical simulations, due to the intermittent nature of photovoltaic output (e.g., zero output at night) or the system being at certain insensitive operating points, the calculated Malliavin covariance... The result may be very close to or even equal to 0. Directly performing division can lead to numerical overflow or huge computational errors; this phenomenon is called degradation. To ensure the robustness of the computation, regularization correction is introduced. During implementation, a non-degeneracy check is first performed to examine... Is it less than a preset threshold, such as... If the value is less than the threshold, it is determined to be degenerate or close to degenerate. Regularization correction is then performed. Tikhonov regularization or a simple truncation strategy is used. In this embodiment, a small perturbation term is added, and the calculation formula is adjusted as follows: In the formula, It is a very small positive number, for example This is used to ensure that the denominator is never zero. Continuing with the previous example, if the calculated... =0.009, and set = Then regularize the inverse covariance The calculation is as follows: This value This will serve as the core scaling factor, determining the weight of the sample in the probability density estimation. For those samples with extremely low sensitivity ( (Very small) sample The regularization parameter will limit it to a controllable range to prevent it from having a destructive impact on the overall estimate.

[0044] In particular, in the complex physical scenario of voltage over-limit risk assessment for photovoltaic microgrids, a unique coupling relationship exists between system state and numerical calculation stability. When the output of the photovoltaic system fluctuates drastically due to low radiation periods (such as early morning or dusk) or rapid cloud cover, the sensitivity of the grid node voltage to underlying random disturbance sources, i.e., the Malliavin derivative, drops sharply. This change in physical state directly reflects the numerical degradation of the Malliavin covariance matrix, with some eigenvalues ​​approaching zero, making the matrix highly singular. If a static, fixed-strength regularization method is used for inversion under such ill-conditioned conditions with extremely low gradient energy levels, the correction strength for those tiny eigenvalues ​​is often insufficient, easily leading to a dramatic increase in the weighting factor after inversion—a phenomenon known as weight explosion—ultimately causing the probability density estimation to completely fail. Conversely, applying unnecessary strong regularization under normal photovoltaic conditions with high output and high sensitivity will introduce significant biased errors, impairing the physical accuracy and precision of the assessment results. Furthermore, traditional inversion methods neglect the specific distribution of the covariance matrix's characteristic spectrum. Especially when analyzing sample paths corresponding to low-probability, high-risk voltage exceedance events, the derivative matrix often exhibits ill-conditioned characteristics. Simple truncation or correction cannot find a proper balance between suppressing numerical noise and preserving true physical information. To address the technical challenges arising from the imbalance between numerical stability and statistical unbiasedness, this application employs an adaptive spectral regularization inversion method based on gradient energy sensing. This method introduces a gradient energy sensing factor that dynamically quantifies the degree of covariance matrix degradation, adaptively adjusting the regularization intensity in real time during the matrix inversion process. This ensures robustness of numerical computation across all operating conditions while maximizing the accuracy of the evaluation results.

[0045] Specifically, in an exemplary preferred operation, step S422, regularizing and inverting the Malliavin covariance to obtain the regularized inverse covariance, includes: The Malliavin covariance is decomposed into spectral values ​​and its gradient energy density is calculated to obtain the eigenvalue set, eigenvector matrix, and gradient energy density. First, the complex covariance matrix is ​​decomposed into its most basic constituent units, namely eigenvalues ​​and eigenvectors, to understand its inherent numerical structure and ill-conditioned risks. The magnitude of the eigenvalues ​​directly reflects the sensitivity of the voltage response in different directions, while the distribution of the eigenvalues ​​determines the condition number and inversion stability of the matrix, thus generating a scalar index, the gradient energy density, used to accurately and sensitively quantify the overall degradation degree of the covariance matrix, providing a decision-making basis for subsequent adaptive adjustments. In practice, the input Malliavin covariance matrix is ​​first... Perform spectral decomposition (or eigenvalue decomposition) to obtain ,in It is an orthogonal matrix composed of eigenvectors, that is, an eigenvector matrix. It is composed of eigenvalues The resulting diagonal matrix. Then, based on the eigenvalue set... Calculate gradient energy density The calculation formula is as follows: In this formula, This represents the dimension of the covariance matrix, i.e., the number of photovoltaic nodes. The determinant of a matrix, which is the product of all its eigenvalues, reflects the volume scaling effect of the gradient space. The trace of the matrix represents the sum of all its eigenvalues, and the total gradient energy represents the total gradient energy. This is a very small overflow prevention bias constant used to prevent calculation errors caused by a zero determinant. Its set value is much smaller than the machine's precision, such as... This formula is essentially the ratio of the geometric mean to the arithmetic mean of the eigenvalues. When all eigenvalues ​​are of uniform magnitude, this value approaches 1, indicating a healthy matrix; when any eigenvalue approaches zero, this value rapidly approaches 0, indicating that the matrix is ​​ill-conditioned. For example, in a two-node system (N=2), if the eigenvalues ​​of the covariance matrix at a certain moment are {0.9, 0.8}, then its... Approximately 0.994, close to 1, indicates a good system state; if the characteristic value changes to {0.9, 0.001} at another moment, then... The value plummeted to approximately 0.067, close to 0, clearly indicating the risk of matrix degradation.

[0046] Based on the gradient energy density and eigenvalue set, a diagonal matrix of inverse eigenvalues ​​is constructed after filtering. That is, static regularization cannot distinguish between eigenvalues ​​of different magnitudes and systems with different health states. This step, however, uses a nonlinear filtering function dependent on gradient energy density to apply appropriate corrections to the inversion process of each eigenvalue. Thus, when the system is healthy, almost no correction is performed, maintaining the unbiasedness of the inversion; when the system tends towards ill-conditioning, the regularization strength for small eigenvalues ​​is automatically increased, effectively suppressing the numerical divergence of their reciprocals, thereby ensuring computational stability. In specific implementation, the diagonal matrix... Each original feature value in The corrected inverse value of each is calculated using the following adaptive damping filter function. : In the formula The reference damping coefficient is a preset positive number that provides basic numerical stability for the inversion process. Its value can be set according to the numerical range of the system model; for example, it can be taken as... This is to ensure that the denominator is not zero in extreme cases. This is the energy decay rate parameter, used to adjust the sensitivity of the filter function to changes in gradient energy density. It can be set empirically, for example, to 0.1. (Exponential term) It is the key adaptive switch: when When the value approaches 1 (system is healthy), the exponent term approaches 0, and the denominator is approximately equal to 1. The entire expression degenerates into That is, standard inverse, which guarantees accuracy; when When the exponent approaches 0 (ill-conditioned system), the exponent approaches 1, and the denominator becomes... Even at this time As the denominator approaches zero, the denominator also... The existence of a variable that is not zero effectively avoids numerical explosion. Continuing the example above, when... When ≈0.067, for a small eigenvalue of 0.001, the calculation of its corrected inverse value will be affected. The significant impact was suppressed to a reasonable range, rather than the enormous value of 1000 obtained by simple inversion. All corrected inverse eigenvalues Together they form a new diagonal matrix That is, the inverse eigenvalue diagonal matrix after filtering.

[0047] The filtered inverse eigenvalue diagonal matrix and eigenvector matrix are reconstructed using a regularized inverse covariance matrix to obtain the regularized inverse covariance. It should be understood that the preceding analysis and corrections are performed in the spectral space (eigenspace), and the results need to be transformed back to the original physical coordinate system to obtain a meaningful and complete inverse covariance matrix. This utilizes the unchanged original eigenvector matrix, which carries the system response direction information. The matrix multiplication operation is performed with the stabilized inverse eigenvalue diagonal matrix, thus providing a solid mathematical foundation for accurate and robust assessment of voltage over-limit risks in photovoltaic microgrids under all operating conditions. (Inverse eigenvalue diagonal matrix) This allows for the recombination of a final inverse matrix that is mathematically well-formed and physically preserves the information of the main gradient directions. In practice, based on the inverse theorem of spectral decomposition, reconstruction is achieved through the following matrix multiplication operations: , here It is the eigenvector matrix The transpose of . Through this reconstruction process, the adaptive correction effect applied to the spectral domain is mapped back to the original matrix space. The final regularized inverse covariance is obtained. It solves the problem of instability in the inversion of the Malliavin covariance matrix caused by the randomness of photovoltaic power output, and avoids the accuracy loss caused by static regularization. It can be directly and safely used for subsequent Malliavin weighted samples.

[0048] Finally, step S423 is implemented, which involves generating the Malliavin weighted sample set based on the regularized inverse covariance and the Malliavin derivative matrix. This step is the final synthesis stage of the Malliavin analytical method, where a special stochastic process is constructed and integral approximation is performed to obtain the scalar weighting factors used for probability density estimation. During implementation, an auxiliary process is first constructed based on the integral by parts formula. This process is defined as the product of the Malliavin derivative and the inverse covariance: This process This represents the sensitivity distribution at a normalized covariance scale. Subsequently, the weighting factors are calculated using a discretized approximation of the Skorhod integral (as a generalization of the Itô integral). The general formula for discrete calculation is as follows: ,Will By expanding and substituting, we obtain a more specific numerical implementation: The formula contains two contributions: the first part It is the sensitivity versus the random noise at the current moment ( The weighted product of approximately white noise reflects the direct contribution of the random source to the weights; Part Two The second derivative term (or the divergence of the derivative) is a correction term for the stochastic integral, used to eliminate bias. In actual programming, the second derivative is usually ignored or handled through simplification assumptions, focusing on calculating the first part. Let's illustrate this numerically with an example of the simplified first part: at time step k, the derivative... =0.5, inverse covariance =111.1, Brownian motion increment =0.02, time step =0.01. Therefore, the contribution of this time step to the weight is: (0.5 × 111.1) × (0.02 / 0.01) × 0.01 = 1.111. Repeat this calculation for all N time steps and sum them up to obtain the single scalar weight value for this Monte Carlo sample path. For example, the final calculated weight of a certain path is... =1.5. Finally, the sample set is generated. For the M independent sample paths generated by the Monte Carlo simulation (such as the 2000 paths set in step S1), all operations from S421 to S423 above are repeated. A corresponding weight factor is generated for each path. Finally, these 2000 weighting factors are aggregated and stored to form a Malliavin weighted sample set of dimension [M,1]. This set corresponds one-to-one with the original voltage amplitude sample set, together forming the basic data for subsequent voltage probability density estimation.

[0049] In step S5, voltage probability density estimation and over-limit risk calculation are performed on the Malliavin weighted sample set and voltage amplitude sample set to obtain the voltage over-limit risk probability. In other words, after the calculations in the preceding steps, a large number of discrete voltage amplitude samples and their corresponding Malliavin weighting factors have been obtained. Although these data contain all the response information of the microgrid under photovoltaic random fluctuations, their existence only as a sample set is insufficient to intuitively and quantitatively reflect the degree of voltage over-limit risk. Traditional histogram methods based on frequency statistics often result in drastic fluctuations in estimation results when dealing with low-probability tail events due to sample sparsity, failing to meet the engineering requirements for high-reliability risk assessment. To transform discrete, weighted sample points into continuous, smooth, and accurate probability distribution curves, and to calculate statistically significant over-limit probabilities based on these curves, this application finally adopts a density estimation method based on Malliavin integral theory. This step aims to utilize the smoothing effect unique to Malliavin weights to reconstruct the voltage probability density function through weighted expectation calculations, thereby achieving accurate quantification of voltage over-limit risk (especially extreme risk) without increasing the sampling scale.

[0050] In an exemplary operation, step S5, which involves estimating the voltage probability density and calculating the over-limit risk on the Malliavin weighted sample set and the voltage amplitude sample set to obtain the voltage over-limit risk probability, includes: S51, estimating the probability density function on the voltage amplitude sample set based on the Malliavin weighted sample set to obtain voltage probability density function curve data; S52, performing numerical integration on the voltage probability density function curve data in a first over-limit interval and a second over-limit interval to obtain the voltage over-limit risk probability, wherein the first over-limit interval is... The second over-limit interval is , As the lower limit of allowable voltage, This represents the upper limit of the permissible voltage.

[0051] In the above operation, the specific process of step S5 is as follows: Step S51, as the core calculation link of risk quantification, is mainly responsible for combining discrete Monte Carlo samples with Malliavin weights and generating a continuous voltage probability density function (PDF) through mathematical transformation. The implementation process first constructs the probability distribution evaluation point set. In order to draw a delicate and smooth probability density curve, a series of dense observation points need to be preset within the possible voltage fluctuation range. According to the operating characteristics of the microgrid, a voltage range covering normal operation and possible over-limit areas is selected, such as [0.90, 1.10] per unit value (pu). Within this range, uniform sampling is performed according to a fixed small step size (such as 0.001pu) to generate a vector containing hundreds of discrete voltage values. For example, setting the starting point to 0.9, the ending point to 1.1, and the step size to 0.001 generates a set of 201 evaluation points [0.900, 0.901, ..., 1.100]. These points y will serve as the independent variables of the probability density function, used to calculate their corresponding probability density values ​​p(y) point by point. Subsequently, smoothing and weighted statistical calculations are performed. This is the essence of the Malliavin method; it uses integration by parts to transform the differentiation of the probability density into the expectation of the weighted function, avoiding the bandwidth selection problem in traditional kernel density estimation and eliminating the discrete error of the histogram method. For each voltage value y in the evaluation point set, the following weighted expectation formula is used for estimation: In this formula, The total number of samples for the Monte Carlo simulation, such as the 2000 set in step S1; This represents the voltage amplitude of the k-th sample path at a specific time (or as an average over the entire time period); The Malliavin weight factor is the one corresponding to the k-th sample. This is the smooth approximation function introduced. Although the theoretical Malliavin formula uses Dirac... The derivative form of a function or indicator function is used, but in numerical implementations, to avoid numerical oscillations and singularities, a smooth function (such as the derivative form of the sigmoid function or Gaussian kernel function) is used for approximation. For example, one could choose... ,in This is a smoothing parameter (bandwidth), with a small value such as 0.005. Taking the calculation of the probability density at y=0.95 as an example: iterate through all 2000 samples. For the first sample, such as its voltage... =0.98, weight =1.2. Calculate the difference 0.98-0.95=0.03, substitute it into the smoothing function, and get... If the value is 0.01, then the contribution of this sample is 0.01 × 1.2 = 0.012. For the second sample, its voltage... =0.951, very close to the evaluation point, with a difference of 0.001. Substituting this into the smoothing function yields a large value, such as 10.0, indicating its weight. =0.8, then the contribution is 10.0 × 0.8 = 8.0. Summing the weighted contributions of all 2000 samples and dividing by the total number of samples M, we obtain the probability density value p(0.95) for the point y = 0.95. Finally, the curve data is generated. The probability distribution is evaluated using the point set. Repeat the above calculation process for each point in the vector, and finally obtain the corresponding probability density value vector. The evaluation point vector With density value vector The combination of these parameters constitutes the complete voltage probability density function curve data. This curve not only exhibits smoothness in the high-probability region (near the mean), but more importantly, it shows significant smoothness in the extremely low-probability tail region (such as below 0.92 or above 1.08) due to the introduction of sensitivity weights. This data effectively suppresses Monte Carlo noise and provides non-zero probability density values ​​with statistical confidence. It accurately describes the probability intensity distribution of microgrid node voltages falling within various possible value ranges.

[0052] Sub-step S52, as the final calculation step in the entire voltage over-limit risk detection process, primarily transforms probability density information into risk indicators. This step receives the voltage probability density function curve data generated in step S51, and based on the voltage safety standards preset in step S1, uses a numerical integration algorithm to calculate the total probability that the voltage amplitude falls within the unsafe region, ultimately outputting the voltage over-limit risk probability. This process discretizes the continuous mathematical function and quantifies the risk through surface accumulation, ensuring the accuracy of the assessment results. The implementation process first requires clarifying the boundary conditions of the integration, namely the allowable lower voltage limit. and the upper limit of permissible voltage These two parameters are determined during the system parameter initialization phase in step S1, based on the power system supply voltage deviation standard (such as GB / T 12325) or the specific operating procedures for microgrids. Allowable lower voltage limit. It is typically set to 95% of the nominal voltage (i.e., 0.95 pu), with an upper limit on the allowable voltage. These thresholds are typically set to 105% of the nominal voltage (i.e., 1.05 pu). These two thresholds define the safe corridor for microgrid voltage compliance; any voltage fluctuations exceeding this range (below 0.95 or above 1.05) are considered out-of-limit events, constituting the first out-of-limit interval for subsequent integration calculations. Second over-limit interval Subsequently, the voltage probability density function curve data output in step S51 is read. This data consists of a series of discrete evaluation point pairs, including the voltage abscissa vector. and the corresponding probability density ordinate vector If the evaluation point set covers the range [0.90, 1.10] with a step size of 0.001 pu, containing a total of 201 points, then this discrete data set needs to be partitioned and numerically integrated. Numerical integration is usually achieved using the trapezoidal rule or Simpson's rule. First, the first boundary interval is calculated. Lower bound probability In the numerical implementation, the voltage horizontal axis vector is traversed. Filter out all that meet the requirements Data points with a value ≤ 0.95. For example, within the interval [0.90, 0.95], the probability density value... It exhibits a non-zero distribution. and Taking two adjacent sampling points as an example (such as 0.940 and 0.941), if the corresponding probability density values ​​are respectively =0.50 and =0.52, then the area of ​​probability within this small interval is approximately From =0.90 to =0.95 The lower bound probability can be obtained by summing the areas of all the tiny trapezoids. If, after summation, we get... ≈0.0150. This indicates that, statistically, the probability of this node voltage being below 0.95 pu is 1.5%. Next, the second over-limit interval is calculated. upper limit probability Similarly, filter out all that meet the criteria. Data points ≥1.05. Because fluctuations in photovoltaic output simulated by the Ornstein-Uhlenbeck process may cause voltage increases at certain times, such as the presence of a non-zero probability density distribution within the interval [1.05, 1.10], the trapezoidal integration process is repeated, summing the areas under the probability density curves within this region. The calculated upper limit probability is... ≈0.0192. This indicates that the probability of the voltage exceeding 1.05 pu due to excessive photovoltaic output or other factors is 1.92%. Finally, the total risk probability is synthesized and output. Based on the principle of mutual exclusion of probabilities, the lower limit probability and the upper limit probability are added together to obtain the total voltage limit exceedance risk probability. The calculation formula is as follows: In this formula, the integral symbol In discrete computing, this corresponds to the trapezoidal area summation operation described above; This is the high-precision probability density function obtained in step S51 using Malliavin weight estimation. Substituting the example values ​​above, the total risk probability is... =0.0150 + 0.0192 = 0.0342. This scalar value of 0.0342, or 3.42%, precisely quantifies the probability of a microgrid facing voltage security threats under current photovoltaic configurations and meteorological conditions. This value not only integrates the randomness of photovoltaic output (reflected by Malliavin weights) and the physical constraints of the grid structure (reflected by power flow calculations), but also, due to the use of Malliavin smoothing technology, provides a more stable estimate with smaller variance than the direct frequency statistics of the traditional Monte Carlo method, even with a limited sample size (e.g., 2000 samples). Finally, this risk probability value is output to the result verification and analysis module or the human-computer interaction interface, serving as the core basis for assessing microgrid security.

[0053] In summary, the photovoltaic microgrid voltage over-limit risk detection method based on the embodiments of this application is clarified. By constructing an Ornstein-Uhlenbeck stochastic process model, the random fluctuations and time correlation of photovoltaic power output are accurately simulated, generating a photovoltaic power output sequence and Brownian motion increment that conform to actual physical characteristics. Based on this, the system voltage response is obtained by combining nonlinear power flow simulation, and Malliavin analytical theory is introduced. By constructing the Malliavin derivative matrix and performing covariance analysis, a Malliavin weighted sample set is generated, and then this weight set is used for voltage probability density estimation and over-limit risk calculation. This effectively solves the problems of slow convergence speed and heavy computational burden of existing Monte Carlo simulation methods in low-probability event assessment, and also overcomes the shortcomings of point estimation methods in terms of insufficient accuracy when dealing with strong nonlinearity and non-Gaussian distributions. By introducing the Malliavin weight factor, a sensitivity transfer relationship from random source disturbance to system voltage response is established. Using weight correction technology, a high-precision approximation of the voltage probability distribution (especially the tail over-limit region) can be achieved with a small sample size, thereby significantly improving the computational efficiency of photovoltaic microgrid voltage over-limit risk detection while ensuring assessment accuracy.

[0054] Figure 4 This is a block diagram of a photovoltaic microgrid voltage over-limit risk detection system according to an embodiment of this application. Figure 4 As shown, the photovoltaic microgrid voltage over-limit risk detection system 100 according to an embodiment of this application includes: a parameter acquisition module 110, used to acquire microgrid network structure parameters, photovoltaic configuration parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters, wherein the Ornstein-Uhlenbeck process parameters include photovoltaic output regression rate parameters, long-term mean, and volatility parameters, and the numerical calculation parameters include time step, number of photovoltaic nodes, total number of samples, and total simulation duration; and a modeling and sample sequence generation module 120, used to process the microgrid network structure parameters and Ornstein-Uhlenbeck process parameters. The system employs a stochastic process modeling and sample sequence generation method with numerical calculation parameters to obtain photovoltaic power output time series samples and Brownian motion increment sequences. A simulation module 130 performs nonlinear power flow simulation on the photovoltaic power output time series samples based on microgrid network structure parameters and photovoltaic configuration parameters to obtain a voltage amplitude sample set. A weighted sample set generation module 140 determines the Malliavin weighted sample set based on the voltage amplitude sample set and Brownian motion increment sequences. A voltage over-limit risk probability generation module 150 estimates the voltage probability density and calculates the over-limit risk from the Malliavin weighted sample set and voltage amplitude sample set to obtain the voltage over-limit risk probability.

[0055] Those skilled in the art will understand that the specific operations of each step in the aforementioned photovoltaic microgrid voltage over-limit risk detection system have been referenced above. Figures 1 to 3 The method for detecting voltage over-limit risks in photovoltaic microgrids has been described in detail, and therefore, its repeated description will be omitted.

Claims

1. A method for detecting voltage over-limit risks in photovoltaic microgrids, characterized in that, include: The microgrid network structure parameters, photovoltaic configuration parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters are obtained. The Ornstein-Uhlenbeck process parameters include photovoltaic output regression rate parameters, long-term mean and volatility parameters, and the numerical calculation parameters include time step, number of photovoltaic nodes, total number of samples and total simulation duration. Stochastic process modeling and sample sequence generation were performed on microgrid network structural parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters to obtain photovoltaic power output time series samples and Brownian motion increment sequences. Based on the microgrid network structure parameters and photovoltaic configuration parameters, nonlinear power flow simulation is performed on the photovoltaic power output time series samples to obtain the voltage amplitude sample set; Based on the voltage amplitude sample set and the Brownian motion increment sequence, the Malliavin weighted sample set is determined, including: Based on the voltage amplitude sample set and the Brownian motion increment sequence, the Malliavin derivative matrix is ​​constructed, including: injecting random source perturbation into the Brownian motion increment sequence to obtain the perturbed Brownian motion increment sequence; performing perturbation response propagation on the perturbed Brownian motion increment sequence to obtain the perturbed voltage amplitude sample set; and performing Malliavin derivative numerical difference calculation on the perturbed voltage amplitude sample set and the voltage amplitude sample set to obtain the Malliavin derivative matrix. Covariance analysis is performed on the Malliavin derivative matrix to obtain the Malliavin weighted sample set, including: determining the Malliavin covariance based on the Malliavin derivative matrix; regularizing and inverting the Malliavin covariance to obtain the regularized inverse covariance; and generating the Malliavin weighted sample set based on the regularized inverse covariance and the Malliavin derivative matrix. Voltage probability density estimation and over-limit risk calculation are performed on the Malliavin weighted sample set and voltage amplitude sample set to obtain the voltage over-limit risk probability.

2. The photovoltaic microgrid voltage over-limit risk detection method according to claim 1, characterized in that, Stochastic process modeling and sample sequence generation were performed on microgrid network structural parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters to obtain photovoltaic power output time series samples and Brownian motion increment sequences, including: Based on numerical calculation parameters and Ornstein-Uhlenbeck process parameters, a simulation initialization state set is determined, which includes an initial output vector and a discrete event grid sequence. Based on numerical calculation parameters, correlation decomposition and correlation transformation are performed on the spatial correlation coefficient matrix between implicit photovoltaic nodes in the microgrid network structure parameters to obtain the Brownian motion increment sequence. Based on the Ornstein-Uhlenbeck process parameters and Brownian motion increment sequence, explicit Euler iteration and power sequence synthesis are performed on the simulated initial state set to obtain photovoltaic power output time series samples.

3. The photovoltaic microgrid voltage over-limit risk detection method according to claim 2, characterized in that, Based on numerically calculated parameters, correlation decomposition and correlation transformation are performed on the spatial correlation coefficient matrix between implicit photovoltaic nodes in the microgrid network structure parameters to obtain the Brownian motion increment sequence, including: Generate an independent random number matrix that follows a standard normal distribution. The dimension of the independent random number matrix is ​​. ,in, This refers to the number of photovoltaic nodes; The spatial correlation coefficient matrix between implicit photovoltaic nodes is decomposed to obtain the lower triangular matrix. A linear transformation of the independent random number matrix based on the lower triangular matrix is ​​used to obtain the Brownian motion increment sequence.

4. The photovoltaic microgrid voltage over-limit risk detection method according to claim 1, characterized in that, Based on microgrid network structure parameters and photovoltaic configuration parameters, nonlinear power flow simulation is performed on photovoltaic power output time series samples to obtain a voltage amplitude sample set, including: Based on photovoltaic configuration parameters, net injected power mapping and equation system construction are performed on microgrid network structure parameters and photovoltaic power output time series samples to obtain instances of nonlinear equation systems to be solved. The Newton-Raphson method is used to iterate through the instance of the nonlinear equation system to obtain the converged system state vector. Voltage response is extracted from the converged system state vector to obtain a voltage amplitude sample set.

5. The photovoltaic microgrid voltage over-limit risk detection method according to claim 1, characterized in that, Regularizing and inverting the Malliavin covariance to obtain the regularized inverse covariance includes: Spectral decomposition and gradient energy density calculation were performed on the Malliavin covariance to obtain the eigenvalue set, eigenvector matrix, and gradient energy density. Based on the gradient energy density and eigenvalue set, construct the filtered inverse eigenvalue diagonal matrix; The regularized inverse covariance matrix is ​​reconstructed from the filtered inverse eigenvalue diagonal matrix and eigenvector matrix to obtain the regularized inverse covariance.

6. The photovoltaic microgrid voltage over-limit risk detection method according to claim 1, characterized in that, Voltage probability density estimation and over-limit risk calculation are performed on the Malliavin weighted sample set and voltage amplitude sample set to obtain the voltage over-limit risk probability, including: Based on the Malliavin weighted sample set, the probability density function is estimated from the voltage amplitude sample set to obtain the voltage probability density function curve data; The voltage probability density function curve data is numerically integrated in the first and second over-limit intervals to obtain the voltage over-limit risk probability, where the first over-limit interval is... The second over-limit interval is , As the lower limit of allowable voltage, This represents the upper limit of the permissible voltage.

7. A photovoltaic microgrid voltage over-limit risk detection system, characterized in that, include: The parameter acquisition module is used to acquire microgrid network structure parameters, photovoltaic configuration parameters, Ornstein-Uhlenbeck process parameters, and numerical calculation parameters. Among them, the Ornstein-Uhlenbeck process parameters include photovoltaic output regression rate parameters, long-term mean and volatility parameters, and the numerical calculation parameters include time step, number of photovoltaic nodes, total number of samples and total simulation duration. The modeling and sample sequence generation module is used to perform stochastic process modeling and sample sequence generation on microgrid network structural parameters, Ornstein-Uhlenbeck process parameters and numerical calculation parameters to obtain photovoltaic power output time series samples and Brownian motion increment sequences. The simulation module is used to perform nonlinear power flow simulation on photovoltaic output time series samples based on microgrid network structure parameters and photovoltaic configuration parameters to obtain voltage amplitude sample sets. The weighted sample set generation module is used to determine the Malliavin weighted sample set based on the voltage amplitude sample set and the Brownian motion increment sequence, including: Based on the voltage amplitude sample set and the Brownian motion increment sequence, the Malliavin derivative matrix is ​​constructed, including: injecting random source perturbation into the Brownian motion increment sequence to obtain the perturbed Brownian motion increment sequence; performing perturbation response propagation on the perturbed Brownian motion increment sequence to obtain the perturbed voltage amplitude sample set; and performing Malliavin derivative numerical difference calculation on the perturbed voltage amplitude sample set and the voltage amplitude sample set to obtain the Malliavin derivative matrix. Covariance analysis is performed on the Malliavin derivative matrix to obtain the Malliavin weighted sample set, including: determining the Malliavin covariance based on the Malliavin derivative matrix; regularizing and inverting the Malliavin covariance to obtain the regularized inverse covariance; and generating the Malliavin weighted sample set based on the regularized inverse covariance and the Malliavin derivative matrix. The voltage over-limit risk probability generation module is used to perform voltage probability density estimation and over-limit risk calculation on the Malliavin weighted sample set and voltage amplitude sample set to obtain the voltage over-limit risk probability.

Images