A power system multi-instability coordination optimization control method based on SHAP security boundary

Abstract

The application discloses a power system multi-instability collaborative optimization control method based on a SHAP safety boundary, and comprises the following steps: S1, constructing an instability evaluation model; S2, constructing a SHAP contribution matrix and screening decision variables; S3, constructing a multi-instability collaborative optimization control model; and S4, solving by adopting a hierarchical cascade strategy.The application explicitly expresses the black-box instability risk as an optimizable safety boundary constraint, and can simultaneously solve multiple coupled instabilities only by adjusting the power distribution of network-constructing type and network-following type converters, has small adjustment cost, high success rate, and an interpretable result, and is suitable for the safe and stable operation control of a high-proportion power electronic power system.

Description

Technical Field

[0001] This invention belongs to the field of power system safe and stable operation and optimization control technology, and in particular relates to a power system multi-instability collaborative optimization control method based on SHAP safety boundary. Background Technology

[0002] With the large-scale integration of new energy sources such as wind and solar power into the power grid via power electronic converters, the degree of power electronics integration in the power system continues to deepen. Grid-mounted (GFM) converters and grid-following (GFL) converters have become the two basic forms of new energy grid connection in the new power system: GFM converters can autonomously establish voltage and frequency references and have a certain degree of active support capability; GFL converters rely on phase-locked loops to track grid voltage, offering fast response but limited support capability. Together with the synchronous generators (SG) retained in the grid, they constitute a multi-source power structure, and their power distribution method profoundly affects the dynamic characteristics and stability level of the system.

[0003] As the proportion of synchronous generating units declines, system inertia and short-circuit capacity continue to weaken, leading to a complex situation where multiple instability modes coexist and are coupled together. Under a single large disturbance, the system may simultaneously face multiple typical instability risks, including synchronous instability, oscillatory instability, frequency exceeding limits, low-voltage ride-through exceeding limits, high-voltage ride-through exceeding limits, voltage collapse, continuous current saturation, and repeated low-voltage ride-throughs. Furthermore, the inducing mechanisms of each instability mode are different and intertwined, meaning a single fault can trigger multiple, almost inevitable, instabilities. Traditional stability analysis methods based on mechanistic models require detailed mathematical models for each type of instability, resulting in a heavy modeling burden, strong parameter dependence, and difficulty in coordinating multiple coupled instabilities within a unified framework. They also lack generalization ability and fail to meet the demands of rapid and coordinated decision-making in actual dispatching.

[0004] Data-driven approaches offer new avenues for instability assessment. Machine learning models, such as gradient boosting trees and neural networks, can predict various instability risks with high accuracy, but they are essentially black-box models: the models only output instability probabilities, unable to explain the source of risk to schedulers, and cannot be directly converted into executable constraints for scheduling. Although interpretable artificial intelligence technologies, such as SHAP (SHapley Additive exPlanations), can quantify the contribution of each operational feature to instability prediction, most existing research stops at post-hoc explanations and risk attribution demonstrations, rarely constructing feature-specific explanatory information into executable numerical constraints in optimization models; existing attempts to embed learning models into scheduling optimization either lack interpretability at the feature-specific level or only target a single instability mode, making it difficult to support coordinated control in multi-instability coupled scenarios.

[0005] There are also significant shortcomings at the solution level. First, the safety constraints of the embedded learning model are complex nonlinear functions of the input variables. Black-box numerical solutions are prone to getting trapped in local optima or even failing to converge. The step-like characteristics of the tree model output can also cause algorithms based on numerical gradients to fail in flat regions. Second, although precisely encoding the tree model as mixed integer programming can obtain global optima and rigorous proof of infeasibility, it can only express the constraints at the level of the total output of the model, and cannot constrain the risk contribution pattern of individual features. This may produce fragile solutions where "positive and negative contributions of multiple features cancel each other out." The total amount meets the target, but individual features still carry high positive risks, resulting in poor robustness to model errors and operating condition drift. Third, when the problem becomes infeasible due to limited control resources, existing methods lack a systematic graded relaxation mechanism and failure attribution methods. Adjustments are highly dependent on converter power allocation, and synchronous machine plans cannot be orderly incorporated into the control echelon. At the same time, it is difficult to intuitively explain to dispatchers the core reasons for the limited unit control, and it is also impossible to clearly provide the optimization direction and execution plan for subsequent control operations.

[0006] In summary, to simultaneously ensure the coordinated suppression of multiple coupled instability risks, the minimization of adjustment costs, and the interpretability of control decisions in power systems with a high proportion of power electronics, it is necessary to design a power system multi-instability collaborative optimization control method based on SHAP safety boundary. Summary of the Invention

