Optimized scheduling method and system for photovoltaic power distribution network

Abstract

The invention relates to the technical field of optimal scheduling of a power distribution network, in particular to an optimal scheduling method and system for a photovoltaic power distribution network. The method comprises the following steps: obtaining a historical scheduling state of the photovoltaic power distribution network, and extracting a transformer load state between regional nodes; then, based on a node transformer load state, output voltage instability of the node is quantified and analyzed, an instability strength increment is deduced, and clustering data of instability strength is obtained; and finally, designing load transfer logic of adjacent nodes based on the instability intensity clustering data, reducing the load pressure of the current node, and realizing optimal scheduling of the photovoltaic power distribution network. According to the method, the power distribution network optimization scheduling technology is improved, so that the power distribution network optimization scheduling technology is more perfect.

Description

Technical Field

[0001] This invention relates to the field of power distribution network optimization scheduling technology, and in particular to an optimization scheduling method and system for photovoltaic power distribution networks. Background Technology

[0002] Photovoltaic power generation is characterized by significant intermittency and volatility, with its output power greatly affected by weather, seasons, and sunshine conditions. This presents new challenges to the operation and dispatch of traditional distribution networks. Traditional distribution network dispatching methods typically assume relatively stable loads and power sources. However, frequent power fluctuations in photovoltaic distribution networks can lead to node voltage fluctuations, transformer load imbalances, and even local voltage instability or system instability. Furthermore, photovoltaic grid connection is often concentrated in specific areas, resulting in excessive load pressure on some distribution nodes while other nodes have relatively light loads, causing uneven distribution of grid resources and potential overload risks. However, traditional optimization dispatching methods for photovoltaic distribution networks suffer from inaccurate analysis of transformer node output instability, leading to large errors in the optimal dispatching of the distribution network. Summary of the Invention

[0003] Therefore, it is necessary to provide an optimized scheduling method and system for photovoltaic power distribution networks to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objectives, an optimized scheduling method for a photovoltaic distribution network is provided, the method comprising the following steps: Step S1: Obtain the historical dispatch status of the photovoltaic distribution network; based on the historical dispatch status, perform transformer load status correlation extraction between distribution nodes in different regions to obtain the node transformer load status; Step S2: Quantify the output voltage instability of the distribution network transformers based on the load status of the node transformers to obtain node output voltage instability data; perform time-dimensional instability intensity increment derivation on the node output voltage instability data to obtain node instability intensity clustering data; Step S3: Design load transfer logic between adjacent nodes based on node instability intensity clustering data to reduce the load pressure on the current node and obtain load transfer logic to perform optimized scheduling of photovoltaic power distribution network.

[0005] Preferably, the present invention also provides an optimized scheduling system for a photovoltaic distribution network, used to execute the optimized scheduling method for a photovoltaic distribution network as described above, the optimized scheduling system for the photovoltaic distribution network comprising: The status association extraction module is used to obtain the historical scheduling status of the photovoltaic distribution network; based on the historical scheduling status, it performs transformer load status association extraction between distribution nodes in different regions to obtain the node transformer load status; The instability intensity incremental derivation module is used to quantify the output voltage instability of the distribution network transformer based on the load status of the node transformer to obtain node output voltage instability data; and to perform time-dimensional instability intensity incremental derivation on the node output voltage instability data to obtain node instability intensity clustering data. The load transfer logic design module is used to design load transfer logic between adjacent nodes based on node instability intensity clustering data, thereby reducing the load pressure on the current node and obtaining load transfer logic to perform optimized scheduling of the photovoltaic power distribution network.

[0006] The beneficial effects of this invention are as follows: By acquiring the historical dispatch status of the photovoltaic distribution network, a comprehensive understanding of the network's operation and transformer load can be obtained. By correlating and extracting the load status between distribution nodes in different regions, the load change trends and potential problems of each node can be identified. This process helps reveal areas of load imbalance, predict which nodes are at risk of overload in advance, and provide data support for subsequent optimized dispatch. Through historical data mining, load status analysis can be performed more accurately, improving the precision and efficiency of distribution network dispatch. By quantifying the output voltage instability status of node transformers, it is possible to intuitively identify which nodes experience voltage fluctuations or instability under specific load conditions. Furthermore, by deriving the instability intensity increment over a time dimension, the variation pattern of voltage instability over different time periods can be obtained, thus providing a basis for cluster analysis of instability intensity. Analysis of node output voltage instability data enables refined monitoring of the distribution network, early detection of potential voltage instability risks, and the implementation of corresponding measures to prevent large-scale voltage fluctuations or collapses, thereby improving the stability and security of the power grid. Based on the node instability intensity clustering data, load transfer logic between adjacent nodes can be designed to achieve optimized load dispatch in the distribution network. By rationally transferring loads, the load pressure on certain nodes can be reduced, avoiding voltage instability or equipment damage caused by overload. The load transfer logic design not only balances the load distribution among nodes but also improves the overall operating efficiency of the distribution network and reduces the probability of grid faults. This step, through optimized scheduling strategies, enhances the reliability, flexibility, and scheduling efficiency of the photovoltaic distribution network, thereby ensuring the stability of power supply and supporting the efficient integration of renewable energy. Therefore, this invention is an improvement on a traditional photovoltaic distribution network optimization scheduling method. It solves the problem of inaccurate transformer node output instability analysis, which leads to large optimization scheduling errors in traditional methods. This invention improves the accuracy of transformer node output instability analysis and reduces the optimization scheduling error of the distribution network. Attached Figure Description

[0007] Figure 1 This is a flowchart illustrating the steps of an optimized scheduling method for a photovoltaic power distribution network. Figure 2 for Figure 1 A detailed flowchart illustrating the implementation steps of step S2. Figure 3 for Figure 2 A detailed flowchart illustrating the implementation steps of step S23. Detailed Implementation

[0008] Please see Figures 1 to 3 An optimized scheduling method for a photovoltaic distribution network, the method comprising the following steps: Step S1: Obtain the historical dispatch status of the photovoltaic distribution network; based on the historical dispatch status, perform transformer load status correlation extraction between distribution nodes in different regions to obtain the node transformer load status; Step S2: Quantify the output voltage instability of the distribution network transformers based on the load status of the node transformers to obtain node output voltage instability data; perform time-dimensional instability intensity increment derivation on the node output voltage instability data to obtain node instability intensity clustering data; Step S3: Design load transfer logic between adjacent nodes based on node instability intensity clustering data to reduce the load pressure on the current node and obtain load transfer logic to perform optimized scheduling of photovoltaic power distribution network.

[0009] In this embodiment of the invention, reference Figure 1 The above is a flowchart illustrating the steps of an optimized scheduling method for a photovoltaic distribution network according to the present invention. In this example, the optimized scheduling method for the photovoltaic distribution network includes the following steps: Step S1: Obtain the historical dispatch status of the photovoltaic distribution network; based on the historical dispatch status, perform transformer load status correlation extraction between distribution nodes in different regions to obtain the node transformer load status; In this embodiment of the invention, active power data, reactive power data, voltage amplitude data, current amplitude data, and switch status data of each distribution node, recorded every 5 minutes for the past 365 days, are exported from the SCADA system of the photovoltaic distribution network to form the original historical dispatch status dataset. This dataset is divided into 8 sub-regions according to geographical area, each sub-region containing 12 distribution nodes, for a total of 96 nodes. A sliding time window analysis is performed on the original data, with the time window length set to 4 hours and the step size to 30 minutes. Within 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, with the rated capacity uniformly set to 1000kVA. At the same time, the standard deviation of power fluctuation of the transmission lines between adjacent nodes is calculated, with the standard deviation threshold set to 15kW. When a node is in a continuous If the load rate exceeds 85% for three consecutive time windows and the standard deviation of power fluctuation with adjacent nodes is greater than 15kW, it is marked as a high-load associated period. The node current sequences corresponding to all high-load associated periods are extracted, with a sampling frequency of 10Hz and a duration of 2 hours. A Pearson correlation coefficient matrix is ​​calculated on the current sequences, with a matrix dimension of 96×96 and a correlation coefficient threshold set to 0.72. Node pairs with a correlation coefficient greater than 0.72 are retained as strong-association groups. A load status synchronization test is performed on the nodes within the strong-association groups, using the Dynamic Time Warping (DTW) algorithm with a distance threshold set to 0.35. Finally, the transformer load status of each of the 96 nodes is output, including five parameters: peak load rate, mean load rate, load fluctuation variance, associated node number, and association strength coefficient.

[0010] Step S2: Quantify the output voltage instability of the distribution network transformers based on the load status of the node transformers to obtain node output voltage instability data; perform time-dimensional instability intensity increment derivation on the node output voltage instability data to obtain node instability intensity clustering data; In this embodiment of the invention, a Fast Fourier Transform (FFT) is performed on the current sequence in the load state of each node transformer, with 7200 sampling points and a frequency domain resolution of 0.0139Hz. 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 fundamental amplitude, the current at that moment is marked as an abnormal state, and a binary abnormality marker sequence is generated, with the sequence length equal to the length of the original time series. For each current sample marked as abnormal, a winding heat increment exponential calculation is performed. The calculation method is to multiply the square of the effective current value by the copper loss coefficient of 0.0012, and then multiply by the duration of 0.2 seconds to obtain... The process involves: obtaining the single-point heat increment index; performing cumulative integration on consecutive abnormal points to form a winding heat accumulation data sequence; performing first-order difference on the heat accumulation data to obtain the heat accumulation change rate; determining that the magnetic circuit has entered the saturation trend range when the change rate is greater than 0.8 W / s for 5 consecutive sampling points; sampling the voltage waveform within the saturation range at a sampling frequency of 10 kHz and a sampling length of 200 ms; performing zero-crossing detection on the sampled waveform, recording the peak voltage values ​​of the positive and negative half-cycles, calculating their absolute difference divided by the theoretical peak voltage of 311 V to obtain the waveform asymmetry; and simultaneously calculating the slope of the waveform at a sampling point interval of 0.1 ms, calculating the voltage difference between adjacent points divided by the time difference to obtain the slope. Instantaneous slope; the variance of the slope at 5 sampling points before and after the peak point is calculated to obtain the slope variance; weighted fusion is performed by combining asymmetry and slope variance, with weight coefficients of 0.6 and 0.4 respectively, to output the output voltage distortion intensity sequence; a moving average filter is applied to the distortion intensity sequence with a window length of 50ms and a step size of 10ms to obtain smoothed node output voltage instability data; a time series is constructed from this data, divided into time periods of 1 hour, and the time-series divergence incremental gradient is calculated for each time period by subtracting the value of the previous time period from the value of the next time period to obtain the gradient sequence; fractional derivatives are performed on the gradient sequence with an order of 0.7, using Grünwald-L... The etnikov discretization formula, with a step size of 10 sampling points, yields a divergent numerical fractional-order sequence. Autoregressive fitting is performed on this sequence, with a fixed order of 3. The coefficients are solved using the Yule-Walker equation to obtain a predicted value sequence for the next 3 hours. The predicted value sequence is subtracted from the current value sequence to obtain incremental data on node instability intensity over time. K-means clustering is performed on this incremental data, with 5 cluster centers. The initial cluster centers are determined using the elbow rule, the maximum number of iterations is set to 100, and the convergence threshold is set to 0.001. The instability intensity category label and cluster center coordinates of each node are output, forming the node instability intensity clustering data.

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

