Optimal scheduling method and system of photovoltaic power distribution network

Abstract

The present application relates to power distribution network optimization scheduling technical field, especially to a kind of photovoltaic power distribution network optimization scheduling method and system.The method comprises the following steps: obtaining the historical scheduling state of photovoltaic power distribution network, and extracting the transformer load state between each regional node;Then, based on the node transformer load state, the output voltage instability of node is quantified and analyzed, the instability intensity increment is deduced, and the clustering data of instability intensity is obtained;Finally, based on instability intensity clustering data, the load transfer logic of adjacent node is designed, the load pressure of current node is reduced, and the optimization scheduling of photovoltaic power distribution network is realized.The improvement processing of the present application to power distribution network optimization scheduling technology makes the power distribution network optimization scheduling technology more perfect.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of power distribution network optimization scheduling, and particularly relates to a photovoltaic power distribution network optimization scheduling method and system. BACKGROUND

[0002] Photovoltaic power generation has obvious intermittency and volatility, and its output power is greatly affected by weather, season and sunlight conditions, which brings new challenges to the operation and scheduling of traditional power distribution networks. The traditional power distribution network scheduling method usually assumes that the load and power supply are relatively stable, while the frequent power fluctuations in the photovoltaic power distribution network will cause node voltage fluctuations, transformer load imbalance, and even cause local voltage instability or system instability. In addition, photovoltaic power generation is often concentrated in a specific area, making the load pressure of some power distribution nodes too large, while the load of other nodes is relatively light, resulting in uneven distribution of power grid resources and potential overload risk. However, the traditional photovoltaic power distribution network optimization scheduling method has the problem of inaccurate analysis of transformer node output instability, resulting in large error in the optimization scheduling of the power distribution network. SUMMARY

[0003] Therefore, it is necessary to provide a photovoltaic power distribution network optimization scheduling method and system to solve at least one of the above technical problems.

[0004] To achieve the above purpose, a photovoltaic power distribution network optimization scheduling method, the method comprising the following steps:

[0005] Step S1: obtaining the historical scheduling state of the photovoltaic power distribution network; extracting the transformer load state between different regional power distribution nodes according to the historical scheduling state to obtain the node transformer load state;

[0006] Step S2: quantifying the output voltage instability of the power distribution network transformer according to the node transformer load state to obtain node output voltage instability data; and deriving the instability intensity increment in the time dimension of the node output voltage instability data to obtain node instability intensity clustering data;

[0007] Step S3: designing the load transfer logic between adjacent nodes based on the node instability intensity clustering data to reduce the load pressure of the current node to obtain the load transfer logic to perform the optimization scheduling of the photovoltaic power distribution network.

[0008] Preferably, the present application also provides a photovoltaic power distribution network optimization scheduling system for performing the photovoltaic power distribution network optimization scheduling method as described above, the photovoltaic power distribution network optimization scheduling system comprising:

[0009] The state correlation extraction module is configured to acquire historical scheduling states of the photovoltaic power distribution network; perform transformer load state correlation extraction between different regional power distribution nodes according to the historical scheduling states; and obtain node transformer load states;

[0010] The instability strength increment derivation module is configured to perform output voltage instability quantification of the transformer of the power distribution network according to the node transformer load states, to obtain node output voltage instability data; and perform time-dimension instability strength increment derivation on the node output voltage instability data, to obtain node instability strength clustering data.

[0011] The load transfer logic design module is configured to perform load transfer logic design between adjacent nodes based on the node instability strength clustering data, to reduce the load pressure of the current node, to obtain the load transfer logic, and to perform the optimal scheduling of the photovoltaic power distribution network.

[0012] The beneficial effects of the present application can comprehensively understand the operation status of the distribution network and the load condition of the transformer by obtaining the historical scheduling state of the photovoltaic power distribution network. By correlating and extracting the load states between different regional distribution nodes, the load change trend and potential problems of each node can be identified. This process helps to reveal the areas of load imbalance, predict which nodes are at risk of overload in advance, and provide data support for subsequent optimization scheduling. Through the mining of historical data, more accurate load state analysis can be performed to improve the accuracy and efficiency of distribution network scheduling. By quantifying the output voltage instability state of the node transformer, it can be directly identified which nodes have voltage fluctuations or instability under certain load conditions. Further, through the time dimension instability intensity increment derivation, the voltage instability change rule in different time periods can be obtained, thereby providing a basis for instability intensity clustering analysis. Through the analysis of node output voltage instability data, fine-grained monitoring of the distribution network can be achieved to detect potential voltage instability risks in advance and take appropriate measures to avoid large-scale voltage fluctuations or collapse of the system, improving the stability and security of the power grid. Based on the clustering data of node instability intensity, the load transfer logic between adjacent nodes is designed to achieve the optimization scheduling of the load in the distribution network. Through reasonable load transfer, the load pressure of some nodes can be reduced to avoid voltage instability or equipment damage caused by overload. The load transfer logic design not only balances the load distribution of each node, but also improves the overall operation efficiency of the distribution network, reducing the probability of power grid failure. This step optimizes the scheduling strategy to improve the reliability, flexibility and scheduling efficiency of the photovoltaic power distribution network, thereby ensuring the stability of power supply and supporting the efficient access of renewable energy. Therefore, the present application is an improvement on the traditional optimization scheduling method for a photovoltaic power distribution network, which solves the problem of inaccurate transformer node output instability analysis in the traditional optimization scheduling method for a photovoltaic power distribution network, thereby reducing the error of the optimization scheduling of the distribution network, improving the accuracy of the transformer node output instability analysis, and reducing the error of the optimization scheduling of the distribution network. BRIEF DESCRIPTION OF DRAWINGS

[0013] Figure 1 A step flowchart for an optimization scheduling method for a photovoltaic power distribution network;

[0014] Figure 2 A Figure 1 A detailed implementation step flowchart for step S2 in the method;

[0015] Figure 3 A Figure 2 A detailed implementation step flowchart for step S23 in the method. DETAILED DESCRIPTION

[0016] Please refer to Figures 1 to 3A photovoltaic power distribution network optimization scheduling method, the method comprises the following steps:

[0017] Step S1: obtain the historical scheduling state of the photovoltaic power distribution network; according to the historical scheduling state, the transformer load state correlation extraction between different regional power distribution nodes is carried out, and the node transformer load state is obtained;

[0018] Step S2: according to the node transformer load state, the output voltage instability of the power distribution network transformer is quantified, and the node output voltage instability data is obtained; the instability intensity increment of the node output voltage instability data in time dimension is derived, and the node instability intensity clustering data is obtained;

[0019] Step S3: based on the node instability intensity clustering data, the load transfer logic design between adjacent nodes is carried out, so as to reduce the load pressure of the current node, so as to obtain the load transfer logic, and the optimization scheduling of the photovoltaic power distribution network is carried out.

[0020] In the embodiment of the application, reference is made to Figure 1 The application is a kind of photovoltaic power distribution network optimization scheduling method, and the steps of the method are as follows:

[0021] Step S1: obtain the historical scheduling state of the photovoltaic power distribution network; according to the historical scheduling state, the transformer load state correlation extraction between different regional power distribution nodes is carried out, and the node transformer load state is obtained;

[0022] In the embodiment of the present application, the active power data, the reactive power data, the voltage amplitude data, the current amplitude data and the switch state data of each power distribution node recorded every 5 minutes in the past 365 days are derived from the SCADA system of the photovoltaic power distribution network to form an original historical scheduling state data set; the data set is divided into 8 sub-regions according to geographical regions, each sub-region contains 12 power distribution nodes, and a total of 96 nodes; a sliding time window analysis is performed on the original data, the time window length is set to 4 hours, and the step length is 30 minutes; in each time window, the transformer load rate of each node is calculated, the load rate is equal to the current active power divided by the rated capacity of the transformer, and the rated capacity is uniformly set to 1000 kVA; at the same time, the power fluctuation standard deviation of the transmission line between adjacent nodes is calculated, and the fluctuation standard deviation threshold is set to 15 kW; when the load rate of a node exceeds 85% in the continuous 3 time windows, and the power fluctuation standard deviation with adjacent nodes is greater than 15 kW, it is marked as a high load correlation period; the node current sequence corresponding to all high load correlation periods is extracted, the current sequence sampling frequency is 10 Hz, and the duration is 2 hours; the Pearson correlation coefficient matrix calculation is performed on the current sequence, the matrix dimension is 96x96, the correlation coefficient threshold is set to 0.72, and the node pairs with correlation coefficient greater than 0.72 are reserved as strong correlation groups; the load state synchronicity test is performed on the nodes in the strong correlation group, the dynamic time warping algorithm DTW is used for synchronicity test, and the distance threshold is set to 0.35; finally, the transformer load state of 96 nodes is output, including 5 parameters of load rate peak value, load rate average value, load fluctuation variance, correlation node number and correlation strength coefficient.

[0023] Step S2: quantifying the output voltage instability of the power distribution network transformer according to the node transformer load state to obtain node output voltage instability data; performing instability intensity increment derivation on the node output voltage instability data in the time dimension, thereby obtaining node instability intensity clustering data;

[0024] In the embodiment of the present application, the fast Fourier transform (FFT) is performed on the current sequence in the output of each node transformer load state, the number of sampling points is 7200, the frequency domain resolution is 0.0139 Hz, and the amplitudes of the fundamental component and the 3rd, 5th, 7th, 9th and 11th harmonic components are extracted; if the amplitude of any odd harmonic exceeds 8% of the amplitude of the fundamental component, the current abnormal state at this time is marked, a binary abnormality marker sequence is generated, and the sequence length is equal to the original time sequence length; the winding heat increment index calculation is performed on each current sample marked as abnormal, the calculation method is that the current effective value is squared, multiplied by the copper loss coefficient 0.0012, and then multiplied by the duration 0.2 seconds to obtain the single-point heat increment index; the cumulative integral is performed on the continuous abnormal points to form a winding heat accumulation data sequence; the first-order difference is performed on the heat accumulation data to obtain the heat accumulation change rate, and when the change rate of 5 consecutive sampling points is greater than 0.8 W / s, it is determined that the magnetic circuit enters the saturation trend interval; the voltage waveform in the saturation interval is sampled, the sampling frequency is 10 kHz, and the sampling length is 200 ms; the zero-crossing detection is performed on the sampled waveform, the positive half-cycle and negative half-cycle peak voltage values are recorded, the absolute difference value is calculated, the theoretical peak voltage 311 V is divided, and the waveform asymmetry degree is obtained; at the same time, the slope calculation is performed on the waveform, the sampling point interval is 0.1 ms, the voltage difference between adjacent points is divided by the time difference to obtain the instantaneous slope; the slope variance of the 5 sampling points before and after the peak point is calculated to obtain the slope variance; the weighted fusion is performed on the asymmetry degree and the slope variance, the weight coefficients are 0.6 and 0.4 respectively, and the output voltage distortion strength sequence is output; the sliding average filtering is performed on the distortion strength sequence, the window length is 50 ms, the step length is 10 ms, and the smoothed node output voltage instability data is obtained; the time sequence is constructed for the data, the time period is divided in units of 1 hour, the time series divergence increment gradient calculation is performed on each time period, the method is that the value of the next time period is subtracted from the value of the previous time period to obtain the gradient sequence; the fractional order differentiation is performed on the gradient sequence, the order is 0.7, the Grünwald-Letnikov discretization formula is used, the step length is 10 sampling points, and the divergence numerical fractional order sequence is obtained; the autoregressive fitting is performed on the sequence, the order is fixed to 3, the Yule-Walker equation is used to solve the coefficients, and the future 3-hour prediction value sequence is obtained; the prediction value sequence is subtracted from the current value sequence to obtain the node instability strength increment data in the time dimension; the K-means clustering is performed on the increment data, the number of cluster centers is set to 5, the initial cluster center is determined by the elbow rule, the maximum iteration number is set to 100, the convergence threshold is set to 0.001, and the instability strength category label and the cluster center coordinates to which each node belongs are output to form the node instability strength clustering data.