[0007] To address the challenges in existing technologies such as the difficulty in coordinating the handling of multiple instability coupling risks, the inability to transform black-box risk assessments into executable constraints, the existence of fragile solutions with contribution offsetting in total constraints, and the lack of tiered relaxation and attribution methods when a problem is infeasible, this invention provides a power system multi-instability coordinated optimization control method based on SHAP safety boundaries. This method explicitly transforms the feature-by-feature attribution information from interpretable artificial intelligence into safety boundary constraints and risk budget constraints for the optimization model. By tiered and cascaded solving to coordinate the control sequence of grid-connected converters, grid-linked converters, and synchronous generators, it achieves minimum adjustment costs and transparent decision-making processes while ensuring the coordinated safety of multiple instabilities in the system.

[0008] A collaborative optimization control method for multiple instabilities in power systems based on SHAP safety boundaries includes the following steps: (1) Generate a sample dataset covering normal operation, single instability and compound instability scenarios through electromagnetic transient simulation; take the system operation characteristics as input and the labels of each instability mode as output, use the CatBoost algorithm to train the instability evaluation model of n typical instability modes respectively, and calibrate the instability probability safety threshold of each model. (2) For the working condition to be evaluated, calculate the SHAP contribution value for each feature using each instability assessment model, and construct a SHAP contribution matrix with instability mode as the row and operating feature as the column; filter out the set of hazardous features that have a positive contribution to any instability mode by column, calculate the comprehensive contribution value of each hazardous feature and sort them in descending order to determine the number one control variable; take the intersection of the set of hazardous features and the set of physically adjustable variables as the hazard control group; (3) The objective function is to minimize the sum of the adjustments of all control variables, taking into account the upper and lower limits of unit output, system power balance constraints, node proportion constraints and node output constraints, and constructing a multi-instability SHAP safety boundary constraint to form a multi-instability collaborative optimization control model. (4) The solution is obtained by adopting a hierarchical cascade strategy. First, it enters the robust level and uses the risk control group selected in step (2) as the decision variable. On the basis of the multi-instability SHAP safety boundary constraint, a single SHAP risk budget constraint is added, that is, the contribution of the number one control variable is constrained. The robust solution in the deep safety domain is obtained by black box numerical method. If the standard is not met, it enters the standard level and the tree structure of each instability assessment model is precisely encoded into a mixed integer linear programming. The global optimal solution is obtained on the control variables of all grid-type and follow-grid-type converters or the infeasibility proof under the variable set is given. If it is still not feasible, it enters the extended level, the share of each node synchronous machine is unfrozen and a weighted objective function is adopted so that the solver prioritizes the adjustment of the converter and the synchronous machine only undertakes the necessary adjustment. The solutions obtained at each level are issued for execution after the instability probability is verified. When all three levels are not feasible, the off-field control suggestions based on SHAP attribution are output.

[0009] Furthermore, in step (1), the system operation characteristics include the power reference value and proportion of each node's grid-type generator unit, the power reference value and proportion of the grid-connected generator unit, the power reference value and proportion of the synchronous generator unit, the total output of each node, and the system load level.

[0010] The typical instability modes of type n are: synchronous instability, oscillatory instability, frequency exceeding limit, low ride-through limit, voltage collapse, high ride-through limit, continuous current saturation, and repeated low-voltage ride-through.

[0011] Furthermore, in step (1), during the sample dataset generation stage, a perturbation is applied to the share of the synchronizing machine and the total output of the nodes so that the effective domain of the instability assessment model covers the decision space of the extended level in step (4).

[0012] Further, in step (2), the feature-wise SHAP contribution value is calculated using each instability assessment model, specifically as follows: The SHAP method is used to decompose the predictions of the instability assessment model into a linear sum of the contributions of each input feature: ; In the formula, For the first Instability assessment model for operating conditions The original output, This serves as the baseline value for the model; For the first Operating characteristics under working conditions Next to the The SHAP contribution value for instability-like features is positive, indicating that the feature drives the system toward instability, and negative, indicating that the feature enhances the system's stability.

[0013] Furthermore, in step (2), in the operating characteristics, the power reference value of the same unit and its corresponding proportion physically belong to the same adjustment knob, and the two are merged into the same control group, and their SHAP contribution is calculated by summing and combining within the group.

[0014] Furthermore, in step (3), the power supply of each node consists of three types: synchronous machines, grid-connected units, and grid-following units. The specific node ratio constraint is as follows: ; ; In the formula, , , They are nodes j The proportion of network type, follow-up network type and synchronous machine.

[0015] Furthermore, in step (3), the safety boundary constraints of the multi-instability SHAP are constructed as follows: The overall SHAP value of each instability mode must not exceed the corresponding SHAP safety threshold; where, the first k SHAP safety threshold for instability-like modes By the k Instability probability safety threshold for instability-like modes The only certainty is that the formula is: ; in, For the first The baseline value of SHAP contribution for the instability assessment model.