[0012] In this embodiment of the invention, a one-dimensional convolution operation is performed on the node instability intensity clustering data output in step S2. The convolution kernel length is set to 7, the weight distribution uses a Gaussian function with a standard deviation of 1.2, a stride of 1, and zero padding. This outputs the node instability convolution intensity values ​​for each of the 96 nodes. The convolution intensity value for each node is normalized, with a normalization range of 0 to 1, and used as the node urgency weight. The average transformer load rate for each node is calculated by statistically analyzing load rate data every 5 minutes over the past 24 hours, totaling 288 points, and taking the arithmetic mean. Average load rate; nodes with an average load rate greater than 0.8 are marked as source nodes to be transferred; the configuration status of neighboring nodes of the source node to be transferred is obtained, including transformer rated capacity, line impedance value, and topology connection structure, represented by an adjacency matrix, where a matrix element of 1 indicates a physical direct connection and 0 indicates no direct connection; the load to be transferred is calculated for each source node by subtracting a threshold of 0.75 from the current load rate and then multiplying by the transformer rated capacity of 1000kVA; the carrying capacity difference is calculated for each adjacent node by multiplying the threshold of 0.8 by its rated capacity. Subtract its current average load rate multiplied by its rated capacity, then multiply by the line transmission capacity, where the line transmission capacity is the line's thermal stability limit value, in kVA; calculate the product of the carrying capacity difference and the line transmission capacity as the upper limit of the transferable capacity; sort the nodes in ascending order according to the electrical distance between the source node to be transferred and its adjacent nodes, which is equal to the line resistance value multiplied by 1.2 plus the reactance value multiplied by 0.8; allocate the load to be transferred to the adjacent nodes in sequence, with the allocation rule being: if the upper limit of the transferable capacity of the current adjacent node is greater than the remaining load to be transferred, then all of it is transferred. Otherwise, the maximum transferable capacity is transferred, the remaining load to be transferred is updated, and this process is repeated until the load to be transferred is zero or there are no adjacent nodes to transfer. The target node number, transfer amount, and transfer path number for each transfer are recorded to form a progressive load transfer strategy table. This table contains five fields: source node number, target node number, transfer amount (kVA), transfer path number, and transfer priority number. Based on this table, a control command sequence is generated and sent to the smart terminal of the corresponding node via a GOOSE message to execute the load switching operation and complete the optimized scheduling of the photovoltaic distribution network.

[0013] Step S1 includes the following steps: Step S11: Obtain the historical dispatch status of the photovoltaic distribution network; Step S12: Analyze the high-load periods between distribution nodes in different regions based on the historical scheduling status to obtain the high-load period scheduling status between distribution nodes in different regions. Step S13: Analyze the transmission load power fluctuation between different regional power distribution nodes during the high load period scheduling state to obtain node load power fluctuation difference data; Step S14: Based on the historical scheduling status, extract the transformer load status correlation between different regional distribution nodes by analyzing the difference data of node load power fluctuations, and obtain the node transformer load status.

[0014] In this embodiment of the invention, complete operational data recorded every 5 minutes from January 1, 2023 to December 31, 2023 is extracted from the real-time database deployed at the regional distribution automation master station. Data fields include node number, timestamp, active power value (kW), reactive power value (kvar), three-phase voltage RMS value (V), three-phase current RMS value (A), circuit breaker opening / closing status (0 or 1), transformer oil temperature (°C), and load rate percentage. A total of 96 distribution nodes cover 8 geographical regions, with each region containing 12 nodes. The total original data volume is 10,091,520 records. Missing value imputation is performed on the original data using linear interpolation. The interpolation method uses a fixed time interval of 300 seconds and an interpolation window length of 3 sampling points before and after the data. Outlier removal is performed on current and voltage data, with the removal criterion being data points exceeding the mean plus or minus 3 times the standard deviation. The standard deviation is calculated based on a sliding window with a window length of 720 sampling points, or 6 hours. The processed data is sorted by node number and timestamp to construct a two-dimensional time series matrix, with 96 rows representing the number of nodes and 105120 columns representing the number of sampling points throughout the year. The output format is a CSV file, encoded in UTF-8, with commas as field separators, and timestamp format YYYY-MM-DDHH:MM:SS. Numerical precision is retained to two decimal places, forming a historical scheduling status.

[0015] A threshold filter is applied to the load percentage field in the historical scheduling status output in step S11. The threshold is set to 85%. When a node's load percentage is greater than 85% for four consecutive sampling points (20 minutes), that time period is marked as a high-load candidate period. Regional aggregation analysis is performed on the candidate periods, grouping them by geographical partition number. The number of nodes simultaneously in a high-load state within each partition is counted. If the number of simultaneously high-load nodes in the same partition exceeds 40% of the total number of nodes in that partition (5 nodes), that time period is marked as a regional high-load period. Duration verification is performed on the regional high-load periods. Periods with a duration less than 30 minutes are removed, and periods with a duration of 30 minutes or more are retained. All nodes within the retained periods are... Extract the corresponding active power, reactive power, current value, voltage value, and oil temperature value to form a subset of high-load period scheduling states. Perform time alignment operation on this subset, using the earliest node time to enter the high-load state as the benchmark, extend it forward by 10 minutes and backward by 10 minutes to form a high-load period window including boundary buffers. The total length of the window is fixed at 50 minutes. Output each high-load period window with 8 parameters: start time, end time, list of involved nodes, region number, average load rate, maximum load rate, minimum load rate, and load rate variance. Store these parameters as a structured JSON array. The total number of array elements is determined based on the actual detection results. In the example, a total of 217 high-load period windows that meet the criteria were detected throughout the year.

[0016] For each high-load period window output in step S12, the moving standard deviation is calculated using the sliding window length of 10 sampling points (50 minutes) and the step size of 1 sampling point. The standard deviation of each node within the window is calculated every 5 minutes, in kW. For all nodes within the same region, pairwise difference calculations are performed, where the difference equals the standard deviation of node A minus the standard deviation of node B, generating a power fluctuation difference matrix between nodes within the region with a dimension of 12×12. The same operation is performed on nodes across regions to generate a cross-region fluctuation difference matrix with a dimension of 96×96. All difference data are normalized using the method of subtracting the global minimum. The value is then divided by the global maximum value minus the minimum value, compressing the output range to the interval between 0 and 1. Threshold segmentation is performed on the normalized difference data, with the threshold set at 0.3. Node pairs with differences greater than 0.3 are marked as significant fluctuation difference pairs. For each significant fluctuation difference pair, the original power sequence within its corresponding time period is extracted, and dynamic time warping (DTW) distance calculation is performed. The distance metric is Euclidean distance, the path constraint is Sakoe-Chiba band, and the bandwidth is set to 5 sampling points. The output of each significant fluctuation difference pair includes 5 parameters: DTW distance value, power standard deviation difference value, region combination, time window number, and node number pair, which constitute the node load power fluctuation difference data.

[0017] An inner join operation is performed on the node load power fluctuation difference data table output in step S13 and the historical scheduling status output in step S11, with the join key being the node number and the time window number; Pearson correlation coefficient is calculated on the joined data, with the calculation object being the active power sequence of each pair of nodes within the corresponding high load time window, and the sequence length being uniformly 10 sampling points; Node pairs with an absolute correlation coefficient greater than 0.7 are retained as strong load association groups; Granger causality tests are performed on the nodes within the strong load association groups, with the lag order fixed at 2 and the significance level set at 0.05, retaining only one-way or two-way causal relationship pairs that pass the test; For each node, the number of strong association groups it participates in, the average correlation coefficient, and the... The system extracts the maximum DTW distance, minimum DTW distance, and causal direction identifier (0 indicates no causality, 1 indicates source node, and 2 indicates target node). It also extracts 10 feature parameters for each node during peak load periods throughout the year: peak load rate, valley load rate, average load rate, standard deviation of load rate, and cumulative high load duration (in minutes). These parameters are arranged by node number, outputting a 96-row, 10-column numerical matrix. Each row corresponds to a distribution node, and each column corresponds to a feature parameter. The data type is floating-point, with precision retained to four decimal places. The file format is HDF5, and the attribute fields include node number, region number, and data generation timestamp, forming the final node transformer load status data.

[0018] Step S2 includes the following steps: Step S21: Mark the abnormal current status of the node transformer load and generate the abnormal load current status; Step S22: Quantify the output voltage instability of the distribution network transformer based on the abnormal load current state to obtain node output voltage instability data; Step S23: Perform time-dimensional instability intensity increment derivation on the node output voltage instability data to generate time-dimensional node instability intensity increment data; Step S24: Perform cluster analysis on the incremental data of node instability intensity to obtain node instability intensity cluster data.

[0019] As an example of the present invention, reference is made to... Figure 2 As shown, in this example, step S2 includes: Step S21: Mark the abnormal current status of the node transformer load and generate the abnormal load current status; In this embodiment of the invention, the effective value sequence of the three-phase current of each node within the high-load period window is extracted from the node transformer load status data output in step S14. The sampling frequency is fixed at 10Hz, and the sequence length is dynamically adjusted according to the window duration, with a minimum length of 3000 sampling points corresponding to 5 minutes. A Fast Fourier Transform (FFT) is performed on each phase current sequence, with the number of transformation points taken as the next power of 2 of the sequence length. The frequency domain resolution is equal to the sampling frequency divided by the number of transformation points. The amplitude of the fundamental component and the amplitudes of the 3rd, 5th, 7th, 9th, and 11th odd harmonic components are extracted. The ratio of each odd harmonic amplitude to the fundamental amplitude is calculated. If any harmonic ratio exceeds 0.08 (8%), the sampling point is marked as a harmonic anomaly point. A sliding peak-to-peak detection is performed on the original current sequence. The sliding window length is set to 100 sampling points (10 seconds), and the step size is 10 sampling points. The peak value is obtained by subtracting the minimum value from the maximum value within the window. If the peak value exceeds 1.5 times the rated current amplitude (for example, the threshold is 865.5A when the rated current is 577A), it is marked as an overcurrent anomaly. A logical OR operation is performed on the harmonic anomaly point and the overcurrent anomaly point to generate a binary anomaly marker sequence. The sequence length is the same as the original current sequence. Anomaly points are marked as 1, and normal points are marked as 0. A morphological closing operation is performed on the anomaly marker sequence. The length of the structuring element is set to 50 sampling points (5 seconds) to merge adjacent anomaly segments. The load anomaly current status corresponding to each node is output.

[0020] In another embodiment, the transformer operation data of each distribution node in the photovoltaic distribution network is collected in time periods. The collected operation data includes primary current, secondary current, load current amplitude, and current phase angle. By comparing the difference between the rated current value and the actual current value, a current anomaly threshold range is set. For example, under the condition of rated current of 800A, when the actual current deviates from the rated current by more than ±10%, it is defined as an abnormal state. The current state of this time period is marked as an abnormal current state by the threshold judgment method. The marking result is recorded in binary data form, where 0 represents normal current and 1 represents abnormal current. In actual operation, each sampling period is set to 1 second. By performing sliding window detection on the collected data for 600 consecutive seconds, a load abnormal current state sequence of the node transformer is formed. This sequence serves as the input basis data for subsequent instability quantification.