[0025] Step S3: based on the node instability strength clustering data, the load transfer logic design between adjacent nodes is performed to reduce the load pressure of the current node, so as to obtain the load transfer logic to perform the optimal scheduling of the photovoltaic power distribution network.

[0026] In the embodiment of the present application, a one-dimensional convolution operation is performed on the node instability strength clustering data output in step S2, the convolution kernel length is set to 7, the weight distribution adopts a Gaussian function, the standard deviation is 1.2, the step is 1, and the boundary padding mode is zero padding. 96 node instability convolution strength values of each node are output; the convolution strength value of each node is normalized, and the normalized range is 0 to 1, which is used as the node emergency weight; the average load rate of each node is calculated, which is obtained by counting the load rate data every 5 minutes in the past 24 hours, a total of 288 points, and taking the arithmetic mean value as the average load rate; the nodes with an average load rate greater than 0.8 are marked as transfer source nodes; the configuration state of the adjacent nodes around the transfer source nodes is obtained, including the transformer rated capacity, the line impedance value, and the topological connection structure, which is represented by an adjacency matrix, and the matrix elements are 1 for physical direct connection and 0 for no direct connection; the transfer load of each transfer source node is calculated, which is obtained by subtracting the threshold value of 0.75 from the current load rate and then multiplying the transformer rated capacity of 1000 kVA; the carrying capacity difference of each adjacent node is calculated, which is obtained by multiplying the threshold value of 0.8 by the rated capacity, subtracting the current average load rate multiplied by the rated capacity, and then multiplying the line transmission capacity, which is taken as the line thermal stability limit value, with a unit of kVA; the product of the carrying capacity difference and the line transmission capacity is calculated as the upper limit of the transferable capacity; the electrical distance between the transfer source nodes and the adjacent nodes is calculated, which is equal to the line resistance value multiplied by 1.2, plus the reactance value multiplied by 0.8, and arranged in ascending order; the transferable load is distributed to the adjacent nodes in sequence after the arrangement, and the distribution rule is: if the upper limit of the transferable capacity of the current adjacent node is greater than the remaining transferable load, the transferable capacity is transferred, otherwise the upper limit of the transferable capacity is transferred, the remaining transferable load is updated, and the process is repeated until the transferable load is zero or there is no transferable adjacent node; the target node number, the transfer amount, and the transfer path number of each transfer are recorded to form a progressive load transfer strategy table, which includes the source node number, the target node number, the transfer amount (kVA), the transfer path number, and the transfer priority sequence number, a total of 5 fields; the control instruction sequence is generated according to the table and sent to the intelligent terminal of the corresponding node through the GOOSE message to perform the load switching operation and complete the optimization scheduling of the photovoltaic power distribution network.

[0027] Step S1 includes the following steps:

[0028] Step S11: Obtain the historical scheduling state of the photovoltaic power distribution network;

[0029] Step S12: Perform high-load period analysis between different regional power distribution nodes on the historical scheduling state to obtain a high-load period scheduling state between different regional power distribution nodes;

[0030] Step S13: Transmission load power fluctuation analysis between different regional power distribution nodes is performed on the high load period scheduling state to obtain node load power fluctuation difference data.

[0031] Step S14: According to the historical scheduling state, transformer load state correlation extraction between different regional power distribution nodes is performed on the node load power fluctuation difference data to obtain node transformer load state.

[0032] In the embodiment of the application, complete operation data recorded every 5 minutes from January 1, 2023 to December 31, 2023 is extracted from a real-time database deployed in a regional power distribution automation master station, and the data fields include node number, time stamp, active power value (unit: kW), reactive power value (unit: kvar), three-phase voltage effective value (unit: V), three-phase current effective value (unit: A), circuit breaker opening and closing state (0 or 1), transformer oil temperature (unit: ℃), and load rate percentage, a total of 96 power distribution nodes covering 8 geographical partitions, each partition containing 12 nodes, and the total amount of original data is 10091520 records; missing value filling is performed on the original data, linear interpolation method is adopted, time interval is fixed at 300 seconds, and interpolation window length is set to 3 sampling points before and after; abnormal value elimination is performed on the current and voltage data, and the elimination standard is data points exceeding the mean value plus or minus 3 times the standard deviation, and the standard deviation is calculated based on a sliding window with a window length of 720 sampling points (6 hours); the processed data is sorted according to node number and time stamp to construct a two-dimensional time sequence matrix, the number of rows is 96 representing the number of nodes, and the number of columns is 105120 representing the total number of sampling points in a year; the output format is CSV file, the encoding is UTF-8, the field separator is comma, the time stamp format is YYYY-MM-DDHH:MM:SS, the numerical accuracy is kept to two decimal places after the decimal point, and the historical scheduling state is formed.

[0033] A threshold screening is performed on the load rate percentage field in the historical scheduling state output in step S11, and the threshold is set to 85%. When the load rate of a node is greater than 85% at four consecutive sampling points, that is, within 20 minutes, the time period is marked as a high-load candidate period. A regional aggregation analysis is performed on the candidate period, which is grouped according to the geographic partition number. The number of nodes in each partition that are in a high-load state at the same time is counted. If the number of nodes in a high-load state at the same time in the same partition exceeds 40% of the total number of nodes in the partition, that is, 5 nodes, the time period is marked as a regional high-load period. The duration of the regional high-load period is verified. The time period with a duration less than 30 minutes is removed, and the time period with a duration greater than or equal to 30 minutes is retained. The active power, reactive power, current value, voltage value, and oil temperature value of all nodes in the retained time period are extracted to form a high-load period scheduling state subset. A time alignment operation is performed on the subset. The time of the node entering the high-load state earliest is taken as the reference. The time is extended by 10 minutes forward and 10 minutes backward to form a high-load period window containing a boundary buffer. The total length of the window is fixed at 50 minutes. The start time, end time, involved node list, belonging region number, average load rate, maximum load rate, minimum load rate, and load rate variance of each high-load period window are output. There are 8 parameters in total, which are stored as a structured JSON array. The total number of array elements is determined according to the actual detection result. In the example, 217 high-load period windows that meet the conditions are detected throughout the year.

[0034] A sliding standard deviation calculation is performed on the active power data in each high-load period window output in step S12. The sliding window length is set to 10 sampling points, that is, 50 minutes, and the step is 1 sampling point. The standard deviation value of each node in the window every 5 minutes is calculated, with the unit being kW. A pairwise difference calculation is performed on the standard deviation values of all nodes in the same region. The difference is equal to the standard deviation of node A minus the standard deviation of node B. A power fluctuation difference matrix between nodes within the region is generated, with the matrix dimension being 12x12. The same operation is performed on the cross-region nodes to generate a cross-region fluctuation difference matrix, with the dimension being 96x96. A normalization processing is performed on all difference data. The normalization method is to subtract the global minimum value and then divide by the global maximum value minus the minimum value. The output range is compressed to the interval of 0 to 1. A threshold segmentation is performed on the normalized difference data. The threshold is set to 0.3. Node pairs with a difference greater than 0.3 are marked as significant fluctuation difference pairs. The original power sequence in the corresponding period of each significant fluctuation difference pair is extracted. A dynamic time warping (DTW) distance calculation is performed. The distance metric uses the Euclidean distance, and the path constraint uses the Sakoe-Chiba band. The bandwidth is set to 5 sampling points. The DTW distance value, power standard deviation difference value, belonging region combination, time window number, and node number pair of each significant fluctuation difference pair are output. There are 5 parameters in total, which constitute the node load power fluctuation difference data.

[0035] The node load power fluctuation difference data table output in step S13 is performed with an inner connection operation with the historical scheduling state output in step S11, the connection key being the node number and the time window number; a Pearson correlation coefficient calculation is performed on the connected data, the calculation object being the active power sequence of each pair of nodes in the corresponding high load period window, the sequence length being uniformly 10 sampling points; the node pairs with a correlation coefficient absolute value greater than 0.7 are retained as strong load correlation groups; a Granger causality test is performed on the nodes in the strong load correlation group, the lag order being fixed as 2, the significance level being set as 0.05, and only one-way or two-way causal relationship pairs that pass the test are retained; the number of strong correlation groups participated by each node, the average correlation coefficient, the maximum DTW distance, the minimum DTW distance, and the causal direction identifier (0 represents no causality, 1 represents a source node, and 2 represents a target node) are counted; meanwhile, the load rate peak value, the load rate valley value, the load rate average value, the load rate standard deviation, and the high load cumulative length (in minutes) of the node in the high load period of the whole year are extracted, totaling 10 characteristic parameters; the above parameters are arranged according to the node number, and output as a 96-row 10-column numerical matrix, each row of the matrix corresponding to a power distribution node, each column corresponding to a characteristic parameter, the data type being a floating point number, the precision being retained to four decimal places, the file format being HDF5, and the attribute field including the node number, the region number, and the data generation timestamp, thereby forming the final node transformer load state data.

[0036] Step S2 includes the following steps:

[0037] Step S21: Perform current abnormal state marking on the node transformer load state to generate a load abnormal current state.

[0038] Step S22: Perform output voltage instability quantification of the power distribution network transformer according to the load abnormal current state to obtain node output voltage instability data.

[0039] Step S23: Perform instability intensity increment derivation in the time dimension on the node output voltage instability data to generate time-dimension node instability intensity increment data.

[0040] Step S24: Perform clustering analysis on the node instability intensity increment data to obtain node instability intensity clustering data.

[0041] As an example of the present application, reference is made to FIG. 1, in which the step S2 in the present example includes: Figure 2

[0042] Step S21: Perform current abnormal state marking on the node transformer load state to generate a load abnormal current state.

[0043] ​In the embodiment of the present application, the three-phase current effective value sequence of each node in the high load period window is extracted from the node transformer load state data output in step S14, the sampling frequency is fixed at 10 Hz, the sequence length is dynamically adjusted according to the window length, the minimum length is 3000 sampling points corresponding to 5 minutes, fast Fourier transform (FFT) is performed on each phase current sequence, the transform point number is the next power of 2 of the sequence length, the frequency domain resolution is equal to the sampling frequency divided by the transform point number, and the fundamental component amplitude and the 3, 5, 7, 9, and 11 odd harmonic component amplitudes are extracted; the ratio of each odd harmonic amplitude to the fundamental amplitude is calculated, if any harmonic ratio exceeds 0.08, i.e. 8%, the sampling point is marked as a harmonic abnormal point; at the same time, sliding peak-to-peak value detection is performed on the original current sequence, the sliding window length is set to 100 sampling points, i.e. 10s, and the step length is 10 sampling points, the peak value is calculated by subtracting the minimum value from the maximum value in the window, if the peak value exceeds 1.5 times the rated current amplitude, i.e. the threshold is 865.5A when the rated current is 577A, it is marked as an overcurrent abnormal point; the harmonic abnormal point and the overcurrent abnormal point are subjected to logical or operation to generate a binary abnormality marking sequence, the sequence length is consistent with the original current sequence, the abnormal point is marked as 1 and the normal point is marked as 0; morphological closing operation is performed on the abnormality marking sequence, the structure element length is set to 50 sampling points, i.e. 5 seconds, for merging adjacent abnormal sections; and the load abnormal current state corresponding to each node is output.