[0016] Furthermore, in step (4), the single-group SHAP risk budget constraint is specifically as follows: ; In the formula, For the first The control group for the first k Group contributions that are active and unstable Contributes to power reference value, Contribution to the corresponding percentage; This represents the set of active instability modes under the current operating conditions. For the set of hazard control groups to be transferred; The upper limit of the risk budget is set to 0. When the upper limit is set to 0, the risk driver must be adjusted to a neutral or stable source.

[0017] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention upgrades SHAP from a post-analysis tool to a pre-decision constraint. Existing applications of SHAP in power system security assessments are limited to post-analysis and cannot be transformed into executable control instructions. This invention proves the strict equivalence between the comprehensive SHAP value boundary and the instability probability threshold. Based on this, it constructs multi-instability SHAP safety boundary constraints, explicitly transforming black-box risks into numerical constraints that the solver can directly execute. Furthermore, relying on the SHAP contribution matrix, it ensures that the safety boundaries of n types of coupled instability coexist in the same variable space and are simultaneously satisfied by the same set of control variables, achieving a closed loop from "post-interpretation" to "pre-constraint—optimization decision—control execution."

[0018] 2. This invention realizes the explicit transformation of black-box instability risk into executable safety constraints and multi-instability collaborative defense. By constructing a SHAP contribution matrix of "instability mode × operational characteristics" and proving the strict equivalence between the integrated SHAP boundary and the instability probability threshold, the risk score, which could only be explained ex post factoly, is transformed into a numerical constraint that the optimizer can directly execute. The boundaries of n types of coupled instability coexist in the same variable space and are simultaneously satisfied by a set of control variables, avoiding the heavy burden and insufficient generalization of mechanism modeling methods that model one by one, and realizing integrated collaborative suppression in the scenario of multiple instability induced by a single fault.

[0019] 3. This invention achieves low-cost, interpretable, and precise control through hazard feature screening and minimum adjustment targets. The intersection of the hazard features screened by column in the SHAP matrix and the physically adjustable set constitutes a compact decision space, and the control actions are concentrated on the risk-driven source. The power reference values ​​and proportions of the same unit are merged into the same control group, avoiding duplicate measurement and constraint redundancy. With the goal of minimizing the adjustment amount and keeping the synchronous machine plan unchanged, the triple instability is resolved simultaneously by only fine-tuning the converter allocation in the single-scenario verification. The equipment movement is small, the deviation from the economic dispatch baseline is small, and each constraint and each adjustment step can be mapped to a clear physical variable, which is completely transparent to the dispatchers.

[0020] 4. The risk budget constraint of this invention eliminates the contribution-offsetting fragile solution under the total constraint, thus improving the robustness of the control scheme. In addition to the total probability constraint, a morphological constraint is added to the adjusted hazard control group, which requires that the contribution of a single group of SHAP does not exceed the budget limit. This forces the risk driving source to be adjusted to a neutral or stable source, rather than relying on the negative contribution offset by other features. The total constraint and the budget constraint form a two-layer structure of "necessary condition + morphological constraint", which is not redundant. The resulting robust solution compresses the instability probability to a deeper safety margin below the threshold, and has stronger tolerance to model errors and operating condition drift.

[0021] 5. The cascaded solution mechanism of this invention comprehensively considers the lower bound of success rate, global optimality, and interpretability of failure. The robust level avoids gradient ray scanning combined with local polishing to overcome numerical gradient failure caused by the stair-step characteristics of the tree model; the standard level provides a globally optimal solution and a rigorous proof of infeasibility with the precise mixed-integer linear programming encoding of the tree model, ensuring that the overall success rate is no lower than that of existing single-level methods; the extended level unfreezes the share of the synchronous machine while maintaining the linear structure, and the weighted objective ensures the control order of converter priority and synchronous machine fallback; when all three levels are infeasible, SHAP attribution quantitatively locates the cause of failure to unadjustable external factors, forming a complete decision-making closed loop of "gradual relaxation, failure attribution, and off-site suggestions". The 94.80% success rate and hierarchical statistics of the batch verification in the 1883 scenario prove the engineering practicality of the method. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 This is an overall flowchart of a power system multi-instability collaborative optimization control method based on SHAP safety boundary according to an embodiment of the present invention.

[0024] Figure 2 This is the topology diagram of the 8-machine simulation system in an embodiment of the present invention.

[0025] Figure 3 This is the SHAP force diagram of the optimized first three types of instability assessment models in this embodiment of the invention.

[0026] Figure 4 This is a SHAP risk contribution matrix diagram under the complex instability scenario in this embodiment of the invention.

