Coordinated optimization method and system for micro-grid source-grid-load-storage

Abstract

This invention provides a collaborative optimization method and system for microgrid generation, grid, load, and storage. The method includes clustering uncertain scenarios based on Mahalanobis distance between scenario parameters, calculating scenario weights based on clustering features, and selecting sub-problems to be solved through roulette wheel sampling; generating new tangent planes by solving sub-problems in parallel, merging new tangent planes within the same cluster into consensus tangent planes weighted by constraint violation depth and adding them to the main problem, freezing and re-verifying low-activity tangent planes based on activity scores; reshaping the directional adaptive trust domain according to the normal vector distribution of the active tangent planes, shrinking the boundary in the direction of dense normal vectors and widening the boundary in the direction of sparse normal vectors, and iteratively solving the main problem within this trust domain until the objective function deviation meets the convergence threshold, and outputting a collaborative optimization scheduling scheme for microgrid generation, grid, load, and storage.

Description

Technical Field

[0001] This application belongs to the field of microgrids, and in particular relates to a collaborative optimization method and system for source-grid-load-storage in microgrids. Background Technology

[0002] Microgrids, as an important energy system integrating renewable energy, distributed energy storage, and diverse loads, play a crucial role in improving energy efficiency, enhancing power supply reliability, and promoting the consumption of distributed energy. With the increasing proportion of distributed power sources such as wind and solar power, their output randomness, intermittency, and volatility are becoming more pronounced. Simultaneously, the integration of flexible loads, electric vehicles, and adjustable loads introduces more complex uncertainties into microgrids at the source, grid, load, and energy storage sides. To ensure the safe, economical, and stable operation of microgrids, source-grid-load-storage coordinated optimization has become a critical technical aspect of microgrid dispatch and control. In practical engineering modeling, scenario-based methods are typically used to describe uncertainties such as renewable energy output, load demand, and electricity prices. To cover varied operating conditions, optimization models often need to introduce numerous uncertain scenarios, thus transforming the microgrid source-grid-load-storage coordinated optimization problem into a mathematical programming problem with a large number of decision variables and constraints. As the number of scenarios increases, the model size expands rapidly, significantly increasing solution complexity and leading to problems such as excessive computation time, high memory consumption, and difficulty in meeting real-time dispatch requirements.

[0003] For the aforementioned large-scale optimization problems, existing technologies typically employ mathematical decomposition algorithms, breaking down the original problem into a main problem and multiple scenario sub-problems, and gradually approximating the optimization result through iterative interactions between the main problem and sub-problems. However, under massive scenario conditions, conventional decomposition algorithms usually need to solve a large number of sub-problems in each iteration, resulting in long iteration times per round. Simultaneously, the tangent plane constraints fed back to the main problem by different sub-problems often exhibit strong similarity and redundancy, easily causing the main problem constraint size to continuously expand, thereby reducing solution efficiency and increasing storage pressure. Furthermore, traditional methods lack adaptive management mechanisms for the main problem search space, easily leading to problems such as unreasonable search step sizes, iterative oscillations, or slow convergence speeds in high-dimensional decision spaces. Therefore, existing microgrid source-grid-load-storage collaborative optimization technologies still struggle to achieve a good balance between computational overhead, constraint size control, and iterative stability, urgently requiring a collaborative optimization method suitable for massive uncertainty scenarios. Summary of the Invention

[0004] To address the problems of high overhead in solving subproblems, redundant tangent plane constraints, expansion of the main problem size, and insufficient stability of iterative search in existing technologies under scenarios with massive uncertainty, this invention proposes a collaborative optimization method and system for microgrid sources, grids, loads, and storage.

[0005] In a first aspect, this invention proposes a collaborative optimization method for microgrid sources, grid, load, and storage, comprising the following steps: The microgrid source-grid-load-storage coordinated optimization problem is decomposed into a main problem and several parallel sub-problems corresponding to uncertain scenarios; all uncertain scenarios are clustered based on the Mahalanobis distance between parameters of uncertain scenarios; the importance weight of each uncertain scenario is calculated according to the cluster characteristics, and the parallel sub-problems to be solved are selected by roulette wheel sampling according to the importance weight; The selected parallel subproblems are solved in parallel to generate new tangent planes; the new tangent planes from the same cluster are weighted and fused according to the constraint violation depth relative to the solution of the decision variable of the current main problem to generate consensus tangent planes representing the clusters and add them to the main problem; an initial activity score is set for each consensus tangent plane added to the main problem, the activity score of the consensus tangent plane is updated after the main problem is solved, and the corresponding consensus tangent plane is temporarily frozen when the activity score is lower than the inactivation threshold; By utilizing the normal vector distribution characteristics of all active tangent planes with activity scores higher than the inactivation threshold, the trust domain boundary of the main problem decision variables is shrunk according to a first preset ratio in the principal component direction where the normal vectors are dense, and the trust domain boundary of the main problem decision variables is widened according to a second preset ratio in the secondary component direction where the normal vectors are sparse, thus reshaping the directional adaptive trust domain. The main problem is solved within the directional adaptive trust domain until the objective function deviation between adjacent iterations is less than the convergence threshold, and the final optimized solution is output.

[0006] On the other hand, this invention also proposes a collaborative optimization system for microgrid generation, grid, load, and storage, comprising the following modules: The computation module is used to decompose the microgrid source-grid-load-storage coordinated optimization problem into a main problem and several parallel sub-problems corresponding to uncertain scenarios; cluster all uncertain scenarios based on the Mahalanobis distance between parameters of uncertain scenarios; calculate the importance weight of each uncertain scenario according to the clustering characteristics, and select the parallel sub-problems to be solved by roulette wheel sampling according to the importance weight; The generation module is used to generate new tangent planes by solving selected parallel subproblems in parallel; new tangent planes from the same cluster are weighted and fused according to the constraint violation depth relative to the solution of the decision variables of the current main problem to generate consensus tangent planes representing the clusters and add them to the main problem; an initial activity score is set for each consensus tangent plane added to the main problem, the activity score of the consensus tangent plane is updated after the main problem is solved, and the corresponding consensus tangent plane is temporarily frozen when the activity score is lower than the inactivation threshold; The output module utilizes the normal vector distribution characteristics of all active tangent planes with activity scores higher than the inactivation threshold to shrink the trust domain boundary of the main problem decision variables according to a first preset ratio in the principal component direction with dense normal vectors, and to widen the trust domain boundary of the main problem decision variables according to a second preset ratio in the secondary component direction with sparse normal vectors, thus reshaping the direction-adaptive trust domain; the main problem is solved within the direction-adaptive trust domain until the objective function deviation between adjacent iterations is less than the convergence threshold, and the final optimized solution is output.

