Harmonic power flow calculation method and system of photovoltaic power distribution system based on scene recursion

Abstract

The invention relates to the technical field of harmonic power flow calculation, in particular to a harmonic power flow calculation method and system of a photovoltaic power distribution system based on scene recursion. The method comprises the steps of generating a plurality of typical operation scenes through a Monte Carlo simulation method based on a probability density function of photovoltaic irradiation intensity and basic load, and constructing a multi-scene sample library; performing clustering processing on the historical operation data and the generated scene sample, and identifying and extracting a plurality of typical operation modes; establishing a harmonic source model for a harmonic source in the power distribution network according to the generated typical operation scene; calculating a time-varying harmonic coupling matrix; and selecting a typical working mode, and calculating a harmonic power flow index of each node of the system under each order of harmonic according to the established harmonic source model and the time-varying harmonic coupling matrix. According to the method, the calculation time delay is remarkably reduced while the analysis precision is ensured, and the method is particularly suitable for a modern power distribution system with high new energy permeability and ubiquitous power electronic equipment.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of harmonic flow calculation, in particular to a harmonic flow calculation method and system for a photovoltaic power distribution system based on scenario recursion. BACKGROUND

[0002] With the large-scale access of distributed photovoltaic systems, the proportion of power electronic equipment in the power grid has rapidly increased, and the harmonic pollution problem has become increasingly serious. Traditional harmonic analysis methods usually perform power flow simulation based on fixed operating conditions, ignoring uncertain factors such as light, load fluctuations, etc., resulting in insufficient accuracy of the analysis results in actual operation. In addition, harmonic flow simulation under complex scenarios has large computational load and long time delay, affecting online applications.

[0003] In recent years, common harmonic flow methods include frequency domain node method, harmonic injection model and multi-frequency coupling algorithm, and means such as matrix sparsity optimization, sensitivity pre-computation and parallel acceleration are gradually introduced to improve simulation efficiency. At present, harmonic flow algorithms are widely used in harmonic state estimation, source identification, filter configuration and monitoring point optimization, etc. However, they still face challenges such as complex modeling, large computational overhead and difficulty in real-time response.

[0004] In view of the above problems, there is an urgent need for a harmonic calculation method that can balance uncertainty scenario modeling and analysis efficiency to support faster and more accurate harmonic state evaluation. SUMMARY

[0005] In view of the problems in the prior art, the present application provides a harmonic flow calculation method and system for a photovoltaic power distribution system based on scenario recursion, and the specific technical solutions are as follows:

[0006] A harmonic flow calculation method for a photovoltaic power distribution system based on scenario recursion, comprising the following steps:

[0007] Step S1, based on the probability density function of photovoltaic irradiance intensity and basic load, a plurality of typical operating scenarios are generated by Monte Carlo simulation method, and a multi-scenario sample library is constructed;

[0008] Step S2, the historical operating data and the generated scenario samples are clustered and processed to identify and extract several typical operating modes;

[0009] Step S3, according to the generated typical operating scenarios, a harmonic source model is established for the harmonic sources in the power distribution network;

[0010] Step S4, a time-varying harmonic coupling matrix is calculated;

[0011] Step S5, a typical operating mode is selected, and according to the established harmonic source model and the time-varying harmonic coupling matrix, the harmonic flow indicators of each node of the system under each order harmonic are calculated.

[0012] Preferably, the probability density function f(r) of the photovoltaic irradiance intensity in the step S1 is:

[0013]

[0014] wherein r is the solar irradiance intensity; r max is the maximum irradiance; a, b are shape parameters, and Γ(.) is the Gamma function, and the photovoltaic power is obtained according to the solar irradiance intensity;

[0015] The probability density function f(P L ) of the basic load is:

[0016]

[0017] wherein P L is the load power; s is the standard deviation of the load power; and m is the expected value of the load power.

[0018] Preferably, the step S1 generates a plurality of typical operation scenarios by the Monte Carlo simulation method, and the multi-scenario sample library is constructed, specifically:

[0019] Based on the probability density function f(r) of the photovoltaic irradiance intensity and the probability density function f(P L ) of the basic load, the Monte Carlo simulation method is used to sample the photovoltaic irradiance intensity and the basic load to generate S typical operation scenarios, wherein: S=N PV ×T;

[0020] wherein N PV is the number of photovoltaic scenarios, and T is the number of time periods divided by the typical day;

[0021] The characteristic vector of the s-th typical operation scenario is defined as:

[0022]

[0023] wherein P s is the photovoltaic power sampled by the Beta distribution under the s-th typical operation scenario, and P s is the conventional load sampled by the normal distribution under the s-th typical operation scenario.

[0024] Preferably, the step S2 clusters the historical operation data and the generated scenario samples, identifies and extracts a plurality of typical operation modes, specifically including the following steps:

[0025] (1) When clustering, first set the number of clusters K, and randomly initialize K cluster centers which can be randomly extracted from the samples, and are set as:

[0026]

[0027] wherein, is the initialized kth cluster center;

[0028] (2) For the feature vector x s of the s-th typical operating scenario, calculate its Euclidean distance to each cluster center: d sk = x s - M k ;

[0029] wherein, d sk is the Euclidean distance of the feature vector x s of the s-th typical operating scenario to the kth cluster center M k ;

[0030] (3) Assign the sample to the nearest cluster using the following function:

[0031]

[0032] Cluster(x s ) represents the cluster label to which the feature vector x s of the s-th typical operating scenario is assigned, and the value is an integer 1, 2, …, K; indicates taking the value of k that minimizes the expression, that is, the corresponding cluster label;

[0033] (4) For each cluster C k , update its centroid as the mean of all points within the cluster:

[0034]

[0035] wherein, |C k | is the sample set of the kth cluster, is the updated centroid of the kth cluster at the t+1th iteration;

[0036] (5) If all clusters no longer change, the algorithm stops; otherwise, return to step (2) to recalculate d sk .

[0037] Preferably, the step S3 of establishing a harmonic source model for the harmonic source in the power distribution network according to the generated typical operating scenario specifically comprises the following steps:

[0038] According to the power injection under each working mode M K , a harmonic source model such as a photovoltaic inverter is established according to its corresponding power state, and the harmonic source model is as follows:

[0039]

[0040] wherein, for the operating mode M K the injected current vector at the hth harmonic, for the operating mode M K the injected current of the ith harmonic source at the hth harmonic, modeled as a random variable, satisfying:

[0041]

[0042] wherein, for the operating mode M K the injected current of the ith harmonic source at the hth harmonic, modeled as a random variable, satisfying: for the operating mode M K the injected current of the ith harmonic source at the hth harmonic, modeled as a random variable, satisfying:

[0043] Preferably, the step S4 specifically comprises the following steps:

[0044] A relationship model between harmonic voltages at different time nodes is constructed, specifically as follows:

[0045] V (h) (t+1) = H (h) (t) V (h) (t) ;

[0046] wherein, V (h) (t+1) and V (h) (t) are harmonic voltage vectors of each node at t+1 time and t time under frequency f and hth harmonic; H (h) (t) is a time-varying harmonic coupling matrix at t time under frequency f;

[0047] To solve the time-varying harmonic coupling matrix H (h) (t), a solving equation is constructed as follows:

[0048]

[0049] wherein, is a vector form of the hth harmonic voltage of the n node at t time, is a coupling coefficient of the hth harmonic of the n node and the n node at t time;

[0050] A linear system is constructed as follows:

[0051]

[0052] wherein, n is the number of nodes, and T-1 is the number of samples that can be used for solving.

[0053] The input sample matrix X is singular value decomposed to obtain:

[0054] X = UΣV * ;

[0055] wherein U is a column-orthogonal matrix;∑ is a singular value diagonal matrix; V * is a conjugate transpose of the row-orthogonal matrix V;

[0056] Finally, the pseudo-inverse matrix X + is calculated:

[0057] X + = V∑ + U * ;