[0027] Figure 5 These are the SHAP force diagrams of the three types of instability assessment models optimized in this embodiment of the invention.

[0028] Figure 6 This is a comparison chart of instability probabilities before and after optimization in an embodiment of the present invention.

[0029] Figure 7 This is a comparison chart of the electromagnetic transient dynamic curves of key state variables before and after optimization in an embodiment of the present invention. Detailed Implementation

[0030] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0031] It should be noted that, unless otherwise specified, the features in the following embodiments and implementation methods can be combined with each other.

[0032] A collaborative optimization control method for multiple instabilities in power systems based on the SHAP safety boundary is proposed. This method embeds the SHAP safety boundary constraints and feature-by-feature risk budget constraints into an optimization control model aiming to minimize adjustment amounts. The model is then solved and relaxed in a cascaded order of "robust level—standard level—extended level." Figure 1 As shown, the method mainly includes the following steps: S1. Constructing Instability Assessment Models: Generate sample datasets covering normal operation, single instability, and combined instability scenarios through electromagnetic transient simulation. Using system operation characteristics as input and labels of each instability mode as output, train instability assessment models for n typical instability modes using the CatBoost algorithm, and calibrate the instability probability safety threshold for each model.

[0033] S2. Construct the SHAP contribution matrix and select decision variables: For the operating condition to be evaluated, calculate the feature-by-feature SHAP contribution value using each instability assessment model, and construct the SHAP contribution matrix with instability modes as rows and operating characteristics as columns; select the set of hazardous features that have a positive contribution to any instability mode by column, calculate the comprehensive contribution value of each hazardous feature and sort them in descending order to determine the number one control variable; take the intersection of the set of hazardous features and the set of physically adjustable variables as the hazardous control group, and merge the power reference value of the same unit and its corresponding proportion into the same control group.

[0034] S3. Construct a multi-instability collaborative optimization control model: The objective function is to minimize the sum of the adjustments of all control variables. It takes into account the upper and lower limits of unit output, system power balance constraints, the proportion of node network type / following network type / synchronous machine, and node output constraints. It also constructs a multi-instability SHAP safety boundary constraint, requiring that the comprehensive SHAP value of each active instability mode does not exceed the corresponding safety threshold. This boundary constraint is strictly equivalent to the instability probability not exceeding the probability threshold.

[0035] S4. Cascaded Solving and Execution: First, the robust level is entered. Using the selected hazard control group as the decision variable, a single-group SHAP risk budget constraint is added to the probability constraint. A black-box numerical method is used to obtain a robust solution that penetrates deep into the safety domain. If the standard is not met, the standard level is entered. The tree structure of each instability assessment model is precisely encoded into a mixed-integer linear programming problem. The global optimal solution is obtained on all converter control variables, or an infeasibility proof is given for the set of variables. If it is still infeasible, the extended level is entered. The share of each node's synchronous machine is unfrozen, and a weighted objective function is used so that the solver prioritizes adjusting the converter, and the synchronous machine only undertakes the necessary adjustment. The solutions obtained at each level must be verified by instability probability. Only after all active instability modes meet the standard can they be issued for execution. If all three levels are proven infeasible, the risk contribution locked by the frozen characteristics is given based on the SHAP attribution output. Off-site control suggestions such as relaxing the synchronous machine plan, putting in energy storage, or shedding load are given.

[0036] The method proposed in this invention explicitly transforms the black-box instability risk into two layers of optimizable constraints: the safety boundary and the risk budget, thereby achieving integrated and coordinated suppression of multiple coupled instabilities. By screening for dangerous features to compress the decision space and by coordinating the control sequence of the converter and the synchronous machine through a cascaded relaxation mechanism, it ensures the lower limit of success rate while taking into account the robustness of understanding, adjustment costs, and interpretability.

[0037] Furthermore, the detailed steps of step S1 include: S101. Define operating characteristics and instability labels. Operating characteristics encompass m-dimensional observables characterizing the system's operating points, including but not limited to: power reference values ​​and proportions of grid-connected units at each node, power reference values ​​and proportions of grid-connected units, power reference values ​​and proportions of synchronous units, total output of each node, and system load levels. Instability labels cover n typical instability modes.

[0038] S102. Generate batch electromagnetic transient simulation samples. In the electromagnetic transient simulation platform, large disturbances such as three-phase short circuits are applied to different nodes under different load levels and different renewable energy penetration rates. The operating characteristics before the disturbance and the labels of whether each instability mode occurred after the disturbance are collected to form a sample dataset covering normal, single instability, and compound instability. Furthermore, to support the solution accuracy of the extended level in step S4, disturbances can be applied to the synchronous machine share and the total output of the nodes during the sample generation stage, so that the effective domain of the evaluation model covers the decision space of the extended level.