[0007] This invention decomposes the microgrid source-grid-load-storage collaborative optimization problem into a main problem and scenario-based sub-problems. It then uses Mahalanobis distance to cluster uncertain scenarios and combines roulette wheel sampling to select representative scenarios for solution, reducing the computational cost of each sub-problem and improving solution efficiency. By evaluating the constraint violation depth of the tangent plane, newly emerging tangent planes of the same type are weighted and merged into a consensus tangent plane. An activity score management mechanism is used to freeze and re-verify low-activity tangent planes, reducing the impact of continuous accumulation of redundant constraints on solving the main problem. Furthermore, based on the normal vector distribution characteristics of the active tangent plane, the directional adaptive trust domain is reshaped. The search range is narrowed in directions with dense constraint information and the search boundary is widened in directions with less constraint information, thereby reducing the risk of iterative oscillation and improving the stability and computational efficiency of the collaborative optimization process. Attached Figure Description

[0008] Figure 1 A flowchart of the first embodiment; Figure 2 This is a schematic diagram of the principal component two-dimensional clustering results for an uncertain scenario; Figure 3 This is a schematic diagram of the activity fraction decay characteristic curve during the consensus tangent plane iteration process; Figure 4 A bar chart comparing the maximum memory usage of the model under different algorithm configurations. Detailed Implementation

[0009] To enable those skilled in the art to better understand the technical solutions in this application, the technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of this application. Given the complexity of practical applications, the data types and numerical values ​​involved in the various embodiments in this disclosure are only for illustrative purposes and are not necessarily derived from actual application data. Moreover, different embodiments may use different data, which will not affect the implementation of the present invention. Each embodiment and the specific features contained therein can be freely combined, replaced, or nested after adjusting the operating parameters, adapting the data format, and other conventional modifications in the art. The superior expression of the present invention cannot list all applicable situations. Those skilled in the art can extend its application to other situations based on the teachings of the present invention without departing from the spirit and scope of the present invention.

[0010] In the first embodiment, the present invention proposes a collaborative optimization method for microgrid generation, grid, load, and storage, such as... Figure 1 As shown, it includes the following steps: S1 decomposes the microgrid source-grid-load-storage coordinated optimization problem into a main problem and several parallel sub-problems corresponding to uncertain scenarios; clusters all uncertain scenarios based on the Mahalanobis distance between parameters of uncertain scenarios; calculates the importance weight of each uncertain scenario according to the clustering characteristics, and selects the parallel sub-problems to be solved by roulette wheel sampling according to the importance weight.

[0011] A mixed-integer linear programming mathematical model for microgrids is established. The Benders decomposition algorithm is used to decompose the microgrid source-grid-load-storage coordinated optimization problem into a main problem responsible for optimizing the start-up and shutdown status of generating units and the configuration of equipment capacity, and several sub-problems responsible for optimizing the continuous variable scheduling of wind and solar power output and load fluctuations under given uncertain scenarios.

[0012] An uncertainty scenario parameter matrix is ​​constructed by extracting wind power output sequences, photovoltaic power output sequences, load demand sequences, and real-time market electricity price volatility from historical operational data. The covariance matrix of this parameter matrix is ​​calculated, the inverse matrix is ​​obtained, and the Mahalanobis distance between parameters of each uncertainty scenario is calculated. Based on these Mahalanobis distances, the DBSCAN density-based clustering algorithm is used to cluster all uncertainty scenarios, generating multiple clusters. The number of scenario samples and the core distance density feature within each cluster are statistically analyzed. The product of the number of samples and the density is used as the baseline weight for that cluster, and the importance weights of each uncertainty scenario within the cluster are assigned based on the Mahalanobis distance from the scenario to its cluster center. A roulette wheel sampling method is then used, with the importance weights of the uncertainty scenarios as probability distribution parameters, to select parallel subproblems to be solved in the next stage of the solution process.

[0013] In an optional embodiment, the clustering of all uncertainty scenarios based on the Mahalanobis distance between the uncertainty scenario parameters includes: Construct a parameter vector for an uncertain scenario that includes random disturbances in the microgrid; Calculate the covariance matrix of the overall uncertainty scenario parameter vector, and obtain the inverse matrix of the covariance matrix as the accuracy matrix; The distance between any two sets of uncertain scenario parameter vectors is calculated using the Mahalanobis distance calculation formula, with the accuracy matrix used as the weight matrix in the calculation formula. A density-based spatial clustering algorithm is applied to divide all uncertain scenarios into multiple non-overlapping clusters based on all calculated Mahalanobis distance values.

[0014] The process of constructing parameter vectors for uncertain scenarios requires collecting historical and predicted data on multidimensional random disturbances in the microgrid system, covering wind power output prediction errors, photovoltaic power output prediction errors, load prediction fluctuations, and real-time market electricity price volatility. Let the total number of microgrid scenario samples be N=1000, and the parameter vector for each scenario be denoted as... The feature dimension d=4. The overall covariance matrix is ​​calculated using a formula, where the diagonal elements represent the variance of a single variable, and the off-diagonal elements represent the covariance between random variables such as wind and solar power output and load. After obtaining the covariance matrix, its inverse is used as the precision matrix. The distance between any two sets of scene vectors is calculated using the Mahalanobis distance formula. During the calculation, the precision matrix is ​​used as the weight matrix to reduce the interference from the multidimensional dimensional differences of the microgrid and the correlation between wind and solar loads. After obtaining a distance matrix of size 1000×1000, density-based spatial clustering algorithms such as DBSCAN are used for partitioning. The neighborhood search radius parameter of the DBSCAN algorithm is set. The range is set to 0.8 to 1.5, and the example uses... =1.2 and the minimum number of sample points (MinPts) required to become a core object, ranging from 5 to 15; in this example, MinPts = 10. Using Mahalanobis distance as the spatial metric, high-density perturbation scene clusters are identified and grouped into one cluster. Discrete low-density anomalous scenes are marked as noise. For uncertain scenes marked as noise, they can be retained separately as low-density perturbation classes, or merged into the nearest cluster according to Mahalanobis distance, and participate in subsequent importance weight calculations to ensure that each uncertain scene has a cluster affiliation suitable for weight calculation. Through this mechanism, the aforementioned 1000 uncertain scenes are divided into, for example, K = 8 non-overlapping cluster sets, such as... Figure 2 As shown, the spatial clustering distribution and outlier sample distribution characteristics of various disturbance scenarios are presented.