[0058] Further, the harmonic coupling matrix H (h) (t) is obtained:

[0059] H (h) (t) = YX + = YV∑ + U * ;

[0060] wherein "+" is a pseudo-inverse symbol.

[0061] Preferably, the step S5 of selecting a typical operating mode and calculating the harmonic current of each node of the system under each order harmonic according to the established harmonic source model and time-varying harmonic coupling matrix comprises the following steps:

[0062] First, a typical operating mode M K is selected according to a known typical operating scenario characteristic vector:

[0063]

[0064] At time t, the definition of the hth harmonic injection under the kth typical operating mode M K is as follows:

[0065]

[0066] For each order harmonic h, the system impedance matrix is defined as follows and is directly obtained from the system parameters:

[0067] Z h ∈ C n×n ;

[0068] The ohm law is used to calculate:

[0069]

[0070] Subsequently, the harmonic voltage value at time t+1 can be obtained according to the coupling matrix:

[0071] V (h) (t+1) = H (h) (t) V (h) (t);

[0072] I (h) (t+1)=V (h) (t+1) / Z h ;

[0073] The harmonic flow index of each node under each order harmonic is calculated as follows:

[0074] P (h) (t+1)=V (h) (t+1)×I (h) (t+1)

[0075] V (h) (t+1)=v a (t+1)+jv b (t+1);

[0076] wherein, I (h) (t+1) is the harmonic current vector at t+1 time under frequency f and h order harmonic, P (h) (t+1) is the harmonic power vector at t+1 time under frequency f and h order harmonic; v a (t+1) and v b (t+1) are the real part and imaginary part of V (h) (t+1) respectively.

[0077] A harmonic flow calculation system of a photovoltaic power distribution system based on scene recursion, applying the method, comprising: a running scene generation module for generating a plurality of typical running scenes based on the probability density function of photovoltaic irradiance and basic load through the Monte Carlo simulation method, and constructing a multi-scene sample library;

[0078] A typical pattern recognition module for clustering processing the historical running data and the generated scene samples, identifying and extracting a plurality of typical running patterns;

[0079] A harmonic source establishment module for establishing a harmonic source model according to the generated typical running scenes, and the harmonic source in the power distribution network;

[0080] A harmonic coupling matrix calculation module for calculating a time-varying harmonic coupling matrix;

[0081] A harmonic flow index calculation module for selecting a typical working mode, calculating the harmonic flow index of each node under each order harmonic according to the established harmonic source model and the time-varying harmonic coupling matrix.

[0082] A computer readable storage medium, comprising a stored program, wherein the program controls the device where the computer readable storage medium is located to execute the harmonic flow calculation method of a photovoltaic power distribution system based on scene recursion when the program is running.

[0083] A processor for running a program, wherein the program performs the harmonic power flow calculation method of the scene-based recursive photovoltaic power distribution system when running.

[0084] Compared with the prior art, the present application has the following beneficial effects:

[0085] The method of the present application significantly reduces the calculation delay while ensuring analysis accuracy, and is particularly suitable for modern power distribution systems with high new energy penetration rate and widespread power electronic equipment. Through the mode of "work mode driving + coupling sensitivity analysis", the power distribution network harmonic simulation is changed from "panoramic simulation" to "mode abstract driving", meeting the rapid diagnosis requirements under multiple scenarios and dynamic operation. BRIEF DESCRIPTION OF DRAWINGS

[0086] In order to more clearly illustrate the specific embodiments of the present application or the technical solutions in the prior art, the drawings needed in the specific embodiments or prior art description will be briefly introduced below. In all the drawings, similar elements or parts are generally identified by similar reference numerals. In the drawings, each element or part is not necessarily drawn according to the actual scale.

[0087] Figure 1 The method flowchart of the present application.

[0088] Figure 2 The test result graph of the present application, wherein Figure 2 (a) is a TVE histogram, Figure 2 (b) is a THD line graph.

[0089] Figure 3 The system principle diagram of the present application. DETAILED DESCRIPTION