[0039] S103. Train the instability assessment model and calibrate the threshold. Train a CatBoost binary classification model independently for each type of instability mode, with the model input being the running feature vector. The output is the instability probability of this mode. Let the first... k Instability assessment model is Its instability probability is And calibrate a safety threshold for the probability of instability for each model. In this embodiment, a unified approach is taken. =0.2, meaning an instability probability below 20% is considered satisfactory. To ensure consistency between the constraint and validation criteria, the constraint objective for each mode during optimization is uniformly set to the smaller of the probability threshold and the model classification threshold.

[0040] Furthermore, the detailed steps of step S2 include: S201. Calculate the feature-wise SHAP contribution. The SHAP method decomposes the model prediction into a linear sum of the contributions of each input feature: Equation (1): ; In the formula, For the first Instability assessment model for operating conditions The original output; This serves as the baseline value for the model; For the first Operating characteristics under working conditions Next to the The SHAP contribution value for instability-like features is positive, indicating that the feature promotes system instability, and negative, indicating that the feature enhances system stability. The original output is mapped to instability probability using the sigmoid function: Equation (2): ; In the formula, It is the sigmoid function; For the first The instability probability of the instability-like mode.

[0041] S202. Construct the SHAP contribution matrix. Arrange the feature-wise contributions of the n types of instability assessment models to the same working condition by row to obtain the SHAP contribution matrix: Equation (3): ; In the formula, each row of the matrix corresponds to a type of instability mode, and each column corresponds to a type of operational characteristic; elements Measurement of the first The feature is related to the first The direction and intensity of the impact of instability-like risks.

[0042] S203. Screen and rank hazardous features. Define the set D of hazardous features as all features that make a positive contribution under any instability mode: Equation (4): ; For the Summing each feature column yields its target. Overall contribution of instability-like factors: Equation (5): ; Arrange the hazard characteristics in descending order of their overall contribution value: Equation (6): ; In the formula, d represents the number of dangerous features. The ranking results also provide the strength of each feature driving instability and the priority of regulation; the feature ranked first is the number one control variable.

[0043] S204. Determine the optimal decision variables. The intersection of the set of hazardous features D and the set of physically adjustable variables A is taken as the set of optimal decision variables, ensuring that safety constraints remain compact and that control actions are concentrated on risk drivers. In this invention, the power reference value and its corresponding proportion of the same unit physically belong to the same control knob. They are grouped into the same control group for screening, budget constraints, and control ranking, and their SHAP contribution is calculated by summing within the group.

[0044] Furthermore, the detailed steps of step S3 include: S301. Construct the objective function. The adjustment range of the control variables determines the equipment's movement distance, the intensity of system disturbances, and the scheduling cost. The objective is to minimize the sum of the adjustments of all controllable variables. Equation (7): ; In the formula, C is the set of all controllable variables of the system; To determine the optimized values ​​for the control variables; This is the initial operating value of the variable under the original operating conditions.

[0045] S302. Establish normal operating constraints. The active power output of each generating unit should remain within the rated range: Equation (8): ; In the formula, For the first The active power output of the generator unit; and These are its lower and upper limits of output, respectively; This refers to the collection of all power generation units in the system. The system's active power is balanced in real time, and the sum of the outputs of each node equals the total system load. Equation (9): ; In the formula, For nodes They always contribute their efforts; This represents the total active power load of the system. It is a set of nodes.

[0046] S303. Establish node proportion constraints. The power supply for each node consists of three types: synchronous generators, grid-connected units, and grid-following units. Synchronous generators bear the system inertia and provide the main voltage support; their output plan remains unchanged and they do not participate in optimization. Grid-connected and grid-following converters have fast response and low regulation costs, serving as the main control components. The proportions of the three types must meet the following requirements: Equation (10): ; ; In the formula, , , They are nodes The proportion of network type, follow-up network type and synchronous machine.

[0047] S304. Construct node output constraints. Multiply the total node output by its respective proportion to obtain the corresponding unit output, ensuring consistency between the total node output and the actual output of the two types of converters: Equation (11): ; ; In the formula, and They are nodes The output of the grid-type and the output of the follow-grid type are defined. Under the condition that the total output of the nodes and the output of the synchronous machine remain unchanged, the output of the follow-grid type is linearly derived from the output of the grid-type, and the proportion of each is also a linear function of the control variables.