[0044] In another embodiment, the transformer operation data of each distribution node in the photovoltaic distribution network is collected in each period, the collected operation data includes primary side current, secondary side current, load current amplitude and current phase angle, the current abnormal threshold range is set by comparing the difference between the rated current value and the actual current value, for example, when the rated current is 800A, if the amplitude of the actual current deviates from the rated current by more than ±10%, it is defined as an abnormal state, the current state of this period is marked as an abnormal current state by threshold judgment method, and the marking result is recorded in the form of binary data, wherein 0 represents normal current and 1 represents abnormal current. In actual operation, each sampling period is set to 1s, the load abnormal current state sequence of the node transformer is formed by sliding window detection on the collected data of 600s, and the sequence is used as the input basic data for subsequent instability quantification.

[0045] Step S22: quantifying the output voltage instability of the distribution network transformer according to the load abnormal current state to obtain node output voltage instability data;

[0046] In the embodiment of the present application, index positioning is performed on each node load abnormal current state sequence output in step S21, and the time stamps corresponding to all sampling points marked as 1 are extracted; the three-phase voltage instantaneous value sequence at the corresponding moment is extracted from the historical scheduling state according to the time stamp, the sampling frequency is 10 kHz, the extraction window length is fixed as 200 ms, that is, 2000 sampling points, and the center is aligned with the abnormal point moment; zero-crossing detection is performed on each voltage instantaneous value sequence, the positive half-cycle peak voltage and the negative half-cycle peak voltage are recorded, the absolute difference value is calculated and then divided by the standard peak voltage 311 V to obtain the waveform asymmetry parameter; slope calculation is performed on the voltage sequence, the sampling interval is 0.1 ms, the slope is equal to the voltage value of the next point minus the voltage value of the previous point and then divided by 0.0001 seconds, the slope variance of each peak point and the five sampling points before and after the peak point is calculated to obtain the peak slope variance parameter; at the same time, harmonic analysis is performed on the voltage sequence, FFT with Hanning window is adopted, the window length is 2000 points, the total harmonic distortion rates of 3, 5, 7, 9 and 11 times are extracted, and the calculation method is to square root the sum of the amplitudes of each harmonic and then divide by the amplitude of the fundamental wave; linear weighted fusion is performed on the waveform asymmetry, peak slope variance and total harmonic distortion rate three parameters, the weight coefficients are 0.4, 0.3 and 0.3 respectively, and the single-point voltage instability strength value is output; the instability strength values of all sampling points in the same abnormal period are subjected to arithmetic average to obtain the representative instability strength of the abnormal event; the representative instability strengths of all abnormal events of each node in a year are arranged in time sequence to form a time sequence, the sampling interval is the interval of abnormal events, the minimum interval is 1 s and the maximum interval is 86400 s; three times of spline interpolation is performed on the time sequence, the interpolation target frequency is 1 Hz, and continuous node output voltage instability data is generated.

[0047] In another embodiment, after obtaining the load abnormal current state of the node transformer, the winding current loss data and the core excitation current data corresponding to the time period are first extracted, which are compared with the node voltage sampling data, when continuous load abnormal current state is detected, the voltage fluctuation amplitude and the voltage distortion component of the time period are calculated, for example, under the condition that the voltage rating is 10 kV, when the voltage fluctuation range is detected to be between 9.2 kV and 10.5 kV, and the harmonic distortion component amplitude in this range is greater than 5% of the rated voltage, the node voltage is defined as instability state, the voltage instability quantization value is obtained by interval integration of the voltage fluctuation amplitude of the continuous period, the quantization value represents the instability degree in numerical form, for example, when the cumulative instability quantization value of a certain node reaches 450 kV·s in a 600 s period, the output voltage instability data of the node in the period is recorded, and finally the complete node output voltage instability data is generated.

[0048] Step S23: performing instability strength increment derivation in time dimension on the node output voltage instability data to generate node instability strength increment data in time dimension;

[0049] In the embodiment of the present application, the node output voltage instability data output in step S22 is subjected to sliding difference calculation, the sliding window length is set to 3600s, i.e. 1 hour, the step length is 60s, the hour-level increment value is calculated by subtracting the start point value from the end point value of each window; the first-order autoregressive processing is performed on the increment value sequence, the Yule-Walker equation is used to estimate the coefficient, the order is fixed to 1, and the residual sequence is calculated; the fractional order differential is performed on the residual sequence, the order is 0.7, the Grünwald-Letnikov discrete formula is used, the step length is 60 sampling points, i.e. 60 seconds, the coefficient table is pre-calculated and cached; the moving standard deviation calculation is performed on the differential result, the window length is 720 sampling points, i.e. 12 hours, the step length is 60 seconds, and the local fluctuation intensity is obtained; the original hour increment value, the residual value, the fractional order differential value and the local fluctuation intensity are spliced in columns to form a four-dimensional time sequence matrix, the number of rows of the matrix is 8760 corresponding to the number of hours in a year, and the number of columns is 4; the principal component analysis PCA is performed on the matrix, the principal components with cumulative variance contribution rate greater than 95% are retained, and usually the first two principal components are retained; the principal component scores are taken as the final instability intensity increment feature vector, and each node outputs a 8760x2 numerical matrix; the L2 norm calculation is performed on each row of the matrix to obtain the scalarized comprehensive instability intensity increment value, and the node instability intensity increment data in the time dimension are formed.

[0050] In another embodiment, after obtaining the node output voltage instability data, the instability intensity increment in the time dimension is derived, and in the operation process, first, the 600s period output voltage instability data sequence is divided into 60 subintervals, each subinterval has a length of 10s, the increment change range of the voltage instability quantized value is calculated in each subinterval, the instability intensity increment data between adjacent intervals is obtained through the difference calculation method, for example, the instability quantized value of a node in the first interval is 7.5kV·s, and the instability quantized value in the second interval is 10.2kV·s, then the instability intensity increment of the interval is 2.7kV·s, and the whole 600s sequence is processed in turn to form the node instability intensity increment sequence in the time dimension, which can reflect the instability change process of the node in different time intervals, thereby providing the time-sequenced instability intensity data for subsequent clustering analysis.

[0051] Step S24: performing clustering analysis on the node instability intensity increment data to obtain node instability intensity clustering data.

[0052] In the embodiment of the present application, the instability strength increment data of all 96 nodes output in step S23 is longitudinally spliced to form a 96x8760 feature matrix, each row of the matrix representing the comprehensive instability strength increment value of a node per hour throughout the year; the matrix is subjected to Z-score standardization, the mean and standard deviation of each column being calculated independently, the mean of each column being 0 and the standard deviation being 1 after conversion; the K-means clustering algorithm is initialized, the number of cluster centers K being set to 5, the initial center point being selected by the K-means++ algorithm, the maximum number of iterations being set to 500, and the convergence threshold being set to 0.0001; the clustering calculation is performed, the Euclidean distance being used as the distance measure, and the cluster center being updated as the mean of the samples in the cluster in each iteration; after clustering is completed, the cluster label of each node is output, the value range being 0 to 4, and the coordinates of the five cluster centers are output, each center being a 8760-dimensional vector; the number of nodes included in each cluster, the average distance of the samples in the cluster to the center, and the minimum distance between clusters are calculated for each cluster to form a clustering statistics table; the original feature matrix and the cluster label are horizontally spliced to add a column of cluster_label and saved as a clustering enhanced dataset; the cluster center matrix is output as an independent file centers_5clusters.csv with a dimension of 5x8760 and a data precision of six decimal places; and the final node instability strength clustering data is composed of the node number, cluster label, region number, annual hourly increment sequence, and cluster center number.

[0053] In another embodiment, after obtaining the node instability strength increment data in the time dimension, the data is subjected to clustering analysis. In the implementation process, the instability strength increment sequence of each node is first normalized to map the instability strength increment values of different nodes to a unified numerical interval of 0 to 1, facilitating clustering calculation. Then, the Euclidean distance measure method is used to calculate the similarity between the instability increment sequences of different nodes to construct a node instability distance matrix. Based on the matrix, a hierarchical clustering-based aggregation method is used for processing, gradually merging nodes with high similarity into a cluster category. For example, in a certain calculation, the similarity of the instability strength increment sequences of node A and node B in the 600s period reaches 0.92, and they are merged into the same category of nodes. Through gradual aggregation, a complete clustering structure is formed, and finally the node instability strength clustering data is obtained, which is used to guide the design of load transfer logic between adjacent nodes.

[0054] Step S22 includes the following steps:

[0055] Step S221: performing winding heat increment index calculation of the distribution network transformer according to the load abnormal current state to obtain winding heat increment index;

[0056] Step S222: performing peak heat accumulation coupling on the winding heat increment index to obtain winding heat accumulation data;

[0057] Step S223: quantifying the magnetic circuit saturation trend based on the winding heat accumulation data;

[0058] Step S224: analyzing the output voltage waveform distortion strength according to the magnetic circuit saturation trend to obtain the output voltage distortion strength;

[0059] Step S225: quantifying the output voltage instability of the power distribution network transformer based on the output voltage distortion strength to obtain node output voltage instability data.

[0060] In the embodiment of the present application, the original three-phase current effective value corresponding to all sampling points marked as 1 in the output load abnormal current state sequence in step S21 is extracted, the sampling frequency is 10 Hz, the unit is A, each phase is processed independently, and the maximum value in the three phases is taken as the calculation reference; copper loss calculation is performed on each abnormal sampling point, and the copper loss is equal to the square of the current effective value multiplied by the winding direct current resistance value, the direct current resistance value is fixed as 0.015 , according to the nominal parameters of the transformer model; the copper loss value is multiplied by the sampling interval 0.1 seconds to obtain the single-point instantaneous heat energy increment, the unit is J; sliding accumulation is performed on the heat energy increment, the sliding window length is set to 100 sampling points, i.e. 10s, the step is 1 sampling point, and the local heat accumulation sequence is generated; exponential weighted moving average (EWMA) is performed on the local heat accumulation sequence, the smoothing coefficient is set to 0.2, which is used to suppress high-frequency fluctuations; the winding heat increment index corresponding to each abnormal sampling point is output, the index value is equal to the EWMA smoothed heat accumulation value divided by the standard heat accumulation value 4500J within 10s under the rated load, and the dimensionless is realized; all abnormal points in a year are sorted according to the time stamp to form the winding heat increment index.

[0061] Local maximum value detection is performed on the winding heat increment index sequence output in step S221, the detection method is to compare the current point with each of the 5 sampling points before and after it, i.e. 11 points in total, if the current point is the maximum value and greater than the threshold value 0.8, it is marked as a peak point; for each peak point, the start point of the continuous increasing segment is traced back, and the end point of the continuous decreasing segment is traced back, to form a heat accumulation event interval; numerical integration is performed on all heat increment indexes in each event interval, the integration method adopts the trapezoidal rule, the step is fixed as 0.1s, the single heat accumulation total amount is obtained, the unit is dimensionless integral value; at the same time, the event duration is recorded, the unit is second, the minimum duration threshold is set as 3s, and the events below the threshold are removed; the remaining heat accumulation events are sorted according to the occurrence time to generate an event list, each event contains the starting sampling index, the ending sampling index, the peak sampling index, the heat accumulation total amount, the duration, and the peak heat index; time resampling is performed on the list, the sampling frequency is 1 Hz, and the method is to uniformly distribute the heat accumulation total amount of each event to each second within the event duration, and the time points not covered by the event are assigned as 0; the resampled winding heat accumulation data is output.