[0015] In an optional embodiment, the step of calculating the importance weights of each uncertainty scenario based on the clustering features, and selecting the parallel sub-problems to be solved by roulette wheel sampling according to the importance weights, includes: The total number of uncertain scenario samples within each cluster is used as the cluster size. Calculate the arithmetic mean of the Mahalanobis distances from each uncertainty scenario to the cluster center within the same cluster, and take the reciprocal of the mean as the cluster density value; The cluster size is multiplied by the cluster density to obtain the basic cluster weight. The basic cluster weight is then multiplied by the reciprocal of the Mahalanobis distance from the target uncertainty scenario to its cluster center to obtain the importance weight of the parallel subproblem corresponding to the target uncertainty scenario. Calculate the proportion of each importance weight in the total weight of all parallel subproblems, set the interval width based on this proportion, generate random numbers using the roulette wheel sampling algorithm based on the interval width, and select the corresponding parallel subproblems for parallel solution.

[0016] For each scene cluster generated by clustering, the statistics module counts the number of samples within it and extracts the cluster size. For example, for the k-th cluster, the total number of scenes contained in that cluster is calculated to be... =150. Determine the multidimensional centroid of the cluster, calculate the Mahalanobis distance from each uncertainty scenario within the cluster to the centroid, and take the arithmetic mean of these 150 distance values. Assuming the calculated average Mahalanobis distance is 2.5, take the reciprocal of this Mahalanobis distance to obtain the cluster density value of 0.4. Multiply the cluster size value of 150 by the cluster density value of 0.4 to obtain the basic cluster weight of 60. This operation gives higher basic sampling weights to typical scenario classes with large scale and dense data distribution, i.e., high probability of occurrence. For the target uncertainty scenario within the cluster, calculate the reciprocal of the Mahalanobis distance from itself to the cluster centroid as the center proximity factor. To prevent the denominator from being 0, add 10 to the denominator. -6Constant compensation. Assuming the distance between the target scene and the centroid is 1.5, the reciprocal of which is approximately 0.667, multiplying the base weight of 60 by 0.667, the importance weight of the parallel subproblem corresponding to this scene is calculated to be 40.02.

[0017] After calculating the weights of all 1000 parallel subproblems across the entire network, the total weight base is obtained by summing the weights. The weight of each subproblem divided by the total weight gives the selection percentage for that subproblem. Based on this percentage, the [0,1] interval is divided into 1000 continuous probability intervals of corresponding width. A uniformly distributed pseudo-random number generator is used to repeatedly generate a preset number of random floating-point numbers within the [0,1] range. By determining the specific interval position where each random number falls, the corresponding parallel subproblem is selected, forming a set of subproblems to be solved in parallel. These subproblems are then submitted to a multi-threaded pool for distributed solving. This ensures that the microgrid prioritizes high-frequency typical disturbances during scheduling, while retaining the probability of other representative disturbance scenarios being sampled and solved.

[0018] S2, Parallel solution of selected parallel subproblems generates new tangent planes; New tangent planes originating from the same cluster are weighted and fused according to the constraint violation depth relative to the solution of the decision variables of the current main problem to generate consensus tangent planes representing the clusters and add them to the main problem; An initial activity score is set for each consensus tangent plane added to the main problem, the activity score of the consensus tangent plane is updated after the main problem is solved, and the corresponding consensus tangent plane is temporarily frozen when the activity score is lower than the inactivation threshold.

[0019] The selected parallel subproblems are assigned to different processes or computing nodes, and the Gurobi mathematical programming solver is invoked for parallel solving. After solving, the dual results corresponding to the coupling constraints of the main problem in each subproblem are extracted; wherein, when the subproblem is feasible and obtains the optimal solution, the dual result includes the Lagrange multipliers of the corresponding coupling constraints; when the subproblem is infeasible, the dual result includes the dual ray used to construct the feasible tangent plane. Based on the dual results, the coefficient vector and constant term of the Benders optimal tangent plane or feasible tangent plane are constructed to generate the new tangent plane.

[0020] Substituting the solution of the current master problem's decision variables into the equations of each newly formed tangent plane, the corresponding inequality constraint violations are calculated. The positive values ​​of these violations are used as the constraint violation depths of each newly formed tangent plane. The constraint violation depths of each newly formed tangent plane within the same cluster are normalized to obtain the corresponding fusion weights. The coefficient vectors and constant terms of multiple newly formed tangent planes within the same cluster are weighted and summed using the fusion weights to generate a consensus tangent plane equation representing that cluster. This equation is then added to the master problem model as a linear inequality constraint.

[0021] Each newly added consensus tangent plane to the main problem model is assigned an initial activity score. After the main problem is solved again using the Gurobi solver, the state of the consensus tangent plane is determined based on the dimensionless constraint residuals of the current main problem decision variable solution relative to each consensus tangent plane. When the constraint residual is greater than a preset violation threshold, the consensus tangent plane is determined to be a violated constraint, and its activity score is increased or maintained. When the absolute value of the constraint residual meets a preset fit condition, the activity score of the corresponding consensus tangent plane is increased. When the constraint residual is less than a preset internal threshold, a decay factor is calculated based on the dimensionless orthogonal distance, and the activity score of the corresponding consensus tangent plane is decreased. After each activity score update, a judgment is made. When the activity score of a consensus tangent plane is lower than the inactivation threshold, and the consensus tangent plane is not in a violated state, the constraint condition corresponding to the consensus tangent plane is temporarily removed from the main problem model in the next round of main problem iteration, achieving a freeze process.