[0090] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0091] It should be understood that when used in the present specification and the appended claims, the terms "comprise" and "include" indicate the presence of described features, whole, steps, operations, elements and / or components, but do not exclude the presence or addition of one or more other features, whole, steps, operations, elements, components and / or sets thereof.

[0092] It should also be understood that the terminology used in this specification is for the purpose of describing particular embodiments only and is not intended to limit the invention. As used in this specification and the appended claims, the singular forms “a,” “an,” and “the” are intended to include the plural forms unless the context clearly indicates otherwise.

[0093] It should also be further understood that the term "and / or" as used in this specification and the appended claims refers to any combination of one or more of the associated listed items and all possible combinations, and includes such combinations.

[0094] Example 1:

[0095] like Figure 1 As shown, this embodiment provides a harmonic power flow calculation method for a photovoltaic power distribution system based on scenario recursion, including the following steps:

[0096] Step S1: Based on the probability density functions of photovoltaic irradiance and base load, several typical operating scenarios are generated using Monte Carlo Simulation (MCS) to construct a multi-scenario sample library.

[0097] The probability density function f(r) of the photovoltaic output power is:

[0098]

[0099] Where r is the solar irradiance; r max α represents the maximum irradiance; α and β are shape parameters, and Γ(·) is the Γ function.

[0100] The probability density function f(P) of the base load L )for:

[0101]

[0102] Among them, P L σ is the load power; μ is the standard deviation of the load power; and μ is the expected value of the load power.

[0103] Based on the probability density function f(r) of photovoltaic output power and the probability density function f(P) of base load. L The Monte Carlo simulation method was used to sample the photovoltaic output power and base load, generating S typical operating scenarios, where S = N. PV ×T;

[0104] Where N PV denoted as the number of photovoltaic scenarios, and T as the number of time periods divided into typical days.

[0105] The feature vector of the s-th typical operating scenario is defined as:

[0106]

[0107] wherein, is the photovoltaic power sampled by the Beta distribution under the s-th typical operating scenario, is the conventional load sampled by the normal distribution under the s-th typical operating scenario.

[0108] In step S2, the historical operating data and the generated scenario samples are clustered, and several typical operating modes are identified and extracted as important condition inputs of the harmonic analysis, so as to improve the adaptability and efficiency of the analysis.

[0109] The S typical operating scenarios S and the historical monitoring data such as harmonic voltage and harmonic current are clustered and classified, and K representative typical operating modes are identified. The clustering uses the K-means clustering method, and the specific steps are as follows:

[0110] (1) When clustering, first set the number of clusters K, and randomly initialize K cluster centers which can be randomly extracted from the samples, and are set as:

[0111]

[0112] wherein, is the initialized k-th cluster center;

[0113] (2) For the feature vector x s of the s-th typical operating scenario, calculate the Euclidean distance d sk of each cluster center: s k

[0114] wherein, d sk is the Euclidean distance of the feature vector x s of the s-th typical operating scenario to the k-th cluster center M k .

[0115] (3) The following function is used to assign the sample to the nearest cluster:

[0116]

[0117] Cluster(x s ) represents the cluster label to which the feature vector x s of the s-th typical operating scenario is assigned, and the value is an integer 1, 2, …, K; indicates the k value that makes the expression minimum, that is, the corresponding cluster label. ​​

[0118] (4) Update the centroid of each cluster C k as the mean of all points within the cluster:

[0119]

[0120] where |C k | is the sample set of the kth cluster, is the updated centroid of the kth cluster at iteration t+1.

[0121] (5) If all clusters no longer change, the algorithm stops; otherwise, return to step (2) to recalculate d sk .

[0122] According to the above steps, a set of typical patterns {M1, M2,..., M K} is obtained, where K is the number of patterns. Each pattern corresponds to a typical photovoltaic or load characteristic as the subsequent simulation input.

[0123] Step S3, according to the generated typical operating scenario, a harmonic source model is established for the harmonic source in the power distribution network, covering multiple typical harmonic frequencies, and the uncertainty modeling is introduced for the harmonic current parameters.

