Virtual Power Plant Flexible Load Aggregation and Dispatch Control Method and System

Abstract

This invention provides a method and system for flexible load aggregation and dispatch control in a virtual power plant, relating to the field of power management technology. The method includes acquiring flexible load resource information, calculating the matching degree based on a multi-dimensional matching space, dynamically constructing temporary aggregation units, building a desired response benchmark and individual difference envelope, generating main dispatch instructions and individual compensation instructions, and updating the matching degree weights and envelopes based on actual response data. This invention can improve the response accuracy of the virtual power plant to grid dispatch tasks, achieving efficient aggregation and precise dispatch control of flexible load resources.

Description

Technical Field

[0001] This invention relates to power management technology, and more particularly to a method and system for flexible load aggregation and dispatch control of virtual power plants. Background Technology

[0002] With the rapid development of distributed energy and smart grid technologies, virtual power plants have gradually gained widespread attention as a new energy management approach. However, existing dispatch control methods lack refined modeling and analysis of the individual differences in flexible loads, often employing uniform models and dispatch strategies. This fails to fully consider the differences in response characteristics of different loads, leading to significant deviations between actual and expected responses and affecting dispatch effectiveness. Summary of the Invention

[0003] The present invention provides a method and system for flexible load aggregation and dispatch control of virtual power plants, which can solve the problems in the prior art.

[0004] A first aspect of this invention provides a method for flexible load aggregation and dispatch control in a virtual power plant, comprising: acquiring real-time adjustability information and historical response characteristic data of distributed flexible load resources; receiving a power grid dispatch task and extracting a task feature vector; calculating the matching degree between each distributed flexible load resource and the task feature vector in a multi-dimensional matching space based on the real-time adjustability information, and dynamically grouping resources with matching degrees exceeding a preset matching degree threshold into temporary aggregation units; constructing an expected response benchmark and individual difference envelope for the temporary aggregation unit based on the historical response characteristic data; generating a main dispatch instruction for the temporary aggregation unit based on the expected response benchmark, estimating the execution deviation based on the individual difference envelope, and generating an individual compensation instruction based on the deviation contribution degree when the execution deviation exceeds the limit; issuing the main dispatch instruction and the individual compensation instruction and collecting actual response data, and updating the matching degree weights and the individual difference envelope in the multi-dimensional matching space based on the actual response data.

[0005] Based on the real-time adjustable capability information, the step of calculating the matching degree between each distributed flexible load resource and the task feature vector in the multi-dimensional matching space, and dynamically grouping resources with matching degrees exceeding a preset matching degree threshold into temporary aggregation units includes:

[0006] Extract the historical call time and actual response capacity of each distributed flexible load resource from the historical response characteristic data, calculate the response capacity decay rate, and correct the real-time adjustable capacity information based on the response capacity decay rate and the time interval from the last call.

[0007] The corrected real-time adjustable capability information is mapped to a multi-dimensional matching space to form resource feature coordinates; dimensional weight coefficients are assigned to each coordinate axis based on the dimensionality ratio of the task feature vector; the weighted feature distance between the resource feature coordinates and the task feature vector is calculated based on the dimensional weight coefficients to determine the matching degree; and resources with a matching degree exceeding a preset matching degree threshold are selected as candidate resources.

[0008] Calculate the electrical distance between each pair of candidate resources, and calculate the spatial clustering degree of each candidate resource based on the average electrical distance; select the candidate resource with the highest spatial clustering degree as the initial member of the temporary aggregation unit, and sequentially determine the average electrical distance between the remaining candidate resources and the members already included. If the average electrical distance is less than the preset distance threshold, it is included in the current temporary aggregation unit; otherwise, it is used as the initial member of the new temporary aggregation unit.

[0009] The steps for calculating the weighted feature distance between the resource feature coordinates and the task feature vector based on the aforementioned dimension weight coefficients to determine the matching degree include:

[0010] Extract the three dimensions of required response rate, required adjustment capacity and required duration from the task feature vector, calculate the ratio of each dimension component to the mean of the corresponding dimension as the normalization weight factor, and normalize the normalization weight factor as the dimension weight coefficient.

[0011] For each distributed flexible load resource, extract the coordinate values ​​of the response rate dimension, adjustment capacity dimension, and duration dimension from its resource feature coordinates, and calculate the weighted distance between each dimension coordinate value and the corresponding dimension component of the task feature vector based on the dimension weight coefficient.

[0012] The reciprocal of the weighted distance is calculated as the initial matching degree. The historical task completion rate is extracted from the historical response characteristic data. When the historical task completion rate is lower than the preset completion rate threshold, an attenuation coefficient is applied to the initial matching degree. The attenuation coefficient is positively correlated with the historical task completion rate. The initial matching degree after applying the attenuation coefficient is taken as the final matching degree of the resource.

[0013] The steps for constructing the expected response baseline and individual difference envelope of the temporary aggregation unit based on the historical response characteristic data include:

[0014] Extract the actual response curves of each distributed flexible load resource in the temporary aggregation unit in the historical scheduling tasks from the historical response characteristic data; group the actual response curves according to task type; align the actual response curves of each distributed flexible load resource according to the time axis for each task type and calculate the median value of the response at each moment; use the time series composed of the median value of the response as the expected response benchmark curve for that task type.

[0015] For each task type, the deviation of each actual response curve from the corresponding expected response baseline curve at each time point is calculated. The deviation distribution of each distributed flexible load resource in the temporary aggregation unit under this task type at each time point is statistically analyzed, and the upper and lower bounds of the deviation under the preset confidence level are extracted. The expected response baseline curve is superimposed with the upper bound of the deviation to form the upper envelope boundary, and the expected response baseline curve is superimposed with the lower bound of the deviation to form the lower envelope boundary. The upper and lower envelope boundaries constitute the individual difference envelope.

[0016] The steps of generating a main scheduling instruction for the temporary aggregation unit based on the expected response benchmark, estimating the execution deviation based on the individual difference envelope, and generating an individual compensation instruction based on the deviation contribution when the execution deviation exceeds the limit include:

[0017] Based on the current power grid dispatching task, the expected response baseline curve of the corresponding task type is matched, and the expected response baseline curve is divided into multiple control periods according to the time window. The expected response amount of each period is extracted as the power setpoint to form the main dispatching instruction.

[0018] Extract the upper and lower boundary values ​​of the individual difference envelope in each control period and calculate the envelope width. Determine the envelope symmetry based on the positive deviation of the upper boundary value from the expected response and the negative deviation of the lower boundary value from the expected response. When the envelope width exceeds a preset threshold or the envelope exhibits asymmetrical characteristics, it is marked as a high-risk period.

[0019] Extract the actual response deviation sequence of each member within the unit from historical response characteristic data, construct the correlation matrix between the individual deviation of each member and the overall deviation of the unit, calculate the amplification factor of the unit deviation on the overall deviation based on the correlation matrix, and multiply the amplification factor by the variance of the member's historical deviation to obtain the deviation contribution.

[0020] For the high-risk period, members within the unit are sorted according to their deviation contribution, and members whose cumulative contribution reaches a preset sharing threshold are selected as compensation targets. Compensation power is allocated according to their deviation contribution to form individual compensation instructions.

[0021] The steps of constructing a correlation matrix between the individual deviation and the overall deviation of each member, and calculating the amplification factor of the unit deviation on the overall deviation based on the correlation matrix, include:

[0022] Extract the actual response deviation time series of each member in the temporary aggregation unit in the same type of historical task from the historical response characteristic data, and calculate the overall deviation time series.

[0023] For each member, the individual deviation time series is used as the independent variable and the overall deviation time series of the unit is used as the dependent variable. Multiple time segments are constructed by using a sliding time window. The partial correlation coefficient between the individual deviation change and the overall deviation change is calculated in each time segment. The partial correlation coefficients of each time segment are used to form the association vector of the member. The association vectors of all members in the temporary aggregation unit are arranged in rows to form an association matrix.

[0024] The singular value decomposition of the correlation matrix is ​​performed to extract the right singular vector corresponding to the principal singular value. The inner product of the member correlation vector and the right singular vector is calculated as the coupling strength of the member to the overall deviation. The ratio of the coupling strength to the standard deviation of the member's historical deviation time series is used as the amplification factor.

[0025] The steps of issuing the main scheduling instruction and individual compensation instruction and collecting actual response data, and updating the matching degree weights in the multidimensional matching space and the individual difference envelope based on the actual response data include:

[0026] The main scheduling instruction is issued to each member in the temporary aggregation unit, and an additional individual compensation instruction is issued to the compensation object. The actual response power and response time of each member are collected to form actual response data.

[0027] The response accuracy index is calculated based on the actual response data and the power setting value of the main scheduling command. The variance of the actual response data in each control period is calculated as the response stability index. The execution quality level of each member is evaluated based on the weighted combination of the response accuracy index and the response stability index. The actual response data of members whose execution quality level exceeds the preset quality threshold are marked as high-quality data.

[0028] Adjust the weight coefficients of the corresponding dimensions in the multidimensional matching space based on the execution quality level of each member; after filtering high-quality data and adding it to the historical response characteristic data according to task type and removing abnormal data, recalculate the deviation distribution corresponding to the individual difference envelope, and reconstruct the individual difference envelope based on the updated deviation distribution.

[0029] The present invention constructs a desired response benchmark and individual difference envelope based on historical response characteristic data, enabling the system to accurately predict the response behavior of aggregation units and generate targeted main scheduling instructions and individual compensation instructions. This effectively improves the coordination and control capability of flexible load resources, reduces execution deviations caused by individual differences, and ensures the high-quality completion of scheduling tasks. Attached Figure Description

[0030] Figure 1 This is a flowchart illustrating the virtual power plant flexible load aggregation and dispatch control method according to an embodiment of the present invention;

[0031] Figure 2A flowchart for generating two-level instructions based on deviation prediction is provided. Detailed Implementation

[0032] Figure 1 This is a flowchart illustrating the virtual power plant flexible load aggregation and dispatch control method according to an embodiment of the present invention, as shown below. Figure 1 As shown, the method includes:

[0033] Acquire real-time adjustability information and historical response characteristic data of distributed flexible load resources; receive power grid dispatching tasks and extract task feature vectors;

[0034] Based on the real-time adjustable capability information, the matching degree between each distributed flexible load resource and the task feature vector is calculated in the multi-dimensional matching space. Resources with matching degrees exceeding a preset matching degree threshold are dynamically grouped into temporary aggregation units. Based on the historical response characteristic data, the expected response benchmark and individual difference envelope of the temporary aggregation units are constructed.

[0035] A main scheduling instruction for the temporary aggregation unit is generated based on the expected response benchmark. Execution deviation is estimated based on the individual difference envelope. When the execution deviation exceeds the limit, an individual compensation instruction is generated based on the deviation contribution.

[0036] The main scheduling instruction and individual compensation instruction are issued, and actual response data is collected. The matching degree weights in the multidimensional matching space and the individual difference envelope are updated based on the actual response data.

[0037] In the actual operation of the virtual power plant, edge computing nodes deployed on each distributed flexible load resource side collect adjustable capacity information in real time. This information includes parameters such as the current power and temperature setpoint margin of the air conditioning load, the state of charge and upper and lower limits of charging and discharging power of energy storage devices, the allowable interruption duration and recovery rate of industrial loads, and the remaining charging demand and adjustable time window of electric vehicle charging piles. Edge nodes periodically upload data packets to the cloud aggregation platform. The data packet fields include resource identifier, timestamp, adjustable power range, response latency, continuous capacity duration, and current operating status flag. Historical response characteristic data is stored in a time-series database, retaining response records for nearly 180 days. Each record includes fields such as the time of dispatch command issuance, command target value, actual response curve sampling point, response latency deviation, peak deviation, steady-state deviation, and response completion degree.