[0021] Step S22: Quantify the output voltage instability of the distribution network transformer based on the abnormal load current state to obtain node output voltage instability data; In this embodiment of the invention, index positioning is performed on the abnormal load current state sequence of each node output in step S21, and the timestamps corresponding to all sampling points marked as 1 are extracted; based on the timestamps, the three-phase voltage instantaneous value sequence at the corresponding moment is extracted from the historical scheduling state, with a sampling frequency of 10kHz and an extraction window length of 200ms (i.e., 2000 sampling points), centered on the abnormal point moment; zero-crossing detection is performed on each voltage instantaneous value sequence, the peak voltage of the positive half-cycle and the peak voltage of the negative half-cycle are recorded, their absolute difference is calculated and then divided by the standard peak voltage of 311V to obtain the waveform asymmetry parameter; the slope is calculated on the voltage sequence, with a sampling interval of 0.1ms, 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, and the variance of the slope of the 5 sampling points before and after each peak point is calculated to obtain the peak slope variance parameter; at the same time, the voltage sequence is... Harmonic analysis was performed on the series using an FFT with a Hanning window and a window length of 2000 points. The total harmonic distortion (THD) of the 3rd, 5th, 7th, 9th, and 11th harmonics was extracted. The calculation method was to take the square root of the sum of the squares of the amplitudes of each harmonic and then divide it by the fundamental amplitude. The waveform asymmetry, peak slope variance, and THD were linearly weighted and fused with weighting coefficients of 0.4, 0.3, and 0.3, respectively, to output the single-point voltage instability strength value. The arithmetic mean of the instability strength values ​​of all sampling points within the same abnormal period was performed to obtain the representative instability strength of the abnormal event. The representative instability strengths of all abnormal events for each node throughout the year were arranged in chronological order to form a time series. The sampling interval was the interval between abnormal events, with a minimum interval of 1 second and a maximum interval of 86400 seconds. Cubic spline interpolation was performed on this time series with a target interpolation frequency of 1 Hz to generate continuous node output voltage instability data.

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

[0023] Step S23: Perform time-dimensional instability intensity increment derivation on the node output voltage instability data to generate time-dimensional node instability intensity increment data; In this embodiment of the invention, the node output voltage instability data output in step S22 is subjected to sliding differential calculation. The sliding window length is set to 3600s (1 hour), and the step size is 60s. The hourly increment value is obtained by subtracting the starting value from the end value of each window. The increment value sequence is subjected to first-order autoregression processing, and the coefficients are estimated using the Yule-Walker equation with a fixed order of 1. The residual sequence is then calculated. Fractional derivatives are performed on the residual sequence with an order of 0.7, using the Grünwald-Letnikov discrete formula with a step size of 60 sampling points (60 seconds). The coefficient table is pre-calculated and cached. The moving standard deviation is calculated on the derivative result with a window length of 720 sampling points. The local fluctuation intensity is obtained by using a 12-hour interval and a 60-second step size. The four sequences of original hourly increment, residual, fractional derivative, and local fluctuation intensity are concatenated column-wise to form a four-dimensional time series matrix with 8760 rows corresponding to the number of hours in the year and 4 columns. Principal component analysis (PCA) is performed on this matrix, retaining the principal components with a cumulative variance contribution rate greater than 95%, typically the first two principal components. The principal component scores are used as the final instability intensity increment feature vector, and each node outputs an 8760×2 numerical matrix. The L2 norm is calculated for each row of the matrix to obtain the scalarized comprehensive instability intensity increment value, forming the node instability intensity increment data in the time dimension.

[0024] In another embodiment, after obtaining the node output voltage instability data, the instability intensity increment is derived in the time dimension. The process first divides the 600s period output voltage instability data sequence into 60 sub-intervals, each 10s in length. Within each sub-interval, the increment of the voltage instability quantization value is calculated. The instability intensity increment between adjacent intervals is obtained using a differential calculation method. For example, if a node's instability quantization value is 7.5kV·s in the first interval and 10.2kV·s in the second interval, then the instability intensity increment for that interval is 2.7kV·s. This process is repeated for the entire 600s sequence to form a node instability intensity increment sequence in the time dimension. This sequence reflects the instability change process of the node in different time intervals, thus providing time-series instability intensity data for subsequent cluster analysis.

[0025] Step S24: Perform cluster analysis on the incremental data of node instability intensity to obtain node instability intensity cluster data.

[0026] In this embodiment of the invention, the instability intensity increment data of all 96 nodes output in step S23 are vertically concatenated to form a 96×8760 feature matrix, where each row of the matrix represents the comprehensive instability intensity increment value of a node per hour throughout the year; Z-score standardization is performed on this matrix, and the mean and standard deviation of each column are calculated independently, resulting in a mean of 0 and a standard deviation of 1 for each column after transformation; the K-means clustering algorithm is initialized, with the number of cluster centers K set to 5, the initial center points selected using the K-means++ algorithm, the maximum number of iterations set to 500, and the convergence threshold set to 0.0001; clustering calculation is performed, using Euclidean distance as the distance metric, and the cluster center is updated to the mean of the samples within the cluster in each iteration; after clustering is completed, the result is output. Each node has a cluster label, ranging from 0 to 4. Five cluster center coordinates are also output, each a 8760-dimensional vector. For each cluster, the number of nodes, the average distance from the cluster's samples to the center, and the minimum distance between clusters are calculated to form a clustering statistics table. The original feature matrix and cluster labels are horizontally concatenated, adding a new column `cluster_label`, and the dataset is re-saved as a clustering augmentation dataset. The cluster center matrix is ​​output as a separate file, `centers_5clusters.csv`, with dimensions 5×8760 and data precision retained to six decimal places. The final node instability intensity clustering data consists of node number, cluster label, region number, annual hourly incremental sequence, and cluster center number.

[0027] In another embodiment, after acquiring the node instability intensity increment data in the time dimension, cluster analysis is performed on the data. The process begins by normalizing the instability intensity increment sequence of each node, mapping the instability intensity increment values ​​of different nodes to a unified numerical range of 0 to 1 to facilitate cluster calculation. Then, the Euclidean distance metric is used to calculate the similarity of the instability increment sequences between different nodes, constructing an instability distance matrix between nodes. Based on this matrix, a hierarchical clustering-based aggregation method is used to process the data, gradually merging nodes with high similarity into a single cluster category. For example, in a certain calculation, the instability intensity increment sequences of nodes A and B within a 600s period have a similarity of 0.92 and are merged into the same node category. Through gradual aggregation, a complete cluster structure is formed, ultimately obtaining the node instability intensity cluster data, which is used to guide the load transfer logic design between adjacent nodes.

[0028] Step S22 includes the following steps: Step S221: Calculate the winding heat increment index of the distribution network transformer based on the abnormal load current state to obtain the winding heat increment index; Step S222: Perform peak thermal accumulation coupling on the winding heat increment index to obtain winding thermal accumulation data; Step S223: Quantify the magnetic circuit saturation trend based on winding thermal accumulation data; Step S224: Analyze the output voltage waveform distortion intensity based on the magnetic circuit saturation trend to obtain the output voltage distortion intensity; Step S225: Quantify the output voltage instability of the distribution network transformer based on the output voltage distortion intensity to obtain node output voltage instability data.

[0029] In this embodiment of the invention, the original three-phase current effective values ​​corresponding to all sampling points marked as 1 are extracted from the load abnormal current state sequence output in step S21. The sampling frequency is 10Hz, and the unit is A. Each phase is processed independently, and the maximum value among the three phases is taken as the calculation benchmark. Copper loss is calculated for each abnormal sampling point. The copper loss is equal to the square of the current effective value multiplied by the winding DC resistance value, which is fixed at 0.015. Based on the transformer model and nominal parameters, the copper loss value is multiplied by the sampling interval of 0.1 seconds to obtain the instantaneous heat energy increment at a single point, in J. A sliding accumulation is performed on this heat energy increment, with a sliding window length of 100 sampling points (10 seconds) and a step size of 1 sampling point, generating a local heat accumulation sequence. An exponentially weighted moving average (EWMA) is then applied to the local heat accumulation sequence, with a smoothing coefficient... Set to 0.2 to suppress high-frequency fluctuations; output the winding heat increment index corresponding to each abnormal sampling point. The index value is equal to the heat accumulation value after EWMA smoothing divided by the standard heat accumulation value of 4500J within 10s under rated load, thus achieving dimensionless measurement; sort all abnormal points throughout the year by timestamp to form the winding heat increment index.

[0030] Local maximum detection is performed on the winding heat increment index sequence output in step S221. The detection method is to compare the current point with the five sampling points before and after it, a total of 11 points. If the current point is the maximum value and is greater than the threshold of 0.8, it is marked as a peak point. For each peak point, the starting point of the continuous increasing segment is traced forward, and the ending point of the continuous decreasing segment is traced backward to form a heat accumulation event interval. Numerical integration is performed on all heat increment indices within each event interval. The integration method adopts the trapezoidal rule, and the step size is fixed at 0.1s to obtain the total heat accumulation in a single event, and the unit is the dimensionless integral value. At the same time, the event duration is recorded in seconds. The minimum duration threshold is set to 3s, and events below this value are discarded. The retained heat accumulation events are sorted by occurrence time to generate an event list. Each event includes the start sampling index, end sampling index, peak sampling index, total heat accumulation, duration, and peak heat index. Time resampling is performed on this list at a sampling frequency of 1Hz. The method is to evenly distribute the total heat accumulation of each event to every second within the event duration, and assign a value of 0 to time points not covered by the event. The resampled winding heat accumulation data is output.

[0031] The winding thermal accumulation data output in step S222 is subjected to sliding maximum value filtering with a sliding window length of 600 seconds and a step size of 1 second to extract the local thermal accumulation upper limit trend line. First-order difference is performed on this trend line to obtain the thermal accumulation change rate sequence. The change rate threshold is set to 0.005. When the change rate of 10 consecutive sampling points is greater than 0.005, the system is marked as entering the magnetic circuit pre-saturation stage. Second-order difference is performed on the data within the pre-saturation stage. If the second-order difference value is greater than 0.0001 for 5 consecutive points, the system is determined to have entered the accelerated saturation interval. The starting point of the accelerated saturation interval is located, and the entry and exit times are recorded. The exit condition is when the thermal accumulation data value falls below 0.3. Furthermore, the rate of change is less than 0; for each saturation interval, four parameters are extracted: duration, peak heat accumulation value, average rate of change, and maximum value of the second difference; simultaneously, the original heat accumulation data is normalized within the saturation interval by subtracting the minimum value of the interval and then dividing by the maximum value of the interval minus the minimum value, to obtain a standardized heat accumulation curve in the range of 0 to 1; slope fitting is performed on this curve, and the least squares method is used to fit a straight line, and the output slope value is used as the magnetic circuit saturation rate index; the eight parameters, namely saturation interval number, start time, end time, duration, peak heat accumulation, average rate of change, maximum value of the second difference, and saturation rate slope, are arranged in chronological order to form the magnetic circuit saturation trend.