[0124] According to the power injection under each working mode M K , a harmonic source model such as a photovoltaic inverter is established according to its corresponding power state, and the harmonic source model is as follows:

[0125]

[0126] where, is the injection current vector under the hth harmonic of the working mode M K , represents the injection current of the ith harmonic source under the hth harmonic, which is modeled as a random variable, satisfying:

[0127]

[0128] where, is the average value of the injection current of the ith harmonic source under the hth harmonic of the working mode M K ; is the standard deviation of the injection current of the ith harmonic source under the hth harmonic of the working mode M K , indicating the device difference or control effect. Step S4, calculate the time-varying harmonic coupling matrix. Specifically, the following steps are included:

[0129] The relationship model between the harmonic voltages at different time nodes is constructed, which is as follows:

[0130] V (h)(t+1) = H (h) (t) V (h) (t) ;

[0131] where V (h) (t+1) and V (h) (t) are the harmonic voltage vectors of the n nodes at time t+1 and t, respectively, at frequency f and hth harmonic; H (h) (t) is the time-varying harmonic coupling matrix at time t at frequency f;

[0132] To solve the time-varying harmonic coupling matrix H (h) (t), the solving equation is constructed as follows:

[0133]

[0134] where, is the vector form of the harmonic voltage of the nth node at time t and hth harmonic, is the coupling coefficient of the nth node at time t and hth harmonic of the nth node;

[0135] The linear system is constructed as follows:

[0136]

[0137] where n is the number of nodes, and T-1 is the number of samples that can be used for solving.

[0138] The input sample matrix X is singular value decomposed to obtain:

[0139] X = UΣV * ;

[0140] where U is a column-orthogonal matrix; Σ is a singular value diagonal matrix; V * is the conjugate transpose of the row-orthogonal matrix V;

[0141] Finally, the pseudo-inverse matrix X + is calculated:

[0142] X + = VΣ + U * ;

[0143] The harmonic coupling matrix H (h) (t) is further obtained:

[0144] H (h) (t) = YX + = YVΣ + U * ;

[0145] where "+" is a pseudo-inverse symbol.

[0146] Step S5, selecting a typical working mode According to the established harmonic source model And the time-varying harmonic coupling matrix H (h) (t), the harmonic flow indicators of each node of the system under each order harmonic are calculated, and the specific steps are as follows.

[0147] First, a typical working mode M is selected according to a known typical operating scenario characteristic vector K :

[0148]

[0149] At time t, the hth harmonic injection under the kth typical working mode M is defined as: K

[0150]

[0151] For each order harmonic h, the system impedance matrix is defined as follows and can be directly obtained from the system parameters:

[0152] Z h ∈C n×n ;

[0153] Using Ohm's law, the following is obtained:

[0154]

[0155] Subsequently, the harmonic voltage value at time t+1 can be obtained according to the coupling matrix:

[0156] V (h) (t+1)=H (h) (t)V (h) (t)

[0157] I (h) (t+1)=V (h) (t+1) / Z h The harmonic flow indicators of each node under each order harmonic are calculated as follows:

[0158]

[0159] Where I (h) (t+1) is the harmonic current vector at time t+1 under frequency f and hth harmonic, P (h) (t+1) is the harmonic power vector at time t+1 under frequency f and hth harmonic; v a (t+1) and v b (t+1) are the real part and imaginary part of V (h) (t+1) respectively. In addition, the present application also calculates the following indicators.

[0160] ​

[0161] and ang(t+1) are the amplitude and phase angle of V (h) (t+1) respectively, V 1,j is the fundamental voltage of node j, H is the highest harmonic number considered; V(t) is the actual harmonic voltage, is the harmonic voltage measurement; n is the number of all nodes in the system, and j is a certain node among them. VTHD j is the harmonic voltage distortion rate at node j; TVE j is the total vector error of harmonic voltage at node j; V h,j is the hth harmonic voltage value at node j.