[0048] S305. Construct a multi-instability SHAP safety boundary constraint. To simultaneously defend against all n types of instability risks, the combined SHAP value of each instability mode must not exceed its safety threshold: Equation (12): ; In the formula, For working conditions Next The overall SHAP risk value of instability-like patterns; For the first The SHAP safety threshold for instability-like modes. From equations (1) and (2), it can be seen that this boundary constraint is strictly equivalent to the instability probability constraint: Equation (13): ; In the formula, For the first The instability probability safety threshold for instability-like modes. Based on this, the safety threshold... It can be determined by probability threshold. The only certainty is that the originally inexplicable and unexecutable black-box risk score is explicitly transformed into a safety boundary constraint that the optimizer can directly execute, and multiple coupling instabilities are satisfied by a set of variables simultaneously within the same variable space.

[0049] Furthermore, the detailed steps of step S4 include: S401, Robust Solution. Using the hazard control groups selected in step S2 as decision variables, based on the probability constraints of equation (13), a single-group SHAP risk budget constraint is added to each adjusted hazard control group: Equation (14): ; In the formula, For the first Each control group (power reference value contribution) Contribution of corresponding percentage (sum of) for the first Contributions from groups exhibiting activity instability; This represents the set of active instability modes under the current operating conditions. For the set of hazard control groups to be transferred; The upper limit of the risk budget is 0. When it is set to 0, the risk driving source is required to be adjusted to a neutral or stable source. The total amount constraint (13) ensures that the target is met, which is a necessary condition. The budget constraint (14) further constrains the distribution pattern of the risk. The two are not redundant: the budget constraint removes the fragile solution where "the positive and negative contributions of multiple features cancel each other out and the total amount meets the target but individual features still carry high positive risks", making the obtained solution more robust to model errors and operating condition drift.

[0050] Since the contribution of feature-wise SHAP is a complex nonlinear function of the decision variables and lacks analytical gradient, and the tree model output is stepped, causing the numerical gradient to fail in flat regions, the robust stage adopts a black-box solution strategy combining gradient-free and local refinement: First, a regression ray is constructed along the dangerous regression direction indicated by the SHAP matrix, and the feasible interval is located by coarse scanning of the regression ratio. Then, the boundary is sharpened using a bisection method to obtain the feasible point with the minimum adjustment on the ray. Then, using this feasible point as the initial point, the Sequential Quadratic Programming (SLSQP) algorithm is used with an appropriately amplified numerical gradient step size for local refinement, supplemented by a multi-initial-point mechanism to prevent getting trapped in local optima. The refined solution that passes the verification and has a better objective is adopted; otherwise, the ray solution is retained. The failure of the robust stage does not affect the overall process, and subsequent stages will handle the fallback.

[0051] S402, Standard Level Solution. If the robustness level is not met, the budget constraint is removed, and the exact model is solved using the power reference values ​​of all grid-type converters and their linear derivatives as decision variables. The symmetric tree structure of each instability assessment model is explicitly encoded as a mixed-integer linear programming problem: a binary indicator variable is introduced for each split condition of each tree, and the feature values ​​are associated using the Big M method. Indicator variables with the same feature and the same split threshold are globally deduplicated and shared. Leaf selection is expressed by a combination of continuous variables constrained by the simplex constraint, thus making the original model output... The exact linear expression of the decision variable is given, and the security constraints are written as: Equation (15): ; After linearizing the objective function (7) by introducing auxiliary variables, the entire problem becomes a mixed-integer linear programming problem. The global optimal solution is obtained by the branch and bound solver, or a rigorous proof that the problem is infeasible under the current set of variables is given. The standard level has the same capabilities as the existing single-level methods, ensuring the lower limit of the success rate of the overall solution; the control direction indicated by the SHAP matrix can also be used as the initial value for hot start of mixed-integer programming, accelerating the solution without affecting the optimality and infeasibility proof.

[0052] S403, Extended Stage Solution. If the standard stage proves infeasible, it means that converter power allocation alone cannot simultaneously satisfy all safety boundaries. Therefore, the synchronous machine share P_{SG,j} of each node is unfrozen as a new decision variable. Since the total node output P_{all,j} remains frozen, the network output is still a linear function of the decision variables. Equation (16): ; The structure of the mixed-integer linear programming problem remains completely preserved, with only continuous variables added. The objective function is also changed to a weighted form: Equation (17): ; In the formula, For the set of converter control variables, weights Take 1; For the set of synchronization machine share variables, weights Greater than 1 (10 in this embodiment). The weighted objective ensures that the solver prioritizes the use of the converter's adjustment margin, and the synchronous machine only undertakes the adjustment amount that cannot be changed, so as to retain the inertia and voltage support carried by the original synchronous machine plan to the maximum extent. Before the extended stage is put into operation, it is advisable to expand the sample containing the synchronous machine share disturbance according to the further scheme of step S102 and retrain the evaluation model, and the scheme obtained from the extended stage should be implemented after verification by electromagnetic transient simulation. As a further alternative implementation method, the total output of the node can be unfrozen after the extended stage and the system power balance constraint of equation (9) can be supplemented to realize cross-node power transfer. At this time, a bilinear term appears in equation (11), which can be solved by a solver that supports non-convex quadratic constraints or a black box numerical method.