[0022] In an optional embodiment, the step of weighting and fusing newly generated tangent planes originating from the same cluster according to the constraint violation depth relative to the solution of the current master problem decision variables to generate a consensus tangent plane representing the cluster and adding it to the master problem includes: Extract the coefficient vector and constant term of all newly generated tangent planes belonging to the same cluster in the current iteration; Substitute the solution of the current main problem decision variables into the equations of each newly generated tangent plane, calculate the constraint violation amount, and set the positive part of the constraint violation amount as the constraint violation depth of each newly generated tangent plane; The constraint violation depth of each newly formed tangent plane within the same cluster is normalized to obtain the proportional weight corresponding to each newly formed tangent plane. The coefficient vectors and constant terms of each newly generated tangent plane within the same cluster are linearly weighted and summed using the aforementioned proportional weights to generate a consensus tangent plane representing the cluster. This consensus tangent plane is then added to the main problem as an inequality constraint.

[0023] In the current collaborative optimization iteration process, after parallel solving of a specific cluster, a set of generated Benders tangent planes, such as M=5, is extracted. The mathematical expression for each new tangent plane is as follows: The system extracts the coefficient vector corresponding to the decision dimension. It includes generator output coefficient, energy storage charge / discharge state coefficient, and bias constant term. The scheduling algorithm will use the real solutions to the microgrid decision variables already solved in the current main problem. Substitute the values ​​into the equations of the five newly generated tangent planes one by one, and calculate the difference between the actual calculated value of the equation and 0, which is the constraint violation amount. When the calculation result is greater than 0, it indicates that the current solution is out of bounds, and this positive violation amount is set as the constraint violation depth; if the result is less than or equal to 0, the constraint violation depth is set to 0. Assuming that the calculated constraint violation amounts of the five tangent planes are 1.2, 0.8, -0.5, 2.0, and 0.1, the effective constraint violation depth sequence after filtering out negative values ​​is 1.2, 0.8, 0, 2.0, and 0.1. After obtaining all positive constraint violation depths, normalize these positive constraint violation depths and calculate their respective proportional weights. The effective constraint violation depths are summed, resulting in a total of 1.2 + 0.8 + 2.0 + 0.1 = 4.1. The proportional weights of each tangent plane are 1.2 / 4.1 ≈ 0.293, 0.8 / 4.1 ≈ 0.195, 2.0 / 4.1 ≈ 0.488, and 0.1 / 4.1 ≈ 0.024, respectively. Tangent planes that do not exceed the limits are assigned a weight of 0. Using these proportional weights, a linear weighted summation is performed on the coefficient vector matrix and constant term vector of these five tangent planes to generate a new consensus tangent plane. The consensus tangent plane, as an inequality constraint representing the entire cluster, is added to the constraint matrix of the main problem model in the next round. This operation preserves the direction information of the newly generated tangent planes that have a strong constraint effect on the solution of the current main problem, and reduces the scale of linear constraints of the main problem in subsequent iterations by fusion, thereby reducing the solver's memory usage.

[0024] In an optional embodiment, the step of setting an initial activity score for each consensus tangent added to the main problem, updating the activity score of the consensus tangent after the main problem is solved, and temporarily freezing the corresponding consensus tangent when the activity score falls below an inactivation threshold includes: Based on the scaling benchmark of the decision variables of the main problem, the coefficient vectors and constant terms of the current decision variable solution of the main problem and each consensus tangent plane are uniformly dimensionless, and the dimensionless constraint residuals of the current decision variable solution of the main problem relative to each consensus tangent plane are calculated. If the constraint residual is greater than a preset violation threshold, the consensus tangent plane is determined to be a violated constraint, a preset violation compensation term is added to the activity score of the consensus tangent plane, and the consensus tangent plane is kept participating in the construction of the main problem constraint. If the absolute value of the constraint residual satisfies the preset fitting condition, the consensus tangent plane is determined to be a valid constraint, and a preset positive constant compensation term is added to the activity fraction of the consensus tangent plane. If the constraint residual is less than a preset internal threshold, the consensus tangent plane is determined to be an inactive constraint. The decay factor is calculated with the natural constant as the base and the negative value of the product of the absolute value of the constraint residual and the preset dimensionless decay coefficient as the exponent. The current activity score is then updated by multiplying the decay factor. When the activity score of the consensus tangent is lower than the inactivation threshold and the consensus tangent is not in a violated state, the consensus tangent will be temporarily removed from the constraint set and frozen in the subsequent iteration of the main problem solution, and will not participate in the construction of the main problem constraints.

[0025] After the mixed-integer programming model of the main problem is solved, the constraint management module evaluates all added consensus tangents. When a tangent is first generated and added to the model, it is assigned an initial activity score, for example, by setting an initial value of [value missing]. =10.0. The latest decision solution is obtained in each round of the main problem. Subsequently, based on the nominal capacity, allowable adjustment interval, or historical standard deviation of the main problem decision variables, a scaling benchmark is determined. Then, the coefficient vector and constant term of the current main problem decision variable solution, consensus tangent plane, and other related terms are uniformly dimensionless. After this uniform dimensionless process, the dimensionless decision variable solution is denoted as... Let the dimensionless coefficient vector be denoted as The dimensionless constant term is denoted as .Will Substitute by and The expression for the dimensionless consensus tangent plane is constructed, and the calculation is performed. and The sum of the sums, and the sum divided by the sum. The L2 norm is used to obtain the constraint residual r of the current main problem decision variable solution with respect to the consensus tangent plane. The absolute value of r is taken as the dimensionless orthogonal distance L, and the larger value between r and 0 is taken as the dimensionless constraint violation V.

[0026] Set the numerical tolerance parameter to 10. -5 When V is greater than a preset violation threshold, it indicates that the current solution to the main problem violates the consensus tangent plane. A preset violation compensation term is added to the activity score of this tangent plane, and the consensus tangent plane continues to participate in the construction of the main problem constraints. When |r|≤10 -5 When the current solution is determined to be on the boundary of the tangent plane, the constraint is considered active. A preset positive constant compensation term, for example, 2.0, is added to the activity score of the tangent plane, updating the current activity score to 12.0. If r is less than a preset internal threshold, it indicates that the current principal problem solution is within the feasible region defined by the tangent plane but far from it, and is therefore considered an inactive constraint. An exponential decay factor γ = exp(-α·L) with the natural constant e as its base is constructed, where α is a preset dimensionless decay coefficient, for example, set to 2.0. If the dimensionless orthogonal distance L is 0.5, then the decay factor γ ≈ 0.367, and the current activity score is updated multiplicatively to 10.0 × 0.367 = 3.67. When the constraint residual does not meet any of the above update conditions, the current activity score of the consensus tangent plane remains unchanged.