[0062] Perform sliding maximum filtering on the winding heat accumulation data output in step S222, sliding window length 600 seconds, step 1s, extract the local heat accumulation upper limit trend line; Perform first-order difference on the trend line to obtain the heat accumulation rate sequence, set the rate threshold to 0.005, when the rate of 10 consecutive sampling points is greater than 0.005, mark the entry into the magnetic circuit pre-saturation stage; Perform second-order difference on the data in the pre-saturation stage, if the second-order difference value of 5 consecutive points is greater than 0.0001, it is determined to enter the accelerated saturation interval; Perform start point positioning on the accelerated saturation interval, record the entry time and exit time, the exit condition is that the heat accumulation data value falls below 0.3 and the rate is less than 0; Extract the duration, peak heat accumulation value, average rate, and second-order difference maximum value of each saturation interval; At the same time, perform normalization on the original heat accumulation data within the saturation interval, the method is to subtract the interval minimum value and then divide by the interval maximum value minus the minimum value, to obtain a standardized heat accumulation curve in the range of 0 to 1; Perform slope fitting on the curve, fit a straight line using the least squares method, and output the slope value as the magnetic circuit saturation rate indicator; Arrange the 8 parameters of saturation interval number, start time, end time, duration, peak heat accumulation, average rate, second-order difference maximum value, and saturation rate slope in chronological order to form the magnetic circuit saturation trend.

[0063] According to the start time and end time in the magnetic circuit saturation trend output in step S223, extract the three-phase voltage instantaneous value sequence in the corresponding period from the historical scheduling state, the sampling frequency is 10 kHz, extract the window aligned saturation event boundary, and expand 5 ms buffer before and after; Perform zero-crossing detection on each phase voltage sequence, record the positive and negative half-cycle peak voltage amplitude, and calculate the asymmetry degree equal to the absolute difference between the positive and negative peak values divided by 311V; Perform harmonic analysis on the voltage sequence, use Blackman window FFT, window length 2048 points, extract 3, 5, 7, 9, 11 odd harmonic amplitudes, and calculate the total harmonic distortion rate THD equal to the square root of the sum of the harmonic amplitudes divided by the fundamental amplitude; Perform slope calculation on the voltage waveform, sampling interval 0.1ms, calculate the instantaneous slope for each peak point and the 10 sampling points before and after it, and then calculate the standard deviation of the slope in the interval as the waveform steepness indicator; Calculate the arithmetic mean of the above three indicators for all sampling points in each saturation event to obtain the event-level asymmetry mean, THD mean, and steepness mean; Perform weighted fusion on the three means, the weight coefficients are 0.5, 0.3, and 0.2 respectively, and output the single event output voltage distortion intensity value; Sort all saturation events in a year by time to form the output voltage distortion intensity.

[0064] The output voltage distortion intensity value corresponding to each magnetic circuit saturation event output in step S224 is mapped to a continuous time sequence in seconds according to its start and end time on the time axis, and the time range covers 8760 hours in a year. Each saturation event is assigned the same distortion intensity value at all second points within its duration period, and the time points without saturation are uniformly assigned a value of 0. A sliding maximum value filtering process is performed on the preliminary filled sequence, with a sliding window length of 300s and a step of 1s, to smooth the instantaneous distortion fluctuations and retain the local peak value characteristics. Then, a Gaussian weighted moving average is performed, with a kernel function standard deviation of 60s and a kernel total length of 360s. The convolution boundary uses mirror padding to ensure that the edge data is not distorted. A hard clipping operation is performed on the filtered sequence, with a lower threshold of 0.1 and an upper threshold of 0.9. All values less than 0.1 are replaced with 0, and all values greater than 0.9 are replaced with 0.9, to achieve dynamic range compression and abnormal value suppression. The processed sequence is segmented by hour, and the arithmetic mean of all second-level data within each hour is calculated to generate an output voltage instability intensity sequence at the hour granularity, with a sequence length of 8760, corresponding to a value for each hour of the year. Trend separation is performed on the hour sequence, and the residual sequence is obtained by subtracting the original value from the 12-hour sliding median. The residual sequence is then transformed into an absolute value to form a fluctuation-enhanced instability intensity sequence. Finally, two parallel sequences are output, the first being the original hour mean sequence and the second being the fluctuation-enhanced sequence, which together constitute the node output voltage instability data.

[0065] In another embodiment, after the abnormal current state of the transformer is collected in the photovoltaic power distribution network, the primary side current and the secondary side current are recorded according to the sampling period of 1s, and the deviation amplitude of the actual current from the rated current is calculated based on the historical rated current. The deviation amplitude is combined with the winding resistance parameter to obtain winding instantaneous heating power data. For example, the winding resistance of a certain transformer is 0.15Ω, when the detected current reaches 850A and the rated current is 800A, the current deviation is 50A, and the calculated heating power increment is 6375W. The heating power increment of each second is accumulated and converted into temperature rise data within 600s of continuous monitoring, and the temperature rise change process is fitted in a logarithmic curve manner, and finally the winding heat increment index is obtained. The index is recorded in dimensionless numerical form, which is used to represent the heat accumulation rate of the winding under the current abnormal state.

[0066] After obtaining the winding heat increment index, peak identification is performed on the index according to the time axis, the maximum heat increment index value in each sampling period is extracted, and the peak values of adjacent time periods are cumulatively calculated to form a peak heat accumulation curve to reflect the winding heat accumulation trend under long-period operation conditions. For example, in a 600s monitoring period, if the peak heat increment index of each 100s is 1.05, 1.12, 1.18, 1.15, 1.20 and 1.25 respectively, the peak heat accumulation data of 600s is obtained by step-by-step addition, which is 7.95. This data reflects the total heat accumulation of the winding under long-time abnormal current operation, and finally forms the winding heat accumulation data, which provides input conditions for subsequent quantification of the magnetic circuit saturation trend.

[0067] After obtaining the winding heat accumulation data, the data is coupled with the core permeability parameter for analysis. The core permeability is set to 2.3T·m / A. By comparing the heat accumulation value with the core saturation point flux threshold, when the winding heat accumulation data rapidly increases in a short time, the core internal magnetic resistance increases, causing the magnetic circuit to tend to be saturated. In order to quantify the magnetic circuit saturation trend, the heat accumulation data is subjected to differential processing to calculate its change rate per unit time, and is mapped in combination with the actual sampling data of the magnetic flux density. For example, when the winding heat accumulation data growth rate reaches 0.15 / 100s, the core magnetic flux density rises to 2.1T, it is determined that the node is in the magnetic circuit saturation critical interval, and the quantification result is recorded in the form of a saturation trend index to provide input data for voltage distortion strength analysis.

[0068] After obtaining the magnetic circuit saturation trend index, the excitation current waveform data of the corresponding period is extracted, and Fourier decomposition is performed to extract the fundamental and harmonic components. When the 3rd harmonic amplitude of the excitation current reaches 12% of the fundamental, and the 5th harmonic amplitude reaches 7% of the fundamental, it indicates that the magnetic circuit saturation causes significant waveform distortion. Further, peak cycle detection is performed on the waveform distortion, and by comparing the interval and amplitude variation of the waveform peak points, the cycle distortion rate is calculated. For example, if the peak amplitude difference of the repeated peaks detected in 10 fundamental cycles exceeds 15% of the fundamental peak value, it is recorded as the distortion intensity increment. Finally, the harmonic amplitude, peak cycle distortion rate and waveform asymmetry parameters are integrated to obtain the output voltage distortion intensity data, which represents the distortion degree in percentage form, for example, the output voltage distortion intensity is 8.5%. After obtaining the output voltage distortion intensity, the intensity is coupled with the actual sampling sequence of the output voltage, and the voltage distortion intensity and voltage deviation in each 10s sampling interval are weighted to form the voltage instability quantization value. For example, in a 10s interval, the voltage rated value is 10kV, the actual average voltage is 9.4kV, the voltage deviation is 6%, and the output voltage distortion intensity is 8.5%. By weighting the deviation and distortion intensity with weighting factors 0.6 and 0.4, the voltage instability quantization value of the interval is 7.1%. The voltage instability quantization values of all intervals in a 600s period are accumulated and sequenced to finally form the node output voltage instability data.

[0069] The output voltage waveform distortion intensity analysis includes:

[0070] According to the magnetic circuit saturation trend, the excitation current nonlinear distortion is analyzed to obtain excitation current distortion data;

[0071] The peak value multiple ratio of the excitation current distortion data is derived to obtain the peak value multiple ratio of the peak cycle;

[0072] Based on the peak value multiple ratio of the peak cycle, the odd harmonic frequency variance of the excitation distortion current is calculated from the excitation current distortion data to obtain the odd harmonic frequency variance;

[0073] According to the peak value multiple ratio of the peak cycle and the odd harmonic frequency variance, the multi-peak positive and negative half-waveform asymmetry of the induced electromotive force is analyzed from the excitation current distortion data to generate positive and negative half-waveform asymmetry data;

[0074] The slope variance of the peak points of the positive and negative half-waveform asymmetry data is calculated to obtain the slope variance of the peak points of the induced electromotive force waveform;

[0075] According to the peak value multiple ratio of the peak cycle, the odd harmonic frequency variance and the slope variance, the output voltage waveform distortion intensity is analyzed to obtain the output voltage distortion intensity.

[0076] In the embodiments of the present application, the excitation current nonlinear distortion analysis is performed according to the magnetic circuit saturation trend to obtain excitation current distortion data; the specific embodiments are as follows: according to the time boundary of the magnetic circuit saturation event, the instantaneous value of the three-phase excitation current on the high-voltage side of the transformer is extracted from the power distribution network real-time acquisition system, the sampling frequency is fixed at 20 kHz, and the extraction period covers the start and end time of each saturation event and is extended by 10 ms before and after, to ensure the waveform integrity; the phase-locked loop technology is used to accurately track the 50 Hz fundamental frequency for each phase current sequence to generate an ideal sinusoidal reference waveform with the same frequency and phase; the original current is subtracted from the reference waveform to obtain a residual sequence containing only nonlinear distortion components; the sliding mean filter is performed on the absolute value of the residual sequence, the window length is 200 sampling points, i.e. 10 ms, and the step length is 1 sampling point, to smooth the high-frequency burrs; the local maximum value is detected in the filtered sequence, the judgment condition is that the amplitude of the current sampling point is greater than that of the previous and subsequent 10 sampling points and exceeds the threshold value of 0.5 A, and the condition is met The distortion peak is marked; the time corresponding to each peak, the amplitude, the rising slope (linearly fitted by the previous 5 sampling points), the falling slope (linearly fitted by the subsequent 5 sampling points), and the half-width (the time span of the amplitude decay to half of the peak value) are recorded; all detected peaks are sorted by occurrence time to form excitation current distortion data, including phase identification, time stamp, amplitude, rising slope, falling slope, and half-width of six parameters; the original residual sequence is also retained as continuous waveform data, the sequence length dynamically changes according to the event duration, the shortest is not less than 2000 points corresponding to 100 ms, and the longest is not more than 60000 points corresponding to 3 s, for subsequent joint analysis of periodic structure and harmonic characteristics.