[0032] Based on the start and end times of the magnetic circuit saturation trend output in step S223, the instantaneous three-phase voltage value sequence within the corresponding time period is extracted from the historical scheduling state. The sampling frequency is 10kHz, the extraction window is aligned with the saturation event boundary, and the buffer is extended by 5ms before and after. Zero-crossing detection is performed on each phase voltage sequence, and the peak voltage amplitudes of the positive and negative half-cycles are recorded. The asymmetry is calculated as the absolute difference between the positive and negative peak values ​​divided by 311V. Harmonic analysis is performed on the voltage sequence using Blackman window FFT with a window length of 2048 points. The amplitudes of the 3rd, 5th, 7th, 9th, and 11th odd harmonics are extracted, and the total harmonic distortion (THD) is calculated. The slope is calculated by dividing the square root of the sum of the squares of the harmonic amplitudes by the fundamental amplitude. For the voltage waveform, a slope calculation is performed with a sampling interval of 0.1 ms. The instantaneous slope is calculated for 10 sampling points before and after each peak point, and the standard deviation of the slope in that interval is used as the waveform kurtosis index. For all sampling points within each saturation event, the arithmetic mean of the above three indices is calculated to obtain the event-level asymmetry mean, THD mean, and kurtosis mean. The three means are then weighted and fused with weighting coefficients of 0.5, 0.3, and 0.2, respectively, to output the single-event output voltage distortion intensity value. All saturation events throughout the year are sorted by time to form the output voltage distortion intensity.

[0033] For each magnetic circuit saturation event output in step S224, the corresponding output voltage distortion intensity value is mapped to a continuous time series in seconds based on its start and end times on the time axis. The time range covers 8760 hours throughout the year. Each saturation event is assigned the same distortion intensity value at all seconds within its duration, and the time points where saturation does not occur are uniformly assigned a value of 0. A sliding maximum value filter is then applied to this preliminary filling sequence, with a sliding window length set to 300s and a step size of 1. The kernel is first set to s to smooth transient distortion fluctuations and preserve local peak features. Then, a Gaussian weighted moving average is applied with a kernel standard deviation of 60s and a total kernel length of 360s. Mirror filling is used at the convolution boundaries to ensure no distortion of edge data. A hard truncation operation is performed on the filtered sequence, setting a lower threshold of 0.1 and an upper threshold of 0.9. All values ​​less than 0.1 are forcibly replaced with 0, and all values ​​greater than 0.9 are forcibly replaced with 0.9, achieving dynamic range compression and outlier suppression. The processed sequence is then segmented by hour, and the arithmetic mean of all second-level data within each hour is calculated to generate an hourly-granular output voltage instability intensity sequence with a length of 8760, corresponding to one value per hour throughout the year. Trend separation is performed on the hourly sequence by subtracting the original value from the 12-hour moving median to obtain a residual sequence. An absolute value transformation is then performed on the residual sequence to form a fluctuation-enhanced instability intensity sequence. Finally, two parallel sequences are output: the first is the original hourly mean sequence, and the second is the fluctuation-enhanced sequence, which together constitute the node output voltage instability data.

[0034] In another embodiment, after collecting the abnormal load current state of the transformer in the photovoltaic distribution network, the primary and secondary currents are recorded according to a sampling period of 1 second. Based on the historical rated current, the deviation between the actual current and the rated current is calculated. The deviation is combined with the winding resistance parameter to obtain the instantaneous heat generation power data of the winding. For example, if the winding resistance of a transformer is 0.15Ω, when the detected current reaches 850A and the rated current is 800A, the current deviation is 50A, and the calculated heat generation power increment is 6375W. The heat generation power increment is accumulated every second within 600s of continuous monitoring and converted into temperature rise data. The temperature rise change process is fitted in the form of a logarithmic curve to finally obtain the winding heat increment index. This index is recorded in dimensionless numerical form and is used to characterize the heat accumulation rate of the winding under abnormal current conditions.

[0035] After obtaining the winding heat increment index, the peak value of the index is identified 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 accumulated to form a peak heat accumulation curve, which reflects the trend of winding heat accumulation under long-term operating 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 gradually adding them together. This data reflects the total heat accumulation of the winding under long-term abnormal current operation, and finally forms the winding heat accumulation data, which provides input conditions for subsequent quantification of magnetic circuit saturation trend.

[0036] After obtaining the winding thermal accumulation data, this data is coupled with the core permeability parameter for analysis. The core permeability is set to 2.3 T·m / A. By comparing the thermal accumulation value with the core saturation point magnetic flux threshold, when the winding thermal accumulation data increases rapidly in a short period of time, the internal magnetic reluctance of the core increases, causing the magnetic circuit to tend towards saturation. In order to quantify the magnetic circuit saturation trend, the thermal accumulation data is differentially processed to calculate its rate of change per unit time, and then mapped with the actual magnetic flux density sampling data. For example, when the winding thermal accumulation data growth rate reaches 0.15 / 100s, the corresponding core magnetic flux density rises to 2.1T, and the node is determined to be in the critical range of magnetic circuit saturation. The quantification result is recorded in the form of a saturation trend index to provide input data for voltage distortion intensity analysis.

[0037] After obtaining the magnetic circuit saturation trend index, the excitation current waveform data for the corresponding time period is extracted and subjected to Fourier decomposition to extract the fundamental and harmonic components. When the amplitude of the third harmonic of the excitation current reaches 12% of the fundamental and the amplitude of the fifth harmonic reaches 7% of the fundamental, it indicates that magnetic circuit saturation has caused significant waveform distortion. Further peak period detection is performed on the waveform distortion. By comparing the interval and amplitude change of the waveform peak points, the period distortion rate is calculated. For example, if the difference in peak amplitude of repeated peaks exceeds 15% of the fundamental peak within 10 fundamental periods, it is recorded as the distortion intensity increment. Finally, by combining the harmonic amplitude, peak period distortion rate, and waveform asymmetry parameters, the output voltage distortion intensity data is obtained. This data is expressed as a percentage of the distortion degree, for example, the output voltage distortion intensity is 8.5%. After obtaining the output voltage distortion intensity, this intensity is coupled with the actual output voltage sampling sequence. The voltage distortion intensity and voltage deviation within each 10s sampling interval are weighted and calculated to form the voltage instability quantization value. For example, in a 10s interval, the rated voltage is 10kV, the actual average voltage is 9.4kV, and the voltage deviation is 6%. At this time, the output voltage distortion intensity is 8.5%. By weighting the deviation and distortion intensity with weighting factors of 0.6 and 0.4, the voltage instability quantization value of this interval is obtained as 7.1%. The voltage instability quantization values ​​of all intervals are accumulated and serialized within a 600s period to finally form the node output voltage instability data.

[0038] The analysis of output voltage waveform distortion intensity includes: Based on the magnetic circuit saturation trend, nonlinear distortion analysis of the excitation current is performed to obtain excitation current distortion data; The peak value ratio between peak repetition cycles is derived from the excitation current distortion data to obtain the peak value ratio of the peak cycle. Based on the peak-to-peak ratio, the variance of the odd harmonic frequency of the excitation current is calculated from the excitation current distortion data to obtain the variance of the odd harmonic frequency. Based on the peak-to-peak ratio and odd harmonic frequency variance, the excitation current distortion data is analyzed by multi-peak positive and negative half-cycle waveform asymmetric analysis of induced electromotive force to generate positive and negative half-cycle waveform asymmetric data. The slope variance of the peak points of the induced electromotive force waveform is obtained by calculating the slope variance of the peak points of the asymmetric data of the positive and negative half-cycle waveforms. The output voltage waveform distortion intensity is analyzed based on the peak-to-peak ratio, the odd harmonic frequency variance, and the slope variance to obtain the output voltage distortion intensity.

[0039] In this embodiment of the invention, nonlinear distortion analysis of the excitation current is performed based on the magnetic circuit saturation trend to obtain excitation current distortion data. Specifically, the instantaneous values ​​of the three-phase excitation current on the high-voltage side of the transformer are extracted from the real-time acquisition system of the distribution network based on the time boundary of the magnetic circuit saturation event. The sampling frequency is fixed at 20kHz, and the extraction time period covers the start and end times of each saturation event and extends before and after by 10ms to ensure waveform integrity. A phase-locked loop (PLL) technique is used to accurately track the 50Hz fundamental frequency for each phase current sequence, generating an ideal sinusoidal reference waveform with the same frequency and phase. The original current is subtracted from this reference waveform to obtain a residual sequence containing only the nonlinear distortion component. The absolute value of this residual sequence is taken, and a moving average filter is applied with a window length of 200 sampling points (10ms) and a step size of 1 sampling point to smooth high-frequency glitches. Local glitches are detected in the filtered sequence. The maximum value is determined by the following criteria: the amplitude of the current sampling point is greater than that of the 10 sampling points before and after it and exceeds the 0.5A threshold. Those that meet the criteria are marked as distortion spikes. The time, amplitude, rising slope (linearly fitted from the first 5 sampling points), falling slope (linearly fitted from the last 5 sampling points), and half-width at half-maximum (the time span in which the amplitude decays to half the peak value) of each spike are recorded. All detected spikes are sorted by occurrence time to form excitation current distortion data, which includes six parameters: phase identifier, timestamp, amplitude, rising slope, falling slope, and half-width at half-maximum. At the same time, the original residual sequence is retained as continuous waveform data. The sequence length changes dynamically according to the duration of the event, with a minimum of 2000 points corresponding to 100ms and a maximum of 60000 points corresponding to 3s, for subsequent joint analysis of periodic structure and harmonic characteristics.

[0040] The peak multiple ratio between peak repetition cycles is derived from the excitation current distortion data to obtain the peak multiple ratio of the peak cycle. Specifically, the peak events in the excitation current distortion data table output in the previous step are processed independently by phase. Taking phase A as an example, the timestamps and amplitudes of all phase A peak events are extracted and sorted in ascending order of time. The time interval between adjacent peak events is calculated. If the interval is within the range of 18ms to 22ms, which is close to the power frequency half-cycle of 20ms, it is determined to be an event within the same cycle group. Peak events that continuously meet this interval condition are clustered into a cycle cluster, with no fewer than 3 events within the cluster. For each cycle cluster, all peak amplitudes contained therein are extracted, and the ratio of the largest amplitude to the smallest amplitude is calculated, called the peak multiple ratio within the cycle. Simultaneously, the amplitude ratios of all adjacent peaks within the cluster are calculated. The arithmetic mean of the two values ​​is called the adjacent peak multiple ratio. The larger of the two values ​​is taken as the final peak multiple ratio of the periodic cluster. The same operation is performed on all periodic clusters throughout the year to generate a peak multiple ratio sequence of the peak period. The total number of sequence elements is equal to the number of periodic clusters. In the example, 147 effective periodic clusters were detected in phase A. The same process is repeated for phases B and C to generate their respective peak multiple ratio sequences. The results of the three phases are merged and sorted uniformly according to the time of event occurrence to form a time series of peak multiple ratios of the peak period. 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 its respective periodic cluster. This sequence serves as the core intermediate parameter for voltage waveform distortion intensity analysis and is used to drive the subsequent calculation of odd harmonic frequency variance.