[0027] After updating all activity scores, a threshold check is performed. A fixed inactivation threshold is set, for example, 1.0. When the score of a consensus tangent is below 1.0, and the consensus tangent is not in a violated state, it indicates that the constraint effect of the consensus tangent on the solution of the current main problem is weak, and the state of the consensus tangent is changed to "frozen". When entering the next outer iteration to construct the underlying mathematical model of the main problem, all tangent constraints in the frozen state will be skipped.

[0028] For consensus tangent planes in a frozen state, their coefficient vectors, constant terms, and cluster identifiers are retained in subsequent iterations. The solution to the current master problem's decision variables is substituted into the frozen tangent plane for dimensionless constraint violation verification according to a preset verification period. When the dimensionless constraint violation of the frozen tangent plane exceeds a preset reactivation threshold, or when its cluster generates a similar new tangent plane in the current iteration, the frozen tangent plane is re-added to the master problem's constraint set, and its activity score is restored or reassigned. This mechanism reduces the continuous occupation of the master problem solution by low-activity constraints, helping to maintain the sparsity of the constraint matrix. The activity score decay characteristic curve during consensus tangent plane iteration is shown below. Figure 3 As shown.

[0029] S3 utilizes the normal vector distribution characteristics of all active tangent planes with activity scores higher than the inactivation threshold. In the principal component direction with dense normal vectors, the trust domain boundary of the main problem decision variables is shrunk according to a first preset ratio. In the secondary component direction with sparse normal vectors, the trust domain boundary of the main problem decision variables is widened according to a second preset ratio, thus reshaping the direction-adaptive trust domain. The main problem is solved within the direction-adaptive trust domain until the objective function deviation between adjacent iterations is less than the convergence threshold, and the final optimized solution is output.

[0030] Iterate through all tangent planes in the current main problem, selecting active tangent planes that are not frozen and have an activity score higher than the inactivation threshold. Extract the normal vectors of these active tangent planes to form a normal vector feature matrix. Call the PCA principal component analysis algorithm to perform dimensionality reduction and direction feature extraction calculations on the normal vector feature matrix, obtaining eigenvalues ​​representing the energy distribution of normal vectors and their corresponding direction vectors. Determine the top few direction vectors whose cumulative eigenvalue ratio meets a preset threshold as the principal component directions with dense normal vectors, and determine the remaining direction vectors as the secondary component directions with sparse normal vectors.

[0031] In the principal component direction, the original trust domain search step size of the main problem decision variables is reduced by a first preset ratio to shrink the boundary constraints in that direction. In the secondary component direction, the original trust domain search step size of the main problem decision variables is relaxed by a second preset ratio to widen the boundary constraints in that direction. Based on the adjusted step sizes in each direction, the upper and lower bound constraints of the directional adaptive trust domain are reshaped. The reshaped upper and lower bound constraints of the directional adaptive trust domain are updated in the main problem model and solved using the Gurobi solver to obtain the system operating cost objective function value for this iteration. The absolute difference between the objective function value of this iteration and the objective function value of the previous iteration is calculated. It is determined whether the absolute difference is less than a preset cost deviation convergence threshold. If the difference is greater than or equal to the convergence threshold, the result of the decision variables of the current main problem is passed down to the subproblems and the next iteration begins. If the difference is less than the convergence threshold, the current iteration process is determined to meet the convergence stopping condition and the iteration stops. The `to_csv` function of the Pandas library is used to export the unit output plan matrix and energy storage charging and discharging power sequence, etc., obtained from the main problem and subproblems, to a comma-separated value file, outputting the final optimization scheme.

[0032] In an optional embodiment, the step of utilizing the normal vector distribution characteristics of all active tangent planes with activity scores higher than the inactivation threshold, shrinking the trust domain boundary of the main problem decision variable according to a first preset ratio in the principal component direction with dense normal vectors, and widening the trust domain boundary of the main problem decision variable according to a second preset ratio in the secondary component direction with sparse normal vectors, and reshaping the direction-adaptive trust domain, includes: Extract the normal vectors of all active tangent planes with activity fractions higher than the inactivation threshold, and then perform dimensionless standardization on the normal vectors according to the dimensions of the main problem decision variables before combining them to form a normal vector matrix; Perform eigenvalue decomposition on the covariance matrix corresponding to the normal vector matrix to extract eigenvalues ​​and their corresponding direction vectors; The first few directional vectors whose cumulative eigenvalue ratio is greater than a preset threshold are designated as principal component directions with dense normal vectors, and the trust domain boundary of the main problem decision variable in the principal component direction is shrunk according to a first preset ratio. The remaining direction vectors are labeled as secondary component directions with sparse normal vectors. The trust domain boundaries of the main problem decision variables in the secondary component directions are widened according to a second preset ratio to reshape the direction adaptive trust domain.

[0033] To implement step-size control in the high-dimensional space of a microgrid, consensus tangent planes with activity scores higher than the inactivation threshold (e.g., greater than 1.0) are extracted from the constraint pool, and their respective d-dimensional normal vectors are extracted. Since microgrid decision variables simultaneously include active power output (MW), reactive power support (MVar), and energy storage capacity (MWh), the allowable upper and lower limits of adjustment for each decision variable are extracted as benchmark references. Multiplicative scaling is then performed on the corresponding dimensional components of all normal vectors to achieve dimensionless standardization.

[0034] The processed K effective normal vectors are arranged row by row and assembled into a normal vector matrix of size K×d. The product of the transpose of the normal vector matrix and itself is calculated to obtain a d×d direction covariance matrix. Eigenvalue decomposition is then performed on the direction covariance matrix, yielding d eigenvalues ​​and their corresponding orthogonal eigenvector groups in descending order of eigenvalues. In the spatial projection direction classification stage, the sum of all eigenvalues ​​is calculated, and a preset proportion threshold of 85% is set. The top-ranked eigenvalues ​​are accumulated one by one. When the cumulative proportion first exceeds 85%, the subspace spanned by the corresponding top p eigenvectors is established as the principal component direction with dense normal vectors. The boundary line of the principal component direction, constrained by the main problem's trust domain (i.e., the maximum allowable movement step of the decision variable between two iterations), is shrunk according to a first preset proportion, for example, a shrinkage coefficient of 0.5. The remaining dp direction vectors are determined as secondary component directions with sparse normal vectors. This secondary component subspace is then expanded outward according to a second preset proportion, for example, a relaxation coefficient of 1.25. Through the aforementioned coordinate system transformation and boundary scaling mechanism, a direction-adaptive trust domain polyhedron is reshaped.