[0162] In order to highlight the detection performance of the algorithm, the Monte Carlo simulation (MCS) method and the 2m+1 point estimation (PME) method are selected as the comparative algorithms. Under the power frequency condition, the total vector error (TVE) and the total harmonic voltage distortion error (VTHD) of the three algorithms are compared by randomly selecting 10 nodes to evaluate the performance of the algorithms.

[0163] The test uses the IEEE-33 node system, and the final results are shown in Figure 2 . In the TVE column chart of Figure 2 (a), the MCS method shows the smallest error value on all nodes, and the TVE is basically controlled between 1% and 5%, verifying the high precision characteristics of the MCS method as the benchmark algorithm. In comparison, the TVE of the PME method on most nodes exceeds 6%, while the Recursive method of the present application is superior to the PME method in overall accuracy, and the TVE can be reduced to 1.78% at the lowest, showing strong robustness and local adaptability.

[0164] As can be seen from the THD line chart of Figure 2 (b), the THD error of the MCS method on all nodes remains below 3%, which is the smallest on average. The THD error of the Recursive method on most nodes is below 4%, which is significantly better than the PME method and shows good tracking ability on complex nodes such as nodes 5, 8, and 9. The THD of the PME method on multiple nodes exceeds 5%, which is significantly different from the other two methods. In summary, the Recursive method of the present application not only maintains high computational efficiency, but also achieves an accuracy level close to the MCS method, which is superior to the PME method. The MCS method is a benchmark algorithm in the harmonic flow field, which is specially used for detecting the accuracy of the algorithm, and has high accuracy but takes a long time.

[0165] Example 2:

[0166] As Figure 3As shown, based on the same inventive concept as example 1, the embodiment provides a harmonic power flow calculation system of a photovoltaic power distribution system based on scene recursion, and the method is applied, comprising:

[0167] A running scene generation module is configured to generate a plurality of typical running scenes based on the probability density function of photovoltaic irradiance and basic load by Monte Carlo simulation method, and construct a multi-scene sample library;

[0168] A typical pattern recognition module is configured to cluster the historical running data and the generated scene samples, recognize and extract a plurality of typical running patterns;

[0169] A harmonic source establishment module is configured to establish a harmonic source model according to the generated typical running scenes and the harmonic sources in the power distribution network;

[0170] A harmonic coupling matrix calculation module is configured to calculate a time-varying harmonic coupling matrix;

[0171] A harmonic power flow index calculation module is configured to select a typical working mode, calculate the harmonic power flow index of each node of the system under each order harmonic according to the established harmonic source model and the time-varying harmonic coupling matrix.

[0172] Example 3:

[0173] Based on the same inventive concept as example 1, the embodiment provides a computer readable storage medium, the computer readable storage medium comprises a stored program, wherein the program controls the device where the computer readable storage medium is located to execute the harmonic power flow calculation method of the photovoltaic power distribution system based on scene recursion when the program is running.

[0174] Example 4:

[0175] Based on the same inventive concept as example 1, the embodiment provides a processor, the processor is used to run a program, wherein the program executes the harmonic power flow calculation method of the photovoltaic power distribution system based on scene recursion when the program is running.

[0176] Those of ordinary skill in the art can realize that the units of the examples described in combination with the embodiments disclosed herein can be realized in electronic hardware, computer software or a combination of both. In order to clearly illustrate the interchangeability of hardware and software, the components of the examples have been described in the above description in general terms. Whether the functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of the present application.

[0177] In the embodiments of the present application, it should be understood that the division of units is only a logical function division, and actual implementation can have another division manner, for example, multiple units can be combined as one unit, one unit can be split into multiple units, or some features can be ignored, etc.

[0178] In addition, each functional unit in each embodiment of the present application can be integrated in one processing unit, or each unit can exist physically, or two or more units can be integrated in one unit. The integrated unit can be realized in the form of hardware or in the form of a software functional unit.