[0041] Based on the peak-to-peak ratio of the excitation current distortion data, the frequency variance of the odd harmonics of the excitation current is calculated to obtain the frequency variance of the odd harmonics. Specifically, for each period cluster, a Fast Fourier Transform (FFT) with a Blackman window is performed on the original excitation current distortion residual sequence segment. The window length is fixed at 4096 points, the sampling frequency is 20kHz, and the frequency domain resolution is 4.8828Hz. The amplitudes of the 3rd, 5th, 7th, 9th, and 11th odd harmonic components are extracted. For each odd harmonic component, its amplitude sequence in all periods within the cluster is calculated. For example, the 3rd harmonic has 10 amplitude samples in 10 periods. The sample variance is calculated for this amplitude sequence with n-1 degrees of freedom to obtain the periodic fluctuation variance of that harmonic. The calculated variances for the 3rd, 5th, 7th, 9th, and 11th harmonics are normalized. The method involves dividing by the square of the historical average amplitude of the corresponding harmonic during the unsaturated period to eliminate the influence of the inherent amplitude scale; arranging the normalized fifth-order variance values ​​in frequency order to form a 5-dimensional variance vector; performing principal component projection on this vector, with the projection matrix pre-calculated from historical training data, retaining the first principal component as the comprehensive odd-order harmonic frequency variance index; simultaneously recording the original fifth-order variance values ​​and projection coefficients for traceability analysis; multiplying the comprehensive variance index with the peak-to-peak ratio of the corresponding peak period in step S224 to strengthen the harmonic instability weight of high-ratio events; and outputting the final odd-order harmonic frequency variance values, with each period cluster corresponding to a scalar value. A total of 441 values ​​(three-phase total) were generated in the annual example, ranging from 0 to 3.5, dimensionless, with timestamps aligned with the center time of the period clusters, forming a complete odd-order harmonic frequency variance.

[0042] Based on the peak-to-peak ratio and odd harmonic frequency variance, asymmetric analysis of the induced electromotive force of the excitation current distortion data is performed on the multi-peak positive and negative half-cycle waveforms to generate asymmetric positive and negative half-cycle waveform data. Specifically, for each cycle cluster corresponding to the original excitation current distortion residual sequence segment, the instantaneous value sequence of the transformer low-voltage side output voltage is extracted synchronously at a sampling frequency of 20kHz, with the time window strictly aligned with the current segment. Zero-crossing detection is performed on this voltage sequence, and the positive and negative half-cycle intervals are divided based on a 50Hz fundamental frequency. Each half-cycle is theoretically 10ms long, i.e., 200 sampling points. Local maxima are searched within each half-cycle interval. If multiple maxima exist and the interval between adjacent maxima is greater than 5 sampling points (0.25ms), it is marked as a multi-peak structure. For each multi-peak half-cycle, the amplitude, position index, maximum slope of the rising edge, and falling edge of all peak points are recorded. Along the maximum slope; calculate the absolute difference between the maximum peak value in the positive half-cycle and the maximum peak value in the negative half-cycle, and then divide it by the theoretical peak value of 311V to obtain the basic asymmetry; at the same time, calculate the absolute value of the difference between the number of peaks in the positive and negative half-cycles. If the difference is greater than 1, an additional asymmetry penalty term is added, which is equal to the difference multiplied by 0.05; add the basic asymmetry and the penalty term, and then multiply by the geometric mean of the peak-to-peak ratio of the corresponding period cluster and the variance of the odd harmonic frequency to complete the nonlinear weighted modulation; perform the same operation on all period clusters throughout the year. In the A-phase example, 89 multi-peak asymmetric half-cycles were detected, 94 in the B-phase example, and 87 in the C-phase example; output seven parameters for each asymmetric event: time center point, phase, number of peaks in the positive half-cycle, number of peaks in the negative half-cycle, basic asymmetry, penalty term, and asymmetric intensity after modulation, and sort them by time to form asymmetric data of the positive and negative half-cycle waveforms.

[0043] The slope variance of the peak points is calculated for the asymmetric data of the positive and negative half-cycle waveforms to obtain the slope variance of the peak points of the induced electromotive force waveform. Specifically, for each asymmetric event output in the previous step, the corresponding instantaneous voltage waveform segment is traced back. For each identified peak point, 10 sampling points before and after it are extracted, i.e., a data window of 0.5ms before and after each point. First-order difference is performed on the 21 sampling points within this window to obtain 20 instantaneous slope values ​​(V / s). The sample variance is calculated for these 20 slope values, with 19 degrees of freedom, to obtain the local slope fluctuation intensity of the peak point. If the half-cycle contains multiple peak points, the arithmetic mean of the slope variances of all peak points is calculated as the representative slope variance of the half-cycle. Simultaneously, the slope variance is recorded... The maximum and minimum single-point slope variances are used to characterize the dispersion of waveform steepness distribution. Multiplying the half-cycle representative slope variance by the ratio of the peak period multiple associated with the event strengthens the contribution of high-multiple events to waveform steepness instability. Then, a weighted summation is performed with the odd-order harmonic frequency variances, with weighting coefficients of 0.6 and 0.4, to generate a comprehensive slope distortion index. The same process is applied to all asymmetric events throughout the year, outputting parameters such as the timestamp, phase number, number of peak points, representative slope variance, maximum slope variance, minimum slope variance, and comprehensive slope distortion index for each event. This parameter set is arranged according to the occurrence time, constituting the slope variance of the peak points of the induced electromotive force waveform. In the example, there are a total of 270 events across three phases, with values ​​ranging from... to (V / s)², time stamp accuracy 0.05ms.

[0044] For each asymmetric event, the three parameters associated with it—peak-to-peak ratio, odd-harmonic frequency variance, and comprehensive slope distortion index—are normalized. The normalization method is to subtract the minimum value of the parameter throughout the year and then divide by the maximum value minus the minimum value, so that all three are uniformly mapped to the interval between 0 and 1. The three normalized parameters are then weighted and fused with weight coefficients of 0.4, 0.3, and 0.3, respectively, and the weighted sum is used as the base distortion intensity for the event. This base distortion intensity is then nonlinearly enhanced by taking its 1.2 power to amplify the distinguishability of the high distortion region. The enhanced value is then multiplied by a scaling factor that represents the duration of the event relative to the power frequency cycle. The scaling factor is equal to the actual half-cycle length divided by 1. 0ms, compensating for waveform truncation effect; outputting the final output voltage distortion intensity value, with each asymmetric event corresponding to a scalar, totaling 270 values ​​for the three phases throughout the year, ranging from 0.05 to 0.98, dimensionless; for ordinary periods without multi-peak asymmetry, a uniform 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 generate a continuous time series; the series is then subjected to a 10ms sliding maximum value filter, followed by Gaussian smoothing with a standard deviation of 5ms to obtain the smoothed second-level output voltage distortion intensity; the sampling interval is 1s, serving as the direct input basis for transformer output voltage instability quantization.

[0045] In another embodiment, specifically in the case of nonlinear distortion analysis of excitation current based on magnetic circuit saturation trend, after collecting the excitation current waveform data of the transformer in the photovoltaic distribution network, the sampling frequency is set to 10000Hz to ensure accurate resolution of the fundamental and higher harmonics. First, when the magnetic circuit saturation trend index exceeds the critical value of 0.85, the excitation current data for the corresponding time period is intercepted. The waveform is then split into frequency domains by fast Fourier decomposition to extract the fundamental component and the amplitude of each harmonic. The distortion data of the excitation current is recorded, where the fundamental frequency is 50Hz. The amplitude of the 3rd harmonic is detected to be 11% of the fundamental, the amplitude of the 5th harmonic is 7% of the fundamental, and the amplitude of the 7th harmonic is 3% of the fundamental. The presence 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 series matrix for subsequent peak periodic characteristic analysis.

[0046] In a specific embodiment of deriving the peak multiple ratio between repetitive peaks in excitation current distortion data, the peak current point within each fundamental period is extracted from the distortion data, the peak amplitude of adjacent periods is recorded, and the multiple relationship between the peak amplitudes of each period is calculated to form a peak-to-peak 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. Then, a multiple ratio sequence is formed between each period, such as 1.12, 0.96, 1.10, 0.97, 1.04, 0.94, 1.08, 0.97, and 1.02. By statistically analyzing this ratio sequence, the peak-to-peak multiple ratio of the peak period is found to be 1.02, which indicates the amplification trend when the current peaks recur in multiple periods. This ratio data serves as the input parameter for subsequent odd harmonic frequency variance calculation.

[0047] In a specific embodiment of calculating the variance of odd harmonic frequencies of excitation current based on the peak-to-peak ratio of the excitation current distortion data, the extracted odd harmonic components are statistically analyzed. The amplitudes of the 3rd, 5th, 7th, 9th, and 11th harmonics are selected as the calculation objects. The variance of each odd harmonic amplitude is calculated after normalization to obtain the variance value of the odd harmonic frequency. For example, the normalized amplitudes are 0.11, 0.07, 0.03, 0.015, and 0.008, respectively. The variance is 0.00162. This variance value is correlated with the peak-to-peak ratio to characterize the degree of unevenness in the distribution of harmonic frequencies. The variance result of the odd harmonic frequency is used as input data to pass to the multi-peak waveform asymmetric analysis step.

[0048] In a specific embodiment of the multi-peak positive and negative half-cycle waveform asymmetric analysis of excitation current distortion data based on the peak-to-peak ratio of the peak period and the variance of the odd harmonic frequency, the peak amplitude and duration of the positive and negative half-cycles are extracted from the distorted current waveform, and their waveform differences are compared. When the average peak amplitude of the positive half-cycle is 480A and the average peak amplitude of the negative half-cycle is 450A, and the duration of the positive half-cycle is 9.9ms and the duration of the negative half-cycle is 10.3ms, positive and negative half-cycle asymmetric data are formed, and the asymmetric amplitude ratio is recorded as 1.07 and the asymmetric time ratio as 0.96. These asymmetric data are formed into a matrix and used as input parameters in the slope variance calculation step.

[0049] In a specific embodiment of calculating the slope variance of the peak point for asymmetric data of positive and negative half-cycle waveforms, the rate of change of current within a 2ms range before and after the peak point is extracted in each half-cycle waveform to form a peak point slope sequence. For example, the slope values ​​detected in the positive half-cycle are 125A / ms, 118A / ms, and 130A / ms, and the slope values ​​detected in the negative half-cycle are 110A / ms, 105A / ms, and 115A / ms. By comparing and calculating the slope variance of the positive and negative half-cycles, the variance of the positive half-cycle is 28, and the variance of the negative half-cycle is 27. The combined slope variance of the peak point of the induced electromotive force waveform is obtained by weighting the two.