[0053] S404, Compliance Verification and Execution. Candidate solutions obtained at each stage are substituted back into the original CatBoost evaluation model for instability probability verification, requiring all active instability modes to satisfy... ,in For the first The classification threshold of the model is used to ensure that the optimization constraint and the verification constraint are strictly consistent; after the verification is passed, the difference between the control variables before and after optimization is used as the power adjustment command to be issued and executed.

[0054] S405. Infeasibility Attribution. If all three levels are proven infeasible by the solver, a diagnostic conclusion is output based on SHAP attribution: For each active instability mode, the total positive SHAP contribution corresponding to unadjustable frozen features such as load level and predetermined synchronous machine plan can be statistically analyzed and the features can be sorted to accurately pinpoint the external influencing factors restricting the closure of the system's safety boundary. On this basis, off-site control decision suggestions such as relaxing the synchronous machine output plan, commissioning energy storage devices, and emergency load shedding can be provided to dispatchers, constructing an "explainable, traceable, and guiding" scheduling closed-loop mechanism under control-constrained scenarios, effectively solving the problem that traditional control modes cannot explain the causes of constraints or clarify subsequent operations.

[0055] Furthermore, to verify the effectiveness of the proposed method, simulation verification was performed on an eight-machine power system containing grid-type and grid-connected converters and synchronous generators.

[0056] Figure 2 The simulation system topology is 8 units. A three-phase short-circuit fault is applied at node 5, forming a complex instability scenario with triple coupling of oscillation instability, frequency over-limit, and high voltage ride-through over-limit. After the fault, although each unit is still within the rated limits and can be checked by conventional constraints, the predicted instability probabilities of the three types of instability modes are as high as 98.5%, 90.9%, and 88.3%, respectively, all far exceeding the 20% safety threshold. Figure 3 To optimize the SHAP force diagrams of the first three types of instability assessment models, Figure 4For the SHAP risk contribution matrix in this scenario: the column sum of variable Pref_2_1 (power reference value of node 2 grid-type unit) is significantly higher than that of the other variables, and it has a strong positive contribution to all three types of instability. It is identified as the number one control variable and given priority for regulation. The other risk variables are entered into the regulation sequence in order of comprehensive contribution value.

[0057] The optimization solution only adjusts the power distribution of grid-type and grid-following type converters, while keeping the output of synchronous machines unchanged. Table 1 shows the comparison of the output of each power source before and after optimization.

[0058] Table 1 Comparison of power output before and after optimization Figure 5 The SHAP plot for the three types of instability assessment models after optimization shows that, compared with the original model, the overall positive contribution of the three models has significantly decreased, the positive weight of the number one variable Pref_2_1 on the three types of instability has decreased significantly, and the risk distribution pattern has been reshaped. Figure 6 To optimize the comparison of instability probabilities before and after: all three types of instability probabilities have fallen from the high-risk range to within the safe threshold of 0.2, and multiple instabilities have been resolved simultaneously. Figure 7 To compare the electromagnetic transient dynamic curves of key state quantities before and after optimization: Before optimization, the fault caused severe frequency oscillations, continuous frequency exceeding the limit and over-frequency operation. After optimization, the transient process was effectively suppressed, and the active power, frequency and voltage quickly recovered to near the rated values.

[0059] Further batch verification was conducted on 5460 operating scenarios. These scenarios were generated by random disturbances to daily load levels and renewable energy penetration rates during typical winter and summer operating days. Three-phase short-circuit faults were applied at six grid-connected nodes, covering a large number of complex instability conditions. Considering the limited effect of regulating the active power output of generating units on the latter four types of instability (n types), this batch simulation verification only targeted the first four types of 1883 operating scenarios. All scenarios were automatically solved using the same process, without individual parameter tuning. Success was only achieved when the probability of all active instability modes in the operating condition dropped below the threshold. Table 2 presents the statistical results for each mode.

[0060] Table 2. Statistical results of batch optimization for 1883 scenarios Batch results show that the average instability probability of each instability mode decreased by 0.62 to 0.96, and the overall success rate reached 94.80% according to the scenario statistics.