[0038] When the power grid dispatch center issues a dispatch task, the task message includes fields such as dispatch type identifier, target regulation power, expected response time, duration, power ramp rate constraint, and frequency regulation dead zone. After parsing the task message, a task feature vector is extracted. This vector consists of six normalized dimensions: the proportion of regulation to total system capacity, response time urgency index, normalized duration value, normalized ramp rate value, task priority weight, and frequency deviation sensitivity coefficient. The regulation proportion is calculated by dividing the target value by the current total online resource capacity, with a value ranging from 0.01 to 0.95. A value exceeding 0.95 triggers a capacity shortage warning. The response time urgency index is obtained by multiplying the reciprocal of the expected response time by a calibration coefficient of 100 seconds, resulting in an index of 10 for a 10-second response requirement and approximately 1.67 for a 60-second response requirement. The duration normalization value is normalized based on 3600 seconds. The ramp rate normalization value is calculated by dividing the system's rated capacity by kilowatts per second. Task priority weights are mapped from the priority field in the scheduling message: 1.5 for urgent tasks, 1.0 for regular tasks, and 0.7 for pre-scheduled tasks. The frequency deviation sensitivity coefficient is obtained by normalizing the inverse of the frequency regulation dead zone. For frequency regulation tasks, the coefficient is mapped to 1.0 to 5.0 based on the dead zone value; for non-frequency regulation tasks, this coefficient is 0.

[0039] In one optional implementation, the step of calculating the matching degree between each distributed flexible load resource and the task feature vector in a multi-dimensional matching space based on the real-time adjustable capability information, and dynamically grouping resources with matching degrees exceeding a preset matching degree threshold into temporary aggregation units includes:

[0040] Extract the historical call time and actual response capacity of each distributed flexible load resource from the historical response characteristic data, calculate the response capacity decay rate, and correct the real-time adjustable capacity information based on the response capacity decay rate and the time interval from the last call.

[0041] The corrected real-time adjustable capability information is mapped to a multi-dimensional matching space to form resource feature coordinates; dimensional weight coefficients are assigned to each coordinate axis based on the dimensionality ratio of the task feature vector; the weighted feature distance between the resource feature coordinates and the task feature vector is calculated based on the dimensional weight coefficients to determine the matching degree; and resources with a matching degree exceeding a preset matching degree threshold are selected as candidate resources.

[0042] Calculate the electrical distance between each pair of candidate resources, and calculate the spatial clustering degree of each candidate resource based on the average electrical distance; select the candidate resource with the highest spatial clustering degree as the initial member of the temporary aggregation unit, and sequentially determine the average electrical distance between the remaining candidate resources and the members already included. If the average electrical distance is less than the preset distance threshold, it is included in the current temporary aggregation unit; otherwise, it is used as the initial member of the new temporary aggregation unit.

[0043] For example, historical call times and actual response capabilities for each resource are extracted from historical response characteristic data. The data extraction range is set to the past 90 days, and at least 5 valid records are retained for each resource. The response capability decay rate is obtained by calculating the ratio of the change in response capability between two adjacent schedulings to the time interval. The median decay rate of multiple schedulings is taken to exclude the impact of abnormal fluctuations. The unit of calculation for the decay rate is kilowatts per day, and the value range is limited to between -5 kilowatts per day and +10 kilowatts per day.

[0044] The time interval between the current and last call is calculated by subtracting the time of the most recent call to the resource from the current time, with the time interval in days and rounded to two decimal places. When correcting the real-time adjustable capability information, the response capability decay rate is multiplied by the time interval to obtain the decay amount, which is then subtracted from the real-time adjustable capability information to obtain the corrected adjustable capability value. A lower limit protection is set during the correction process; the corrected adjustable capability must not be lower than 20% of the resource's rated capacity. If it falls below this threshold, it is forcibly rolled back to 20% of the rated capacity, triggering a resource performance degradation alarm. For resources that have not been called for more than 30 days, an additional long-term idle correction factor of 0.92 is applied, which multiplies the corrected adjustable capability by 0.92 to reflect the performance uncertainty caused by long-term idleness. After correction, the adjustable capability information is updated in the resource status cache table. The cache table uses a hash index structure with the resource identifier as the key and supports millisecond-level query response.

[0045] When the corrected real-time adjustable capacity information is mapped to the multi-dimensional matching space, the adjustable capacity information of each resource is converted into a six-dimensional coordinate vector. The first dimension is the ratio of the corrected adjustable capacity to the rated capacity of the resource, ranging from 0.20 to 1.00. The second dimension is the reciprocal of the historical average response delay of the resource multiplied by a calibration coefficient of 60 seconds; resources with a delay of less than 10 seconds are mapped to 6.0, and resources with a delay of 60 seconds are mapped to 1.0. The third dimension is the ratio of the resource's sustainable adjustment time to the standard time of 3600 seconds, ranging from 0.10 to 2.00. The fourth dimension is the historical average ramp rate of the resource divided by the system's unit capacity baseline ramp rate, with the baseline value set at 0.05 per second per kilowatt, ranging from 0.50 to 3.00. The fifth dimension is the resource's scheduling success rate over the past 30 days, defined as the proportion of scheduling times where the actual response volume reaches 85% or more of the target value out of the total number of scheduling times, ranging from 0.30 to 1.00. The sixth dimension is the resource's geographic location code, achieved by mapping the distribution area number where the resource is located to a range of 0 to 1. Resources in the same area have the same code value. The six-dimensional coordinate vector is stored in a multi-dimensional index structure, using a KD-tree partitioning method to accelerate neighborhood queries. The tree splitting threshold is set so that each leaf node contains no more than 8 resources.

[0046] The task feature vector is also constructed as a six-dimensional vector, with each dimension corresponding to the dimension definition of resource feature coordinates. The first dimension is the proportion of the task's adjustment capacity requirement to the total capacity; the second dimension is the urgency index of the response time requirement; the third dimension is the normalized value of the duration; the fourth dimension is the normalized value of the ramp rate; the fifth dimension is the task priority weight; and the sixth dimension is the geographic concentration coefficient of the task requirements. The geographic concentration coefficient is determined based on whether the task requires resources to be realized in a specific region: 1.0 for regional tasks and 0 for tasks across the entire network. The numerical ratios between the dimensions of the task feature vector are calculated. Dimensions with ratios greater than 2 are considered key constraint dimensions and assigned higher dimension weight coefficients. The dimension weight coefficient allocation strategy is as follows: key constraint dimensions are weighted at 0.25, secondary constraint dimensions at 0.15, and general dimensions at 0.10, with the total weight of the six dimensions normalized to 1.0. In actual allocation, the constraint level is determined by comparing the deviation of each dimension's value from its historical statistical mean; dimensions with deviations exceeding 1.5 standard deviations are judged as key constraints.

[0047] When calculating the weighted feature distance between resource feature coordinates and task feature vectors, the absolute value of the difference between the coordinate value and the vector value is calculated for each dimension. This difference is then multiplied by the corresponding dimension's weight coefficient and summed to obtain the weighted Manhattan distance. The weighted feature distance ranges from 0 to 6; a smaller distance indicates a better match between the resource and the task. The matching degree is calculated by inversely mapping the weighted feature distance using a linear transformation function. A distance of 0 results in a matching degree of 100, a distance of 6 results in a matching degree of 0, and intermediate values ​​are linearly interpolated proportionally. The matching degree calculation result is preserved to one decimal place, ranging from 0.0 to 100.0. The preset matching degree threshold is dynamically adjusted based on the total number of available resources and task capacity requirements. When available resources are sufficient, the threshold is set to 70.0; when resources are scarce, the threshold decreases to 55.0; and for urgent tasks, the threshold increases to 80.0. Resources with matching degrees exceeding the preset matching degree threshold are included in the candidate resource set. Candidate resources are sorted by matching degree from highest to lowest and stored in a priority queue with a maximum queue length of 200 to control subsequent computational overhead.

[0048] The electrical distance between each pair of candidate resources is obtained by querying the impedance matrix in the power grid topology database. The impedance matrix pre-stores the per-unit equivalent impedance values ​​between each node. The electrical distance is defined as the ohm value obtained by multiplying the per-unit equivalent impedance between two resource access nodes by the system reference impedance, ranging from 0.001 ohms to 50 ohms. For resources accessing the same node, the electrical distance is set to 0.001 ohms to avoid calculation errors caused by zero values. The electrical distance matrix is ​​constructed as a symmetric matrix and stored in memory. The matrix size is the square of the number of candidate resources. When there are more than 100 candidate resources, sparse matrix compression storage is used to reduce memory usage.

[0049] The spatial clustering degree of each candidate resource is obtained by calculating the average electrical distance between that resource and all other candidate resources. The smaller the average distance, the closer the resource is to the electrical center of the candidate resource group. The clustering degree value is the reciprocal of the average electrical distance multiplied by a calibration factor of 1000 ohms. The clustering degree value ranges from 20 to 1000. Resources with high clustering degree are preferentially selected as the initial members of the temporary aggregation unit. The resource with the highest clustering degree is found by traversing the candidate resource set, marked as the initial member of the first temporary aggregation unit, and removed from the candidate resource set.

[0050] The remaining candidate resources are evaluated one by one to determine whether to include them in the current temporary aggregation unit. During the evaluation process, the electrical distance between the resource to be evaluated and all existing members in the current temporary aggregation unit is calculated, and the arithmetic mean of these electrical distances is obtained. The preset distance threshold is set according to the system voltage level and network topology density: 5 ohms for 10 kV distribution network, 15 ohms for 35 kV distribution network, and 30 ohms for 110 kV transmission network. When the average electrical distance is less than the preset distance threshold, the resource to be evaluated meets the electrical compactness requirement, is included in the current temporary aggregation unit, and the unit member list is updated. When the average electrical distance is greater than or equal to the preset distance threshold, the resource to be evaluated is too far from the current unit, and the resource is used as the initial member of a new temporary aggregation unit to create a new unit. After the new unit is created, the remaining candidate resources are evaluated again. The average electrical distance of each candidate resource is calculated with all existing temporary aggregation units in turn, and the unit with the smallest average electrical distance that is less than the preset distance threshold is selected for inclusion. If the average electrical distance with all units exceeds the threshold, it becomes a new unit on its own.

[0051] The entire resource pooling process employs a greedy, iterative strategy, processing one candidate resource per iteration. The iteration terminates when the candidate resource set is empty or the preset maximum number of units is reached. The maximum number of units is set to 20; exceeding this limit stops creating new units, and remaining candidate resources are forcibly allocated to the nearest existing unit. After the temporary aggregation units are assembled, the system calculates the number of members and total adjustable capacity for each unit. Units with fewer than 3 members or a total capacity less than 5% of the task requirements are marked as low-quality units and removed. Resources within these removed units are returned to the candidate resource pool for reallocation.

[0052] This invention improves the accuracy of resource assessment by correcting real-time adjustable capability information through response capability decay rate, and reduces network loss and voltage fluctuation risks during scheduling by using electrical distance constraints to form spatially compact temporary aggregation units.

[0053] In an optional implementation, the step of calculating the weighted feature distance between the resource feature coordinates and the task feature vector based on the dimension weight coefficients to determine the matching degree includes:

[0054] Extract the three dimensions of required response rate, required adjustment capacity and required duration from the task feature vector, calculate the ratio of each dimension component to the mean of the corresponding dimension as the normalization weight factor, and normalize the normalization weight factor as the dimension weight coefficient.