[0050] In a specific embodiment of analyzing the output voltage waveform distortion intensity based on the peak-to-peak ratio, odd harmonic frequency variance, and slope variance, the three types of data are used as multi-dimensional inputs and weighted by coefficients of 0.4, 0.35, and 0.25 respectively to form a comprehensive distortion intensity index. For example, if the peak-to-peak ratio is 1.02, the odd harmonic frequency variance is 0.00162, and the comprehensive slope variance is 27.5, the output voltage waveform distortion intensity is calculated to be 8.7% through weighted calculation. This intensity value is recorded as a quantitative index in the node output voltage distortion data sequence for subsequent instability quantification analysis.

[0051] Step S23 includes the following steps: Step S231: Based on the voltage instability data output by the node, plot the voltage instability curve to construct the voltage instability curve; Step S232: Perform time-series divergent incremental gradient derivation on the voltage instability curve to obtain the time-series divergent incremental gradient of the voltage instability state; Step S233: Perform fractional numerical differentiation on the time-series divergent incremental gradient to obtain the divergent numerical fractional order; Step S234: Use an autoregressive model to predict and summarize the divergence trend of the fractional order of the divergent values, and obtain the divergence trend prediction and summary data; Step S235: Based on the divergence trend prediction and induction data, derive the time dimension instability intensity increment to generate the time dimension node instability intensity increment data.

[0052] As an example of the present invention, reference is made to... Figure 3 As shown, step S23 in this example includes: Step S231: Based on the voltage instability data output by the node, plot the voltage instability curve to construct the voltage instability curve; In this embodiment of the invention, a voltage instability curve is plotted based on the node output voltage instability data to construct a voltage instability curve. Specifically, the input is node output voltage instability data recorded throughout the year, with a sampling interval of 1 second and a value range of 0 to 0.9. The sequence is then resampled along the time axis, with the target frequency adjusted to a 1-minute granularity. This is achieved by calculating the arithmetic mean of every 60 consecutive second-level data points, generating a new sequence of length 52560, corresponding to one instability intensity value per minute throughout the year. The timestamps are then uniformly formatted as YYYY-MM-DD. HH:MM:00, ensuring integer alignment; perform cubic spline interpolation on the resampled sequence, with an interpolation target of one point every 10 seconds, generating an intermediate transition sequence for smooth visual presentation, and setting the interpolation boundary condition to a natural boundary, i.e., the second derivative is zero; superimpose the original minute mean points as anchor constraints on the interpolated sequence to avoid overfitting; perform moving median filtering on the interpolated sequence, with a window length of 5 minutes (30 10-second points) and a step size of 1 point to suppress local impulse noise; segment the final sequence by day. Each day, 1440 10-second points are collected, totaling 365 segments. Each segment is independently normalized to the interval between 0 and 1 by subtracting the minimum value of the day and then dividing by the maximum value of the day minus the minimum value, preserving the relative trend of change within the day. 365 sets of normalized 10-second granular sequences are output, each set containing 1440 values, forming a set of voltage instability curves. Each curve represents the voltage instability evolution trajectory of the entire day. The horizontal axis is the time scale from 00:00:00 to 23:59:50, with a step size of 10 seconds, and the vertical axis is the normalized instability intensity.

[0053] In another embodiment, the amplitude of output voltage instability data at different time points is collected. The data sampling frequency is set to 1000Hz, the collection duration is 600s, and the data is stored as a time-series vector set. A voltage instability curve is plotted using a Cartesian coordinate system, with the horizontal axis representing time t and the vertical axis representing the voltage instability amplitude. When plotting, the voltage instability amplitude at each moment is connected sequentially to form a continuous curve, thus constructing the voltage instability curve. For example, in the time period from 0s to 600s, the node voltage instability amplitude gradually increases from 0.5V to 12.6V, forming a nonlinear rising curve. The voltage instability curve is used to provide a data basis for the subsequent derivation of the time-series divergent incremental gradient.

[0054] Step S232: Perform time-series divergent incremental gradient derivation on the voltage instability curve to obtain the time-series divergent incremental gradient of the voltage instability state; In this embodiment of the invention, for each daily voltage instability curve output in step S231, a 10-second granular sequence of 1440 points is taken, and a first-order forward difference is performed. The calculation method is to subtract the value of the nth point from the value of the (n+1)th point, generating a difference sequence of length 1439, with the unit being dimensionless per 10s. Sign separation is performed on this difference sequence; the positive values ​​are retained and marked as positive divergence, while the negative values ​​are taken as absolute values ​​and marked as negative convergence. The sliding standard deviation is calculated for the positive divergence part, with a window length of 12 points (2 minutes) and a step size of 1 point, to obtain the local fluctuation intensity. The same operation is performed on the negative convergence part, forming a positive divergence fluctuation sequence and a negative convergence fluctuation sequence respectively. The original difference sequence, positive fluctuation sequence, and negative fluctuation sequence are then combined. The three elements are concatenated according to their element positions to form a three-dimensional gradient feature vector sequence with a length of 1439. Principal component analysis is performed on this three-dimensional sequence to extract the first principal component as the comprehensive time-series divergent incremental gradient sequence, retaining information with an original variance contribution rate greater than 92%. At the same time, the original difference value, positive fluctuation value, negative fluctuation value, and principal component score corresponding to each time point are recorded. The same processing is repeated for each curve for 365 days of the year, outputting 365 sets of time-series divergent incremental gradients, each set containing 1439 scalar values, with time labels aligned to the start time of every 10 seconds, such as 00:00:10, 00:00:20, etc. This reflects the acceleration or deceleration trend of voltage instability in a short time scale, and the values ​​can be positive or negative, ranging from -0.15 to +0.18.

[0055] In another embodiment, the voltage instability curve data is discretized by using a sampling interval of 1ms as the time step, transforming the continuous curve into 600,000 discrete data points, and then calculating the difference between adjacent data points. To obtain the instantaneous rate of change of the voltage instability curve, a sliding window method was further used with a window width of 50ms and a step size of 10ms to smooth the differential data and eliminate high-frequency noise, thereby obtaining the time-series divergent incremental 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 shows an accelerated divergence trend in the later stage.

[0056] Step S233: Perform fractional numerical differentiation on the time-series divergent incremental gradient to obtain the divergent numerical fractional order; In this embodiment of the invention, for each set of time-series divergent incremental gradients output in step S232, with a length of 1439 and a sampling interval of 10s, fractional derivatives are performed using the Grünwald-Letnikov discretization formula, with the order α fixed at 0.7; the number of coefficients is equal to the sequence length, and the recursive formula is as follows: , ;in The gamma function is used, and k increases from 1 to 1438. Convolution is performed on the sequence, with zero-padding used for boundary handling; the padding length is equal to the sequence length minus 1. The output convolution sequence is still 1439 points long, with each point representing the 0.7th derivative at that moment. A moving maximum filter is applied to this derivative sequence, with a window length of 6 points (60 seconds) and a step size of 1 point, preserving local extrema. An exponentially weighted moving average is then applied, with a smoothing coefficient of 0.3 to suppress high-frequency oscillations. Z-score standardization is applied to the filtered sequence, subtracting the sequence mean and dividing by the standard deviation to achieve an output distribution with a mean of 0 and a standard deviation of 1. The standardized sequence is multiplied by a fixed scaling factor of 0.85 to control the dynamic range. The final divergent numerical fractional order is output, with 1439 points per group, ranging from -2.1 to +1.9, dimensionless, and a time resolution of 10 seconds. All 365 groups of sequences for the year are merged to form the complete annual divergent numerical fractional order.

[0057] In another embodiment, the Grünwald–Letnikov fractional difference method is used, with the order set to 0.75 and the numerical deviation step length set to 1 ms. Fractional operations are performed on the time-series divergent incremental gradient sequence. During the operation, the gradient value at each time step is weighted and superimposed with its historical time points. The weighting coefficients are controlled by binomial expansion coefficients, thereby achieving the approximate calculation of the fractional derivative. For example, at 400 s, the result of the traditional integer first derivative is 0.045 V / ms, while the result after the 0.75th order fractional numerical differentiation is 0.038 V / ms. This processing preserves the nonlocality of the divergence process, and the generated divergent numerical fractional sequence is used for subsequent trend prediction and induction.

[0058] Step S234: Use an autoregressive model to predict and summarize the divergence trend of the fractional order of the divergent values, and obtain the divergence trend prediction and summary data; In this embodiment of the invention, the input is the fractional order of the divergent numerical values ​​for the whole year output in step S233, divided into 365 groups by day, with 1439 points in each group and a sampling interval of 10 seconds. A stationarity test is performed on each group of daily sequences using the augmented Dickey-Fowler test (ADF), with a fixed lag order of 12. If the p-value is less than 0.05, the sequence is considered stationary; otherwise, first-order differencing is performed until the stationarity condition is met. An autoregressive model (AR) is fitted to the stationary sequence. The order p is determined by the truncation position of the partial autocorrelation function (PACF), and the maximum search order is set to 24, corresponding to a 4-minute historical dependency. Finally, the order that minimizes the Akaike Information Criterion (AIC) is selected as the optimal p. The Yule-Walker equation is used to solve the autoregressive coefficients, without introducing a moving average term to ensure a pure autoregressive structure. An AR(p) model is trained independently for each group of sequences. For example, after testing, a daily sequence has p=8, and the coefficients... to The calculation is performed precisely using matrix inversion. Rolling predictions are then made for the last 60 points (the last 10 minutes) of the sequence using the trained coefficients, with a prediction step size of 1. Each prediction uses the previous p true values ​​to predict the next time step, generating a total of 60 predicted values. The 60 predicted values ​​are then subtracted point-by-point from their corresponding true values ​​to obtain the prediction residual sequence. Three statistical measures are calculated on the residual sequence: Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Maximum Absolute Error (MAXAE). Simultaneously, 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 daily AR order, the 60-step prediction sequence, the three error indicators, and the trend representative value—a total of 64 parameters—are packaged to form a daily divergence trend prediction summary record. This process is repeated for 365 days throughout the year, outputting 365 records, constituting the divergence trend prediction summary data.

[0059] In another embodiment, an autoregressive (AR) method with an order of 5 is used to divide the divergent numerical fractional order sequence into a training segment and a validation segment. The training segment is from 0s to 480s, and the validation segment is from 480s to 600s. During the training process, the least squares method is used to estimate the autoregressive coefficients, and the numerical fractional order value at each time step is predicted and compared with the actual value. The mean square error is used as the convergence criterion. When the mean square error is less than 0.01, the fitting is considered to have converged. The final prediction and induction results show that the numerical fractional order gradually increases from 0.041 to 0.073 in the interval from 540s to 600s. The trend prediction and induction data accurately reflect the divergent acceleration process of voltage instability.

[0060] Step S235: Based on the divergence trend prediction and induction data, derive the time dimension instability intensity increment to generate the time dimension node instability intensity increment data.