[0061] The embodiments described above provide a detailed explanation of the technical solutions and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the present invention. Any modifications, additions, and equivalent substitutions made within the scope of the principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A collaborative optimization control method for multiple instabilities in power systems based on SHAP safety boundaries, characterized in that, Includes the following steps: (1) Generate a sample dataset covering normal operation, single instability and compound instability scenarios through electromagnetic transient simulation; Using system operating characteristics as input and labels of each instability mode as output, the CatBoost algorithm is used to train instability assessment models for n typical instability modes, and the instability probability safety threshold of each model is calibrated. (2) For the working condition to be evaluated, calculate the SHAP contribution value for each feature using each instability assessment model, and construct a SHAP contribution matrix with instability mode as the row and operating feature as the column; filter out the set of hazardous features that have a positive contribution to any instability mode by column, calculate the comprehensive contribution value of each hazardous feature and sort them in descending order to determine the number one control variable; take the intersection of the set of hazardous features and the set of physically adjustable variables as the hazard control group; (3) The objective function is to minimize the sum of the adjustments of all control variables, taking into account the upper and lower limits of unit output, system power balance constraints, node proportion constraints and node output constraints, and constructing a multi-instability SHAP safety boundary constraint to form a multi-instability collaborative optimization control model. (4) The solution is obtained by adopting a hierarchical cascade strategy. First, it enters the robust level and uses the risk control group selected in step (2) as the decision variable. On the basis of the multi-instability SHAP safety boundary constraint, a single SHAP risk budget constraint is added, that is, the contribution of the number one control variable is constrained. The robust solution in the deep safety domain is obtained by black box numerical method. If the standard is not met, it enters the standard level and the tree structure of each instability assessment model is precisely encoded into a mixed integer linear programming. The global optimal solution is obtained on the control variables of all grid-type and follow-grid-type converters or the infeasibility proof under the variable set is given. If it is still not feasible, it enters the extended level, the share of each node synchronous machine is unfrozen and a weighted objective function is adopted so that the solver prioritizes the adjustment of the converter and the synchronous machine only undertakes the necessary adjustment. The solutions obtained at each level are issued for execution after the instability probability is verified. When all three levels are not feasible, the off-field control suggestions based on SHAP attribution are output.

2. The power system multi-instability collaborative optimization control method based on SHAP safety boundary as described in claim 1, characterized in that, In step (1), the system operation characteristics include the power reference value and proportion of each node's grid-type units, the power reference value and proportion of the grid-connected units, the power reference value and proportion of the synchronous units, the total output of each node, and the system load level.

3. The power system multi-instability collaborative optimization control method based on SHAP safety boundary according to claim 1, characterized in that, In step (1), during the sample dataset generation stage, a perturbation is applied to the share of the synchronizing machine and the total output of the nodes so that the effective domain of the instability assessment model covers the decision space of the extended level in step (4).

4. The power system multi-instability collaborative optimization control method based on SHAP safety boundary according to claim 1, characterized in that, In step (2), the feature-wise SHAP contribution value is calculated using each instability assessment model, specifically as follows: The SHAP method is used to decompose the predictions of the instability assessment model into a linear sum of the contributions of each input feature: ； In the formula, For the first Instability assessment model for operating conditions The original output, This serves as the baseline value for the model; For the first Operating characteristics under working conditions Next to the The SHAP contribution value for instability-like features is positive, indicating that the feature drives the system toward instability, and negative, indicating that the feature enhances the system's stability.

5. The power system multi-instability collaborative optimization control method based on SHAP safety boundary according to claim 1, characterized in that, In step (2), in the operating characteristics, the power reference value of the same unit and its corresponding proportion physically belong to the same adjustment knob. The two are merged into the same control group, and their SHAP contribution is calculated by summing and combining within the group.

6. The power system multi-instability collaborative optimization control method based on SHAP safety boundary according to claim 1, characterized in that, In step (3), the power supply of each node consists of three types: synchronous machines, grid-connected units, and grid-following units. The specific node ratio constraint is as follows: ； ； In the formula, , , They are nodes j The proportion of network type, follow-up network type and synchronous machine.

7. The power system multi-instability collaborative optimization control method based on SHAP safety boundary according to claim 1, characterized in that, In step (3), the safety boundary constraints of the multi-instability SHAP are constructed, specifically as follows: The overall SHAP value of each instability mode must not exceed the corresponding SHAP safety threshold; where, the first k SHAP safety threshold for instability-like modes By the k Instability probability safety threshold of instability-like modes The only certainty is that the formula is: ； in, For the first The baseline value of SHAP contribution for the instability assessment model.

8. The power system multi-instability collaborative optimization control method based on SHAP safety boundary according to claim 1, characterized in that, In step (4), the single-group SHAP risk budget constraint is specifically as follows: ； In the formula, For the first The control group for the first k Group contributions from active and unstable groups Contributes to power reference value, Contribution to the corresponding percentage; This represents the set of active instability modes under the current operating conditions. For the set of hazard control groups to be transferred; The upper limit of the risk budget is set to 0. When the upper limit is set to 0, the risk driver must be adjusted to a neutral or stable source.

Images