[0055] For each distributed flexible load resource, extract the coordinate values ​​of the response rate dimension, adjustment capacity dimension, and duration dimension from its resource feature coordinates, and calculate the weighted distance between each dimension coordinate value and the corresponding dimension component of the task feature vector based on the dimension weight coefficient.

[0056] The reciprocal of the weighted distance is calculated as the initial matching degree. The historical task completion rate is extracted from the historical response characteristic data. When the historical task completion rate is lower than the preset completion rate threshold, an attenuation coefficient is applied to the initial matching degree. The attenuation coefficient is positively correlated with the historical task completion rate. The initial matching degree after applying the attenuation coefficient is taken as the final matching degree of the resource.

[0057] For example, three dimensions—required response rate, required regulation capacity, and required duration—are extracted from the task feature vector data structure. The task feature vector is stored as a floating-point array with six elements. The second element corresponds to the required response rate in kilowatts per second (kW), ranging from 0.01 to 5.00; the first element corresponds to the required regulation capacity in kilowatts (kW), ranging from 10 to 50000; and the third element corresponds to the required duration in seconds (seconds), ranging from 60 to 7200. The extraction operation is implemented through array index access: the response rate is read from index 1, the regulation capacity from index 0, and the duration from index 2. The extraction results are cached in a temporary variable for subsequent calculations.

[0058] The ratios of each dimensional component to its corresponding dimensional mean are calculated using dimensional mean parameters stored in the historical task statistics database. These dimensional mean parameters are updated every 24 hours and are obtained by calculating the arithmetic mean of the dimensional components of all scheduled tasks over the past 30 days. The default dimensional mean for response rate is 0.8 kW / s, for regulation capacity it is 1200 kW, and for duration it is 1800 seconds. The ratios are calculated by dividing the current task's dimensional component by its corresponding dimensional mean. The response rate ratio is obtained by dividing the required response rate by the mean response rate; the regulation capacity ratio is obtained by dividing the required regulation capacity by the mean regulation capacity; and the duration ratio is obtained by dividing the required duration by the mean duration. These three ratios are denoted as response rate ratio, regulation capacity ratio, and duration ratio, respectively. Their values ​​typically range from 0.01 to 10.00. Outliers exceeding 10.00 are truncated to 10.00, and outliers below 0.01 are truncated to 0.01.

[0059] The normalized weighting factor is achieved by scaling the three ratios proportionally until their sum equals 1.0. During the calculation, the three ratios are summed to obtain the total value. The response rate ratio divided by the total value yields the response rate normalized weighting factor, the accommodation capacity ratio divided by the total value yields the accommodation capacity normalized weighting factor, and the duration ratio divided by the total value yields the duration normalized weighting factor. The normalized weighting factor is numerically accurate to four decimal places, the sum of the three factors is strictly equal to 1.0000, and the error tolerance is ±0.0001. The normalized weighting factors are directly used as dimension weight coefficients. The response rate dimension weight coefficient, accommodation capacity dimension weight coefficient, and duration dimension weight coefficient each correspond to one of the three normalized weighting factors. These weight coefficients are stored in the weight vector structure for use in matching degree calculations.

[0060] Resource characteristic coordinates are extracted from the resource status database, with one resource characteristic record maintained for each distributed flexible load resource. The resource characteristic record includes fields such as response rate dimension coordinates, adjustment capacity dimension coordinates, and duration dimension coordinates. These coordinate values ​​are calculated by fusing real-time monitoring data with historical performance data. The response rate dimension coordinate reflects the amount of power the resource can change per unit time after receiving a scheduling command, measured in kilowatts per second, ranging from 0.01 to 8.00. This value is obtained from the average ramp rate of the resource's last 10 scheduling responses. The adjustment capacity dimension coordinate reflects the maximum adjustable power the resource can currently provide, measured in kilowatts, ranging from 20% to 100% of the resource's rated capacity. This value is calculated from the resource's current operating status and remaining adjustment margin. The duration dimension coordinate reflects the maximum time the resource can maintain its adjusted state, measured in seconds, ranging from 60 to 14400. This value is determined based on the resource type and energy storage capacity; for energy storage resources, it is calculated based on remaining electricity and discharge power; for load resources, it is set based on the longest interruption time allowed by the user.

[0061] The weighted distance calculation is performed separately for each of the three dimensions. The weighted distance for the response rate dimension is obtained by multiplying the absolute value of the difference between the resource response rate dimension coordinate value and the required response rate of the task by the response rate dimension weight coefficient. The weighted distance for the adjustment capacity dimension is obtained by multiplying the absolute value of the difference between the resource adjustment capacity dimension coordinate value and the required adjustment capacity of the task by the adjustment capacity dimension weight coefficient. The weighted distance for the duration dimension is obtained by multiplying the absolute value of the difference between the resource duration dimension coordinate value and the required duration of the task by the duration dimension weight coefficient. The sum of the weighted distances of the three dimensions yields the total weighted distance, which is a dimensionless numerical value ranging from 0 to 50000. The smaller the value, the closer the resource characteristics are to the task requirements. During the calculation, the absolute values ​​of the differences are rounded to two decimal places, the weight coefficients to four decimal places, the weighted distances to two decimal places, and the total weighted distance to two decimal places.

[0062] The initial matching degree is obtained by calculating the reciprocal of the total weighted distance. Before calculating the reciprocal, an offset of 1.0 is applied to the total weighted distance to avoid an infinite result caused by zero distance. The corrected calculation method is to add 1.0 to the total weighted distance and then take the reciprocal. The initial matching degree ranges from 0.00002 to 1.00000, with a larger value indicating a higher degree of matching. The initial matching degree is multiplied by a calibration coefficient of 100 to obtain a percentage-based initial matching degree score, ranging from 0.002 to 100.000, retaining three decimal places of precision. The percentage-based initial matching degree score is stored in a temporary calculation cache as a baseline value for subsequent attenuation processing.

[0063] Historical task completion rates are extracted from the resource historical response characteristic data table, which records all scheduled tasks participated in by each resource within the past 90 days and their completion status. Task completion status is determined by comparing the actual response volume of the resource with the required response volume of the task. A task is considered complete if the actual response volume reaches 90% or higher of the required response volume; otherwise, it is considered incomplete. The historical task completion rate is calculated by dividing the number of completed tasks by the total number of tasks, with a value ranging from 0.00 to 1.00, rounded to two decimal places. The completion rate statistics window is 90 days. Resources with fewer than 5 tasks within the window are marked as having insufficient data for completion rate calculations, and a default completion rate of 0.75 is used for subsequent calculations. Completion rate data is updated in real time after each scheduled task ends, using a sliding window mechanism to maintain the latest 90-day record. Historical records older than 90 days are archived to offline storage and do not participate in real-time calculations.

[0064] The preset completion rate threshold is set to 0.85, which serves as the dividing line for determining whether to apply attenuation to the historical performance of a resource. When the historical task completion rate is lower than the preset completion rate threshold of 0.85, an attenuation coefficient is applied to the initial matching degree, reducing the final matching degree score of the resource. The attenuation coefficient is positively correlated with the historical task completion rate, and this positive correlation is achieved through a linear mapping function. A completion rate of 0.00 corresponds to an attenuation coefficient of 0.50, a completion rate of 0.85 corresponds to an attenuation coefficient of 1.00, and intermediate values ​​are calculated using linear interpolation. The specific calculation method is: attenuation coefficient = 0.50 + (historical task completion rate ÷ 0.85) × 0.50. This formula ensures that the lower the completion rate, the smaller the attenuation coefficient, thus weakening the initial matching degree more significantly. The attenuation coefficient ranges from 0.50 to 1.00, retaining four decimal places of precision. When the completion rate is equal to or higher than 0.85, the attenuation coefficient is constant at 1.00, and no attenuation is applied.

[0065] The initial matching degree after applying the attenuation coefficient is calculated by multiplying the initial matching degree score (out of 100) by the attenuation coefficient. The result serves as the final matching degree for the resource. The final matching degree ranges from 0.001 to 100.000, retaining three decimal places. A higher value indicates a better overall resource matching degree. The final matching degree calculation result is written to the resource matching degree result table. This table uses the resource identifier as the primary key, the matching degree as the numeric field, and also records auxiliary fields such as the calculation timestamp, task identifier, initial matching degree, and attenuation coefficient for traceability and auditing. The matching degree result table supports fast queries sorted in descending order of matching degree, using a B-tree index to accelerate the sorting operation, with query response time controlled within 10 milliseconds.

[0066] This invention dynamically determines the dimension weight coefficient by the ratio of the task dimension component to the historical mean, so that the matching degree calculation adapts to the differences in task features. Based on the inverse of the weighted distance, the matching degree between resources and tasks is quantified, and a decay coefficient driven by the historical task completion rate is introduced to punish resources with poor performance.

[0067] In one optional implementation, the step of constructing the expected response baseline and individual difference envelope of the temporary aggregation unit based on the historical response characteristic data includes:

[0068] Extract the actual response curves of each distributed flexible load resource in the temporary aggregation unit in the historical scheduling tasks from the historical response characteristic data; group the actual response curves according to task type; align the actual response curves of each distributed flexible load resource according to the time axis for each task type and calculate the median value of the response at each moment; use the time series composed of the median value of the response as the expected response benchmark curve for that task type.

[0069] For each task type, the deviation of each actual response curve from the corresponding expected response baseline curve at each time point is calculated. The deviation distribution of each distributed flexible load resource in the temporary aggregation unit under this task type at each time point is statistically analyzed, and the upper and lower bounds of the deviation under the preset confidence level are extracted. The expected response baseline curve is superimposed with the upper bound of the deviation to form the upper envelope boundary, and the expected response baseline curve is superimposed with the lower bound of the deviation to form the lower envelope boundary. The upper and lower envelope boundaries constitute the individual difference envelope.

[0070] For example, the actual response curves of each resource within a temporary aggregation unit in historical scheduling tasks are extracted from the historical response characteristic database. The actual response curves are stored in time-series format, including a time field and a response quantity field. The time field records the relative time since the scheduling instruction was issued, in seconds, and the response quantity field records the actual power output or reduction of the resource at that time, in kilowatts. The historical time window is set to the past 180 days by default, and the window length can be adjusted according to data sufficiency, with a minimum of 60 days and a maximum of 365 days. During data extraction, invalid records with a duration of less than 30 seconds or fewer than 10 data sampling points are filtered out. At least 5 valid curves are required for each resource.

[0071] Actual response curves are grouped by task type. Task types are read from the scheduling task metadata table, which records the type identifier field for each scheduled task. The task type identifier uses enumeration coding: peak shaving tasks are coded as 1, valley filling tasks as 2, frequency adjustment tasks as 3, reserve capacity call as 4, and demand response as 5. The grouping operation clusters all actual response curves according to their task type identifier, grouping curves of the same type into the same set. Each set contains response curves from all historical scheduling of all resources under that type, typically ranging from 50 to 500 curves. After grouping, the number of curves for each task type is counted. Task types with fewer than 20 curves are marked as insufficient samples; for these types, no expected response baseline curve or individual difference envelope is constructed, and a system-preset general baseline curve is used instead.