[0077] The peak multiple ratio of the peak repetition period of the excitation current distortion data is derived to obtain the peak multiple ratio of the peak period. The specific embodiment is: the peak events in the excitation current distortion data table output in the previous step are independently processed according to the phase. Taking the A phase as an example, the time stamp and amplitude of all A-phase peak events are extracted and arranged in ascending order of time. The time interval of adjacent peak events is calculated. If the interval is in the range of 18 ms to 22 ms, that is, close to the half cycle of power frequency 20 ms, it is determined as an event in the same cycle group. The peak events that continuously satisfy the interval condition are clustered into a cycle cluster, and the number of events in the cluster should not be less than 3. For each cycle cluster, all peak amplitudes contained therein are extracted, and the ratio of the maximum amplitude to the minimum amplitude is calculated, which is called the peak multiple ratio in the cycle. At the same time, the arithmetic mean value of the ratio of all adjacent peak amplitudes in the cluster is calculated, which is called the adjacent peak multiple ratio. The larger value of the two is taken as the final peak multiple ratio of the cycle cluster. The same operation is performed on all cycle clusters in a year to generate a peak period peak multiple ratio sequence. The total number of sequence elements is equal to the number of cycle clusters. In the example, 147 valid cycle clusters are detected for the A phase. The same process is repeated for the B phase and the C phase to generate their respective peak multiple ratio sequences. The three-phase results are combined and sorted according to the event occurrence time to form a time sequence of the peak period peak multiple ratio. Each element in the sequence is associated with the start time, end time, phase number, number of peaks in the cluster, maximum amplitude, minimum amplitude, and peak multiple ratio value of the cycle cluster to which it belongs. The sequence is used as the core intermediate parameter for voltage waveform distortion strength analysis to drive the subsequent odd harmonic frequency variance calculation.

[0078] The field current distortion data is compared with the peak multiple of the peak value of the spike cycle to obtain the odd harmonic frequency variance of the field distortion current. The specific embodiment is: for each original field current distortion residual sequence segment corresponding to each cycle cluster, perform fast Fourier transform (FFT) with Blackman window, the window length is fixed at 4096 points, the sampling frequency is 20 kHz, the frequency domain resolution is 4.8828 Hz, and the amplitude of 3, 5, 7, 9 and 11 odd harmonics is extracted; for each odd harmonic component, calculate the amplitude sequence in all cycles in the cluster, for example, the 3rd harmonic has 10 amplitude samples in 10 cycles; perform sample variance calculation on the amplitude sequence, and the degrees of freedom are n-1, to obtain the inter-cycle fluctuation variance of the harmonic; the variance values calculated for the 3rd, 5th, 7th, 9th and 11th harmonics are normalized by dividing the square of the historical average amplitude of the corresponding harmonic in the unsaturated period, to eliminate the influence of the inherent amplitude scale; the normalized five variance values are arranged in order of frequency to form a 5-dimensional variance vector; perform principal component projection on the vector, and the projection matrix is calculated in advance from historical training data, and the first principal component is retained as the comprehensive odd harmonic frequency variance index; at the same time, record the original five variance values and the projection coefficients for traceability analysis; perform product operation on the comprehensive variance index and the corresponding peak multiple of the peak value of the spike cycle in step S224 to strengthen the harmonic instability weight of high multiple events; output the final odd harmonic frequency variance value, each cycle cluster corresponds to a scalar value, a total of 441 values (three phases in total) are generated in the whole year example, the value range is 0 to 3.5, dimensionless, the time stamp is aligned with the center moment of the cycle cluster, to form a complete odd harmonic frequency variance.

[0079] According to the peak multiple ratio of the spike cycle and the variance of the odd harmonic frequency, the multi-peak positive and negative half-wave form asymmetric data is generated by analyzing the induced electromotive force of the excitation current distortion data. The specific embodiment is: for each corresponding original excitation current distortion residual sequence segment of the cycle cluster, the instantaneous value sequence of the output voltage of the low-voltage side of the transformer is extracted synchronously, the sampling frequency is 20 kHz, and the time window is strictly aligned with the current segment; zero-crossing detection is performed on the voltage sequence, the positive half cycle and the negative half cycle are divided based on the 50 Hz fundamental frequency, and the length of each half cycle is 10 ms, that is, 200 sampling points; the local maximum value points are searched in each half cycle interval, if there are multiple maximum values and the interval between adjacent maximum values is greater than 5 sampling points, that is, 0.25 ms, then it is marked as a multi-peak structure; for each multi-peak half cycle, record all peak point amplitudes, position indexes, rising edge maximum slope, and falling edge maximum slope; calculate the absolute difference value between the maximum peak value in the positive half cycle and the maximum peak value in the negative half cycle, and then divide by the theoretical peak value 311 V to obtain the basic asymmetry degree; at the same time, the absolute value of the difference between the number of peaks in the positive half cycle and the number of peaks in the negative half cycle is calculated, if the difference is greater than 1, an asymmetric penalty term is additionally added, and the penalty term is equal to the difference multiplied by 0.05; add the basic asymmetry degree and the penalty term, and then multiply by the geometric mean value of the spike cycle peak multiple ratio and the variance of the odd harmonic frequency corresponding to the cycle cluster to complete the nonlinear weighted modulation; the same operation is performed on all cycle clusters in a year, 89 multi-peak asymmetric half cycles are detected in the A phase example, 94 in the B phase, and 87 in the C phase; the time center point, the phase to which it belongs, the number of positive half cycle peaks, the number of negative half cycle peaks, the basic asymmetry degree, the penalty term, and the asymmetric intensity after modulation are output as seven parameters, which are sorted by time to form the positive and negative half-wave form asymmetric data.

[0080] The slope variance of the peak point of the positive and negative half-wave asymmetric data is calculated to obtain the slope variance of the peak point of the induced electromotive force waveform. The specific embodiment is: for each asymmetric event output in the previous step, backtrack the corresponding voltage instantaneous value waveform segment, and extract the data window of 10 sampling points before and after each identified peak point, i.e. 0.5 ms before and after; perform first-order difference on the 21 sampling points in the window to obtain 20 instantaneous slope values, with the unit of V / s; calculate the sample variance of the 20 slope values, with the degree of freedom being 19, to obtain the local slope fluctuation intensity of the peak point; if the half-wave contains multiple peak points, calculate the arithmetic mean of the slope variances of all peak points as the representative slope variance of the half-wave; at the same time, record the maximum single-point slope variance and the minimum single-point slope variance to depict the discrete nature of the waveform steepness distribution; multiply the half-wave representative slope variance by the peak cycle peak value multiple ratio associated with the event to strengthen the waveform steepness instability contribution of high multiple events; then perform weighted summation with the odd harmonic frequency variance, with the weight coefficients being 0.6 and 0.4, to generate a comprehensive slope distortion index; perform the same process on all asymmetric events throughout the year to output the timestamp, phase number, peak point number, representative slope variance, maximum slope variance, minimum slope variance, and comprehensive slope distortion index of each event; the parameter set is arranged according to the occurrence time to form the slope variance of the peak point of the induced electromotive force waveform, and the number of events in the example is 270 in total, with the value range to (V / s)², with the time label accuracy being 0.05 ms.

[0081] For each asymmetric event, the peak multiple ratio of the associated spike cycle, the odd harmonic frequency variance, and the comprehensive slope distortion index are normalized by subtracting the minimum value of the parameter throughout the year and dividing by the maximum value minus the minimum value, so that the three are mapped to the 0 to 1 interval; the three normalized parameters are fused by weighted fusion, and the weight coefficients are 0.4, 0.3, and 0.3, respectively, and the weighted sum is taken as the basic distortion intensity of the event; the nonlinear enhancement is performed on the basic distortion intensity, and the method is to take its 1.2 power to enlarge the discrimination of the high distortion region; the enhanced value is multiplied by the proportion factor of the event duration to the power cycle, and the proportion factor is equal to the actual half cycle length divided by 10ms, which compensates for the waveform truncation effect; the final output voltage distortion intensity value is output, each asymmetric event corresponds to a scalar, and there are 270 values in total for three phases throughout the year, with a value range of 0.05 to 0.98, and dimensionless; for ordinary cycles without multiple peak asymmetry, a value of 0.02 is assigned as the background distortion baseline; all events are expanded along the time axis, with the event center time as the anchor point, and linear interpolation is performed within ±5ms to fill in, generating a continuous time sequence; the sequence is filtered by 10ms sliding maximum value, and then smoothed by Gaussian smoothing standard deviation of 5ms, to obtain the smoothed second-level output voltage distortion intensity; the sampling interval is 1s, which is used as the direct input basis for transformer output voltage instability quantification.

[0082] In another embodiment, according to the magnetic circuit saturation trend, the specific embodiment of the excitation current nonlinear distortion analysis is as follows: after collecting the excitation current waveform data of the transformer in the photovoltaic power distribution network, the sampling frequency is set to 10000Hz to ensure accurate resolution of the fundamental wave and high-order harmonics. First, when the magnetic circuit saturation trend index exceeds the critical value of 0.85, the excitation current data of the corresponding period is intercepted. The waveform is frequency domain split by fast Fourier decomposition to extract the fundamental component and harmonic amplitude. The distortion data of the excitation current is recorded, wherein the fundamental frequency is 50Hz, the 3rd harmonic amplitude is 11% of the fundamental wave, the 5th harmonic amplitude is 7% of the fundamental wave, and the 7th harmonic amplitude is 3% of the fundamental wave. The existence of the above harmonic amplitudes indicates that the excitation current exhibits obvious nonlinear distortion under the magnetic circuit saturation trend. The distortion data is saved in the form of a time sequence matrix for subsequent spike cycle characteristic analysis.

[0083] In the embodiment of deriving the peak value ratio of the peak repetition period of the excitation current distortion data, the peak value points in each fundamental period are extracted from the distortion data, the peak amplitudes of adjacent periods are recorded, and the ratio between the peak amplitudes of each period is calculated to form a peak period peak value ratio sequence. For example, in 10 consecutive fundamental periods, the peak amplitudes are 420A, 470A, 450A, 495A, 480A, 500A, 470A, 510A, 495A, and 505A, respectively. The ratio sequence is 1.12, 0.96, 1.10, 0.97, 1.04, 0.94, 1.08, 0.97, 1.02, respectively. Through statistics of the ratio sequence, the peak period peak value ratio is 1.02, indicating the amplification trend of the current peak repetition. The ratio data is used as an input parameter for subsequent odd harmonic frequency variance calculation.

[0084] In the embodiment of calculating the odd harmonic frequency variance of the excitation current distortion data based on the peak period peak value ratio, the extracted odd harmonic components are frequency counted, and the 3rd, 5th, 7th, 9th, and 11th harmonic amplitudes are selected as the calculation objects. After normalization of the odd harmonic amplitudes, the variance is calculated to obtain the odd harmonic frequency variance value. For example, the normalized amplitudes are 0.11, 0.07, 0.03, 0.015, and 0.008, respectively. The result of 0.00162 is obtained by statistical variance. The variance value is associated with the peak period peak value ratio to represent the unevenness of the harmonic frequency distribution. The odd harmonic frequency variance result is transmitted as input data to the multi-peak waveform asymmetry analysis step.

[0085] In the embodiment of analyzing the multi-peak positive and negative half-waveform asymmetry of the induced electromotive force based on the peak period peak value ratio and the odd harmonic frequency variance, the peak amplitudes and durations of the positive and negative half-waves are extracted from the distortion current waveform, and the waveform differences are compared. When the average peak amplitude of the positive half-wave is 480A and the average peak amplitude of the negative half-wave is 450A, and the duration of the positive half-wave is 9.9ms and the duration of the negative half-wave is 10.3ms, the positive and negative half-wave asymmetry data are formed, and the asymmetry amplitude ratio is 1.07 and the asymmetry time ratio is 0.96. These asymmetry data are formed into a matrix as input parameters for the slope variance calculation step.