[0179] When the integrated unit is realized in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer readable storage medium. Based on such understanding, the technical solutions of the present application essentially or the part that contributes to the prior art or the whole or part of the technical solutions can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in each embodiment of the present application. The foregoing storage medium includes: a U disk, a read-only memory (ROM, Read-Only Memory), a random access memory (RAM, Random Access Memory), a mobile hard disk, a magnetic disk or an optical disk, and various program code storage media.

[0180] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application, and they should be covered in the scope of the claims and the specification of the present application.

Claims

1. A method for harmonic power flow calculation of a photovoltaic power distribution system based on scenario recursion, characterized in that, The method comprises the following steps: Step S1, generating a plurality of typical operation scenarios based on the probability density function of photovoltaic irradiance intensity and basic load through a Monte Carlo simulation method to construct a multi-scenario sample library; Step S2, clustering the historical operation data and the generated scenario samples to identify and extract a plurality of typical operation modes; Step S3, establishing a harmonic source model for the harmonic source in the distribution network according to the generated typical operation scenarios; Step S4, calculating a time-varying harmonic coupling matrix; Step S5, selecting a typical working mode, and calculating the harmonic current index of each node at each order harmonic according to the established harmonic source model and the time-varying harmonic coupling matrix.

2. The method of claim 1, wherein, The probability density function f(r) of the photovoltaic irradiance intensity in the step S1 is: where r is the solar irradiance; r max is the maximum irradiance; a, b are shape parameters, and G(.) is the Gamma function, the photovoltaic power is obtained according to the solar irradiance; The probability density function f(P L ) of the base load is: where P L is the load power; σ is the standard deviation of the load power; μ is the expected value of the load power.

3. The method of claim 1, wherein, The step S1 of generating a plurality of typical operation scenarios through the Monte Carlo simulation method to construct a multi-scenario sample library specifically comprises: based on the probability density function f(r) of the photovoltaic irradiance intensity, the probability density function f(P L ), the photovoltaic irradiance intensity and the basic load are sampled by using the Monte Carlo simulation method to generate S typical operation scenarios, wherein: S = N PV x T; where N PV is the number of photovoltaic scenarios, and T is the number of time periods in a typical day. The s-th typical operation scenario characteristic vector is defined as: wherein, Psis the photovoltaic power sampled from the Beta distribution for the s-th typical operating scenario, Psis the regular load sampled from the normal distribution for the s-th typical operating scenario.

4. The method of claim 1, wherein, The step S2 of clustering the historical operation data and the generated scenario samples to identify and extract a plurality of typical operation modes specifically comprises the following steps: (1) Clustering, first set the number of clusters K, randomly initialize K cluster centers Randomly selected from the sample, set to: wherein, is the initialized kth cluster center; (2) the feature vector x for the s-th typical operating scenario s , compute its Euclidean distance to each cluster center: d sk = x s - M k ; wherein d sk is the Euclidean distance of the feature vector x s to the k-th cluster center M k . (3) the following function is used to assign the samples to the nearest cluster: Cluster(x s ) represents the feature vector x of the s-th typical running scenario s The cluster label assigned to it, which is an integer value of 1, 2, …, K; represents the value of k that minimizes the expression, that is, the corresponding cluster label; (4) For each cluster C k , update its centroid to be the mean of all points within the cluster: where |C k is the sample set of the kth cluster, is the updated centroid of the kth cluster at the t+1th iteration. (5) If all clusters no longer change, the algorithm stops; otherwise, return to step (2) to recalculate d sk .

5. The method of claim 1, wherein, The step S3 of establishing a harmonic source model for the harmonic source in the distribution network according to the generated typical operation scenarios specifically comprises the following steps: According to the power injection under each working mode M K According to the corresponding power state, the harmonic source model of the photovoltaic inverter is established, and the harmonic source model is as follows: wherein, for operating mode M K the injected current vector at the hth harmonic, denotes operating mode M K the injected current at the hth harmonic from the ith harmonic source, modeled as a random variable, satisfying: wherein, is the operating mode M K is the average value of the injected current of the i-th harmonic source at the h-th harmonic; is the operating mode M K is the standard deviation of the injected current of the i-th harmonic source at the h-th harmonic.