[0072] For each task type, the actual response curves of each distributed flexible load resource are aligned along the time axis. The alignment operation needs to address the issue of inconsistent sampling times for different curves. A unified time grid is defined, with a starting point of 0 seconds and an ending point equal to the longest duration in the historical records for that task type. The grid step size is set to 10 seconds. For curve data whose sampling time is not at a grid point, a linear interpolation method is used to calculate the estimated response value at the grid point. During interpolation, the two closest actual sampling points to the grid point are taken, and a weighted average is performed based on time distance weights to obtain the grid point response value. The interpolation accuracy is maintained at 0.01 kW. If curve data is missing or interrupted during interpolation, the grid point response value for the missing period is marked as null and not included in subsequent statistics. After alignment, each curve is converted into a fixed-length vector, with a vector length equal to the number of grid points, and vector elements representing the response value or null value marker at each grid point.

[0073] When calculating the median response value at each time step, valid response values ​​of all curves at each grid point are collected, and null values ​​are excluded to form a response sample set. The number of elements in the sample set equals the number of valid curves at that time step, and the sample set is sorted in ascending order of numerical value. The median value calculation is handled differently depending on the number of elements in the sample set: if the number of elements is odd, the value of the middle element after sorting is taken; if the number of elements is even, the arithmetic mean of the two middle elements after sorting is taken. The median value is kept with a precision of 0.01 kW, and the calculation result is stored in a median value vector, where the vector index corresponds to the grid point number, and the vector element is the median response value at that time step. For some grid points where the sample set has fewer than 5 elements, the median value at that time step is marked as unreliable, and the expected response baseline curve at that time step is filled by linear interpolation of the medians of the reliable times before and after.

[0074] A time series based on the median response value is used as the expected response baseline curve for this task type. The expected response baseline curve is stored as a key-value pair sequence of time and median value. The time key corresponds to the grid point time of a unified time grid, and the median value corresponds to the median response value at that grid point. The expected response baseline curve data structure supports fast median value lookup by time, using a hash table to achieve millisecond-level access performance. Each task type maintains an independent expected response baseline curve, and the curve data is persistently stored in a baseline curve configuration table. The table structure includes fields such as task type identifier, grid point time, median response value, sample size, and update timestamp. The baseline curve is updated weekly, with the most recent 180 days of historical data re-extracted and recalculated during updates. The update operation is performed during off-peak business hours to avoid impacting real-time scheduling performance.

[0075] The deviation of each actual response curve from the corresponding expected response baseline curve at each time step is defined as: the response value of the actual response curve at a given time step minus the median value of the expected response baseline curve at that time step. A positive deviation indicates that the actual response is higher than the expected baseline, while a negative deviation indicates that the actual response is lower than the expected baseline. The deviation is calculated for each grid point of each curve individually. The response value at that time step is read from the aligned vector of the actual response curve, and the median value at that time step is read from the expected response baseline curve. The difference between the two yields the deviation. The deviation accuracy is maintained at 0.01 kW. The deviation vector has the same dimension as the actual response curve vector, and the vector elements are the deviation values ​​at each time step. If the actual response curve has a null value at a given time step, the deviation value at that time step is not calculated and is directly marked as null.

[0076] This task involves statistically analyzing the deviation distribution of distributed flexible load resources within a temporary aggregation unit at various time points. The statistical operation is performed independently for each grid point. The deviation distribution at a given grid point is obtained by collecting the effective deviation values ​​of all curves at that grid point, forming a deviation sample set. The number of elements in the sample set equals the number of effective deviations at that time. The deviation sample set is sorted in ascending order of numerical value, and the quantile characteristics of the sample set are calculated. The preset confidence level is set to 95%, corresponding to the 2.5% quantile and the 97.5% quantile. The lower bound of the deviation is extracted as the deviation value corresponding to the 2.5% quantile, and the upper bound is extracted as the deviation value corresponding to the 97.5% quantile. Quantile calculation uses a linear interpolation method. The position corresponding to the 2.5% quantile is the sample set length multiplied by 0.025. When the position is not an integer, linear interpolation is performed between adjacent integer positions to obtain the quantile value. The upper and lower bounds of the deviation are kept with an accuracy of 0.01 kW. Each grid point corresponds to a pair of upper and lower bound values. The upper bound values ​​of all grid points constitute the upper bound vector of the deviation, and the lower bound values ​​of all grid points constitute the lower bound vector of the deviation.

[0077] The upper envelope boundary is formed by superimposing the upper bound of the expected response baseline curve with the upper bound of the deviation. The superposition operation is a time-by-time addition operation. For each grid point, the median value at that time is read from the expected response baseline curve, and the upper bound of the deviation at that time is read from the upper bound vector of the deviation. The two are added together to obtain the upper envelope boundary value at that time. The accuracy of the upper envelope boundary value is retained to 0.01 kW. The upper envelope boundary values ​​at all times constitute the upper envelope boundary curve. The lower envelope boundary is formed by superimposing the lower bound of the expected response baseline curve with the lower bound of the deviation. The superposition operation is also a time-by-time addition operation. For each grid point, the median value at that time is read from the expected response baseline curve, and the lower bound of the deviation at that time is read from the lower bound vector of the deviation. The two are added together to obtain the lower envelope boundary value at that time. The accuracy of the lower envelope boundary value is retained to 0.01 kW. The lower envelope boundary values ​​at all times constitute the lower envelope boundary curve.

[0078] The upper and lower envelope boundaries constitute the individual-variable envelope, which is stored in a hyperbolic data structure, containing two sub-curves: the upper boundary curve and the lower boundary curve. Each sub-curve is a key-value pair sequence of time and boundary value, with the time corresponding to a grid point in a unified time grid, and the boundary value corresponding to the envelope boundary response at that time. The individual-variable envelope data is persistently stored in an envelope configuration table, whose structure includes fields such as task type identifier, aggregation unit identifier, grid point time, upper boundary value, lower boundary value, confidence level, sample size, and update timestamp. The envelope configuration table supports joint queries by task type and aggregation unit identifier, with a query response time controlled within 50 milliseconds. The individual-variable envelope is updated weekly in sync with the expected response baseline curve. During updates, a transaction mechanism is used to ensure consistency between the baseline curve and the envelope data; if an update fails, it rolls back to the previous version.

[0079] This invention provides a reliable basis for assessing response capabilities and risk control boundaries for scheduling decisions, significantly improving the accuracy and robustness of virtual power plant aggregation scheduling.

[0080] In one optional implementation, the steps of generating a main scheduling instruction for the temporary aggregation unit based on the expected response benchmark, estimating the execution deviation based on the individual difference envelope, and generating an individual compensation instruction based on the deviation contribution when the execution deviation exceeds the limit include:

[0081] Based on the current power grid dispatching task, the expected response baseline curve of the corresponding task type is matched, and the expected response baseline curve is divided into multiple control periods according to the time window. The expected response amount of each period is extracted as the power setpoint to form the main dispatching instruction.

[0082] Extract the upper and lower boundary values ​​of the individual difference envelope in each control period and calculate the envelope width. Determine the envelope symmetry based on the positive deviation of the upper boundary value from the expected response and the negative deviation of the lower boundary value from the expected response. When the envelope width exceeds a preset threshold or the envelope exhibits asymmetrical characteristics, it is marked as a high-risk period.

[0083] Extract the actual response deviation sequence of each member within the unit from historical response characteristic data, construct the correlation matrix between the individual deviation of each member and the overall deviation of the unit, calculate the amplification factor of the unit deviation on the overall deviation based on the correlation matrix, and multiply the amplification factor by the variance of the member's historical deviation to obtain the deviation contribution.

[0084] For the high-risk period, members within the unit are sorted according to their deviation contribution, and members whose cumulative contribution reaches a preset sharing threshold are selected as compensation targets. Compensation power is allocated according to their deviation contribution to form individual compensation instructions.

[0085] Combination Figure 2The flowchart for generating a two-layer instruction based on deviation prediction is explained below. Based on the current power grid dispatch task, the task type identifier field is read from the task metadata table. The task type identifier is stored in an enumerated encoding format and a matching relationship is established with the task type identifier field in the expected response baseline curve configuration table. The matching operation is achieved by executing a database join query, with the query condition being that the task type identifiers are equal. The query result returns the complete dataset of the expected response baseline curve for the corresponding task type. The expected response baseline curve dataset includes a time field and a median response value field. The time field records the relative time from the start of dispatch, in seconds, with a precision of 10 seconds, and a value range from 0 to the preset duration of the task. The median response value field records the expected response power at that time, in kilowatts, with a precision of 0.01 kilowatts. If matching fails, a task type not configured alarm is triggered, and the dispatch task is suspended pending manual intervention or switching to a general baseline curve.

[0086] The expected response baseline curve is divided into multiple control periods by time windows. The time window length is set according to the characteristics of the task type. The default time window for peak shaving and valley filling tasks is 300 seconds, the default time window for frequency regulation tasks is 60 seconds, and the default time window for reserve capacity mobilization is 180 seconds. The time window length can be configured from 30 seconds to 600 seconds, and the configuration parameter is stored in the time window field of the task type configuration table. The division operation starts from the scheduling start time and divides the time axis of the expected response baseline curve into continuous control periods with the time window length as the step size. Each control period contains several grid points. The control period numbers are incremented from 1 in chronological order, and the number of control periods is equal to the task duration divided by the time window length and rounded up. The actual length of the last control period is less than the time window length, and the effective duration of this period is based on the task end time.

[0087] The expected response value for each time period is extracted as the power setpoint. This extraction operation is performed on all grid points within each control time period. The expected response value for a given control time period is calculated as the arithmetic mean of the median response values ​​of all grid points within that time period. The average value is calculated by summing the median response values ​​of all grid points within the time period, and then dividing the sum by the number of grid points in the time period to obtain the expected response value for that time period. The accuracy of the expected response value for each time period is maintained at 0.01 kW. The calculation results are stored in a control time period array, with the array index being the time period number and the array element being the power setpoint for that time period. The power setpoint ranges from 0 to the total rated capacity of the temporary aggregation unit. Values ​​exceeding this range trigger a capacity over-limit alarm and are clipped to the capacity boundary.

[0088] The main scheduling instruction is generated from a control time period array. The main scheduling instruction data structure includes fields such as instruction identifier, temporary aggregation unit identifier, task type identifier, scheduling start time, and control time period sequence. The control time period sequence is a nested structure, with each time period element containing sub-fields such as time period number, time period start time, time period end time, and power setpoint. The time period start time and end time are absolute timestamps with second-level precision. The start time = scheduling start time + (time period number - 1) × time window length, and the end time equals the start time plus the time window length. The main scheduling instruction is sent to the temporary aggregation unit controller using a message queue mechanism. The message format uses key-value pair encoding, and message priority is set according to task type: frequency regulation tasks have the highest priority of 1, peak shaving and valley filling tasks have a priority of 2, and other tasks have a priority of 3. The instruction issuance timeout is set to 5 seconds. If no acknowledgment response is received within the timeout period, a retransmission mechanism is triggered, with a maximum of 3 retransmissions and a retransmission interval of 2 seconds.

[0089] Extract the upper and lower boundary values ​​of the individual difference envelope for each control period. The extraction operation reads the upper and lower boundary curve data from the envelope configuration table. For each control period, collect the upper and lower boundary values ​​of all grid points within the period. The upper boundary value of the period is calculated as the arithmetic mean of the upper boundary values ​​of all grid points within the period, and the lower boundary value of the period is calculated as the arithmetic mean of the lower boundary values ​​of all grid points within the period. The accuracy of the average value is retained to 0.01 kW. The calculation results are stored in a period envelope array, and the array elements contain fields such as period number, upper boundary value of the period, and lower boundary value of the period.