[0086] In a specific embodiment of slope variance calculation of peak points of positive and negative half-wave form asymmetric data, the current change rate in the range of 2 ms before and after the extraction of peak points in each half-wave is calculated to form a peak point slope sequence, for example, the slope values detected in the positive half-wave are 125 A / ms, 118 A / ms, and 130 A / ms, and the slope values detected in the negative half-wave are 110 A / ms, 105 A / ms, and 115 A / ms, the slope variances of the positive half-wave and the negative half-wave are calculated by comparison, the positive half-wave variance is 28, the negative half-wave variance is 27, and the comprehensive slope variance result of the peak points of the induced electromotive force waveform is obtained by weighting the two.

[0087] In a specific embodiment of output voltage waveform distortion strength analysis according to the peak value multiple ratio of the peak period, the odd harmonic frequency variance, and the slope variance, the three types of data are taken as multi-dimensional input quantities, and weighting coefficients 0.4, 0.35, and 0.25 are respectively assigned for weighted superposition to form a comprehensive distortion strength index, for example, the peak value multiple ratio of the peak period is 1.02, the odd harmonic frequency variance is 0.00162, and the comprehensive slope variance is 27.5, the output voltage waveform distortion strength is 8.7% obtained by weighted calculation, and the strength value is recorded as a quantitative index in the node output voltage distortion data sequence for subsequent instability quantitative analysis.

[0088] Step S23 includes the following steps:

[0089] Step S231: voltage instability curve drawing is performed based on the node output voltage instability data to construct a voltage instability curve.

[0090] Step S232: time sequence divergence incremental gradient derivation is performed on the voltage instability curve to obtain a time sequence divergence incremental gradient of the voltage instability state.

[0091] Step S233: fractional order numerical differentiation is performed on the time sequence divergence incremental gradient to obtain a divergence numerical fractional order.

[0092] Step S234: divergence trend prediction induction is performed on the divergence numerical fractional order by a self-recurrence model to obtain divergence trend prediction induction data.

[0093] Step S235: instability strength increment data in the time dimension of the node is generated by performing time dimension instability strength increment derivation according to the divergence trend prediction induction data.

[0094] As an example of the present application, referring to FIG. 1, in the present example, the step S23 includes: Figure 3

[0095] Step S231: voltage instability curve drawing is performed based on the node output voltage instability data to construct a voltage instability curve. ​

[0096] In the embodiment of the present application, the voltage instability curve is plotted based on the node output voltage instability data to construct the voltage instability curve; in a specific embodiment, the input is the recorded node output voltage instability data for the whole year, the sampling interval is 1 second, and the numerical range is 0 to 0.9; time axis resampling is performed on the sequence, the target frequency is adjusted to 1 minute granularity, and the method is to calculate the arithmetic mean value for every 60 consecutive second-level data points to generate a new sequence with a length of 52560, corresponding to one instability intensity value per minute for the whole year; the time stamp is uniformly formatted as YYYY-MM-DD HH:MM:00 to ensure that the integral is aligned; three times of spline interpolation are performed on the resampled sequence, the interpolation target is one point every 10s, an intermediate transition sequence is generated for smooth visual presentation, and the interpolation boundary condition is set as a natural boundary, i.e., the second derivative is zero; the original minute average points are superimposed on the interpolated sequence as anchor points to avoid overfitting; the sliding median filter is performed on the interpolated sequence, the window length is 5 minutes, i.e., 30 10s points, and the step length is 1 point to suppress local impulse noise; the final sequence is divided by day, there are 1440 10s points per day, a total of 365 segments, each segment is independently normalized to the interval of 0 to 1, the method is to subtract the daily minimum value and then divide by the daily maximum value minus the minimum value to retain the relative change trend within the day; 365 groups of normalized 10s granularity sequences are output, each group has 1440 values, and a voltage instability curve set is constructed; each curve represents the voltage instability evolution trajectory of a single day, the horizontal axis is the time scale from 00:00:00 to 23:59:50 with a step length of 10 seconds, and the vertical axis is the normalized instability intensity.

[0097] In another embodiment, the amplitude in the output voltage instability data at different time nodes is collected, the data sampling frequency is set to 1000 Hz, the collection duration is 600s, the data is stored as a time sequence vector set, the voltage instability curve is plotted in a Cartesian coordinate system, the abscissa is time t, and the ordinate is the voltage instability amplitude , and the voltage instability amplitude at each time is sequentially connected to form a continuous curve to construct the voltage instability curve, for example, in the time period of 0s to 600s, the node voltage instability amplitude gradually increases from 0.5V to 12.6V to form a nonlinear rising curve, and the voltage instability curve is used to provide a data basis for subsequent time sequence divergence incremental gradient derivation.

[0098] Step S232: performing time sequence divergence incremental gradient derivation on the voltage instability curve to obtain the time sequence divergence incremental gradient of the voltage instability state;

[0099] In the embodiment of the present application, for each daily voltage instability curve output in step S231, 1440 points of 10-second granularity sequence are taken, first-order forward difference is performed, the calculation method is the value of the nth+1 point minus the value of the nth point, a difference sequence with a length of 1439 and a unit of dimensionless per 10s is generated, sign separation is performed on the difference sequence, the positive value part is marked as positive divergence with the original value, and the negative value part is marked as negative convergence after taking the absolute value; sliding standard deviation calculation is performed on the positive divergence part, the window length is 12 points, i.e. 2 minutes, and the step length is 1 point, to obtain the local fluctuation intensity; the same operation is performed on the negative convergence part, to form a positive divergence fluctuation sequence and a negative convergence fluctuation sequence respectively; the original difference sequence, the positive fluctuation sequence and the negative fluctuation sequence are spliced according to the element position to form a three-dimensional gradient feature vector sequence with a length of 1439; principal component analysis is performed on the three-dimensional sequence, the first principal component is extracted as a comprehensive time sequence divergence increment gradient sequence, and information with an original variance contribution rate greater than 92% is retained; the original difference value, the positive fluctuation value, the negative fluctuation value and the principal component score corresponding to each time point are recorded; the same processing is repeated for each curve of 365 days in a year, and 365 groups of time sequence divergence increment gradients are output, each group has 1439 scalar values, and the time label is aligned to the starting time of every 10s, such as 00:00:10, 00:00:20 and the like; the acceleration or deceleration change trend of the voltage instability state in a short time scale is reflected, the value can be positive or negative, and the range is-0.15 to +0.18.

[0100] In another embodiment, the voltage instability curve data is discretized to take 1ms as a time step, to convert the continuous curve into 600000 discrete data points, and then the difference between adjacent data points is calculated to obtain the instantaneous change rate of the voltage instability curve, and the sliding window method is further used to set the window width to 50ms and the step length to 10ms, to smooth the difference data to eliminate high-frequency noise, so as to obtain the time sequence divergence increment gradient data of the voltage instability state, for example, the gradient at 300s is 0.026V / ms, and the gradient at 500s reaches 0.087V / ms, indicating that the voltage instability presents an accelerating divergence trend in the later period.

[0101] Step S233: performing fractional order numerical differentiation on the time sequence divergence increment gradient to obtain a divergence numerical fractional order;

[0102] In the embodiment of the present application, for each group of time sequence divergence increment gradients output in step S232, the length is 1439 and the sampling interval is 10s, Grünwald-Letnikov discretization formula is used to perform fractional order differentiation, the order α is fixed to 0.7, the number of coefficients is equal to the sequence length, and the recursive formula is , ; wherein For the gamma function, k increases from 1 to 1438; a convolution operation is performed on the sequence, and the boundary processing adopts a pre-zero padding method, and the padding length is equal to the sequence length minus 1; the convolution result sequence is output, and the length is still 1439, and each point represents the 0.7 order derivative value at the time; sliding maximum value filtering is performed on the derivative sequence, the window length is 6 points, that is, 60 seconds, the step length is 1 point, and the local extreme value characteristics are retained; exponential weighted moving average is further performed, the smoothing coefficient is set to 0.3, and high-frequency oscillation is suppressed; Z-score standardization is performed on the filtered sequence, the sequence mean is subtracted and then divided by the standard deviation, so that the output distribution mean is 0 and the standard deviation is 1; the standardized sequence is multiplied by a fixed scaling factor 0.85 to control the dynamic range; the final divergence numerical score is output, and there are 1439 points in each group, the numerical range is-2.1 to +1.9, and the time resolution is 10 seconds; 365 groups of sequences in a year are combined to form a complete annual divergence numerical score.

[0103] In another embodiment, the Grunwald-Letnikov fractional order difference method is used, the order is set to 0.75 order, the numerical discrete step length is set to 1ms, the time sequence divergence increment gradient sequence is processed by fractional order operation, and the gradient value at each time and its historical time sequence point are weighted and superimposed during operation. The weight coefficient is controlled by the binomial expansion coefficient, so that the approximation calculation of the fractional order derivative is realized, for example, the traditional integer order first derivative result is 0.045V / ms at 400s, and the result after 0.75 order fractional order numerical differentiation is 0.038V / ms. The processing makes the non-local characteristics of the divergence process be retained, and the generated divergence numerical fractional order sequence is used for subsequent trend prediction induction.

[0104] Step S234: performing divergence trend prediction induction on the divergence numerical fractional order by using an autoregressive model to obtain divergence trend prediction induction data;

[0105] In the embodiment of the application, the input is the annual divergence numerical fractional order output by step S233, which is divided into 365 groups per day, each group has 1439 points, and the sampling interval is 10s; stationary test is performed on each daily sequence, augmented Dickey-Fuller test ADF is used, the lag order is fixed at 12, if the p value is less than 0.05, the sequence is determined to be stationary, otherwise, first-order difference is performed until the stationary condition is met; the autoregressive model AR is fitted to the stationary sequence, the order p is determined by the PACF truncation position, the maximum search order is set to 24 corresponding to 4 minutes of historical dependence, and finally the order that makes the Akaike information criterion AIC minimum is selected as the optimal p; the autoregressive coefficients are solved by using the Yule-Walker equation, and no moving average term is introduced to ensure the pure autoregressive structure; the AR(p) model is independently trained for each sequence, for example, the p of a certain daily sequence is 8 after the test, the coefficients are to The last 60 points of the sequence, i.e., the last 10 minutes, are executed by rolling prediction using the coefficients obtained by training, the prediction step is 1, and the first p real values are used to predict the value at the next time point each time, and 60 predicted values are generated; the 60 predicted values are subtracted from the corresponding real values point by point to obtain a prediction residual sequence; the root mean square error RMSE, the mean absolute error MAE, and the maximum absolute error MAXAE of the residual sequence are calculated; the last predicted value of the AR model is extracted as the representative value of the divergence trend at the end of the day; the AR order, the 60-step prediction sequence, the three error indicators, and the trend representative value of each day, a total of 64 parameters, are packaged to form a single-day divergence trend prediction induction record; the same process is repeated for 365 days a year, and 365 records are output to form the divergence trend prediction induction data.

[0106] In another embodiment, the autoregressive AR method is used with an order of 5, the divergence numerical score sequence is divided into a training section and a verification section, the training section is from 0s to 480s, and the verification section is from 480s to 600s, the autoregressive coefficients are estimated by the least square method during the training process, the numerical score value at each time point is predicted and compared with the actual value, the mean square error is used as the convergence standard, when the mean square error is less than 0.01, the fitting is considered to be converged, the final prediction induction result is that the numerical score value gradually increases from 0.041 to 0.073 in the interval of 540s to 600s, and the trend prediction induction data accurately reflects the divergence acceleration process of voltage instability.

[0107] Step S235: According to the divergence trend prediction induction data, the instability strength increment data in the time dimension is derived, and the node instability strength increment data in the time dimension is generated.