[0035] In an optional embodiment, solving the master problem within the orientation-adaptive trust domain until the objective function deviation between adjacent iterations is less than a convergence threshold, and outputting the final optimized solution, includes: At the end of each iteration, the objective function value of the current main problem is used as the current overall estimated cost; The absolute difference between the current total estimated cost of the current iteration and the total estimated cost of the previous iteration is calculated as the deviation of the objective function; When the deviation of the objective function is less than the set convergence threshold for a preset number of consecutive times, the algorithm is determined to have converged, and a stop command is triggered to terminate the iteration. The optimal decision variables at the convergence of the current main problem are extracted as the final control parameters of the microgrid source-grid-load-storage collaborative optimization scheme, and the output is the final optimization scheme.

[0036] After solving the mixed-integer master problem model with added constraints on the orientation-adaptive trust domain boundary and consensus tangent plane, the returned objective function value is obtained at the end of the t-th iteration. This real scalar represents the current estimated total cost of the microgrid, including economic factors such as generator fuel costs, depreciation costs for energy storage equipment, costs of purchasing electricity from outside the grid, and penalties for default. For example, the estimated total cost obtained in the t-th iteration is... =23,546.8 yuan, while the estimated cost obtained in the previous (t-1)th iteration was =23,551.2 yuan. The absolute difference between the two is calculated as the evaluation index, i.e., the objective function deviation ΔF = |23,546.8 - 23,551.2| = 4.4 yuan. A fixed convergence threshold is set, for example, an absolute threshold of 5.0 yuan, and a preset number of consecutive judgments is set, for example, 3 times. In the collaborative solution process, the main control module maintains a sliding verification window to monitor the historical sequence of the objective function deviation in real time. If the absolute difference calculated in the t-th, t+1-th, and t+2-th consecutive times all fall below 5.0 yuan, it indicates that the optimization process meets the stopping condition, and the algorithm is determined to have converged. A stop command is triggered to terminate the parallel iteration loop, and all decision variables in the converged state are extracted from the memory of the main control node. These variables include the active power output of each diesel generator set at each time period, such as the output of Unit 3 being 500kW, the charging and discharging power and SOC curve of the energy storage battery at each time period, and the timing instructions for the activation and deactivation of flexible loads. These variables will be packaged and converted into executable dispatch instruction messages, which will be sent to the energy management system at the bottom of the microgrid as the final optimization scheme for implementation.

[0037] The experimental platform utilizes a computing device equipped with 32GB of memory and employs a commercial mathematical programming solver within the algorithm development environment to conduct microgrid collaborative optimization verification. Based on historical operation and forecast data of real microgrids, the total number of scenario samples with uncertain wind and solar loads was set to 1000, the feature dimension of the scenario parameter vector was four-dimensional, the neighborhood radius of the density clustering algorithm was set to 1.2, the minimum number of samples was set to 10, and the tangent plane inactivation threshold was set to 1.0. The convergence stopping condition was set to an absolute deviation of the objective function between adjacent iterations of less than 5.0 yuan for three consecutive iterations. The comparative evaluation indicators of the ablation experiment covered four core dimensions: daily comprehensive estimated cost, total number of iterations, average time per iteration, and maximum memory usage.

[0038] The baseline traditional decomposition algorithm yields a daily estimated cost of 24,510.3 yuan, requiring a total of 82 iterations, with an average iteration time of 15.4 seconds and a maximum memory usage of 4.2 GB. After eliminating the Mahalanobis distance-based clustering and weighted sampling mechanism, the daily estimated cost increases to 24,035.6 yuan compared to the complete solution because the algorithm cannot select representative perturbation scenarios based on the scene clustering structure. Eliminating the consensus tangent plane fusion and liveness freezing mechanism, while retaining all newly generated constraints, leads to a geometric expansion of the main problem constraint size, increasing the iteration time to 11.8 seconds and the maximum memory usage to 3.5 GB. Eliminating the direction-adaptive trust domain mechanism causes oscillations in the solution process due to the fixed optimization step size, increasing the total number of iterations to 58. In contrast, applying the complete solution yields a lower daily estimated cost of 23,546.8 yuan, reducing the total number of iterations to 35, with an average iteration time of only 4.6 seconds and a maximum memory usage of 1.3 GB. Maximum memory usage of the model under different algorithm configurations, for example Figure 4 As shown.

[0039] Clustering and center proximity factor sampling mechanisms reallocate computational resources, prioritizing high-frequency typical disturbances while retaining the probability of other representative disturbance scenarios being sampled and solved. This reduces potential over-limit penalty costs for system operation and improves scheduling economy. Tangential plane fusion based on constraint violation depth and activity freezing operations based on dimensionless constraint residuals and dimensionless orthogonal distance identify and temporarily remove redundant low-activity constraints within the feasible region, mitigating increased memory consumption and computational latency as iteration deepens. The direction-adaptive trust domain reshaping mechanism relies on feature decomposition to identify the spatial distribution patterns of constraint boundaries. It shrinks the search step size in directions with dense and sensitive constraint information to reduce over-limit risks, and appropriately expands the search range in directions with less constraint information to improve exploration efficiency. This enhances the optimization stability and computational efficiency of the microgrid source-grid-load-storage collaborative optimization scheme.