[0090] The envelope width is calculated by subtracting the lower boundary value from the upper boundary value of the time period. The envelope width reflects the uncertainty range of the resource response within that time period, measured in kilowatts (kW) with a precision of 0.01 kW. The envelope width is stored in the envelope width field of the time period envelope array. The positive deviation magnitude is calculated by subtracting the expected response amount from the upper boundary value of the time period, while the negative deviation magnitude is calculated by subtracting the lower boundary value from the expected response amount. Both positive and negative deviation magnitudes have a precision of 0.01 kW and are stored in the positive and negative deviation fields of the time period envelope array, respectively.

[0091] Envelope symmetry is determined based on the ratio of positive to negative deviation. The symmetry index is calculated by dividing the positive deviation by the negative deviation. A symmetry index close to 1 indicates a symmetrical envelope, while a significant deviation from 1 indicates an asymmetrical envelope. The symmetry threshold is set between 0.7 and 1.3. A symmetry index within this range is considered a symmetrical envelope, while an index less than 0.7 or greater than 1.3 is considered an asymmetrical envelope. The symmetry determination result is stored in the symmetry flag field of the time-segment envelope array. The symmetry flag is a Boolean type, with a true value indicating symmetry and a false value indicating asymmetry.

[0092] Envelope width exceeding limits is determined based on a preset threshold, which is set according to the task type and the rated capacity of the aggregation unit. The default envelope width threshold for peak shaving and valley filling tasks is 15% of the total rated capacity of the aggregation unit; for frequency modulation tasks, it's 8%; and for standby capacity call-up, it's 12%. The threshold can be configured from 5% to 25% of the total rated capacity, and the configuration parameters are stored in the envelope width threshold field of the task type configuration table. Envelope width exceeding limits is determined by comparing the time-period envelope width with the preset threshold; if the envelope width exceeds the preset threshold, it is considered exceeding the limit.

[0093] High-risk periods are marked as high-risk periods if either the envelope width exceeds the limit or the envelope exhibits asymmetric characteristics. High-risk periods are stored in the risk flag field of the period envelope array. The risk flag is a Boolean type; a true value indicates high risk, and a false value indicates low risk. High-risk periods trigger the individual compensation instruction generation process; low-risk periods do not generate compensation instructions but only execute the main scheduling instruction.

[0094] The actual response deviation sequence of each member within the unit is extracted from the historical response characteristic database. The deviation sequence is time series data, recording the deviation between the actual response and the expected response at each moment in the historical scheduling task. The deviation value equals the actual response minus the median value of the expected response baseline curve at the corresponding moment. A positive deviation value indicates overshoot, and a negative deviation value indicates undershoot. The deviation sequence extraction range is historical tasks of the same type as the current task within the past 90 days. Each member must extract at least 10 deviation sequences. Members that do not meet the condition are excluded from the compensation target candidate pool.

[0095] For each member, an association matrix is ​​constructed between their individual deviation and the overall deviation of the unit. The dimension of the association matrix is ​​the number of historical tasks multiplied by the number of time grid points. The elements of the individual deviation matrix are the deviation values ​​of a member at a specific time in a given task, while the elements of the overall unit deviation matrix are the sum of the deviation values ​​of all members within the unit at that time in the same task. The association matrix is ​​constructed element-wise, with the matrix row index corresponding to the historical task number, the matrix column index corresponding to the time grid point number, and the matrix elements being binary tuples containing both individual and overall deviations.

[0096] The amplification factor is calculated based on linear regression analysis using the correlation matrix, with the overall unit deviation as the dependent variable and the individual deviation as the independent variable. The slope parameter obtained from the regression fit is the amplification factor of the unit deviation on the overall deviation. The amplification factor reflects the strength of the influence of the member deviation on the overall unit deviation. The amplification factor is reserved with a precision of 0.001, and its value range is usually between 0.01 and 2.00. An amplification factor greater than 1 indicates that the member deviation has an amplifying effect on the overall deviation, while an amplification factor less than 1 indicates that the member deviation has a weak influence on the overall deviation. The amplification factor is stored in the amplification factor field of the member characteristic table, which stores the amplification factor under different task types in separate columns.

[0097] The historical deviation variance of a member is calculated as the sample variance of all historical deviation values ​​for that member, with a variance precision of 0.01, and the unit is kilowatt squared. A larger variance indicates poorer member response stability, while a smaller variance indicates better member response stability. The historical deviation variance is stored in the deviation variance field of the member characteristics table.

[0098] The deviation contribution is calculated by multiplying the amplification factor by the historical deviation variance, with a precision of 0.01, and the unit is kilowatt squared. The deviation contribution comprehensively considers the intensity and fluctuation range of the member's deviation; a higher contribution indicates a more significant contribution of that member to the overall deviation of the unit. The deviation contribution is stored in the contribution field of the member characteristic table, and this field is updated in real time before each compensation instruction is generated to ensure data timeliness.

[0099] For high-risk periods, members within a unit are sorted by their deviation contribution, and a member sorting list is generated by sorting the members in descending order of their numerical values. Each element in the sorting list contains a member identifier and a deviation contribution.

[0100] The preset sharing threshold is set to 80% of the total deviation contribution of all members within the unit. The sharing threshold can be configured from 60% to 95%, and the configuration parameter is stored in the sharing threshold field of the compensation strategy configuration table. Compensation objects are selected starting from the first member of the sorted list, accumulating the members' deviation contributions until the cumulative contribution reaches or exceeds the sharing threshold. Members selected during the accumulation process constitute the compensation object set. The size of the compensation object set is typically 20% to 50% of the total number of members in the unit, containing a minimum of one member and a maximum of 70% of the total number of members.

[0101] The compensation power is allocated proportionally based on the deviation contribution. The compensation power for a given member is equal to the envelope width during the high-risk period multiplied by that member's deviation contribution, divided by the sum of the contributions of the entire set of compensation objects. The accuracy of the compensation power is retained at 0.01 kW. The compensation direction is determined based on the envelope symmetry index. When the symmetry index is less than 1, the compensation direction is positive, increasing power output; when the symmetry index is greater than 1, the compensation direction is negative, decreasing power output; when the symmetry index is equal to 1, the compensation direction is dynamically adjusted based on real-time deviation monitoring.

[0102] The individual compensation instruction data structure includes fields such as instruction identifier, member identifier, control period number, compensation power amount, compensation direction, and execution priority. Individual compensation instructions and main scheduling instructions are issued in parallel through a message queue. The compensation instruction has a higher priority than the main scheduling instruction, and the execution sequence is that the compensation instruction arrives at the member controller before the main scheduling instruction. The compensation instruction issuance timeout is set to 3 seconds. If the timeout occurs, an alarm is triggered, and the member is removed from the current compensation object set. The compensation power amount is then redistributed to the remaining compensation objects.

[0103] This invention generates master dispatch instructions based on the expected response baseline curve to provide clear power setting targets, identifies high-risk periods by predicting execution deviations through individual difference envelopes, quantifies the influence weight of each member on the overall deviation based on deviation contribution metric, and generates individual compensation instructions for high-contribution members to achieve precise control, significantly reducing aggregate response deviations and improving the dispatch execution accuracy of temporary aggregation units and the quality of power grid dispatch task completion.

[0104] In one optional implementation, the step of constructing a correlation matrix between the individual deviation and the overall deviation of the unit for each member, and calculating the amplification factor of the unit deviation on the overall deviation based on the correlation matrix, includes:

[0105] Extract the actual response deviation time series of each member in the temporary aggregation unit in the same type of historical task from the historical response characteristic data, and calculate the overall deviation time series.

[0106] For each member, the individual deviation time series is used as the independent variable and the overall deviation time series of the unit is used as the dependent variable. Multiple time segments are constructed by using a sliding time window. The partial correlation coefficient between the individual deviation change and the overall deviation change is calculated in each time segment. The partial correlation coefficients of each time segment are used to form the association vector of the member. The association vectors of all members in the temporary aggregation unit are arranged in rows to form an association matrix.

[0107] The singular value decomposition of the correlation matrix is ​​performed to extract the right singular vector corresponding to the principal singular value. The inner product of the member correlation vector and the right singular vector is calculated as the coupling strength of the member to the overall deviation. The ratio of the coupling strength to the standard deviation of the member's historical deviation time series is used as the amplification factor.

[0108] For example, the actual response deviation time series of each member within a temporary aggregation unit in historical tasks of the same type is extracted from the historical response characteristic database. This extraction is achieved by executing a database query. The query conditions include fields such as a list of temporary aggregation unit member identifiers, task type identifiers, and the start and end times of the historical time window. The task type identifier is consistent with the currently scheduled task type. The historical time window is set to the past 90 days by default, and the window length can be adjusted according to data sufficiency, with a minimum of 30 days and a maximum of 180 days. The actual response deviation time series records the deviation value of each member at each moment in each historical task. The deviation value is defined as the actual response power at that moment minus the median value of the expected response baseline curve at the corresponding moment, in kilowatts, with a precision of 0.01 kilowatts. The deviation time series is stored as key-value pairs of timestamps and deviation values. The timestamps are expressed in absolute time with a precision of seconds, and the deviation values ​​range from the negative to the positive interval of the member's rated capacity. The extracted results are grouped by member identifier. Each member corresponds to a deviation time series set, which contains the deviation sequences of that member in all historical tasks of the same type. The set size requires that each member contain at least 10 valid deviation sequences. Members that do not meet the conditions are excluded from subsequent calculations and a data shortage alarm is triggered.

[0109] The overall deviation time series calculation is achieved by summing the deviation values ​​of all members within the unit at the same time. The calculation process iterates through each historical task time-by-time, reading the deviation values ​​of all members at that time, and summing them to obtain the overall deviation value for that time. The overall deviation value retains a precision of 0.01 kW, and its range is from negative to positive of the total rated capacity of the aggregated unit. The time axis of the overall deviation time series is aligned with that of the individual deviation time series, with a one-to-one correspondence of timestamps. If some members' deviation values ​​are missing at a certain time, only the valid deviation values ​​are summed at that time; missing members are considered to have a deviation of zero and are not included in the summation. The overall deviation time series is stored in the overall deviation table of the historical characteristics database. The table structure includes fields such as task identifier, timestamp, overall deviation value, and number of valid members, supporting fast querying by task identifier and timestamp range.

[0110] For each member, their individual deviation time series is used as the independent variable, and the overall deviation time series of the unit is used as the dependent variable for correlation analysis. The individual deviation time series is read line by line from the member deviation sequence set, and the overall deviation time series is read from the record corresponding to the task identifier in the overall deviation table. The two are linked through a timestamp field. The independent variable vector consists of the deviation values ​​of the individual deviation sequences arranged in chronological order, and the dependent variable vector consists of the deviation values ​​of the overall deviation sequence arranged in chronological order. The vector length is equal to the number of sampling points for that task, typically between 100 and 3600.

[0111] A sliding time window constructs multiple time segments. The length of the sliding time window is set according to the task duration: 60 seconds for tasks lasting less than 600 seconds, 120 seconds for tasks lasting 600 to 1800 seconds, and 300 seconds for tasks lasting more than 1800 seconds. The window length is configurable from 30 to 600 seconds, and the configuration parameter is stored in the "Time Window Length" field of the association analysis configuration table. The sliding step size is set to 50% of the time window length, and the step size is configurable from 25% to 75% of the time window length. The configuration parameter is stored in the "Sliding Step Size" field of the association analysis configuration table. The sliding operation starts from the beginning of the time series and sequentially extracts segments of the time window length at intervals of the sliding step size until the right boundary of the window reaches the end of the time series. The number of time segments equals the total length of the time series minus the time window length, divided by the sliding step size, rounded down, and then incremented by 1. The actual length of the last segment is less than the time window length, and the effective range of this segment is based on the end of the series.