[0108] In the embodiment of the present application, the 365 divergence trend prediction induction records output in step S234 are extracted, the last prediction value in each record is constituted to form a trend end point sequence with a length of 365, each value representing the predicted evolution direction of the voltage instability divergence state at the end time; a sliding difference is performed on the sequence, the window span is 7 days, the week increment value is calculated by subtracting the prediction value of the nth day from the prediction value of the n+7th day, and a week increment sequence with a length of 358 is generated; an exponential weighted moving average (EWMA) is performed on the week increment sequence, the smoothing coefficient is set to 0.4, and short-term fluctuations are suppressed; a second difference is performed on the smoothed sequence to obtain an acceleration change sequence for identifying trend turning points; the original week increment, the EWMA smoothed increment, and the acceleration change are spliced according to elements to form a three-dimensional feature vector sequence; principal component analysis (PCA) is performed on the three-dimensional sequence, the first principal component is retained, the cumulative variance contribution rate needs to be greater than 90%, and a scalarized comprehensive increment sequence is obtained after projection; the comprehensive increment sequence is normalized by subtracting the minimum value and dividing by the range, which is mapped to the interval of 0 to 1; the normalized sequence is expanded back to the original 365-day scale according to the day, the missing days are filled by linear interpolation, and the first 6 days and the last 1 day are filled by boundary replication; and the node instability strength increment data in the final time dimension is output.

[0109] In another embodiment, the prediction induction data is subjected to time integration processing to obtain a cumulative divergence strength sequence, the trapezoidal integration method is used for time integration, the integration step is 1 ms, and then the difference value of the cumulative divergence strength in adjacent time periods is calculated, which is defined as the instability strength increment data in the time dimension, for example, the cumulative divergence strength increases from 22.7 to 27.9 in the 550s-560s time period, the increment is 5.2, and the cumulative divergence strength increases from 35.4 to 46.1 in the 580s-590s time period, the increment is 10.7, indicating that the node instability strength presents an accelerated increment trend over time, and the finally generated time dimension instability strength increment data provides an input parameter for subsequent node instability clustering and scheduling logic design.

[0110] Step S3 includes the following steps:

[0111] Step S31: performing convolution processing on the node instability strength clustering data to obtain node instability convolution strength;

[0112] Step S32: performing load transfer logic design between adjacent nodes based on the node instability convolution strength, thereby reducing the load pressure of the current node, to obtain the load transfer logic for performing the optimal scheduling of the photovoltaic power distribution network.

[0113] In the embodiment of the present application, the node instability strength clustering data is convoluted to obtain node instability convolution strength. Specifically, the input is the 96 power distribution node annual hourly node instability strength clustering data output by step S24, each node sequence length is 8760, the sampling interval is 1 hour, and the value has been normalized to a mean of 0 and a standard deviation of 1. The 96 sequences are arranged according to the geographical topological adjacency relationship to construct a two-dimensional grid structure, the number of rows is 12 and the number of columns is 8, corresponding to 8 regions each with 12 nodes, physically adjacent nodes are adjacent in the grid, and boundary nodes are filled with zeros to form a 14x10 extended grid. One-dimensional convolution is independently performed on the 8760-hour sequence of each node, the convolution kernel adopts a discrete form of a Gaussian function, the kernel length is 15 hours, the standard deviation σ is 3 hours, the center weight is maximum and exponentially decays on both sides, and the kernel coefficients are pre-calculated and normalized to make the sum equal to 1. The convolution operation adopts an effective mode, the step is 1 hour, the boundary is not filled, and the output sequence length is 8746. The sliding maximum value filtering is performed on the sequence after convolution, the window length is 5 hours, and the local extreme value characteristics are retained. Then, Z-score standardization is performed to make each node output sequence meet the mean of 0 and the standard deviation of 1 again. The standardized sequence and the original clustering label are weighted and fused, the weight is the reciprocal normalized value of the clustering center distance, the weight is higher when the distance is closer, and the consistency of the same type of node trend is strengthened. 96 groups of node instability convolution strength sequences are output, each group has 8746 points, the time stamp is aligned from the 8th hour to the 8753rd hour, and covers the whole year of effective analysis period. The annual maximum value, the 95th percentile value, the daily average value and the weekly fluctuation standard deviation of each node sequence are extracted as static convolution strength features. The feature set and the dynamic sequence jointly constitute the complete expression of the node instability convolution strength, which is used to drive the adjacent node load transfer logic design.

[0114] The load transfer logic between adjacent nodes is designed based on the instability convolution intensity of the nodes, so as to reduce the load pressure of the current node and obtain the load transfer logic to perform the optimal scheduling of the photovoltaic power distribution network. Specifically, for each node output in step S31, the maximum value of the annual static convolution intensity feature is extracted as an emergency index, and if the value is greater than 1.8, the node is marked as a high-risk node. The current load rate of the high-risk node is calculated, and the method is to take the average load rate in the last 24 hours. If it exceeds 0.8, the load transfer process is started. Get the list of all physically connected adjacent nodes of the node, which is determined according to the power grid topology adjacency matrix, and the maximum number of adjacent nodes is not more than 6. The remaining carrying capacity of each adjacent node is calculated, which is equal to 0.8 times the rated capacity of the transformer 1000kVA minus the current average load rate times 1000kVA. At the same time, the available transmission capacity of the connecting line is calculated, which is the minimum value of the line thermal stability limit and the protection setting value, with the unit of kVA. Multiply the remaining carrying capacity by the available transmission capacity to obtain the upper limit of the transferable capacity. The adjacent nodes are arranged in ascending order according to the electrical distance, and the electrical distance is equal to the line resistance multiplied by 1.2 plus the reactance multiplied by 0.8. The sorted list is executed in a progressive distribution, starting from the first order node. If the upper limit of the transferable capacity is greater than the amount to be transferred, the full amount is transferred, otherwise the upper limit value is transferred and the remaining amount to be transferred is updated. The initial value of the amount to be transferred is equal to the current load rate of the high-risk node minus 0.75 times 1000kVA. Repeat the distribution process until the amount to be transferred is zero or there is no available adjacent node. Record the target node number, transfer amount, line number, and priority sequence of each transfer. Generate a load transfer instruction sequence, each instruction containing the source node ID, target node ID, load switching amount kVA, execution time, and duration in minutes. The instructions are sorted according to the priority and electrical distance, and the shortest path has the highest priority and is executed first. The instructions are packaged through the IEC61850 GOOSE message protocol and sent to the corresponding intelligent terminal controller to trigger the breaker opening and closing action. After completing a round of load transfer, update the real-time load rate and instability convolution intensity of each node, and enter the next round of evaluation cycle. The scheduling period is fixed at 15 minutes to ensure system dynamic balance.

[0115] Step S32 includes the following steps:

[0116] Step S321: Calculate the average load rate of the node transformer according to the instability convolution intensity of the node; obtain the configuration state and topological connection structure of the peripheral adjacent nodes of the current node;

[0117] Step S322: Identify the amount of load to be transferred for the transformer average load rate to obtain the amount of load to be transferred for the current node;

[0118] Step S323: Calculate the product of the carrying capacity difference between adjacent nodes and the line transmission capacity according to the topological connection structure and configuration state;

[0119] Step S324: According to the distance between the current node and the adjacent node, the product of the difference between the carrying capacity of the adjacent nodes and the line transmission capacity is progressively transferred by the load to be transferred, and a progressive load transfer strategy is obtained.

[0120] Step S325: Based on the progressive load transfer strategy, the load transfer logic design between adjacent nodes is carried out, so as to reduce the load pressure of the current node, and the load transfer logic is obtained.

[0121] In the embodiment of the application, for each power distribution node, the last 24 points in the 8746 data points corresponding to the last 24 hours in the node instability convolution intensity dynamic sequence output from step S31 are extracted, corresponding to one value per hour in the past 24 hours; the arithmetic mean of the 24 points is obtained to obtain the time average convolution intensity value; at the same time, the load rate data recorded every 5 minutes of the node in the same period is extracted from the historical scheduling database, a total of 288 points, and the arithmetic mean is obtained to obtain the node transformer average load rate, the value range is 0 to 1; if the average load rate is greater than 0.8 and the time average convolution intensity value is greater than 1.5, the node is determined as the object to be processed; the grid physical topology map of the node is queried to obtain the adjacent node number of all direct electrical connections, the maximum number of which is not more than 6; the configuration parameters of each adjacent node are extracted, including the transformer rated capacity unified to 1000kVA, the line positive sequence resistance unit Ω / km multiplied by the length km to obtain the total resistance value, the line positive sequence reactance unit Ω / km multiplied by the length km to obtain the total reactance value, the line thermal stability limit current unit A multiplied by the nominal voltage 10kV and then divided by the square root of 3 to obtain the thermal stability limit capacity unit kVA; at the same time, the current average load rate of the adjacent node is extracted, and the calculation method is the same as the first half of step S321; the adjacency structure table is constructed, each row of the table corresponds to an adjacent node, and contains five parameters of adjacent node number, resistance value, reactance value, thermal stability limit capacity and current average load rate; the table is arranged in ascending order of electrical distance, and the electrical distance is equal to the resistance value multiplied by 1.2 plus the reactance value multiplied by 0.8; the average load rate of the current node and the adjacency structure table are output as the basis input for subsequent transfer amount calculation and strategy generation.

[0122] ​For the node determined as the processing object in step S321, the average load rate is recorded as L_avg. If L_avg is greater than 0.8, the transition amount calculation process is started. The safety load threshold is set to 0.75, and the transition load amount ΔP is equal to L_avg minus 0.75 multiplied by the rated capacity of the transformer 1000 kVA, with the unit being kVA. For example, if a node L_avg = 0.86, then ΔP = (0.86-0.75) x 1000 = 110 kVA. Hard cut-off is performed on the transition load amount. If ΔP is less than 10 kVA, it is forced to be zero and no transition is performed. If ΔP is greater than 300 kVA, it is forced to be limited to 300 kVA to prevent excessive transition. At the same time, the maximum cut-off load capacity of the photovoltaic inverter connected to the current node is checked. If the total installed capacity of the inverter minus the current minimum technical output is less than ΔP, then ΔP is corrected to the difference. For example, if the total capacity of the photovoltaic inverter is 800 kVA and the minimum technical output is 550 kVA, then the maximum cut-off amount is 250 kVA. If the original ΔP = 300 kVA, it is corrected to 250 kVA. The final transition load amount ΔP is output, with the value range being 10 kVA to 300 kVA, the precision being retained to an integer, and the unit being kVA. This value is used as the total amount input of the progressive distribution algorithm, and the original calculation value and the correction reason are recorded for audit tracing. The transition load amount is recalculated at the start of each dispatching period, with the dispatching period being fixed at 900 seconds, i.e. 15 minutes, to ensure real-time response to load changes.