[0061] In this embodiment of the invention, the last predicted value of the AR model in each of the 365 divergence trend prediction summarization records output in step S234 is extracted to form a trend endpoint sequence of length 365. Each value represents the predicted evolution direction of the voltage instability divergence state at the end of the day. A sliding difference is performed on this sequence with a window span of 7 days. The weekly increment value is obtained by subtracting the predicted value of day n from the predicted value of day n+7, generating a weekly increment sequence of length 358. An exponentially weighted moving average (EWMA) is applied to the weekly increment sequence with a smoothing coefficient of 0.4 to suppress short-term fluctuations. A second difference is then performed on the smoothed sequence to obtain the acceleration change sequence. The algorithm is used to identify trend turning points. The original weekly increment, EWMA smoothed increment, and acceleration change are concatenated element-wise to form a three-dimensional feature vector sequence. Principal component analysis (PCA) is performed on this three-dimensional sequence, retaining the first principal component, with a cumulative variance contribution rate greater than 90%. After projection, a scalarized comprehensive increment sequence is obtained. This comprehensive increment sequence is normalized by subtracting the minimum value and then dividing by the range, mapping it to the 0-1 interval. The normalized sequence is then extended back to the original 365-day scale by day, using linear interpolation to fill missing days, with boundary replication filling for the first 6 days and the last 1 day. The final time-dimensional node instability intensity increment data is output.

[0062] In another embodiment, the predicted inductive data is processed by time integration to obtain a cumulative divergence intensity sequence. The time integration adopts the trapezoidal integration method with an integration step size of 1ms. Then, the difference in cumulative divergence intensity in adjacent time periods is calculated and defined as the instability intensity increment data in the time dimension. For example, the cumulative divergence intensity increases from 22.7 to 27.9 in the time period from 550s to 560s, with an increment of 5.2. In the time period from 580s to 590s, the cumulative divergence intensity increases from 35.4 to 46.1, with an increment of 10.7. This shows that the node instability intensity shows an accelerating increment trend as time goes by. The finally generated time dimension instability intensity increment data provides input parameters for subsequent node instability clustering and scheduling logic design.

[0063] Step S3 includes the following steps: Step S31: Perform convolution processing on the node instability intensity clustering data to obtain the node instability convolution intensity; Step S32: Design load transfer logic between adjacent nodes based on node instability convolution strength, thereby reducing the load pressure on the current node and obtaining load transfer logic to perform optimized scheduling of photovoltaic power distribution network.

[0064] In this embodiment of the invention, convolution processing is performed on the node instability intensity clustering data to obtain the node instability convolution intensity. Specifically, the input is the annual hourly node instability intensity clustering data of 96 distribution nodes output in step S24. Each node sequence is 8760 bytes long, with a sampling interval of 1 hour, and the values ​​have been normalized to a mean of 0 and a standard deviation of 1. The 96 sequences are arranged according to geographical topological adjacency relationships to construct a two-dimensional grid structure with 12 rows and 8 columns, corresponding to 8 regions, each with 12 nodes. Physically adjacent nodes are adjacent in the grid, and boundary nodes are padded with zeros to form a 14×10 extended grid. One-dimensional convolution is performed independently on the 8760-hour sequence of each node. The convolution kernel uses a Gaussian function discretization form, with a kernel length of 15 hours, a standard deviation σ = 3 hours, and a maximum weight at the center that decays exponentially towards both sides. The kernel coefficients are pre-calculated and normalized to a sum of 1. The convolution operation uses an efficient mode, step... The output sequence is 1 hour long with no boundary padding, resulting in a length of 8746. A sliding maximum filter is applied to the convolutional sequence with a window length of 5 hours to preserve local extrema. Z-score standardization is then performed to ensure that each node's output sequence again satisfies a mean of 0 and a standard deviation of 1. The standardized sequence is then weighted and fused with the original cluster labels, with the weights being the inverse normalized value of the distance to the cluster centers; closer clusters have higher weights to enhance the trend consistency of similar nodes. 96 sets of node instability convolution intensity sequences are output, each with 8746 points, with timestamps aligned to hours 8 to 8753, covering the entire year's effective analysis period. Simultaneously, four statistics are extracted for each node's sequence: the annual maximum value, the 95th percentile value, the daily average value, and the weekly standard deviation of fluctuation, serving as static convolution intensity features. This feature set, together with the dynamic sequence, constitutes a complete expression of node instability convolution intensity, used to drive the load transfer logic design between adjacent nodes.

[0065] The load transfer logic between adjacent nodes is designed based on the node instability convolution intensity to alleviate the load pressure on the current node, thereby obtaining the load transfer logic to perform optimized scheduling of the photovoltaic distribution network. Specifically, for each node output in step S31, the maximum annual value of its static convolution intensity feature is extracted as an urgency index. If this value is greater than 1.8, it is marked as a high-risk node. The current load rate of the high-risk node is calculated by taking the average load rate of the most recent 24 hours. If it exceeds 0.8, the load transfer process is initiated. A list of all physically directly connected adjacent nodes of the node is obtained, determined according to the grid topology adjacency matrix, with a maximum of 6 adjacent nodes. The remaining carrying capacity of each adjacent node is calculated, equal to 0.8 multiplied by its transformer rated capacity of 1000kVA minus its current average load rate multiplied by 1000kVA. Simultaneously, the available transmission capacity of the connecting line is calculated, taking the minimum value between the line thermal stability limit and the protection setting, in kVA. The remaining carrying capacity is multiplied by the available transmission capacity to obtain the upper limit of the transferable capacity. Adjacent nodes are arranged in ascending order of electrical distance. The distance between the line and the electrical distance is equal to the line resistance multiplied by 1.2 plus the reactance multiplied by 0.8. A progressive allocation is performed on the sorted list, starting from the first-priority node. If its maximum transferable capacity is greater than the amount to be transferred, the entire amount is transferred; otherwise, its maximum 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 and then multiplied by 1000kVA. The allocation process is repeated until the amount to be transferred is zero or there are no available adjacent nodes. The target node number, transfer amount, line number, and priority number are recorded for each transfer. A load transfer instruction sequence is generated, with each instruction containing the source node ID, target node ID, load amount to be switched (kVA), execution time, and duration (minutes). Instructions are sorted by priority and electrical distance, with the shortest path and highest priority executed first. Instructions are encapsulated using the IEC61850 GOOSE message protocol and sent to the corresponding intelligent terminal controller, triggering the circuit breaker's opening and closing actions. After completing one round of load transfer, the real-time load rate and instability convolution intensity of each node are updated, and the next evaluation cycle begins. The scheduling cycle is fixed at 15 minutes to ensure dynamic system balance.

[0066] Step S32 includes the following steps: Step S321: Calculate the average load rate of the node transformer based on the node instability convolution intensity; obtain the configuration status and topology of the surrounding neighboring nodes of the current node; Step S322: Identify the load to be transferred based on the average load rate of the transformer to obtain the load to be transferred at the current node; Step S323: Calculate the product of the carrying capacity difference between adjacent nodes and the line transmission capacity based on the topology and configuration status; Step S324: Based on the distance between the current node and its neighboring nodes, a progressive load transfer strategy is obtained by multiplying the load to be transferred by the product of the carrying capacity difference between the neighboring nodes and the line transmission capacity. Step S325: Design the load transfer logic between adjacent nodes based on the progressive load transfer strategy, thereby reducing the load pressure on the current node and obtaining the load transfer logic.

[0067] In this embodiment of the invention, for each distribution node, the last 24 data points out of the 8746 data points corresponding to the most recent 24 hours are extracted from the node instability convolution intensity dynamic sequence output in step S31, corresponding to one value per hour in the past 24 hours; the arithmetic mean of these 24 points is performed to obtain the time-averaged convolution intensity value; simultaneously, the load rate data recorded every 5 minutes for this node in the same time period, totaling 288 points, are extracted from the historical scheduling database, and the arithmetic mean of these data points is performed to obtain the node transformer average load rate, with a value range of 0 to 1; if the average load rate is greater than 0.8 and the time-averaged convolution intensity value is greater than 1.5, then the node is determined to be an object to be processed; the power grid physical topology diagram is queried for this node to obtain the numbers of all its directly electrically connected adjacent nodes, with a maximum number not exceeding 6; the configuration parameters of each adjacent node are extracted, including the transformer rated capacity being uniformly 1000kVA and the line positive sequence resistance unit. Multiply by the length in km to obtain the total resistance value; multiply the line positive sequence reactance (Ω / km) by the length in km to obtain the total reactance value; multiply the line thermal stability limit current (A) by the nominal voltage (10kV) and divide by the square root of 3 to obtain the thermal stability limit capacity (kVA); simultaneously extract the current average load rate of adjacent nodes, calculated using the same method as the first half of step S321; construct an adjacency structure table, where each row corresponds to an adjacent node and includes five parameters: adjacent node number, resistance value, reactance value, thermal stability limit capacity, and current average load rate; this table is sorted in ascending order by electrical distance, where electrical distance equals resistance value multiplied by 1.2 plus reactance value multiplied by 0.8; output the current node's average load rate and the adjacency structure table as the basic input for subsequent transfer quantity calculations and strategy generation.

[0068] For nodes identified as pending targets in step S321, their average load rate is denoted as L_avg. If L_avg is greater than 0.8, the transfer calculation process is initiated. The safe load threshold is set to 0.75. The load to be transferred, ΔP, is equal to L_avg minus 0.75 multiplied by the transformer's rated capacity of 1000 kVA (unit: kVA). For example, if L_avg = 0.86 for a node, then ΔP = (0.86 - 0.75) × 1000 = 110 kVA. A hard cutoff is performed on the load to be transferred. If ΔP is less than 10 kVA, it is forcibly set to zero and no transfer is performed. If ΔP is greater than 300 kVA, it is forcibly limited to 300 kVA to prevent over-transfer. At the same time, the photovoltaic inverters connected to the current node are checked. The maximum load shelving capacity is calculated as follows: if the total installed capacity of the inverter minus the current minimum technical output is less than ΔP, then ΔP is corrected to this difference. For example, if the total capacity of the photovoltaic inverter is 800kVA and the minimum technical output is 550kVA, then the maximum shelving capacity is 250kVA. If the original ΔP = 300kVA, then it is corrected to 250kVA. The final load to be transferred, ΔP, is output, with a value ranging from 10kVA to 300kVA, with integer precision, in kVA units. This value serves as the total input for the progressive allocation algorithm, and the original calculated value and the reason for correction are recorded for auditing and traceability. The load to be transferred is recalculated at the beginning of each scheduling cycle, which is fixed at 900 seconds (15 minutes) to ensure real-time response to load changes.

[0069] The product of the carrying capacity difference between adjacent nodes and the line transmission capacity is calculated based on the topology and configuration status. Specifically, for each adjacent node in the adjacency table output in step S321, its remaining carrying capacity R_cap is calculated, which is equal to 0.8 multiplied by 1000kVA minus its current average load rate multiplied by 1000kVA, in kVA. If R_cap is less than 0, it is forcibly set to zero, indicating no acceptance capacity. Simultaneously, the thermal stability limit capacity T_lim of the corresponding connecting line of the adjacent node is extracted, in kVA. The product of the carrying capacity difference and the line transmission capacity C_prod is calculated, which is equal to R_cap multiplied by T_lim, in kVA². For example, if the current average load rate of an adjacent node is 0.65, then R_cap = 0.8 × 1000 - 0.65 × 1000 = 150kVA. If the corresponding line T_lim = 800kVA, then C_prod = 150 × 800 = 120000kVA². This product value is normalized by dividing it by the largest C_prod value among all adjacent nodes, compressing the output range to 0 to 1. The original C_prod value is retained for physical constraint verification. Two new columns are added to the adjacency structure table: the original product value and the normalized product value. The adjacency table is rearranged in descending order of the normalized product value, with higher priority nodes at the top. If the normalized product values ​​are the same, they are sorted in ascending order of electrical distance. The updated adjacency structure table is output, with each row containing the adjacent node number, original product value, normalized product value, electrical distance, and current average load rate. This table serves as the capacity constraint basis for progressive load transfer, ensuring that the transfer path meets both node margin and line thermal limits.