[0112] Within each time segment, the individual deviation change and the overall deviation change are calculated. The change is defined as the first-order difference sequence of the deviation values ​​within the segment. The individual deviation change sequence is obtained by subtracting adjacent values ​​from the individual deviation value sequence within the segment. The change at a given time moment equals the deviation value at that time moment minus the deviation value at the previous time moment. The length of the change sequence equals the segment length minus 1. The overall deviation change sequence is calculated in the same way, by subtracting adjacent values ​​from the overall deviation value sequence within the segment. The accuracy of the change is retained to 0.01 kW, and the value range is the maximum change amplitude under the constraint of the rated capacity change rate.

[0113] Partial correlation coefficients measure the strength of the linear correlation between individual deviation changes and overall deviation changes, while controlling for the influence of deviation changes in other members. The partial correlation coefficient uses a partial correlation extension of the Pearson correlation coefficient, and its calculation requires constructing a covariance matrix and an inverse covariance matrix. The covariance matrix has a dimension equal to the number of members plus one, and its elements are the covariance values ​​between the individual deviation changes of each member and the overall deviation changes. Covariance = Sum of the deviation products of the two variation sequences ÷ (Sequence length - 1), where the deviation product is the product of the two sequences at each time step after subtracting their respective means. The covariance precision is retained to 0.001, and the covariance matrix is ​​a symmetric positive definite matrix. The inverse covariance matrix is ​​obtained through matrix inversion, using Gaussian elimination or Cholleysky decomposition. Numerical stability is guaranteed by a condition number test; if the condition number exceeds 1000, a numerical instability alarm is triggered, and regularization is applied. The partial correlation coefficient is extracted from the inverse covariance matrix. The partial correlation coefficient of a member is calculated as: -(the element in the row corresponding to the member and the column corresponding to the overall deviation in the inverse covariance matrix) ÷ [(the square root of the diagonal element corresponding to the member × the diagonal element corresponding to the overall deviation)]. The partial correlation coefficient is kept to a precision of 0.001 and ranges from -1 to +1. An absolute value close to 1 indicates a strong correlation, while an absolute value close to 0 indicates a weak correlation.

[0114] The partial correlation coefficients of each time segment constitute the association vector of a member. The length of the association vector is equal to the number of time segments, and the vector elements are the partial correlation coefficients of each segment. The association vector is stored in the member association characteristic table, which includes fields such as member identifier, task identifier, segment number, partial correlation coefficient, and segment start and end time. For segments whose partial correlation coefficient calculation fails or has abnormal values, the partial correlation coefficient of that segment is set to 0 to indicate no association. The failure count of segments is accumulated, and when the proportion of failed segments exceeds 20%, the association vector of that member is marked as unreliable and excluded from subsequent calculations.

[0115] The association vectors of all members within a temporary aggregation unit are arranged in rows to form an association matrix. The number of rows in the association matrix equals the number of members within the unit, the number of columns equals the number of time segments, and the matrix elements are the partial correlation coefficients of a specific member within a specific time segment. The association matrix is ​​stored using a two-dimensional array data structure, where row indices correspond to member identifiers, column indices correspond to segment numbers, and element values ​​are partial correlation coefficients. The matrix size is typically 8 to 50 rows multiplied by 10 to 100 columns.

[0116] Singular value decomposition (SVD) of the incidence matrix extracts the right singular vectors corresponding to the principal singular values. SVD decomposes the incidence matrix into a product of a left singular matrix, a singular value diagonal matrix, and a right singular matrix. The decomposition operation satisfies the following conditions: the incidence matrix equals the left singular matrix multiplied by the singular value diagonal matrix, then multiplied by the transpose of the right singular matrix. The number of rows in the left singular matrix equals the number of rows (members) in the incidence matrix, and the number of columns equals the matrix rank. The number of rows in the right singular matrix equals the number of columns (time segments) in the incidence matrix, and the number of columns equals the matrix rank. The SVD algorithm is implemented using either the Jacobian rotation or a divide-and-conquer approach. The computation library calls the SVD interface from the linear algebra computation library. Interface parameters include the input matrix, decomposition mode, and numerical precision. The decomposition mode is set to full decomposition or economic decomposition. Full decomposition returns the full-size singular matrix, while economic decomposition only returns valid columns within the rank range. The default is economic decomposition mode to reduce computational load. The numerical precision is set to double-precision floating-point arithmetic, the iteration convergence threshold is 0.00001, and the maximum number of iterations is 500. A decomposition failure alarm is triggered if convergence is not achieved after exceeding the iteration limit. The diagonal elements of the singular value diagonal matrix are singular values, arranged in descending order. The singular value precision is retained at 0.0001, and the value range is non-negative real numbers. The principal singular value is defined as the element with the largest value in the singular value sequence, and the right singular vector corresponding to the principal singular value is the first column vector of the right singular matrix. The length of the right singular vector is equal to the number of columns in the incidence matrix, i.e., the number of time segments. The vector elements are real numbers, and the vector magnitude is normalized to 1. Normalization is achieved by dividing each element by the square root of the sum of the squares of all elements, retaining an element precision of 0.0001.

[0117] The inner product of a member's association vector and its right singular vector represents the coupling strength of the member's deviation from the overall pattern. This inner product is calculated by summing element-wise multiplications. A member's association vector is represented by a row in the association matrix, and the right singular vector is the extracted principal singular vector. The inner product calculation iterates through each position of the vectors, multiplying the element of the member's association vector at that position with the element of the right singular vector at that position. All these multiplications are summed to obtain the inner product value. The inner product value is kept to a precision of 0.001 and ranges from the negative vector magnitude to the positive vector magnitude. A positive inner product indicates that the member's deviation changes in the same direction as the dominant pattern, while a negative inner product indicates that the member's deviation changes in the opposite direction to the dominant pattern. The absolute value of the inner product reflects the coupling strength. The coupling strength is stored in the coupling strength field of the member characteristics table, which records the coupling strength under different task types.

[0118] The standard deviation of a member's historical deviation time series is calculated as the sample standard deviation of all historical deviation values ​​for that member. The standard deviation calculation involves three steps: subtracting the mean deviation from each deviation value in the deviation series to obtain the deviation residuals; summing the squared residuals to obtain the sum of squares; dividing the sum of squares by (sample size - 1) to obtain the variance; and taking the square root of the variance to obtain the standard deviation. The standard deviation is rounded to 0.01 kW. A larger standard deviation indicates greater fluctuation in the member's response, while a smaller standard deviation indicates smaller fluctuations. The historical deviation standard deviation is stored in the deviation standard deviation field of the member characteristics table.

[0119] The amplification factor is calculated by dividing the coupling strength by the historical deviation standard deviation, with the division operation precision retained to 0.001. The amplification factor is dimensionless. It reflects the contribution weight of the member unit deviation standard deviation to the overall deviation dominance pattern. An amplification factor greater than 1 indicates that the member deviation has an amplifying effect on the overall deviation; a factor less than 1 indicates that the member deviation has a weak impact on the overall deviation; and a negative amplification factor indicates that the member deviation has an inverse effect on the overall deviation. The amplification factor typically ranges from -2 to +2. Values ​​outside this range trigger outlier checks. If the check passes, the calculation result is retained; if the check fails, the amplification factor is set to 0 and marked as invalid. The amplification factor is stored in the amplification factor field of the member characteristic table. This field is updated in real-time before each compensation instruction is generated to ensure data timeliness.

[0120] This invention constructs multiple time segments using a sliding time window to calculate partial correlation coefficients, capturing the dynamic correlation characteristics between individual deviations and overall deviations at different time periods. The coupling strength is normalized by combining the historical deviation standard deviation to obtain the amplification coefficient, accurately assessing the contribution weight of each member to the overall deviation of the aggregation unit, and providing a scientific basis for the selection of compensation objects and the allocation of compensation power.

[0121] In one optional implementation, the steps of issuing the main scheduling instruction and individual compensation instruction and collecting actual response data, and updating the matching degree weights in the multidimensional matching space and the individual difference envelope based on the actual response data include:

[0122] The main scheduling instruction is issued to each member in the temporary aggregation unit, and an additional individual compensation instruction is issued to the compensation object. The actual response power and response time of each member are collected to form actual response data.

[0123] The response accuracy index is calculated based on the actual response data and the power setting value of the main scheduling command. The variance of the actual response data in each control period is calculated as the response stability index. The execution quality level of each member is evaluated based on the weighted combination of the response accuracy index and the response stability index. The actual response data of members whose execution quality level exceeds the preset quality threshold are marked as high-quality data.

[0124] Adjust the weight coefficients of the corresponding dimensions in the multidimensional matching space based on the execution quality level of each member; after filtering high-quality data and adding it to the historical response characteristic data according to task type and removing abnormal data, recalculate the deviation distribution corresponding to the individual difference envelope, and reconstruct the individual difference envelope based on the updated deviation distribution.

[0125] For example, a main scheduling instruction is issued to each member within the temporary aggregation unit, and this instruction issuance is achieved through a message queue mechanism. The main scheduling instruction message format includes fields such as instruction identifier, member identifier, task type identifier, and control time period sequence. The control time period sequence includes subfields such as time period number, time period start time, time period end time, and power setting value. Message priority is set according to the task type. After receiving the message, the member controller returns an acknowledgment response. The scheduling control center sets the response timeout limit to 5 seconds. If no response is received within the timeout period, a retransmission mechanism is triggered, with a maximum of 3 retransmissions.

[0126] Additional individual compensation instructions are issued to the compensation targets. The compensation targets read the member identifier list from the compensation target set. The individual compensation instruction message format includes fields such as instruction identifier, member identifier, control period number, compensation power amount, compensation direction, and execution priority. The compensation power amount is in kilowatts with a precision of 0.01 kilowatts. The compensation direction uses enumeration encoding: a positive increase in power is encoded as 1, and a negative decrease in power is encoded as -1. The execution priority is higher than the main scheduling instruction; the smaller the priority value, the higher the priority. The compensation instruction priority is set to 1, and the main scheduling instruction priority is set to 2. The compensation instruction is issued before the main scheduling instruction arrives at the member controller, with a time interval set to 10 seconds, ensuring that the compensation instruction executes before the main scheduling instruction. The compensation instruction and the main scheduling instruction are associated through the instruction identifier field. When the member controller executes, it adds the compensation power amount to the power setting value of the main scheduling instruction, with the addition direction determined by the compensation direction field. The compensation instruction response timeout is set to 3 seconds. Upon timeout, an alarm is triggered, and the member is removed from the compensation target set. The compensation power amount is then redistributed to the remaining compensation targets.

[0127] Actual response data is generated by collecting the actual response power and response time of each member device. Data acquisition is achieved through the power sampling module of the member controller. The power sampling module reads the power sensor output of the member device, with a sampling frequency set to 1 Hz (1 sample per second) and a sampling accuracy of 0.01 kW. The sampling time is recorded as an absolute timestamp with second-level accuracy. The actual response power is the actual output power or reduced power of the member device at the sampling time, ranging from 0 to the member's rated capacity. The collected actual response data is stored in time-series format, including timestamp and actual power fields. The actual response data is uploaded to the dispatch control center periodically via a data reporting interface, with a reporting cycle set to 10 seconds. After receiving the data, the dispatch control center stores it in the actual response data table, which includes fields such as member identifier, task identifier, timestamp, actual power, and sampling quality flag.