[0040] In the second embodiment, the present invention also proposes a collaborative optimization system for microgrid generation, grid, load, and storage, comprising the following modules: The computation module is used to decompose the microgrid source-grid-load-storage coordinated optimization problem into a main problem and several parallel sub-problems corresponding to uncertain scenarios; cluster all uncertain scenarios based on the Mahalanobis distance between parameters of uncertain scenarios; calculate the importance weight of each uncertain scenario according to the clustering characteristics, and select the parallel sub-problems to be solved by roulette wheel sampling according to the importance weight; The generation module is used to generate new tangent planes by solving selected parallel subproblems in parallel; new tangent planes from the same cluster are weighted and fused according to the constraint violation depth relative to the solution of the decision variables of the current main problem to generate consensus tangent planes representing the clusters and add them to the main problem; an initial activity score is set for each consensus tangent plane added to the main problem, the activity score of the consensus tangent plane is updated after the main problem is solved, and the corresponding consensus tangent plane is temporarily frozen when the activity score is lower than the inactivation threshold; The output module utilizes the normal vector distribution characteristics of all active tangent planes with activity scores higher than the inactivation threshold to shrink the trust domain boundary of the main problem decision variables according to a first preset ratio in the principal component direction with dense normal vectors, and to widen the trust domain boundary of the main problem decision variables according to a second preset ratio in the secondary component direction with sparse normal vectors, thus reshaping the direction-adaptive trust domain; the main problem is solved within the direction-adaptive trust domain until the objective function deviation between adjacent iterations is less than the convergence threshold, and the final optimized solution is output.

[0041] The above detailed embodiments further illustrate the purpose, technical solution, and beneficial effects of the embodiments of this application. It should be understood that the above are merely specific implementations of the embodiments of this application and are not intended to limit the scope of protection of the embodiments of this application. Any modifications, equivalent substitutions, improvements, etc., made based on the technical solutions of the embodiments of this application should be included within the scope of protection of the embodiments of this application. All data interaction, information collection, and analysis steps involved in this specification should be understood as being performed within a legal framework that fully complies with national laws, social ethics, and public interests.

Claims

1. A collaborative optimization method for power generation, grid, load, and storage in microgrids, characterized in that, Includes the following steps: The microgrid source-grid-load-storage coordinated optimization problem is decomposed into a main problem and several parallel sub-problems corresponding to uncertain scenarios; all uncertain scenarios are clustered based on the Mahalanobis distance between parameters of uncertain scenarios; the importance weight of each uncertain scenario is calculated according to the cluster characteristics, and the parallel sub-problems to be solved are selected by roulette wheel sampling according to the importance weight; The selected parallel subproblems are solved in parallel to generate new tangent planes; the new tangent planes from the same cluster are weighted and fused according to the constraint violation depth relative to the solution of the decision variable of the current main problem to generate consensus tangent planes representing the clusters and add them to the main problem; an initial activity score is set for each consensus tangent plane added to the main problem, the activity score of the consensus tangent plane is updated after the main problem is solved, and the corresponding consensus tangent plane is temporarily frozen when the activity score is lower than the inactivation threshold; By utilizing the normal vector distribution characteristics of all active tangent planes with activity scores higher than the inactivation threshold, the trust domain boundary of the main problem decision variables is shrunk in the principal component direction with dense normal vectors according to a first preset ratio, and the trust domain boundary of the main problem decision variables is widened in the secondary component direction with sparse normal vectors according to a second preset ratio, thus reshaping the direction-adaptive trust domain. Solve the master problem within the direction-adaptive trust domain until the objective function deviation between adjacent iterations is less than the convergence threshold, and output the final optimized solution.

2. The method according to claim 1, characterized in that, The clustering of all uncertain scenarios based on the Mahalanobis distance between parameters of the uncertain scenarios includes: Construct a parameter vector for an uncertain scenario that includes random disturbances in the microgrid; Calculate the covariance matrix of the overall uncertainty scenario parameter vector, and obtain the inverse matrix of the covariance matrix as the accuracy matrix; The distance between any two sets of uncertain scenario parameter vectors is calculated using the Mahalanobis distance calculation formula, with the accuracy matrix used as the weight matrix in the calculation formula. A density-based spatial clustering algorithm is applied to divide all uncertain scenarios into multiple non-overlapping clusters based on all calculated Mahalanobis distance values.

3. The method according to claim 1, characterized in that, The process of calculating the importance weights of each uncertainty scenario based on the characteristics of clustering, and then using roulette wheel sampling to select the parallel sub-problems to be solved based on these importance weights, includes: The total number of uncertain scenario samples within each cluster is used as the cluster size. Calculate the arithmetic mean of the Mahalanobis distances from each uncertainty scenario to the cluster center within the same cluster, and take the reciprocal of the mean as the cluster density value; The cluster size is multiplied by the cluster density to obtain the basic cluster weight. The basic cluster weight is then multiplied by the reciprocal of the Mahalanobis distance from the target uncertainty scenario to its cluster center to obtain the importance weight of the parallel subproblem corresponding to the target uncertainty scenario. Calculate the proportion of each importance weight in the total weight of all parallel subproblems, set the interval width based on this proportion, generate random numbers using the roulette wheel sampling algorithm based on the interval width, and select the corresponding parallel subproblems for parallel solution.

4. The method according to claim 1, characterized in that, The step of weighting and fusing newly generated tangent planes originating from the same cluster according to the constraint violation depth relative to the solution of the current master problem's decision variables to generate a consensus tangent plane representing the cluster and adding it to the master problem includes: Extract the coefficient vector and constant term of all newly generated tangent planes belonging to the same cluster in the current iteration; Substitute the solution of the current main problem decision variables into the equations of each newly generated tangent plane, calculate the constraint violation amount, and set the positive part of the constraint violation amount as the constraint violation depth of each newly generated tangent plane; The constraint violation depth of each newly formed tangent plane within the same cluster is normalized to obtain the proportional weight corresponding to each newly formed tangent plane. The coefficient vectors and constant terms of each newly generated tangent plane within the same cluster are linearly weighted and summed using the aforementioned proportional weights to generate a consensus tangent plane representing the cluster. This consensus tangent plane is then added to the main problem as an inequality constraint.