6. The harmonic power flow calculation method of a scene recursion based photovoltaic power distribution system according to claim 1, wherein, The step S4 specifically The step S4 specifically The relationship model between the harmonic voltages at different time nodes is constructed, and specifically as follows: V (h) (t+1) = H (h) (t) V (h) (t); where V (h) (t+1) and V (h) (t) are the harmonic voltage vectors of the nodes at time t+1 and t at frequency f and hth harmonic; H (h) (t) is the time-varying harmonic coupling matrix at time t at frequency f. To solve the time-varying harmonic coupling matrix H (h) (t), the solving equation is constructed as follows: wherein, is the vector form of the harmonic voltage of the nth time at the hth node t, is the coupling coefficient of the nth time at the hth harmonic of the n node and the n node t; The linear system is constructed as follows: Wherein, n is the number of nodes, and T-1 is the number of samples that can be used for solving. The singular value decomposition is performed on the input sample matrix X to obtain: X = U∑V * ; wherein U is a column-orthogonal matrix;∑ is a singular-value diagonal matrix; V * is the conjugate transpose of the row-orthogonal matrix V. Finally, the pseudo-inverse matrix X is calculated + : X + = V∑ + U * ; Further, the harmonic coupling matrix H is obtained (h) (t): H (h) (t) = YX + = YV∑ + U * ; Wherein, "+" is a pseudo-inverse symbol.

7. The method of claim 1, wherein, The step S5 of selecting a typical working mode, and calculating the harmonic current of each node at each order harmonic according to the established harmonic source model and the time-varying harmonic coupling matrix specifically comprises the following steps: First, a typical operating mode M is selected according to a known typical operating scenario characteristic vector K : At time t, for the k-th typical operating mode M K is defined as: For each harmonic h, the system impedance matrix is defined as follows, which is directly obtained from the system parameters: Z h ∈C n×n ; The harmonic voltage value at the t+1 time is obtained according to the coupling matrix: The harmonic current index of each node at each order harmonic is calculated as follows: V (h) (t+1) = H (h) (t) V (h) (t); I (h) (t+1) = V (h) (t+1) / Z h ; The method of any one of claims 1 to 7 comprises: P (h) (t+1) = V (h) (t+1) x I (h) (t+1) V (h) (t+1) = v a (t+1) + jv b (t+1); where I (h) (t+1) is the harmonic current vector at time t+1 at frequency f and hth harmonic, P (h) (t+1) is the harmonic power vector at time t+1 at frequency f and hth harmonic; v a (t+1), v b (t+1) are the real and imaginary parts of V (h) (t+1), respectively.

8. A system for harmonic power flow calculation of a photovoltaic power distribution system based on scenario recursion, characterized by, An operation scenario generation module configured to generate a plurality of typical operation scenarios based on the probability density function of photovoltaic irradiance intensity and basic load through a Monte Carlo simulation method to construct a multi-scenario sample library; A typical mode identification module configured to cluster the historical operation data and the generated scenario samples to identify and extract a plurality of typical operation modes; A harmonic source establishment module configured to establish a harmonic source model for the harmonic source in the distribution network according to the generated typical operation scenarios; A harmonic coupling matrix calculation module configured to calculate a time-varying harmonic coupling matrix; A harmonic current index calculation module configured to select a typical working mode, and calculate the harmonic current index of each node at each order harmonic according to the established harmonic source model and the time-varying harmonic coupling matrix. ​ 9. A computer-readable storage medium, characterized in that, The computer readable storage medium comprises a stored program, wherein the program controls the device where the computer readable storage medium is located to execute the harmonic power flow calculation method of the scene recursion-based photovoltaic power distribution system according to any one of claims 1 to 7 when the program is running.

10. A processor, comprising: The processor is configured to run a program, wherein the program executes the harmonic power flow calculation method of the scene recursion-based photovoltaic power distribution system according to any one of claims 1 to 7 when the program is running.

Images