[0128] During actual response data acquisition, the system synchronously records the measured values ​​of response delay and ramp rate for each member. Response delay is defined as the time difference between the member controller receiving the main scheduling command and the start of power change. The measured ramp rate is the power change per unit time during the power change process. Response delay is calculated by comparing the command issuance timestamp with the timestamp of the first power deviation from the baseline value exceeding 2% of the rated capacity. The delay accuracy is in the second range, with a value between 0 and 300 seconds. The measured ramp rate is obtained by performing a differential operation on the power time series. The power difference between adjacent sampling points divided by the sampling time interval yields the instantaneous ramp rate. The median value of all instantaneous ramp rates within a time period is used as the measured ramp rate for that period, in kilowatts per second, with an accuracy of 0.01 kilowatts per second. The measured values ​​of response delay and ramp rate are used to verify the accuracy of the member's corresponding dimensional coordinate values ​​in the multidimensional matching space. When the deviation between the measured value and the coordinate value exceeds 20%, a coordinate value correction process is triggered. The coordinate value correction adopts an exponential moving average strategy. The new coordinate value = measured value × 0.3 + original coordinate value × 0.7. The corrected coordinate value is updated to the resource feature database and takes effect when the next aggregation unit is built.

[0129] The response accuracy index is calculated based on the actual response data and the power setpoint of the main dispatch command. The calculation process is performed one by one for each control period. The response accuracy for a certain control period is calculated by inverting and normalizing the average of the absolute values ​​of the deviations between the actual response power and the power setpoint within that period. The absolute value of the deviation is calculated by subtracting the power setpoint from the actual response power and taking the absolute value, in kilowatts, with a precision of 0.01 kilowatts. The sum of the absolute values ​​of the deviations at each sampling point within the period is divided by the number of sampling points in the period to obtain the average absolute value of the deviation. The average absolute value of the deviation is divided by the power setpoint to obtain the relative deviation, which ranges from 0 to positive infinity. The smaller the relative deviation, the higher the accuracy. The response accuracy index is defined as 1 minus the relative deviation, ranging from negative infinity to 1. An accuracy index close to 1 indicates accurate response, while an accuracy index close to 0 or a negative value indicates a large response deviation. The accuracy index precision is retained at 0.001 and stored in the accuracy field of the member execution quality table.

[0130] When calculating the response accuracy index, differentiated deviation tolerances are applied for different task types. For frequency regulation tasks, due to their high accuracy requirements, the relative deviation threshold is set at 5%. When this threshold is exceeded, the accuracy index is calculated with double the tolerance, resulting in an accuracy index of 1 - (relative deviation × 2). The relative deviation threshold is set at 10% for peak shaving and valley filling tasks, 8% for reserve capacity call-up tasks, and 12% for demand response tasks. This differentiated tolerance setting ensures that the quality assessment better aligns with the actual needs of different dispatch scenarios in the virtual power plant. The threshold parameters are stored in the deviation threshold field of the task type configuration table. For response performance with a relative deviation below the threshold, the accuracy index is calculated using the standard formula. For response performance with a relative deviation exceeding the threshold, the accuracy index for task types other than frequency regulation tasks is calculated with a 1.5-fold tolerance, resulting in an accuracy index of 1 - (relative deviation × 1.5). This differentiated attenuation mechanism ensures that the execution quality level assessment results accurately reflect the response adaptability of members in specific task scenarios.

[0131] The response stability index is calculated as the variance of the actual response data across each control period. Variance calculation is performed on the actual power sampling points within each control period. The variance calculation for a given control period involves three steps: subtracting the mean power value of each actual power value within the period to obtain the power residual; summing the squared power residuals to obtain the sum of squares; and dividing the sum of squares by (number of sampling points within the period - 1) to obtain the variance. The mean power value is the arithmetic mean of all actual power values ​​within the period. The power residual is in kilowatts (kW) with a precision of 0.01 kW. The variance is in the square of kilowatts (kW) with a precision of 0.01 kW. A smaller variance indicates better response stability, while a larger variance indicates more severe response fluctuations. The response stability index is defined as the variance after inversion and normalization. Normalization is achieved by dividing by the square of the rated capacity: Stability Index = -(Variance ÷ Rated Capacity) 2The stability index ranges from -1 to 0. A stability index close to 0 indicates stability, while a stability index close to -1 indicates instability. The precision of the stability index is reserved to 0.001 and is stored in the stability field of the member execution quality table.

[0132] The execution quality level assessment is based on a weighted combination of response accuracy and response stability indicators. The weighted combination is calculated as: Accuracy Indicator × Accuracy Weight + Stability Indicator × Stability Weight. The default weight for accuracy is 0.6, and the default weight for stability is 0.4, with a total weight of 1. The weights can be configured from 0.3 to 0.7, and the configuration parameters are stored in the weight field of the quality assessment configuration table. The weighted combination result is a comprehensive quality score, with a precision of 0.001, ranging from -1 to 1. A higher score indicates better execution quality. Execution quality levels are determined based on the comprehensive quality score: a score greater than or equal to 0.8 is excellent, 0.6 to 0.8 is good, 0.4 to 0.6 is acceptable, and a score below 0.4 is unacceptable. Level codes are represented by integers: excellent is 4, good is 3, acceptable is 2, and unacceptable is 1. The execution quality level is stored in the quality level field of the member execution quality table, which records the quality assessment result for each task according to the task identifier.

[0133] The preset quality threshold is set to the lower limit of the overall quality score corresponding to the "good" level, which is 0.6. The threshold can be configured from 0.5 to 0.8, and the configuration parameter is stored in the "Quality Threshold" field of the quality assessment configuration table. The condition for determining whether the execution quality level exceeds the preset quality threshold is that the overall quality score is greater than or equal to 0.6. The actual response data of members meeting this condition are marked as high-quality data. The high-quality data marker is stored in the "Quality Marker" field of the actual response data table. The quality marker is a Boolean type; a true value indicates high quality, and a false value indicates low-quality or anomalous data. High-quality data is used for updating historical response characteristic data and reconstructing individual difference envelopes, while low-quality data is not updated and is only used for anomaly analysis.

[0134] The adjustment of the corresponding dimension weight coefficients in the multidimensional matching space is based on the execution quality level of each member. The adjustment operation is performed on the dimension weight of each member in the multidimensional matching space. The weight adjustment factor of a member in a certain dimension is calculated according to the execution quality level: 1.1 for excellent level, 1.05 for good level, 1.0 for qualified level, and 0.9 for unqualified level. The adjustment factor is multiplied by the current weight of the member in that dimension to obtain the adjusted weight. The precision of the adjusted weight is retained to 0.001, and the value range is from 0.001 to 2.0. Values ​​outside the range are pruned to the boundary value. Weight adjustment is performed after each task is completed. The adjusted weight is stored in the dimension weight field of the member matching characteristic table, which records the real-time weight value of each dimension. The weight adjustment adopts an exponential moving average smoothing strategy. New weight = adjusted weight × smoothing coefficient + historical weight × (1 - smoothing coefficient). The smoothing coefficient is set to 0.3 by default, and the configurable range is 0.1 to 0.5. The configuration parameter is stored in the smoothing coefficient field of the weight update configuration table.

[0135] During the adjustment of dimensional weight coefficients, the system prioritizes adjusting the weights of dimensions most relevant to the current task type. For frequency regulation tasks, the weight adjustment magnitude for the response rate and ramp rate dimensions is 1.5 times that of other dimensions, meaning the adjustment factor is multiplied by 1.5. For peak shaving and valley filling tasks, the weight adjustment magnitude for the duration and capacity dimensions is 1.5 times that of other dimensions. For standby capacity call-up tasks, the weight adjustment magnitude for the capacity and reliability dimensions is 1.5 times that of other dimensions. For demand response tasks, the weight adjustment magnitude for the response rate and reliability dimensions is 1.5 times that of other dimensions. This differentiated adjustment strategy ensures that the weight evolution direction of the multi-dimensional matching space matches the task characteristic distribution in the actual operation of the virtual power plant, improving the targeting of subsequent aggregation unit construction. The identification of priority adjustment dimensions is achieved by reading the key dimension fields of the task type configuration table. This field stores the key dimension identifiers corresponding to each task type in array form. During adjustment, the array is traversed to determine whether the current dimension is a key dimension.

[0136] High-quality data is appended to the historical response characteristic database according to task type. The append operation filters records marked as true from the actual response data table, reading fields such as member identifier, task type identifier, timestamp, and actual power. Appended data is stored in the historical response characteristic data table, which includes fields such as member identifier, task type identifier, task identifier, timestamp, actual power, deviation value, and data version. The deviation value is calculated by subtracting the median value of the expected response baseline curve at the corresponding time from the actual power, with a precision of 0.01 kilowatts. The data version uses an incrementing integer encoding; the version number is incremented by 1 each time data is appended, for data backtracking and version management. The append operation uses a transaction mechanism to ensure data consistency. A transaction consists of two steps: data insertion and index update; if either step fails, the entire transaction is rolled back.

[0137] Anomaly removal identifies and deletes abnormal records from the historical response characteristic data table, using a statistical threshold method. For a specific member and task type, the historical deviation sequence is used to calculate the mean and standard deviation. Records with deviation values ​​less than (mean - 3 times the standard deviation) or greater than (mean + 3 times the standard deviation) are marked as anomalies. Anomaly markers are stored in the anomaly marker field of the historical response characteristic data table. The anomaly marker is a Boolean type; a true value indicates anomaly, and a false value indicates normal behavior. The removal operation deletes records marked as true from the data table, backing them up to an anomaly data archive table before deletion for subsequent analysis. The removal operation is performed after data is appended, and the removal ratio is controlled to within 5% of the total historical data. If it exceeds 5%, a data quality alarm is triggered, and the removal operation is paused pending manual review.

[0138] The updated deviation distribution is recalculated based on historical response characteristic data after removing outliers, and the calculation process is the same as the initial individual difference envelope construction. Deviation values ​​of members in similar historical tasks are extracted from the historical response characteristic data table and aligned by time grid points to form a deviation matrix. For each grid point, deviation values ​​from all members across all historical tasks are collected to form a deviation sample set, which is sorted in ascending order of numerical value. The confidence level is preset to 95%, and the 2.5% and 97.5% quantiles are extracted as the lower and upper bounds of the deviation, respectively. Quantile calculation uses linear interpolation; the position corresponding to the 2.5% quantile is equal to the sample set length × 0.025. When the position is not an integer, linear interpolation is performed between adjacent integer positions to obtain the quantile value, with a quantile precision of 0.01 kW. The updated deviation lower and upper bound vectors are stored in the envelope update cache table, which includes fields such as task type identifier, aggregation unit identifier, grid point time, deviation lower bound, deviation upper bound, and update timestamp.

[0139] Individual difference envelope reconstruction is based on the updated deviation distribution. The reconstruction operation reads the lower and upper bound vectors of deviation from the envelope update cache table and the median vector from the expected response baseline curve configuration table. The upper envelope boundary curve is calculated by adding the upper bound of deviation to the expected response baseline curve at each time step, and the lower envelope boundary curve is calculated by adding the lower bound of deviation to the expected response baseline curve at each time step. The precision of the upper and lower envelope boundaries is retained at 0.01 kW, and the boundary values ​​are stored in the envelope configuration table, overwriting the original envelope data. Envelope updates use a version control mechanism. Before updating, the current envelope data is backed up to the historical version table, with the version number incremented. If the update fails, it is rolled back to the previous version from the historical version table. After the envelope update is completed, envelope validity verification is triggered. The verification checks whether the upper envelope boundary is greater than the lower envelope boundary at all times and whether the envelope width is within a reasonable range. If the verification fails, the update is rolled back and an alarm is triggered.