5. The method according to claim 1, characterized in that, The process of setting an initial activity score for each consensus tangent plane added to the main problem, updating the activity score of the consensus tangent plane after the main problem is solved, and temporarily freezing the corresponding consensus tangent plane when the activity score falls below the inactivation threshold includes: Based on the scaling benchmark of the decision variables of the main problem, the coefficient vectors and constant terms of the current decision variable solution of the main problem and each consensus tangent plane are uniformly dimensionless, and the dimensionless constraint residuals of the current decision variable solution of the main problem relative to each consensus tangent plane are calculated. If the constraint residual is greater than a preset violation threshold, the consensus tangent plane is determined to be a violated constraint, a preset violation compensation term is added to the activity score of the consensus tangent plane, and the consensus tangent plane is kept participating in the construction of the main problem constraint. If the absolute value of the constraint residual satisfies the preset fitting condition, the consensus tangent plane is determined to be a valid constraint, and a preset positive constant compensation term is added to the activity fraction of the consensus tangent plane. If the constraint residual is less than a preset internal threshold, the consensus tangent plane is determined to be an inactive constraint. The decay factor is calculated with the natural constant as the base and the negative value of the product of the absolute value of the constraint residual and the preset dimensionless decay coefficient as the exponent. The current activity score is then updated by multiplying the decay factor. When the activity score of the consensus tangent is lower than the inactivation threshold and the consensus tangent is not in a violated state, the consensus tangent will be temporarily removed from the constraint set and frozen in the subsequent iteration of the main problem solution, and will not participate in the construction of the main problem constraints.

6. The method according to claim 1, characterized in that, The method of utilizing the normal vector distribution characteristics of all active tangent planes with activity fractions higher than the inactivation threshold to shrink the trust domain boundary of the main problem decision variables according to a first preset ratio in the principal component direction where the normal vectors are dense, and to widen the trust domain boundary of the main problem decision variables according to a second preset ratio in the secondary component direction where the normal vectors are sparse, thereby reshaping the direction-adaptive trust domain, includes: Extract the normal vectors of all active tangent planes with activity fractions higher than the inactivation threshold, and then perform dimensionless standardization on the normal vectors according to the dimensions of the main problem decision variables before combining them to form a normal vector matrix; Perform eigenvalue decomposition on the covariance matrix corresponding to the normal vector matrix to extract eigenvalues ​​and their corresponding direction vectors; The first few directional vectors whose cumulative eigenvalue ratio is greater than a preset threshold are designated as principal component directions with dense normal vectors, and the trust domain boundary of the main problem decision variable in the principal component direction is shrunk according to a first preset ratio. The remaining direction vectors are labeled as secondary component directions with sparse normal vectors. The trust domain boundaries of the main problem decision variables in the secondary component directions are widened according to a second preset ratio to reshape the direction adaptive trust domain.

7. The method according to claim 1, characterized in that, The process of solving the master problem within the direction-adaptive trust domain until the objective function deviation between adjacent iterations is less than the convergence threshold, and outputting the final optimized solution, includes: At the end of each iteration, the objective function value of the current main problem is used as the current overall estimated cost; The absolute difference between the current total estimated cost of the current iteration and the total estimated cost of the previous iteration is calculated as the deviation of the objective function; When the deviation of the objective function is less than the set convergence threshold for a preset number of consecutive times, the algorithm is determined to have converged, and a stop command is triggered to terminate the iteration. The optimal decision variables at the convergence of the current main problem are extracted as the final control parameters of the microgrid source-grid-load-storage collaborative optimization scheme, and the output is the final optimization scheme.

8. A collaborative optimization system for source-grid-load-storage in microgrids, characterized in that, Includes the following modules: The computation module is used to decompose the microgrid source-grid-load-storage coordinated optimization problem into a main problem and several parallel sub-problems corresponding to uncertain scenarios; cluster all uncertain scenarios based on the Mahalanobis distance between parameters of uncertain scenarios; calculate the importance weight of each uncertain scenario according to the clustering characteristics, and select the parallel sub-problems to be solved by roulette wheel sampling according to the importance weight; The generation module is used to generate new tangent planes by solving selected parallel subproblems in parallel; new tangent planes from the same cluster are weighted and fused according to the constraint violation depth relative to the solution of the decision variables of the current main problem to generate consensus tangent planes representing the clusters and add them to the main problem; an initial activity score is set for each consensus tangent plane added to the main problem, the activity score of the consensus tangent plane is updated after the main problem is solved, and the corresponding consensus tangent plane is temporarily frozen when the activity score is lower than the inactivation threshold; The output module is used to utilize the normal vector distribution characteristics of all active tangent planes with activity scores higher than the inactivation threshold to shrink the trust domain boundary of the main problem decision variables in the principal component direction with dense normal vectors according to a first preset ratio, and to widen the trust domain boundary of the main problem decision variables in the secondary component direction with sparse normal vectors according to a second preset ratio, thereby reshaping the direction-adaptive trust domain. Solve the master problem within the direction-adaptive trust domain until the objective function deviation between adjacent iterations is less than the convergence threshold, and output the final optimized solution.

9. The system according to claim 8, characterized in that, The clustering of all uncertain scenarios based on the Mahalanobis distance between parameters of the uncertain scenarios includes: Construct a parameter vector for an uncertain scenario that includes random disturbances in the microgrid; Calculate the covariance matrix of the overall uncertainty scenario parameter vector, and obtain the inverse matrix of the covariance matrix as the accuracy matrix; The distance between any two sets of uncertain scenario parameter vectors is calculated using the Mahalanobis distance calculation formula, with the accuracy matrix used as the weight matrix in the calculation formula. A density-based spatial clustering algorithm is applied to divide all uncertain scenarios into multiple non-overlapping clusters based on all calculated Mahalanobis distance values.

10. The system according to claim 8, characterized in that, The process of calculating the importance weights of each uncertainty scenario based on the characteristics of clustering, and then using roulette wheel sampling to select the parallel sub-problems to be solved based on these importance weights, includes: The total number of uncertain scenario samples within each cluster is used as the cluster size. Calculate the arithmetic mean of the Mahalanobis distances from each uncertainty scenario to the cluster center within the same cluster, and take the reciprocal of the mean as the cluster density value; The cluster size is multiplied by the cluster density to obtain the basic cluster weight. The basic cluster weight is then multiplied by the reciprocal of the Mahalanobis distance from the target uncertainty scenario to its cluster center to obtain the importance weight of the parallel subproblem corresponding to the target uncertainty scenario. Calculate the proportion of each importance weight in the total weight of all parallel subproblems, set the interval width based on this proportion, generate random numbers using the roulette wheel sampling algorithm based on the interval width, and select the corresponding parallel subproblems for parallel solution.

Images