[0123] According to the topology connection structure and configuration state, a product of a difference in carrying capacity between adjacent nodes and a line transmission capacity is calculated; a specific embodiment is: for each adjacent node in the adjacency structure table output in step S321, a residual carrying capacity R_cap is calculated, which is equal to 0.8 times 1000 kVA minus a current average load rate times 1000 kVA, in units of kVA; if R_cap is less than 0, it is forced to be zero to represent no receiving capacity; at the same time, a thermal stability limit capacity T_lim of a connection line corresponding to the adjacent node is extracted, in units of kVA; a product C_prod of the difference in carrying capacity and the line transmission capacity is calculated, which is equal to R_cap times T_lim, in units of kVA²; for example, if the current average load rate of a certain adjacent node is 0.65, then R_cap = 0.8 x 1000 - 0.65 x 1000 = 150 kVA, and if the corresponding line T_lim = 800 kVA, then C_prod = 150 x 800 = 120000 kVA²; the product value is normalized, the method is to divide by the maximum C_prod value in all adjacent nodes, so that the output range is compressed to 0 to 1; at the same time, the original C_prod value is retained for physical constraint verification; two columns are added to the adjacency structure table, which are the original product value and the normalized product value; the adjacency table is rearranged in descending order of the normalized product value, and the nodes with high priority are arranged in the front row; if the normalized product values are the same, they are sorted in ascending order of the electrical distance; an updated adjacency structure table is output, each row of which contains the adjacent node number, the original product value, the normalized product value, the electrical distance, and the current average load rate; the table is used as a capacity constraint basis for progressive load transfer, to ensure that the transfer path meets both the node margin and the line thermal limit.

[0124] According to the distance between the current node and the adjacent node, the product of the difference of the carrying capacity and the line transmission capacity between the adjacent nodes is progressively transferred by the load to be transferred, and a progressive load transfer strategy is obtained. The specific embodiment is: input is the load to be transferred ΔP output by step S322 and the sorted adjacency list table output by step S323; initialize the remaining load to be transferred R_remain equal to ΔP; start processing from the first row of the adjacency list, for the adjacent node of the current row, take its original product value C_prod as the upper limit of the acceptable maximum transfer amount; if C_prod is greater than R_remain, then allocate the full amount of R_remain to the node, record the transfer amount equal to R_remain, and the target node is the current row node number; otherwise, allocate the full amount of C_prod to the node, record the transfer amount equal to C_prod, and update R_remain equal to R_remain minus C_prod; continue to process the next row until R_remain is zero or the adjacency list is traversed; if R_remain is still greater than 0 after the adjacency list is traversed, then start the second transfer, which allows transfer through intermediate nodes, the transfer path length is not more than 2 hops, and the transfer capacity is constrained by the minimum value of the two-stage line thermal limit; for each successful allocation record, add the transfer priority sequence number, which starts from 1 and increments according to the allocation order; at the same time, calculate the actual electrical path length, which is equal to the weighted sum of all segment resistance and reactance between the source node and the target node; output the transfer strategy table, each row containing the target node number, the transfer amount kVA, the priority sequence number, the path length Ω, and the direct connection identifier; for example, a certain node ΔP=110kVA, the first order node C_prod=90kVA, and the second order C_prod=60kVA, then the first node allocates 90kVA, the second node allocates 20kVA, and the priority is 1 and 2 respectively; the strategy table is the only basis for generating the final control instruction.

[0125] The load transfer logic between adjacent nodes is designed based on a progressive load transfer strategy, so as to reduce the load pressure of the current node, to obtain the load transfer logic; specific embodiments are: the transfer strategy table is arranged in ascending order of priority, priority 1 to 6 respectively corresponds to execution delay 0 seconds, 5 seconds, 10 seconds, 15 seconds, 20 seconds, 25 seconds, to avoid concurrent impact; each strategy generates a GOOSE control instruction, the content contains source node number, target node number, transfer amount kVA, execution delay, duration 900 seconds; the instruction encapsulation complies with the IEC61850-8-1 standard, the APPID is increased in priority from 0x4001, the multicast address is fixed as 01-0C-CD-01-00-01; the intelligent terminal drives the circuit breaker to act within 50ms after receiving the instruction, and the specified load is switched; after the action is completed, the remote signaling variable position is uploaded, the master station confirms and updates the real-time load rate; if a certain instruction fails, the next order node in the adjacency table is immediately redistributed, and the redistribution amount is equal to the original failure amount; after all operations are completed, the execution log is recorded, including instruction serial number, timestamp, actual transfer amount, target node load change value; the log is used to update the convolution intensity and clustering state before the next round of scheduling period, to realize closed-loop optimization.

[0126] The application further provides an optimal scheduling system of a photovoltaic power distribution network, which is used for executing the optimal scheduling method of the photovoltaic power distribution network.

[0127] The state correlation extraction module is used for acquiring a historical scheduling state of the photovoltaic power distribution network; the transformer load states between different regional power distribution nodes are correlated and extracted according to the historical scheduling state, to obtain node transformer load states;

[0128] The instability intensity increment derivation module is used for quantifying the output voltage instability of the transformer of the power distribution network according to the node transformer load state, to obtain node output voltage instability data; and the instability intensity increment of the node output voltage instability data is derived in the time dimension, to obtain node instability intensity clustering data.

[0129] The load transfer logic design module is used for designing the load transfer logic between adjacent nodes based on the node instability intensity clustering data, to reduce the load pressure of the current node, to obtain the load transfer logic, and to execute the optimal scheduling of the photovoltaic power distribution network.

[0130] The above is only a specific embodiment of the application, which enables those skilled in the art to understand or implement the application. Various modifications of these embodiments will be apparent to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the application. Therefore, the application will not be limited to these embodiments shown herein, but will conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for optimal dispatch of a photovoltaic power distribution network, characterized in that, The method comprises the following steps: Step S1: obtaining historical scheduling states of the photovoltaic power distribution network; extracting transformer load states between different regional power distribution nodes according to the historical scheduling states to obtain node transformer load states; Step S2: quantifying output voltage instability of the power distribution network transformer according to the node transformer load states to obtain node output voltage instability data; and deriving instability intensity increments in the time dimension of the node output voltage instability data to obtain node instability intensity clustering data; Step S3: designing load transfer logic between adjacent nodes based on the node instability intensity clustering data to reduce the load pressure of the current node to obtain the load transfer logic to perform the optimal scheduling of the photovoltaic power distribution network; wherein step S3 comprises the following steps: Step S31: performing convolution processing on the node instability intensity clustering data to obtain node instability convolution intensity; Step S32: designing load transfer logic between adjacent nodes based on the node instability convolution intensity to reduce the load pressure of the current node to obtain the load transfer logic to perform the optimal scheduling of the photovoltaic power distribution network; wherein step S32 comprises the following steps: Step S321: calculating the node transformer average load rate according to the node instability convolution intensity; obtaining the configuration state and topological connection structure of the surrounding adjacent nodes of the current node; Step S322: identifying the amount of load to be transferred according to the transformer average load rate to obtain the amount of load to be transferred of the current node; Step S323: calculating the product of the load carrying capacity difference and the line transmission capacity between adjacent nodes according to the topological connection structure and the configuration state; Step S324: performing progressive load transfer on the product of the load carrying capacity difference and the line transmission capacity between adjacent nodes according to the distance between the current node and the adjacent nodes through the amount of load to be transferred to obtain a progressive load transfer strategy; Step S325: designing load transfer logic between adjacent nodes based on the progressive load transfer strategy to reduce the load pressure of the current node to obtain the load transfer logic.

2. The optimal dispatching method of photovoltaic power distribution network according to claim 1, characterized in that, Step S1 comprises the following steps: Step S11: obtaining historical scheduling states of the photovoltaic power distribution network; Step S12: analyzing high load periods between different regional power distribution nodes for the historical scheduling states to obtain high load period scheduling states between different regional power distribution nodes; Step S13: analyzing transmission load power fluctuation between different regional power distribution nodes for the high load period scheduling states to obtain node load power fluctuation difference data; Step S14: extracting transformer load states between different regional power distribution nodes for the node load power fluctuation difference data according to the historical scheduling states to obtain node transformer load states.

3. The optimal dispatching method of photovoltaic power distribution network according to claim 1, characterized in that, Step S2 comprises the following steps: Step S21: marking current abnormal states for the node transformer load states to generate abnormal current states of the load; Step S22: quantifying output voltage instability of the power distribution network transformer according to the abnormal current states of the load to obtain node output voltage instability data; Step S23: time-dimension instability strength increment derivation is carried out on the node output voltage instability data, and node instability strength increment data in the time dimension is generated; Step S24: clustering analysis is carried out on the node instability strength increment data, so as to obtain node instability strength clustering data.

4. The optimal dispatching method of photovoltaic power distribution network according to claim 3, characterized in that, Step S22 includes the following steps: Step S221: winding heat increment index calculation of the distribution network transformer is carried out according to the abnormal current state of the load, and winding heat increment index is obtained; Step S222: peak heat accumulation coupling is carried out on the winding heat increment index, and winding heat accumulation data is obtained; Step S223: magnetic circuit saturation trend is quantified based on the winding heat accumulation data; Step S224: output voltage waveform distortion strength analysis is carried out according to the magnetic circuit saturation trend, so as to obtain output voltage distortion strength; Step S225: output voltage instability quantification of the distribution network transformer is carried out based on the output voltage distortion strength, and node output voltage instability data is obtained.

5. The optimal dispatching method of photovoltaic power distribution network according to claim 4, characterized in that, The output voltage waveform distortion strength analysis includes: According to the magnetic circuit saturation trend, excitation current nonlinear distortion analysis is carried out, and excitation current distortion data is obtained; Peak multiple ratio derivation is carried out on the excitation current distortion data between sharp peak repetition periods, and a sharp peak period peak multiple ratio is obtained; Based on the sharp peak period peak multiple ratio, odd harmonic frequency variance calculation of excitation distortion current is carried out on the excitation current distortion data, and odd harmonic frequency variance is obtained; According to the sharp peak period peak multiple ratio and the odd harmonic frequency variance, multi-peak positive and negative half-waveform asymmetry analysis of induced electromotive force is carried out on the excitation current distortion data, and positive and negative half-waveform asymmetry data is generated; Slope variance calculation of peak points is carried out on the positive and negative half-waveform asymmetry data, and slope variance of the induced electromotive force waveform peak point is obtained; According to the sharp peak period peak multiple ratio, the odd harmonic frequency variance and the slope variance, output voltage waveform distortion strength analysis is carried out, so as to obtain output voltage distortion strength.

6. The optimal dispatching method of photovoltaic power distribution network according to claim 3, characterized in that, Step S23 includes the following steps: Step S231: voltage instability curve drawing is carried out based on the node output voltage instability data, so as to construct a voltage instability curve; Step S232: time sequence divergence increment gradient derivation is carried out on the voltage instability curve, so as to obtain a time sequence divergence increment gradient of the voltage instability state; Step S233: fractional order numerical differentiation is carried out on the time sequence divergence increment gradient, so as to obtain a divergence numerical fractional order; Step S234: divergence trend prediction induction is carried out on the divergence numerical fractional order through an autoregressive model, so as to obtain divergence trend prediction induction data; Step S235: time-dimension instability strength increment derivation is carried out according to the divergence trend prediction induction data, and node instability strength increment data in the time dimension is generated.

7. An optimal dispatching system of a photovoltaic power distribution network, characterized in that, The photovoltaic power distribution network optimization scheduling method is used to execute the photovoltaic power distribution network optimization scheduling system as claimed in claim 1, and the photovoltaic power distribution network optimization scheduling system includes: A state association extraction module is configured to acquire historical scheduling states of the photovoltaic power distribution network; and perform transformer load state association extraction between different regional distribution nodes according to the historical scheduling states, so as to obtain node transformer load states. The instability strength increment derivation module is configured to quantify the output voltage instability of the transformer of the power distribution network according to the node transformer load state, to obtain node output voltage instability data; and to derive the instability strength increment in the time dimension for the node output voltage instability data, to obtain node instability strength clustering data. The load transfer logic design module is configured to design the load transfer logic between adjacent nodes based on the node instability strength clustering data, to reduce the load pressure of the current node, to obtain the load transfer logic, and to perform the optimal dispatching of the photovoltaic power distribution network.

Images