[0140] After the individual envelope update is completed, the system calculates the width change rate of the new and old envelopes. The width change rate is defined as the new envelope width minus the old envelope width, divided by the old envelope width. The envelope width is the average difference between the upper and lower envelope boundaries at each time point, in kilowatts, with a precision of 0.01 kilowatts. The width change rate is kept to a precision of 0.001, ranging from -1 to positive infinity. A negative width change rate indicates improved unit response consistency, while a positive rate indicates increased unit response dispersion. When the absolute value of the width change rate exceeds 15%, a stability assessment of the aggregated unit is triggered. The assessment process calculates the standard deviation of the execution quality level of each member within the unit for the most recent 5 tasks. A standard deviation greater than 0.5 indicates an unstable unit, while a standard deviation less than or equal to 0.5 indicates a stable unit. If the assessment result is unstable, the system reassembles a temporary aggregated unit in the next task. The reassembly operation marks the current unit as failed. When receiving the next scheduling task, the matching degree is recalculated and electrical distance clustering is performed again from the candidate resource pool to ensure the dynamic adaptability of the virtual power plant's aggregated scheduling. When the evaluation result is stable, the current unit configuration is retained. The next time a task of the same type is performed, the combination of members of this unit will be reused first. The reuse strategy is implemented by setting priority flags in the unit configuration table. Stable units are given high priority, newly formed units are given medium priority, and failed units are given low priority.

[0141] This invention evaluates the execution quality level through actual response data and dynamically adjusts the weights of the multidimensional matching space. It filters high-quality data to update the historical characteristic database and reconstructs the individual difference envelope, thereby achieving closed-loop optimization and continuously improving the accuracy of temporary aggregation unit construction and scheduling execution quality.

[0142] A second aspect of this invention provides a virtual power plant flexible load aggregation and dispatch control system, comprising: a data acquisition module, configured to acquire real-time adjustability information and historical response characteristic data of distributed flexible load resources; receive power grid dispatch tasks and extract task feature vectors; a resource aggregation module, configured to calculate the matching degree between each distributed flexible load resource and the task feature vector in a multi-dimensional matching space based on the real-time adjustability information, and dynamically group resources with matching degrees exceeding a preset matching degree threshold into temporary aggregation units; construct the expected response benchmark and individual difference envelope of the temporary aggregation unit based on the historical response characteristic data; an instruction generation module, configured to generate a main dispatch instruction for the temporary aggregation unit based on the expected response benchmark, estimate the execution deviation based on the individual difference envelope, and generate an individual compensation instruction based on the deviation contribution degree when the execution deviation exceeds the limit; and a feedback optimization module, configured to issue the main dispatch instruction and individual compensation instruction and collect actual response data, and update the matching degree weights and individual difference envelopes in the multi-dimensional matching space based on the actual response data.

[0143] A third aspect of the present invention provides an electronic device, comprising: a processor; and a memory for storing processor-executable instructions; wherein the processor is configured to invoke the instructions stored in the memory to execute the aforementioned method.

[0144] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0145] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

Claims

1. A method for flexible load aggregation and dispatch control in virtual power plants, characterized in that, include: Obtain real-time adjustability information and historical response characteristic data of distributed flexible load resources; Receive power grid dispatching tasks and extract task feature vectors; Based on the real-time adjustable capability information, the matching degree between each distributed flexible load resource and the task feature vector is calculated in the multi-dimensional matching space. Resources with matching degrees exceeding a preset matching degree threshold are dynamically grouped into temporary aggregation units. Specifically, this includes: extracting the historical call time and actual response capability of each distributed flexible load resource from the historical response characteristic data, calculating the response capability decay rate, and correcting the real-time adjustable capability information based on the response capability decay rate and the time interval from the current call to the last call; mapping the corrected real-time adjustable capability information to the multi-dimensional matching space to form resource feature coordinates; assigning dimension weight coefficients to each coordinate axis according to the dimension ratio of the task feature vector, calculating the weighted feature distance between the resource feature coordinates and the task feature vector based on the dimension weight coefficients to determine the matching degree, and screening resources with matching degrees exceeding a preset matching degree threshold as candidate resources; calculating the electrical distance between each pair of candidate resources, and calculating the spatial clustering degree of each candidate resource based on the average electrical distance; selecting the candidate resource with the highest spatial clustering degree as the initial member of the temporary aggregation unit, and sequentially judging the average electrical distance between the remaining candidate resources and the already included members. When the average electrical distance is less than a preset distance threshold, it is included in the current temporary aggregation unit; otherwise, it is used as the initial member of a new temporary aggregation unit. Based on the historical response characteristic data, construct the expected response baseline and individual difference envelope of the temporary aggregation unit; A main scheduling instruction for the temporary aggregation unit is generated based on the expected response benchmark. Execution deviation is estimated based on the individual difference envelope. When the execution deviation exceeds the limit, an individual compensation instruction is generated based on the deviation contribution. The main scheduling instruction and individual compensation instruction are issued, and actual response data is collected. The matching degree weights in the multidimensional matching space and the individual difference envelope are updated based on the actual response data.

2. The method according to claim 1, characterized in that, The steps for calculating the weighted feature distance between the resource feature coordinates and the task feature vector based on the aforementioned dimension weight coefficients to determine the matching degree include: Extract the three dimensions of required response rate, required adjustment capacity and required duration from the task feature vector, calculate the ratio of each dimension component to the mean of the corresponding dimension as the normalization weight factor, and normalize the normalization weight factor as the dimension weight coefficient. For each distributed flexible load resource, extract the coordinate values ​​of the response rate dimension, adjustment capacity dimension, and duration dimension from its resource feature coordinates, and calculate the weighted distance between each dimension coordinate value and the corresponding dimension component of the task feature vector based on the dimension weight coefficient. The reciprocal of the weighted distance is calculated as the initial matching degree. The historical task completion rate is extracted from the historical response characteristic data. When the historical task completion rate is lower than the preset completion rate threshold, an attenuation coefficient is applied to the initial matching degree. The attenuation coefficient is positively correlated with the historical task completion rate. The initial matching degree after applying the attenuation coefficient is taken as the final matching degree of the resource.

3. The method according to claim 1, characterized in that, The steps for constructing the expected response baseline and individual difference envelope of the temporary aggregation unit based on the historical response characteristic data include: Extract the actual response curves of each distributed flexible load resource within the temporary aggregation unit in historical scheduling tasks from the historical response characteristic data; The actual response curves are grouped according to task type. For each task type, the actual response curves of each distributed flexible load resource are aligned by the time axis and the median value of the response at each time point is calculated. The time series composed of the median value of the response is used as the expected response benchmark curve for that task type. For each task type, calculate the deviation of each actual response curve from the corresponding expected response baseline curve at each time point, statistically analyze the distribution of deviation of each distributed flexible load resource in the temporary aggregation unit at each time point under this task type, and extract the upper and lower bounds of the deviation under the preset confidence level. The upper envelope boundary is formed by superimposing the expected response baseline curve with the upper bound of the deviation, and the lower envelope boundary is formed by superimposing the expected response baseline curve with the lower bound of the deviation. The upper envelope boundary and the lower envelope boundary constitute the individual difference envelope.

4. The method according to claim 1, characterized in that, The steps of generating a main scheduling instruction for the temporary aggregation unit based on the expected response benchmark, estimating the execution deviation based on the individual difference envelope, and generating an individual compensation instruction based on the deviation contribution when the execution deviation exceeds the limit include: Based on the current power grid dispatching task, the expected response baseline curve of the corresponding task type is matched, and the expected response baseline curve is divided into multiple control periods according to the time window. The expected response amount of each period is extracted as the power setpoint to form the main dispatching instruction. Extract the upper and lower boundary values ​​of the individual difference envelope in each control period and calculate the envelope width. Determine the envelope symmetry based on the positive deviation of the upper boundary value from the expected response and the negative deviation of the lower boundary value from the expected response. When the envelope width exceeds a preset threshold or the envelope exhibits asymmetric characteristics, it is marked as a high-risk period. Extract the actual response deviation sequence of each member within the unit from historical response characteristic data, construct the correlation matrix between the individual deviation of each member and the overall deviation of the unit, calculate the amplification factor of the unit deviation on the overall deviation based on the correlation matrix, and multiply the amplification factor by the variance of the member's historical deviation to obtain the deviation contribution. For the high-risk period, members within the unit are sorted according to their deviation contribution, and members whose cumulative contribution reaches a preset sharing threshold are selected as compensation targets. Compensation power is allocated according to their deviation contribution to form individual compensation instructions.

5. The method according to claim 4, characterized in that, The steps of constructing a correlation matrix between the individual deviation and the overall deviation of each member, and calculating the amplification factor of the unit deviation on the overall deviation based on the correlation matrix, include: Extract the actual response deviation time series of each member in the temporary aggregation unit in the same type of historical task from the historical response characteristic data, and calculate the overall deviation time series. For each member, the individual deviation time series is used as the independent variable and the overall deviation time series of the unit is used as the dependent variable. Multiple time segments are constructed by using a sliding time window. The partial correlation coefficient between the individual deviation change and the overall deviation change is calculated in each time segment. The partial correlation coefficients of each time segment are used to form the association vector of the member. The association vectors of all members in the temporary aggregation unit are arranged in rows to form an association matrix. The singular value decomposition of the correlation matrix is ​​performed to extract the right singular vector corresponding to the principal singular value. The inner product of the member correlation vector and the right singular vector is calculated as the coupling strength of the member to the overall deviation. The ratio of the coupling strength to the standard deviation of the member's historical deviation time series is used as the amplification factor.

6. The method according to claim 1, characterized in that, The steps of issuing the main scheduling instruction and individual compensation instruction and collecting actual response data, and updating the matching degree weights in the multidimensional matching space and the individual difference envelope based on the actual response data include: The main scheduling instruction is issued to each member in the temporary aggregation unit, and an additional individual compensation instruction is issued to the compensation object. The actual response power and response time of each member are collected to form actual response data. The response accuracy index is calculated based on the actual response data and the power setting value of the main scheduling command, and the variance of the actual response data in each control period is calculated as the response stability index. The execution quality level of each member is evaluated based on a weighted combination of the response accuracy index and the response stability index. The actual response data of members whose execution quality level exceeds the preset quality threshold are marked as high-quality data. Adjust the weight coefficients of the corresponding dimensions in the multidimensional matching space based on the execution quality level of each member; After filtering high-quality data and adding it to the historical response characteristic data according to task type, and removing abnormal data, the deviation distribution corresponding to the individual difference envelope is recalculated, and the individual difference envelope is reconstructed based on the updated deviation distribution.

7. A virtual power plant flexible load aggregation and dispatch control system, used to implement the method of any one of claims 1-6, characterized in that, include: The data acquisition module is used to acquire real-time adjustability information and historical response characteristic data of distributed flexible load resources; Receive power grid dispatching tasks and extract task feature vectors; The resource aggregation module is used to calculate the matching degree between each distributed flexible load resource and the task feature vector in a multi-dimensional matching space based on the real-time adjustable capability information, and dynamically organize resources with matching degrees exceeding a preset matching degree threshold into temporary aggregation units. Based on the historical response characteristic data, construct the expected response benchmark and individual difference envelope of the temporary aggregation unit; The instruction generation module is used to generate a main scheduling instruction for the temporary aggregation unit based on the expected response benchmark, estimate the execution deviation based on the individual difference envelope, and generate an individual compensation instruction based on the deviation contribution degree when the execution deviation exceeds the limit. The feedback optimization module is used to issue the main scheduling instruction and the individual compensation instruction and collect actual response data, and update the matching degree weight in the multidimensional matching space and the individual difference envelope according to the actual response data.

8. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 6.

9. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 6.

Images