[0070] Based on the distance between the current node and its neighboring nodes, a progressive load transfer strategy is obtained by multiplying the load to be transferred by the product of the carrying capacity difference between the neighboring nodes and the line transmission capacity. Specifically, the inputs are: the load to be transferred ΔP output in step S322 and the sorted adjacency structure table output in step S323; the remaining load to be transferred R_remain is initialized to equal ΔP; processing is performed row by row starting from the first row of the adjacency table; for the neighboring nodes in the current row, their original product value C_prod is taken as the upper limit of the acceptable maximum transfer amount; if C_prod is greater than R_remain, R_remain is fully allocated to the node, the transfer amount is recorded as equal to R_remain, and the target node is the node number of the current row; otherwise, C_prod is fully allocated to the node, the transfer amount is recorded as equal to C_prod, and R_remain is updated to equal R_remain minus C_prod. Continue processing the next line until R_remain reaches zero or the adjacency list is traversed. If R_remain is still greater than 0 after the adjacency list is traversed, a secondary transfer is initiated, allowing transfer through intermediate nodes. The transfer path length cannot exceed 2 hops, and the transfer capacity is constrained by the minimum thermal limit of the two-level lines. For each successful allocation record, a transfer priority number is appended, incrementing from 1 according to the allocation order. Simultaneously, the actual electrical path length is calculated, which is equal to the weighted sum of the resistance and reactance of all segments between the source node and the target node. The transfer strategy table is output, with each line containing the target node number, transfer amount (kVA), priority number, path length (Ω), and whether it is a direct connection. For example, if a node ΔP = 110kVA, the first priority node C_prod = 90kVA, and the second priority node C_prod = 60kVA, then the first node is allocated 90kVA, and the second node is allocated 20kVA, with priorities of 1 and 2 respectively. This strategy table serves as the sole basis for generating the final control command.

[0071] The load transfer logic between adjacent nodes is designed based on a progressive load transfer strategy to reduce the load pressure on the current node and achieve load transfer logic. Specifically, the transfer strategy table is arranged in ascending order of priority, with priorities 1 to 6 corresponding to execution delays of 0 seconds, 5 seconds, 10 seconds, 15 seconds, 20 seconds, and 25 seconds respectively, to avoid concurrent impact. Each strategy generates a GOOSE control instruction containing the source node number, target node number, transfer amount (kVA), execution delay, and duration (900 seconds). The instruction encapsulation conforms to the IEC61850-8-1 standard, with APPID starting from 0. Starting with x4001, priorities increase sequentially, with the multicast address fixed at 01-0C-CD-01-00-01. Within 50ms of receiving the instruction, the intelligent terminal drives the circuit breaker to operate, switching the specified load. After the operation is completed, a remote signaling change is sent up, and the master station confirms and updates the real-time load rate. If an instruction fails, the next-ranked node in the adjacency list is immediately used for reassignment, with the reassignment amount equal to the original failed amount. After all operations are completed, an execution log is recorded, including the instruction sequence number, timestamp, actual transfer amount, and target node load change value. The log is used to update the convolution strength and clustering state before the next scheduling cycle, achieving closed-loop optimization.

[0072] The present invention also provides an optimized scheduling system for a photovoltaic distribution network, used to execute the optimized scheduling method for a photovoltaic distribution network as described above, the optimized scheduling system for the photovoltaic distribution network comprising: The status association extraction module is used to obtain the historical scheduling status of the photovoltaic distribution network; based on the historical scheduling status, it performs transformer load status association extraction between distribution nodes in different regions to obtain the node transformer load status; The instability intensity incremental derivation module is used to quantify the output voltage instability of the distribution network transformer based on the load status of the node transformer to obtain node output voltage instability data; and to perform time-dimensional instability intensity incremental derivation on the node output voltage instability data to obtain node instability intensity clustering data. The load transfer logic design module is used to design load transfer logic between adjacent nodes based on node instability intensity clustering data, thereby reducing the load pressure on the current node and obtaining load transfer logic to perform optimized scheduling of the photovoltaic power distribution network.

[0073] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. An optimized scheduling method for a photovoltaic distribution network, characterized in that, Includes the following steps: Step S1: Obtain the historical dispatch status of the photovoltaic distribution network; based on the historical dispatch status, perform transformer load status correlation extraction between distribution nodes in different regions to obtain the node transformer load status; Step S2: Quantify the output voltage instability of the distribution network transformers based on the load status of the node transformers to obtain node output voltage instability data; perform time-dimensional instability intensity increment derivation on the node output voltage instability data to obtain node instability intensity clustering data; Step S3: Design load transfer logic between adjacent nodes based on node instability intensity clustering data to reduce the load pressure on the current node and obtain load transfer logic to perform optimized scheduling of photovoltaic power distribution network.

2. The optimized scheduling method for photovoltaic distribution networks according to claim 1, characterized in that, Step S1 includes the following steps: Step S11: Obtain the historical dispatch status of the photovoltaic distribution network; Step S12: Analyze the high-load periods between distribution nodes in different regions based on the historical scheduling status to obtain the high-load period scheduling status between distribution nodes in different regions. Step S13: Analyze the transmission load power fluctuation between different regional power distribution nodes during the high load period scheduling state to obtain node load power fluctuation difference data; Step S14: Based on the historical scheduling status, extract the transformer load status correlation between different regional distribution nodes by analyzing the difference data of node load power fluctuations, and obtain the node transformer load status.

3. The optimized scheduling method for photovoltaic distribution networks according to claim 1, characterized in that, Step S2 includes the following steps: Step S21: Mark the abnormal current status of the node transformer load and generate the abnormal load current status; Step S22: Quantify the output voltage instability of the distribution network transformer based on the abnormal load current state to obtain node output voltage instability data; Step S23: Perform time-dimensional instability intensity increment derivation on the node output voltage instability data to generate time-dimensional node instability intensity increment data; Step S24: Perform cluster analysis on the incremental data of node instability intensity to obtain node instability intensity cluster data.

4. The optimized scheduling method for photovoltaic distribution networks according to claim 3, characterized in that, Step S22 includes the following steps: Step S221: Calculate the winding heat increment index of the distribution network transformer based on the abnormal load current state to obtain the winding heat increment index; Step S222: Perform peak thermal accumulation coupling on the winding heat increment index to obtain winding thermal accumulation data; Step S223: Quantify the magnetic circuit saturation trend based on winding thermal accumulation data; Step S224: Analyze the output voltage waveform distortion intensity based on the magnetic circuit saturation trend to obtain the output voltage distortion intensity; Step S225: Quantify the output voltage instability of the distribution network transformer based on the output voltage distortion intensity to obtain node output voltage instability data.

5. The optimized scheduling method for photovoltaic distribution networks according to claim 4, characterized in that, The analysis of output voltage waveform distortion intensity includes: Based on the magnetic circuit saturation trend, nonlinear distortion analysis of the excitation current is performed to obtain excitation current distortion data; The peak value ratio between peak repetition cycles is derived from the excitation current distortion data to obtain the peak value ratio of the peak cycle. Based on the peak-to-peak ratio, the variance of the odd harmonic frequency of the excitation current is calculated from the excitation current distortion data to obtain the variance of the odd harmonic frequency. Based on the peak-to-peak ratio and odd harmonic frequency variance, the excitation current distortion data is analyzed by multi-peak positive and negative half-cycle waveform asymmetric analysis of induced electromotive force to generate positive and negative half-cycle waveform asymmetric data. The slope variance of the peak points of the induced electromotive force waveform is obtained by calculating the slope variance of the peak points of the asymmetric data of the positive and negative half-cycle waveforms. The output voltage waveform distortion intensity is analyzed based on the peak-to-peak ratio, the odd harmonic frequency variance, and the slope variance to obtain the output voltage distortion intensity.

6. The optimized scheduling method for photovoltaic distribution networks according to claim 3, characterized in that, Step S23 includes the following steps: Step S231: Based on the voltage instability data output by the node, plot the voltage instability curve to construct the voltage instability curve; Step S232: Perform time-series divergent incremental gradient derivation on the voltage instability curve to obtain the time-series divergent incremental gradient of the voltage instability state; Step S233: Perform fractional numerical differentiation on the time-series divergent incremental gradient to obtain the divergent numerical fractional order; Step S234: Use an autoregressive model to predict and summarize the divergence trend of the fractional order of the divergent values, and obtain the divergence trend prediction and summary data; Step S235: Based on the divergence trend prediction and induction data, derive the time dimension instability intensity increment to generate the time dimension node instability intensity increment data.

7. The optimized scheduling method for photovoltaic distribution networks according to claim 1, characterized in that, Step S3 includes the following steps: Step S31: Perform convolution processing on the node instability intensity clustering data to obtain the node instability convolution intensity; Step S32: Design load transfer logic between adjacent nodes based on node instability convolution strength, thereby reducing the load pressure on the current node and obtaining load transfer logic to perform optimized scheduling of photovoltaic power distribution network.

8. The optimized scheduling method for photovoltaic distribution networks according to claim 7, characterized in that, Step S32 includes the following steps: Step S321: Calculate the average load rate of the node transformer based on the node instability convolution intensity; obtain the configuration status and topology of the surrounding neighboring nodes of the current node; Step S322: Identify the load to be transferred based on the average load rate of the transformer to obtain the load to be transferred at the current node; Step S323: Calculate the product of the carrying capacity difference between adjacent nodes and the line transmission capacity based on the topology and configuration status; Step S324: Based on the distance between the current node and its neighboring nodes, a progressive load transfer strategy is obtained by multiplying the load to be transferred by the product of the carrying capacity difference between the neighboring nodes and the line transmission capacity. Step S325: Design the load transfer logic between adjacent nodes based on the progressive load transfer strategy, thereby reducing the load pressure on the current node and obtaining the load transfer logic.

9. An optimized dispatching system for a photovoltaic distribution network, characterized in that, For executing the optimized scheduling method for a photovoltaic distribution network as described in claim 1, the optimized scheduling system for the photovoltaic distribution network includes: The status association extraction module is used to obtain the historical scheduling status of the photovoltaic distribution network; based on the historical scheduling status, it performs transformer load status association extraction between distribution nodes in different regions to obtain the node transformer load status; The instability intensity incremental derivation module is used to quantify the output voltage instability of the distribution network transformer based on the load status of the node transformer to obtain node output voltage instability data; and to perform time-dimensional instability intensity incremental derivation on the node output voltage instability data to obtain node instability intensity clustering data. The load transfer logic design module is used to design load transfer logic between adjacent nodes based on node instability intensity clustering data, thereby reducing the load pressure on the current node and obtaining load transfer logic to perform optimized scheduling of the photovoltaic power distribution network.

Images