A park power restoration optimization method considering ice disaster icing prediction and multi-source collaboration

Abstract

The application provides a park power restoration optimization method considering ice disaster icing prediction and multi-source cooperation, and belongs to the technical field of electric power. The application predicts icing thickness time series and calculates ice load degradation coefficient by physically embedding an icing agent model, then corrects the edge weight of the power restoration path algorithm with the degradation coefficient, locates the fault position and the root cause by using an artificial intelligence causal double-flow model, outputs a power flow feasible power restoration scheme by using a Benders decomposition method, quantifies the secondary fault probability under different closing timing by combining Monte Carlo tree search and an icing shedding stochastic process model, and coordinates the distributed photovoltaic, energy storage and diesel generator multi-source output to complete park multi-island power restoration within an acceptable risk threshold, thereby solving the technical problems that the power restoration path planning cannot dynamically integrate the icing thickness time series prediction result in the ice disaster scenario, and that the line carrying capacity evaluation is inaccurate after power restoration, resulting in frequent secondary tripping.

Description

Technical Field

[0001] This invention belongs to the field of power technology, and specifically relates to an optimized method for power restoration in industrial parks that considers ice storm icing prediction and multi-source coordination. Background Technology

[0002] Power grid restoration optimization in industrial parks is an important research direction in the field of power system post-disaster recovery. In ice storm scenarios, traditional restoration methods typically rely on static line parameters, employing minimum spanning tree algorithms or integer programming to generate restoration paths, and combining this with manual inspection results to determine line availability. Currently, these methods are widely used in power grid restoration scenarios following natural disasters such as typhoons and floods. Dispatchers execute closing operations based on preset line rated current carrying capacity and topology constraints, and some systems introduce power flow calculations to verify the feasibility of restoration schemes.

[0003] However, the aforementioned traditional methods have inherent flaws under ice storm conditions. Icing causes an additional thermal resistance layer to form on the conductor surface, reducing the actual heat dissipation capacity of the line and thus compressing the actual current carrying capacity. However, traditional methods still use the rated current carrying capacity as the constraint boundary during the route planning stage, which cannot reflect the dynamic capacity degradation process caused by icing.

[0004] Under conditions where ice storms continue to develop and ice thickness increases over time, traditional methods lack the ability to predict the temporal evolution of ice thickness. This makes it impossible to identify the degree of capacity degradation of heavily iced lines during power restoration path selection. Consequently, the actual current carrying capacity of the lines after reactivation exceeds the safety threshold under icing conditions, leading to repeated tripping and causing repeated power restoration failures on multiple isolated areas within the industrial park. In other words, existing technologies suffer from the technical problem of failing to dynamically integrate the temporal prediction results of ice thickness in power restoration path planning during ice storm scenarios, resulting in inaccurate assessments of line capacity after power restoration and frequent re-tripping. Summary of the Invention

[0005] In view of this, the present invention provides a power restoration optimization method for industrial parks that considers ice accretion prediction and multi-source collaboration. This method can solve the technical problem in the prior art that the power restoration path planning of the distribution network in the ice disaster scenario cannot dynamically integrate the time-series prediction results of ice thickness, which leads to frequent re-tripping after power restoration due to inaccurate line carrying capacity assessment.

[0006] This invention is implemented as follows: This invention provides an optimized method for power restoration in industrial parks that considers ice storm icing prediction and multi-source coordination, comprising the following steps:

[0007] Collect conductor surface temperature, ambient temperature, wind speed, liquid water content and median volume diameter, input them into the physical embedded icing proxy model, and output the time series prediction results of icing thickness and ice load degradation coefficient.

[0008] The protection device action time stamp sequence, fault waveform electrical quantities, distribution network topology adjacency matrix, and node icing degree characteristics are input into the causal dual-current fault tracing and localization model, and the fault location set and fault root cause identification results are output.

[0009] Using the ice load degradation coefficient as the edge weight correction, the minimum cost power restoration path algorithm after ice disaster is called to generate the optimal power restoration path set that satisfies the radial topology constraint and the line current carrying capacity constraint.

[0010] Taking the optimal power restoration path set as input, the Benders decomposition method is used to divide the power restoration optimization into the main problem of topology restoration binary decision and the sub-problem of power flow feasibility verification, and outputs feasible power restoration schemes;

[0011] Using the ice adhesion strength decay curve and current meteorological conditions as input, a joint decision-making framework is constructed based on Monte Carlo tree search and ice shedding stochastic process model, and the probability of re-failure under different waiting times is output.

[0012] Based on feasible power restoration plans and the probability of recurrence, the optimal timing for closing the circuit is determined within an acceptable risk threshold. The output of distributed photovoltaic, energy storage and diesel generators is coordinated, and the closing command is executed in order of priority to complete the power restoration of multiple isolated areas in the park.

[0013] The physical embedded icing proxy model is based on deep learning. It uses the partial differential relation of the Makkonen equation as the physical loss term to embed into the neural network training process, and uses the neural network to replace numerical iteration to solve the transient thermal balance equation.

[0014] The physical loss term is composed of a weighted sum of the residuals of the heat conduction equation, the convection heat transfer equation, and the latent heat conservation equation for phase change. The weighting coefficients are determined through a grid search experiment on the validation set that jointly minimizes the prediction error and the physical violation.

[0015] The ice load degradation coefficient is defined as the ratio of the actual current carrying capacity to the rated current carrying capacity of the line under the current icing conditions. The value ranges from 0 to 1 and is calculated by substituting the time-series prediction results of icing thickness into the line heat balance equation.

[0016] The input of the causal dual-flow fault tracing and localization model includes a temporal flow and a spatial flow. The temporal flow is segmented and encoded using a patch embedding encoding method, while the spatial flow is encoded using a graph Transformer. The temporal flow and spatial flow are fused in the middle layer of the model through a cross-attention mechanism.

[0017] The causal dual-flow fault tracing and localization model is connected to a differentiable causal graph module based on a structural causal model after the cross-attention mechanism layer. The differentiable causal graph module maps the features extracted by the network to a causal variable chain consisting of icing thickness, mechanical stress, wire breakage and protection action. The intervention effect is calculated by do-calculus differentiable approximation to identify the root cause of the fault.

[0018] The causal dual-flow fault tracing and localization model internally sets up an iterative refinement loop. In each loop, the fault probability vector corresponding to the current output fault location set is scaled by a temperature coefficient to form a soft label mask, which is then superimposed on the feature vector of the input node and the inference is re-executed. The loop is terminated early when the confidence level meets the confidence level threshold.

[0019] The minimum cost power restoration path algorithm after an ice disaster maps the power restoration problem to a constrained Steiner tree problem. It uses the Dreyfus-Wagner dynamic programming algorithm to solve the problem accurately on the key subgraph. After using the minimum spanning tree approximation as the initial solution for the whole graph, it improves the solution through local search iteration. The lower the ice load degradation coefficient, the higher the edge weight of the corresponding line.

[0020] In the Benders decomposition method, the topology recovery binary decision master problem uses a graph neural network to generate a hot-start initial solution, and the graph neural network shares the spatial flow encoder parameters of the causal dual-flow fault tracing and localization model.

[0021] The power flow feasibility verification subproblem adopts the linearized DistFlow equation, with branch active power, reactive power and node voltage squares as state variables. After ignoring second-order minor quantities, a linear equation system is formed to complete the power flow verification.

[0022] The stochastic process model for ice shedding is based on the ice adhesion strength decay curve, which is established by laboratory measurement experiments of ice adhesion strength changing with temperature and time. Monte Carlo tree search enumerates the closing timing nodes corresponding to different waiting times in the decision tree and estimates the probability of re-failure through Monte Carlo stochastic simulation.

[0023] The training dataset for the causal dual-flow fault tracing and localization model is generated by synthesizing fault waveform data through electromagnetic transient simulation software, and then collecting actual historical fault waveform data. Samples are generated by supplementing rare fault types through variational autoencoders. The mixing ratio of synthetic data and measured data is determined by ablation experiments.

[0024] The training of the causal dual-flow fault tracing and localization model uses the fault localization cross-entropy loss, the weighted sum of the causal graph acyclic constraint penalty term and the physical consistency regularization term as the total loss function. The optimizer adopts the AdamW algorithm, and the initial learning rate is dynamically adjusted through learning rate warm-up and cosine annealing strategies.

[0025] The do-calculus differentiable approximation replaces the hard cutting of parent node edges with the introduction of a differentiable soft shielding function, making the calculation process of the intervention effect differentiable with respect to network parameters; the causal graph topology of the differentiable causal graph module of the causal dual-flow fault tracing and localization model is initialized by expert knowledge and the edge weights are automatically corrected during training using the NOTEARS algorithm.

[0026] In the grid search experiment, the search range of each weight coefficient is 0.01 to 10, with values ​​taken at logarithmic intervals, and a total of no less than 200 training iterations are performed; the segment length of the patch embedding encoding is selected by comparing the positioning accuracy corresponding to 8, 16, 32, and 64 sampling points; the temperature coefficient is determined between 0.1 and 1.0 through ablation experiments; the confidence threshold is determined through receiver operation characteristic curve analysis experiments with a false negative rate of no more than 2% as a constraint; and the acceptable risk threshold is adjusted in accordance with safety management regulations within the range of 2% to 10%.

[0027] This invention employs a physically embedded icing proxy model to predict icing thickness over time and converts the prediction results into an ice load degradation coefficient. This degradation coefficient is used as the edge weight correction amount of the power restoration path algorithm, enabling the path planning stage to automatically avoid heavily iced lines. This solves the technical problem of frequent re-tripping due to inaccurate line carrying capacity assessment after power restoration.

[0028] This invention deeply couples icing thickness prediction results with power restoration path planning. Lines with lower ice load degradation coefficients correspond to higher edge weights. The path optimization process automatically avoids lines with severely degraded carrying capacity while satisfying radial topology constraints and line current carrying capacity constraints, ensuring that the line's operating state after energization remains within the safe current carrying range under icing conditions. Furthermore, the joint decision-making framework of Monte Carlo tree search and the ice shedding stochastic process model further quantifies the risk of recurrence in the time dimension, providing a probabilistic safety guarantee for the selection of energization timing.

[0029] In summary, this invention solves the technical problem mentioned in the background art: the inability to dynamically integrate the time-series prediction results of ice thickness in the power restoration path planning under ice disaster scenarios leads to inaccurate assessment of line carrying capacity after power restoration and frequent re-tripping. Attached Figure Description

[0030] Figure 1 This is a flowchart of the method of the present invention.

[0031] Figure 2 This is a time-series prediction result of the icing thickness at various route nodes in the park.

[0032] Figure 3 A distribution diagram of the activation path and node causal contribution weights of the causal variable chain in the causal dual-flow fault tracing and localization model.

[0033] Figure 4 This is a diagram showing the results of the minimum cost power restoration route planning after an ice storm.

[0034] Figure 5 This is a distribution of the probability of recurrence under different waiting times.

[0035] Figure 6 A timeline diagram for coordinating the multi-source output of each isolated island. Detailed Implementation

[0036] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0037] like Figure 1 The diagram shown is a flowchart of an optimized power restoration method for industrial parks that considers ice storm icing prediction and multi-source coordination, provided by the present invention. This method includes the following steps:

[0038] S01. Perform matrix completion on the missing post-disaster measurements, construct a Bayesian network, and output the posterior probability distribution of the faulty segment.

[0039] S02. Using the posterior probability distribution of the fault section as input, construct a two-layer heterogeneous coupling model of the transportation network and the distribution network. Use a hierarchical decoupling iterative architecture to decompose the dual-network coordination problem into an upper-layer scheduling optimization sub-problem and a lower-layer power flow verification sub-problem.

[0040] S03. Offline pre-calculation of the Lyapunov stability domain boundary under typical island topology, training of hierarchical interpretable recovery decision model, and online invocation of the hierarchical interpretable recovery decision model to complete the real-time determination of island transient voltage stability.

[0041] S04. Call the adaptive large neighborhood search algorithm based on electrical distance weighting to solve the upper-level scheduling optimization subproblem and generate candidate recovery sequences;

[0042] S05. Call the simulated annealing cut algorithm based on Kirchhoff's thermodynamic analogy to perform global optimization on the candidate recovery sequence and output the optimal recovery scheme that satisfies the physical constraints.

[0043] S06. Pass the optimal recovery scheme to the lower-level power flow verification subproblem, call the hierarchical interpretable recovery decision model for feasibility verification. If the verification fails, return to the upper layer to re-iterate the cutting plane constraint until the convergence condition is met, and then output the final cooperative recovery sequence.

[0044] In S01, the matrix completion method is a graph signal processing measurement reconstruction method. This method utilizes the low-rank characteristic of the Laplace matrix of the distribution network topology to transform the missing measurement problem into a matrix completion problem. By minimizing the kernel norm, the missing current direction measurement is recovered. The low-rank characteristic means that the effective rank of the distribution network Laplace matrix is ​​much smaller than the matrix dimension, ensuring that the matrix completion still has theoretically guaranteed reconstruction accuracy when the measurement missing rate is less than 40%. The Bayesian network uses the reconstruction result of the graph signal processing measurement reconstruction method as the observed variable and the fault section disconnection state as the latent variable. It calculates the posterior probability distribution of the fault section through a variational inference algorithm. The posterior probability distribution of the fault section quantifies the topological uncertainty and transmits the topological uncertainty to the upper-level scheduling optimization subproblem in S02 for explicit processing in the form of robust constraints.

[0045] In the dual-layer heterogeneous coupling model described in S02, the transportation network uses vehicle segment flow as the state variable, while the distribution network uses node voltage and branch current as state variables. The two networks are bidirectionally coupled through the load demand of charging stations. The hierarchical decoupling iterative architecture introduces the alternating direction multiplier method to transfer dual variables between the upper-layer scheduling optimization subproblem and the lower-layer power flow verification subproblem, compressing the scale of each subproblem to a level that can be solved in parallel. The hierarchical decoupling iterative architecture also introduces graph sparsification technology to prune weakly coupled node pairs. The graph sparsification technology calculates the sensitivity coefficient matrix between each node pair and prunes node pairs whose absolute sensitivity coefficient is lower than the sparsification threshold. The sparsity threshold is determined by iterative experiments on multiple standard test systems, gradually tightening the threshold and recording the computational accuracy loss and computational time compression ratio. The maximum threshold corresponding to a computational accuracy loss of no more than 1% is taken as the determination.

[0046] In S03, the Lyapunov stability domain boundary is obtained through offline time-domain simulation of a typical island topology, with the simulation step size being... s to For each topology configuration, several operating points are sampled within the output disturbance range of the distributed power source. The Lyapunov function value corresponding to each operating point is calculated, and the critical Lyapunov function value corresponding to the stability domain boundary is used as the core content of the offline database. The hierarchical interpretable recovery decision model takes the network topology feature vector as input and the island transient voltage stability judgment result as output. When called online, the single judgment time is compressed to the millisecond level, replacing the time domain simulation to complete the real-time judgment. The network topology feature vector is formed by flattening the adjacency matrix of the current island topology and concatenating it with the output state vector of the distributed power source.

[0047] The specific process of the adaptive large neighborhood search algorithm based on electrical distance weighting described in S04 is as follows: First, the electrical distance matrix between all pairs of nodes in the distribution network is calculated. The electrical distance matrix is ​​constructed based on the difference between the diagonal and off-diagonal elements of the node impedance matrix to quantify the electrical topological physical proximity between nodes. The destruction operator preferentially selects clusters of mutually coupled nodes with smaller electrical distances as destruction targets. The repair operator adopts a greedy priority sorting based on DC power flow approximation and restores nodes in order of decreasing marginal benefit to the voltage support of the entire network after restoration. After each repair step, a lightweight traffic flow allocation model is called to verify the accessibility constraints of emergency repair vehicles. The adaptive weighting mechanism dynamically adjusts the selection probability of each operator based on the improvement in solution quality of each combination of destruction and repair operators in the past several iterations. The algorithm terminates when the computation time budget is exhausted. If the number of iterations with no improvement exceeds the termination threshold, a candidate recovery sequence is output. The computation time budget is set to 3 to 5 minutes based on the recovery decision time window requirements. The termination threshold with no improvement is determined by recording the algorithm convergence curves on different scale test systems and taking the average number of iterations corresponding to a solution quality improvement of less than 0.1%, plus 20% to 30%. The lightweight traffic flow assignment model is based on the Wardrop user equilibrium principle and uses the linearized approximation of the road segment impedance function to replace the original nonlinear function, transforming the traffic flow assignment problem into a linear programming problem. The adaptive large neighborhood search algorithm based on electrical distance weighting uses the posterior probability distribution of the fault segment output by S01 as the initial prior for the selection weight of the destruction operator, and prioritizes the destruction operation for the neighborhood of the fault segment with a higher posterior probability.

[0048] The technical benefits of the adaptive large neighborhood search algorithm based on electrical distance weighting are reflected in the following aspects: Traditional large neighborhood search algorithms randomly select variables to be disturbed by the destruction operator. In the distribution network scenario, a large number of random disturbances do not affect weakly coupled variables that have power flow, resulting in low search efficiency. After introducing the electrical distance matrix to define the destruction neighborhood in a physically heuristic way, each destruction operation is concentrated on the cluster of electrically strongly coupled nodes, ensuring that the disturbance operation has an effective impact on the objective function. The synergistic use of the greedy repair operator and the lightweight traffic flow allocation model ensures that each repair operation simultaneously satisfies the power flow constraints and the accessibility constraints of the repair vehicles, avoiding the ineffective computation of generating a large number of traffic infeasibility candidate solutions. The adaptive weighting mechanism enables the algorithm to automatically tilt towards combinations of operators with excellent historical performance during the search process, increasing the probability of finding high-quality feasible solutions within a limited computation time budget.

[0049] The specific process of the simulated annealing graph cut algorithm based on Kirchhoff's thermodynamic analogy described in S05 is as follows: The distribution network topology is mapped to a weighted graph, where edge weights are the line thermal stability margins and node weights are the recovery priorities weighted by the critical load coefficients. The system energy function is defined as the sum of weighted outage losses and topology switching costs. Candidate new states are generated by graph cut operations, randomly selecting an edge to perform a cut or closure to generate a new power supply path. The state transition probability adopts the improved Metropolis criterion, multiplying the standard Boltzmann factor by a physical correction term. This physical correction term uses the degree to which Kirchhoff's current law is satisfied at the cut point as the numerator, making the state transition with better current balance have a higher acceptance probability. The cooling strategy adopts adaptive geometric cooling, which is applied when multiple consecutive acceptances occur. When the rate of acceptance is below the lower limit threshold for cooling, the cooling rate is slowed down; both demand response constraints and traffic network constraints are added to the system energy function in the form of penalty functions, imposing a large energy penalty on switching operations where the demand response adjustment exceeds the upper and lower limits or where emergency repair vehicles cannot reach; the critical load factor is given in advance by the power supply authority according to the load type classification standard; the line thermal stability margin is calculated by the difference between the line's rated current carrying capacity and the current load; the lower limit threshold for cooling is determined by recording the solution quality under different combinations of initial temperature and cooling coefficient on multiple test systems, and taking the temperature value corresponding to the acceptance rate being below 5%; the demand response constraints use the upper limit of load adjustment, the lower limit of load adjustment, the adjustment response delay time, and the upper limit of cumulative adjustment times as parameters, and these parameters are obtained by statistical analysis of historical demand response project measured data.

[0050] The technical benefits of the simulated annealing graph cut algorithm based on Kirchhoff's thermodynamic analogy are reflected in the following aspects: The standard simulated annealing algorithm does not distinguish the physical rationality of state transitions in distribution network problems, leading to a large number of candidate states that violate Kirchhoff's current law being accepted but subsequently rejected, resulting in invalid searches. Introducing a physical correction term based on the degree of Kirchhoff's current law satisfaction increases the acceptance probability of physically more rational state transitions, causing the search trajectory to shift along the physically feasible direction in the solution space, reducing the proportion of invalid state transitions. The graph cut operation directly perturbs the topology, highly matching the physical operation form of distribution network topology switching, ensuring a one-to-one correspondence between the algorithm's search steps and actual recovery operations, and directly mapping the generated optimal recovery scheme to an executable sequence of switching operations. The introduction of the demand response penalty function and the traffic network penalty function allows the algorithm to simultaneously consider load-side regulation capacity and the accessibility of emergency repair vehicles during distribution network topology optimization, achieving integrated processing of the two-network collaborative constraints and demand response constraints.

[0051] The hierarchical interpretable recovery decision model is based on deep learning, with the following structure: the network is divided into two levels: a macro-policy layer and a micro-policy layer. The macro-policy layer adopts an Actor-Critic structure, with the Actor network having a fully connected layer as the backbone. The input is a global state encoding vector, which includes the topology encoding after flattening the entire network's topology adjacency matrix, the power outage degree encoding for each segment, and the traffic network accessibility encoding. The Actor network has three hidden layers, with each layer containing 256 to 512 neurons. The neurons are sparsely connected, with a sparsity rate between 0.3 and 0.7. This sparsity rate is verified by gridding on a standard testing system using the power recovery rate as an indicator. The search experiments determined that the Actor network outputs discrete actions for selecting the target sub-region in the current recovery phase; the Critic network is based on a graph neural network, with the node feature dimension set to 64 to 128 dimensions and the number of layers set to 2 to 4. Feature transmission between layers adopts a message passing mechanism, and the memory between layers is dynamically allocated according to the actual size of the feature matrix of each layer's nodes. The Critic network embeds the entire network state graph as a global value estimation scalar; the micro-policy layer uses a sequence decision network based on a Transformer encoder for the target sub-region selected by the macro-policy layer, with the number of attention heads in the Transformer encoder set to 4 to 8 and the number of encoder layers set to... The network is configured with 2 to 6 layers. Computational tasks between attention heads are executed in parallel via a CUDA stream allocation mechanism. Each CUDA stream corresponds to matrix multiplication operations for different attention heads. The number of CUDA thread blocks is adaptively set based on the size of the attention matrix. The input is the current state sequence of all switches within the target sub-region, and the output is the priority order and execution probability of switch operations within the target sub-region. The global state encoding vector is transferred between the macro-policy layer and the micro-policy layer via a shared memory buffer. The target sub-region index output by the macro-policy layer is written to the shared memory, and the micro-policy layer reads from the shared memory and constructs the local state sequence. The memory allocation is pre-calculated and fixed based on the maximum number of nodes in the target sub-region. The two layers communicate via an internal reward mechanism. The incentive mechanism connects the micro-policy layer, which feeds back the intrinsic reward to the macro-policy layer after completing the sub-objective. The macro-policy layer takes the external power recovery reward as the ultimate goal. An action masking mechanism with physical rule constraints is designed within the network. After the policy network outputs the action probability, a rule engine based on the radial constraints and thermal stability constraints of the distribution network is introduced to forcibly reset the action probability that violates the physical constraints to zero and then re-normalize it. The interpretability module extracts the attention weights of the Transformer encoder and outputs the information of the neighboring nodes that have the greatest impact on the current decision in the form of a weight matrix. The specific steps for establishing the training dataset of the hierarchical interpretable recovery decision model include: based on historical disaster records and Monte Carlo fault scenario generation methods, constructing a dataset containing no less than [number missing] [items missing]. A simulation sample library for each fault scenario is provided. Each sample includes the post-fault topology state, distributed power supply output state, traffic network access state, and corresponding optimal recovery action sequence labeling. The optimal recovery action sequence labeling is obtained by offline execution of an adaptive large neighborhood search algorithm based on electrical distance weighting and a simulated annealing graph cut algorithm based on Kirchhoff's thermodynamic analogy. The samples are divided into training, validation, and test sets in an 8:1:1 ratio. The training steps of the hierarchical interpretable recovery decision model specifically include: the macro-policy layer and micro-policy layer adopt a phased training strategy. First, the macro-policy layer is trained with fixed micro-policy layer parameters. After the macro-policy layer converges, the two layers are jointly trained. The optimizer adopts an adaptive moment estimation algorithm, and the initial learning rate is set to... to The intervals are dynamically adjusted by the network convergence state adjustment function. The batch size is selected between 32 and 128 based on the memory capacity. During the training process, the power recovery rate index is evaluated on the validation set every few rounds. The condition for early stopping is that the power recovery rate of the validation set does not increase continuously.

[0052] The hierarchical interpretable recovery decision model brings the following technical benefits to the entire scheme: The macro-policy layer decomposes the exponentially growing joint action space into hierarchical subspaces by encoding the global topology state and outputting the target sub-region selection. This allows the micro-policy layer to search only within the local switch operation space, reducing the search size from exponential to polynomial as the number of nodes increases. The action masking mechanism of physical rule constraints directly excludes inactive actions at the policy output, ensuring that the exploration is always limited to the physically feasible domain, significantly reducing the convergence sluggishness caused by invalid exploration. The attention weight matrix output by the interpretability module enables the dispatcher to trace the basis nodes of the model's decision, improving the credibility of the automatic decision-making scheme in actual power system operation scenarios. The overall architecture decomposes the long-term recovery task into the gradual completion of several sub-objectives through an inherent reward mechanism, alleviating the gradient vanishing problem caused by reward sparsity in long-term reinforcement learning.

[0053] The network convergence state adjustment function is used to dynamically adjust the learning rate parameter of the hierarchical interpretable recovery decision model; the network convergence state adjustment function uses the variance of the power collection recovery rate in the most recent training rounds as the basis for its calculation. The mean training loss over the most recent training rounds and the ratio of the current training rounds to the total training rounds. Using the input, calculate the convergence exponent. The calculation formula is: ,in (unit ), (Unit is) / loss unit) (Dimensionless) is the dimensional compatibility coefficient, which is determined through preliminary experiments on a standard testing system to ensure... If the value falls within the preset range, adjust the values ​​of the three coefficients accordingly. , The dimensionless quantity; when At that time, the learning rate is multiplied by the contraction factor. To suppress oscillations; when When the learning rate remains constant; when At that time, the learning rate is multiplied by the expansion factor. To accelerate escape from the flat area; the and By repeatedly training and recording the convergence curve shape on multiple test systems of different scales, the corresponding value at the time of oscillation is taken. The average value is increased by 10% to 20% as Take the corresponding gradient vanishing event. The average value is adjusted downwards by 10% to 20% as... .

[0054] Optionally, the present invention also provides a computer-based approach to form a demand-response-integrated transportation network and distribution network collaborative recovery optimization system, wherein the computer is equipped with a readable storage medium storing program instructions, and the program instructions execute the above-described method when the computer is run.

[0055] The specific implementation of step S01 is as follows: After a disaster, damage to the distribution network leads to the failure of some measurement terminals, resulting in missing current direction measurements and topological uncertainty in fault location. This step employs a graph signal processing measurement reconstruction method. Based on the Laplace matrix corresponding to the distribution network topology, and utilizing its low-rank characteristic (effective rank much smaller than matrix dimension), the problem of reconstructing missing measurements is transformed into a matrix completion problem. Minimizing the nuclear norm is used as the objective function, and the missing current direction measurements are recovered through iterative solving using a convex optimization solver. When the measurement loss rate is below 40%, this method has theoretically guaranteed reconstruction accuracy. Using the reconstructed measurements as observed variables and the on / off state of each fault section as latent variables, a Bayesian network is constructed. A variational inference algorithm is used to approximate the posterior distribution of the latent variables, outputting the posterior probability distribution of each section. This posterior probability distribution is robustly constrained and passed to subsequent upper-level scheduling optimization sub-problems, allowing topological uncertainty to be explicitly handled within the optimization framework, avoiding the scheme errors caused by deterministic assumptions about fault sections in traditional methods.

[0056] The specific implementation of step S02 is as follows: Using the traffic flow of vehicles on road segments and the voltage and branch current of distribution network nodes as their respective state variables, the two networks establish a bidirectional coupling relationship through the load demand of charging stations, forming a two-layer heterogeneous coupling model. To avoid scale expansion caused by joint solutions, a layered decoupling iterative architecture is introduced. The alternating direction multiplier method is used to alternately transfer dual variables between the upper-layer scheduling optimization subproblem and the lower-layer power flow verification subproblem. The two subproblems are iterated alternately until the dual residuals satisfy the convergence criterion. Simultaneously, graph sparsification technology is introduced to calculate the sensitivity coefficient matrix between each pair of nodes, pruning node pairs whose absolute sensitivity coefficient values ​​are lower than the sparsification threshold. The sparsification threshold is determined by gradually tightening the threshold on multiple standard test systems and recording the computational accuracy loss and computational time compression ratio. The maximum threshold corresponding to a computational accuracy loss of no more than 1% is taken, thereby maximizing the compression of the subproblem scale while ensuring accuracy, enabling each layer of subproblems to be solved in parallel.

[0057] The specific implementation of step S03 is as follows: If the determination of islanded transient voltage stability relies on real-time time-domain simulation, the time consumed by a single evaluation far exceeds the recovery decision time window requirement. Therefore, this step adopts a combination of offline pre-calculation and online fast invocation. In the offline stage, time-domain simulation is performed on a typical islanded topology, and the simulation step size is taken as... s to For each topology configuration, several operating points are sampled within the output disturbance range of the distributed generation source. The Lyapunov function value corresponding to each operating point is calculated, and an offline database is constructed using the critical Lyapunov function values ​​corresponding to the stability domain boundary. Based on the offline database, a hierarchical interpretable recovery decision model is trained. The input is a network topology feature vector formed by concatenating the flattened adjacency matrix of the current islanded topology with the output state vector of the distributed generation source. The output is the islanded transient voltage stability determination result. When called online, the single determination time is compressed to the millisecond level, replacing time-domain simulation to complete real-time determination and meeting the recovery decision time window requirements.

[0058] The specific implementation of step S04 is as follows: This step calls an adaptive large neighborhood search algorithm based on electrical distance weighting to solve the upper-level scheduling optimization sub-problem. First, an electrical distance matrix is ​​constructed based on the difference between the diagonal and off-diagonal elements of the node impedance matrix to quantify the electrical topological physical proximity between each pair of nodes. The destruction operator prioritizes strongly coupled node clusters with smaller electrical distances as disturbance targets to ensure that each destruction operation has an effective impact on the objective function. The repair operator adopts a greedy priority sorting based on DC power flow approximation, and restores nodes in order of decreasing marginal benefit to the overall network voltage support after restoration. After each repair step, a lightweight traffic flow allocation model is called to verify the accessibility constraints of emergency repair vehicles. The lightweight traffic flow allocation model is based on the Wardrop user equilibrium principle and uses a linearized approximation of the road segment impedance function to transform the traffic flow allocation problem into a linear programming problem. The adaptive weighting mechanism dynamically adjusts the selection probability according to the improvement of the solution quality of each operator combination in historical iterations. The algorithm terminates when the computation time budget is exhausted (3 min to 5 min) or the number of consecutive iterations without improvement exceeds the consecutive no-improvement termination threshold, and outputs a candidate restoration sequence.

[0059] The specific implementation of step S05 is as follows: This step maps the distribution network topology to a weighted graph. The edge weights are the line thermal stability margins (calculated from the difference between the line's rated current carrying capacity and the current load), and the node weights are the recovery priorities weighted by the critical load coefficient. The system energy function is defined as the sum of the weighted outage loss and the topology switching cost. Demand response constraints and traffic network constraints are added in the form of penalty functions, imposing a large energy penalty on operations that violate the constraints. Candidate new states are generated by graph cut operations, randomly selecting an edge to perform a cut or closure to generate a new power supply path. The state transition probability adopts the improved Metropolis criterion, multiplying the standard Boltzmann factor by a physical correction term. The physical correction term uses the degree of satisfaction of Kirchhoff's current law at the cut point as the numerator, making the state transition with better current balance have a higher acceptance probability. The cooling strategy adopts adaptive geometric cooling. When the acceptance rate is lower than the cooling lower limit threshold for multiple consecutive times (determined by recording the solution quality under different temperature combinations on multiple test systems and taking the temperature value corresponding to the acceptance rate being lower than 5%), the cooling rate is slowed down, and finally the optimal recovery scheme that satisfies the physical constraints is output.

[0060] The specific implementation of step S06 is as follows: The optimal recovery scheme output in S05 is passed to the lower-level power flow verification subproblem, and the hierarchical interpretable recovery decision model is invoked for feasibility verification. The feasibility verification checks whether the optimal recovery scheme meets physical constraints such as node voltage, branch thermal stability, and islanded transient voltage stability. If the verification fails, cut plane constraints are generated based on the reason for the violation, and these cut plane constraints are returned to the upper-level scheduling optimization subproblem as additional constraints, triggering a new round of iteration. The generation of cut plane constraints is based on the type and degree of constraint violation by the current scheme, mapped to the constraint space of the upper-level problem through dual variables. The above iterative process continues until the lower-level power flow verification passes or the number of iterations reaches the upper limit of the convergence condition, outputting the final cooperative recovery sequence.

[0061] It should be noted that the key technologies of this invention include: a graph signal processing measurement and reconstruction method based on the low-rank characteristics of the Laplace matrix, the physical basis of which is that the sparse connection of the distribution network topology makes the graph signal have low-frequency concentrated characteristics in the frequency domain, so that uncertainty can be quantified rather than enumerated; a hierarchical decoupling iterative architecture relaxes the coupling constraints of the two networks into dual variables that are alternately passed between the upper and lower layers through the alternating direction multiplier method, so that the solution complexity of the joint optimization problem is reduced from the product of the number of nodes in the two networks to the order of magnitude of the sum of the number of nodes in the two networks; a hierarchical interpretable recovery decision model limits the exploration process to the polynomial-level subspace of the physically feasible domain through macro-micro two-level action space decomposition and action masking mechanism; and a simulated annealing graph cut algorithm based on Kirchhoff's thermodynamic analogy embeds the degree of satisfaction of physical laws into the state transition probability, so that the search trajectory is offset along the physically feasible direction. The four technologies mentioned above support each other: measurement reconstruction provides a confident topological input for subsequent optimization; hierarchical decoupling compresses the high-dimensional joint problem into parallelizable subproblems; hierarchical interpretable models further compress the action space and provide fast feasibility determination; graph cut algorithm completes efficient global optimization in the compressed space; and the synergistic effect of the four technologies enables multi-source heterogeneous constraints to be processed in an integrated manner within a finite time budget.

[0062] It should be noted that in islanded operation scenarios with high penetration of distributed power sources, there is a technical challenge in determining transient voltage stability where real-time performance and accuracy are mutually exclusive. When the output of a distributed power source within the island experiences a disturbance, the system voltage trajectory may cross the stability boundary within tens of milliseconds. While traditional time-domain simulation methods can accurately track the voltage trajectory, a single simulation typically takes minutes, far exceeding the requirements of the recovery decision-making time window, making it impossible to invoke in real-time within the online decision-making process. The reason for this technical problem is that transient stability is essentially a high-dimensional nonlinear determination problem dependent on the initial state of the system and the amplitude of the disturbance. Its analytical solution does not exist; it must be solved step by step using numerical integration methods to solve a system of differential equations. The computational load is proportional to the simulation time, and it cannot be replaced by simple algebraic approximations. Common solutions to this technical problem include using a direct method based on energy functions or a simplified equivalent model for approximate determination. However, the direct method is difficult to construct energy functions in islanded systems containing inverter-type distributed power sources, and the simplified equivalent model suffers significant accuracy loss in multi-machine islanded scenarios. Neither method can simultaneously guarantee millisecond-level determination speed and sufficient determination accuracy. This invention solves this technical problem by separating offline pre-computation from online rapid inference. In the offline phase, for typical islanded topology configurations, operating points are densely sampled within the output disturbance range of distributed power sources. The Lyapunov function values ​​corresponding to each operating point are calculated through precise time-domain simulation. An offline database covering typical scenarios is constructed using the critical values ​​corresponding to the stability domain boundaries, and the high-precision simulation results are solidified in the database in a structured form. In the online phase, the hierarchical interpretable recovery decision model takes the current islanded topology feature vector as input and directly outputs the stability determination result through forward inference. The inference process only involves matrix multiplication and activation function operations, compressing the single determination time to the millisecond level, fully meeting the real-time requirements of online decision-making. Simultaneously, since the model training data comes from precise time-domain simulation, its determination accuracy inherits the physical consistency of the offline simulation, overcoming the shortcomings of insufficient accuracy in simplified equivalent models.

[0063] Specifically, the principle of this invention is:

[0064] The present invention can solve the above-mentioned core technical problems. The fundamental reason is that it reduces the dimensionality of the originally indivisible high-dimensional heterogeneous constraint optimization problem layer by layer through a multi-level decomposition mechanism, so that the size of each sub-problem is always within the range that can be solved within the time budget.

[0065] First, topological uncertainty is the initial source of constraint dimensionality explosion. Post-disaster measurement gaps leave the state of faulty sections unknown; traditional methods require enumerating all possible topologies, leading to an exponential increase in the number of scenarios. This invention recovers missing measurements using graph signal processing measurement reconstruction, and then transforms the topological uncertainty into a posterior probability distribution via a Bayesian network. This distribution is then robustly passed to upper-level optimization, incorporating uncertainty into a deterministic framework and fundamentally avoiding the dimensionality explosion caused by scenario enumeration.

[0066] Secondly, the two-layer heterogeneous coupling model decouples the coupling constraints of the transportation network and the distribution network into two sub-problems that are solved alternately by passing dual variables between the upper and lower layers using the alternating direction multiplier method. Furthermore, graph sparsity techniques are used to prune weakly coupled node pairs, further compressing the effective size of each layer's sub-problems. This transforms the original dual-network coordination problem, which required searching in the joint space, into two single-network problems that can be processed in parallel, significantly reducing computational complexity.

[0067] Furthermore, the hierarchical interpretable recovery decision model decomposes the exponential joint action space into a polynomial-level local search space through hierarchical action space decomposition of macro-policy layers and micro-policy layers. The action masking mechanism of physical rule constraints directly sets the probability of actions that violate radial constraints and thermal stability constraints to zero at the policy output, ensuring that the exploration process of reinforcement learning is always confined to the physically feasible region, avoiding a large amount of ineffective exploration that consumes time budget.

[0068] Finally, the electrically distance-weighted adaptive large neighborhood search algorithm, through physical heuristics, breaks the neighborhood definition, concentrating each perturbation on electrically strongly coupled node clusters, ensuring the effective impact of perturbation operations on the objective function. The simulated annealing graph cut algorithm, based on Kirchhoff's thermodynamic analogy, introduces a physical correction term centered on the degree of satisfaction of Kirchhoff's current law into the state transition probability, causing the search trajectory to shift along the physically feasible direction and reducing the proportion of invalid state transitions. The two algorithms work together, converging with a higher probability to the optimal recovery scheme that satisfies all constraints within a finite time budget. These mechanisms support each other, forming the technical foundation for integrated collaborative optimization capabilities across four levels: uncertainty handling, problem decomposition, space compression, and search guidance.

[0069] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0070] The specific implementation method of step S01 is as follows.

[0071] Post-disaster distribution networks suffer from measurement gaps. Utilizing the low-rank property of the Laplace matrix in the distribution network topology, the problem of reconstructing missing measurements is modeled as a matrix completion problem. Let the measurement matrix be... ,in To measure the number of nodes, The number of sampling points within the time window, denoted as . The set of observed elements is denoted as . The problem of minimizing the nuclear norm is formulated as follows:

[0072]

[0073] In the formula, To reconstruct the complete measurement matrix, for The nuclear norm is the sum of all singular values. Here is the regularization weight coefficient, with an empirical value of [value missing]. to , The node index and time index of the observed element. These are the measured values ​​at the corresponding locations. These are the reconstructed values ​​at the corresponding locations. After reconstruction, the direction of the reconstructed current is measured as the observed variable. The fault section's open state is used as a hidden variable. A Bayesian network is constructed, and the posterior probability distribution is calculated using a variational inference algorithm. And it is passed on to the upper-level scheduling optimization subproblem in the form of robust constraints.

[0074] The specific implementation method of step S02 is as follows.

[0075] The transportation network is based on vehicle traffic flow on road segments. For state variables, Indicates road segment Vehicle traffic flow on the distribution network is measured by node voltage. With branch current Let be the state variable, where For nodes voltage amplitude, branch road The current amplitude, the two networks through the charging station load Forming a two-way coupling, Let the active power of the charging station be denoted as . The alternating direction multiplier method is introduced to decompose the cooperative problem into an upper-level scheduling optimization subproblem and a lower-level power flow verification subproblem, with dual variables passed between the two layers. , Let Lagrange multiplier vectors be the coupling constraints. The iterative update formula is expressed as follows:

[0076]

[0077] In the formula, For the first The dual variable vector at the next iteration The index for the iteration number. To augment the Lagrange penalty coefficient, the empirical value is... to , and These are the coupling constraint coefficient matrices for the upper and lower layers, respectively. For the first The upper-level decision variable vector at the next iteration specifically includes the state variables of each branch switch and the demand response adjustment. For the first The lower-level state variable vector at the next iteration specifically includes the voltages of each node. With branch current Graph sparsification techniques calculate the sensitivity coefficient matrix. Pruning weakly coupled node pairs, matrix elements Represents a node Power disturbances on nodes The influence of voltage, the cutting condition is ,in The sparsity threshold is determined by taking the maximum value corresponding to a computational accuracy loss of no more than 1% through iterative experiments.

[0078] The specific implementation method of step S03 is as follows.

[0079] Offline time-domain simulation of a typical island topology was performed, with the simulation step size set to... s to For each topology configuration, sample several operating points within the output disturbance range of the distributed power source, and calculate the Lyapunov function value corresponding to each operating point. ,in For the state vector of the islanded system, Let the dimension of the state variables be . The statement is as follows:

[0080]

[0081] In the formula, The matrix is ​​positive definite and symmetric. It was obtained by solving the Lyapunov equations for each topological configuration through offline time-domain simulation, and the critical value was determined. The corresponding state when the system changes from stable to unstable. value, The core content constituting the offline database. Network topology feature vectors of a hierarchical interpretable and recoverable decision model. It is composed of the flattened adjacency matrix vector and the distributed power source output state vector, and is described as follows:

[0082]

[0083] In the formula, This is the current island topology adjacency matrix. The number of nodes within the isolated island. For matrix flattening operators, This represents the output state vector of the distributed power source. For the number of distributed power sources, This is the concatenated feature vector. When calling the model online, it will... Inputting into a hierarchical interpretable recovery decision model compresses the single decision time to the millisecond level, replacing time-domain simulation to complete real-time determination of transient voltage stability.

[0084] The specific implementation method of step S04 is as follows.

[0085] First, construct the electrical distance matrix based on the node impedance matrix. The elements are described as follows:

[0086]

[0087] In the formula, The node impedance matrix is ​​the first Line number Column elements, and They are nodes With nodes Self-impedance, The smaller the value, the more likely it is to be a node. With nodes The tighter the electrical coupling, the more likely the operator will fail. To prioritize the perturbation of strongly coupled node clusters for weighting, the initial prior is the posterior probability distribution output in step S01. Assignment. The marginal benefit of the repair operator on the overall network voltage support after node-by-node recovery. The descending order is restored sequentially. For nodes The sum of voltage increases across all network nodes after restoration is obtained through DC power flow approximation. The lightweight traffic flow assignment model is based on the Wardrop user equalization principle and uses a linearized approximation of the road segment impedance function. The linearized road segment impedance is expressed as follows:

[0088]

[0089] In the formula, For road section With traffic Travel time at that time For road section Free-flow travel time, For road section The linearized slope coefficient is obtained by fitting the road segment capacity parameters. The adaptive weighting mechanism is based on the past performance of each operator combination. Improvement in solution quality in the next iteration Dynamically adjust the selection probability. The historical window length is empirically between 10 and 30 iterations. The operator selection probability update formula is as follows:

[0090]

[0091] In the formula, For the first Operator combination in the next iteration The probability of choosing, Operator Combination In the Improvement of the objective function in the next iteration Take the initial objective function value as the reference value of the objective function. Let be the traversal index for all operator combinations. and All are dimensionless quantities, making The value is a dimensionless probability. The algorithm terminates when the computation time budget is exhausted (3 to 5 minutes) or the number of consecutive iterations without improvement exceeds the termination threshold.

[0092] The specific implementation method of step S05 is as follows.

[0093] Mapping the distribution network topology to a weighted graph, system energy function Defined as the sum of weighted outage losses and topology switching costs, and incorporating demand response penalty functions and traffic network penalty functions. The dimensions are unified as The expression is as follows:

[0094]

[0095] In the formula, For load The critical load factor is dimensionless and its value range is [value range missing]. The value is predetermined by the power supply authority based on the load type classification standard, with important loads taking a value close to 1. For load The power supply status indicator variable has a value of 1 indicating that power has been restored and a value of 0 indicating that power has not been restored. For load The active power, with dimensions of , branch road The topology switching cost, with dimensions of The value is calculated based on the equivalent power outage loss corresponding to a single switching operation. For load Demand response adjustment quantity, with dimensions of , For load The upper bound of the demand response adjustment is given by the dimension of , The coefficients of the demand response penalty function are dimensionless, and their empirical value is... to , For road section The accessibility indicator variable for emergency repair vehicles, with a value of 1 indicating accessibility and 0 indicating inaccessibility, is dimensionless. Let be the coefficients of the traffic network penalty function, with dimensions . Experience value to , For load set, For branch road collection, This is the set of road sections requiring emergency repairs. The state transition probabilities adopt the improved Metropolis criterion, as follows:

[0096]

[0097] In the formula, Let be the energy difference between the candidate new state and the current state, with dimensions . , The current annealing temperature has dimensions of . , making the exponential term Dimensionless This is a quantification of the degree to which the candidate new state satisfies Kirchhoff's current law at the cut point. The ratio is a quantification of the degree to which the current state satisfies Kirchhoff's current law at the cutting point. Both have the same dimensions. As a dimensionless physical correction term, it increases the acceptance probability of states with better current balance. Let be the dimensionless acceptance probability. The adaptive geometric cooling strategy is described as follows:

[0098]

[0099] In the formula, For the first The annealing temperature of the step, For the first Average acceptance rate of steps, dimensionless. The lower limit threshold for cooling is dimensionless and is determined by the temperature at which the acceptance rate is below 5%. To mitigate the cooling coefficient, dimensionless, This is the normal cooling coefficient, dimensionless.

[0100] The specific implementation method of step S06 is as follows.

[0101] The optimal recovery scheme output from step S05 is passed to the lower-level power flow verification subproblem, and the hierarchical interpretable recovery decision model is invoked for feasibility verification. If the verification fails, cutting plane constraints are generated. And return to the upper level to iterate again, where The constraint gradient vector is calculated from the constraint violation direction when the power flow check fails. To constrain the right-hand side scalar, the boundary values ​​that currently violate the constraints are determined until the convergence conditions are met, at which point the final cooperative recovery sequence is output.

[0102] Hierarchical interpretable recovery decision model network convergence state adjustment function to verify the variance of power recovery rate Mean training loss and training progress ratio Using the input as input, calculate the convergence state exponent, as described below:

[0103]

[0104] In the formula, This is the dimensionless compatibility factor, with dimensions of... , This is the dimensional compatibility factor, with dimensions of 1 / loss unit. This is a dimensionless coefficient of conformity. The variance of the power collection recovery rate over the most recent training rounds is defined by the dimension of... , This represents the average training loss over the most recent training rounds, measured in units of loss. This represents the ratio of the current training round to the total number of training rounds. Dimensionless This is the dimensionless convergent state exponent. When... At that time, the learning rate is multiplied by the contraction factor. ;when When the learning rate remains constant; when At that time, the learning rate is multiplied by the expansion factor. , and By repeatedly training and recording the convergence curve, the corresponding value at the occurrence of oscillations is taken. The average value is increased by 10% to 20% as Take the corresponding gradient vanishing event. The average value is reduced by 10% to 20% as The action masking mechanism outputs an action probability vector in the policy network. Then, the probabilities of actions that violate the radial constraints and thermal stability constraints are forcibly set to zero and renormalized, as follows:

[0105]

[0106] In the formula, For action The original output probability is dimensionless. For action Feasibility mask, For the action space, Used as the traversal index for the action space. The normalized probability of the action is dimensionless.

[0107] To better understand and implement this invention, the following is a specific application scenario of this invention, Example 2:

[0108] To verify the effectiveness of this invention, technicians constructed a test environment, selecting a typical power distribution network in a certain area as the test object. The test system included 118 distribution nodes, 153 feeder branches, and 34 distributed power source access points. The regional road network comprised 86 road network nodes and 124 road sections. Twelve charging stations were set up, covering major commercial and residential areas. The test scenario simulated the power grid damage state after a typhoon: 19 feeders failed, with measurement terminals simultaneously failing in 7 of the fault sections, resulting in a measurement loss rate of 37%; 11 roads in the road network were damaged, causing localized road network interruptions; and there were 8 distributed power source islands in the area, with a total installed capacity of 42 MW.

[0109] In step S01, technicians first used graph signal processing measurement reconstruction methods to complete the matrix of the seven missing current direction measurements. Based on the low-rank characteristics of the 118-node Laplace matrix of the distribution network, the measurement reconstruction was completed by minimizing the nuclear norm. Figure 2 The figure shows a comparison of the posterior probabilities of faults in each section before and after measurement reconstruction. Subsequently, a Bayesian network was constructed using the reconstructed measurements as the observed variables, and the posterior probability distributions of the 19 fault sections were output through variational inference algorithms. The posterior probability distributions of each section are shown in Table 1.

[0110] Table 1. Posterior probability distribution of typical fault sections

[0111]

[0112] In step S02, technicians constructed a two-layer heterogeneous coupling model of the transportation network and the distribution network. Using the load demand of charging stations as the bidirectional coupling interface, they introduced the alternating direction multiplier method for hierarchical decoupling and performed graph sparsification on the 118-node distribution network. After pruning node pairs whose absolute values ​​of sensitivity coefficients were lower than the sparsification threshold, the number of effective coupled node pairs was reduced from the original 6903 pairs to 874 pairs, and the size of the subproblem was significantly reduced.

[0113] In step S03, technicians conducted offline analysis of eight typical island topologies within the region, to... s is the simulation step size for time-domain simulation. 2000 operating points are sampled for each configuration, constructing a time-domain simulation containing... An offline database of samples was used as the basis for training a hierarchical interpretable recovery decision model, reducing the online single-test stability determination time to less than 3 ms, thus meeting real-time requirements.

[0114] In step S04, technicians invoked an adaptive large neighborhood search algorithm based on electrical distance weighting. Using the posterior probability distribution of the faulty sections as the initial prior for the destruction operator, the algorithm prioritized destruction operations on the neighborhoods of high-posterior-probability sections such as L-07 and L-13. After a 3-minute time budget, the algorithm converged, generating a set of candidate recovery sequences, such as... Figure 3 The figure shows the convergence curve of the objective function value as a function of the number of iterations during the search process. The adaptive weight adjustment mechanism enables the algorithm to focus on historical high-quality operator combinations in the later iterations, and the convergence speed is significantly faster than that of the control group with fixed weights.

[0115] In step S05, technicians invoke the simulated annealing graph cut algorithm based on Kirchhoff's thermodynamic analogy, using a weighted graph of the distribution network topology as input. The system energy function includes a weighted outage loss term, a topology switching cost term, a demand response penalty function term, and a traffic network penalty function term. The physical correction term uses the degree of Kirchhoff's current law satisfaction as the numerator, prioritizing the search trajectory to access physically feasible states. The lower limit threshold for cooling is determined based on pre-experiments, corresponding to a cooling rate of approximately 5%, slowing down the cooling rate, and finally outputting the optimal recovery scheme. In the optimal recovery scheme, the demand response adjustment amount, the switching operation sequence, and the repair vehicle path are shown in Table 2 (some typical operations).

[0116] Table 2 Typical Operation Sequence of Optimal Recovery Scheme

[0117]

[0118] In step S06, the optimal recovery scheme is passed to the lower-level power flow verification sub-problem. The hierarchical interpretable recovery decision model performs feasibility checks on all operation sequences. The initial check finds that the current in branch L-35 of section L-35 exceeds the thermal stability margin after SW-35 closure, resulting in a failed check. After generating cut-plane constraints, the process returns to the upper level for a second iteration. After the second iteration, the operation sequence is adjusted and the demand response adjustment is increased. All physical constraints pass the check, and the final cooperative recovery sequence is output. Figure 4 The figure shows the curve of the overall network power recovery rate as a function of the number of operation steps corresponding to the final collaborative recovery sequence. It can be seen that the hierarchical iteration mechanism enables the recovery rate to increase rapidly in the early operation steps and tend to converge to a stable level in the later stages.

[0119] Compared with traditional methods, the advancements of this invention are reflected in the following aspects: Traditional methods replace topological uncertainties with deterministic assumptions, leading to high sensitivity of optimization results to missing measurements; this invention, through robust constraints of Bayesian posterior probabilities, proactively handles topological uncertainties within the optimization framework in a probabilistic quantification manner, fundamentally eliminating the physical root cause of scheme deviations due to missing measurements. Traditional single-layer optimization suffers from exponentially expanding solution size under the superposition of multi-source heterogeneous constraints; this invention, through a three-layer compression mechanism of two-layer hierarchical decoupling, graph sparsification, and hierarchical interpretable models, progressively reduces the joint search space to a size solvable within the time budget, making integrated collaborative optimization possible within the actual decision-making time window. Traditional simulated annealing algorithms do not distinguish between physical rationality and generate a large number of invalid state transitions; this invention, by introducing Kirchhoff's physical correction term, automatically shifts the search trajectory towards the physically feasible direction, significantly improving the effective information density of each iteration, thereby achieving a higher quality feasible solution with fewer iterations.

[0120] It should be noted that the variables involved in this invention are explained in detail in Tables 3 and 4.

[0121] Table 3. Variable Explanation Table (Part 1)

[0122]

[0123] Table 4. Variable Explanation Table (Part Two)

[0124]

[0125] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for optimizing power restoration in industrial parks that considers ice accretion prediction and multi-source coordination, characterized in that, Includes the following steps: Collect conductor surface temperature, ambient temperature, wind speed, liquid water content and median volume diameter, input them into the physical embedded icing proxy model, and output the time series prediction results of icing thickness and ice load degradation coefficient. The protection device action time stamp sequence, fault waveform electrical quantities, distribution network topology adjacency matrix, and node icing degree characteristics are input into the causal dual-current fault tracing and localization model, and the fault location set and fault root cause identification results are output. Using the ice load degradation coefficient as the edge weight correction, the minimum cost power restoration path algorithm after ice disaster is called to generate the optimal power restoration path set that satisfies the radial topology constraint and the line current carrying capacity constraint. Taking the optimal power restoration path set as input, the Benders decomposition method is used to divide the power restoration optimization into the main problem of topology restoration binary decision and the sub-problem of power flow feasibility verification, and outputs feasible power restoration schemes; Using the ice adhesion strength decay curve and current meteorological conditions as input, a joint decision-making framework is constructed based on Monte Carlo tree search and ice shedding stochastic process model, and the probability of re-failure under different waiting times is output. Based on feasible power restoration plans and the probability of recurrence, the optimal timing for closing the circuit is determined within an acceptable risk threshold. The output of distributed photovoltaic, energy storage and diesel generators is coordinated, and the closing command is executed in order of priority to complete the power restoration of multiple isolated areas in the park.

2. The method for optimizing power restoration in industrial parks considering ice storm icing prediction and multi-source synergy as described in claim 1, characterized in that, The physical embedded icing proxy model is based on deep learning. It uses the partial differential relation of the Makkonen equation as the physical loss term to embed into the neural network training process, and uses the neural network to replace numerical iteration to solve the transient thermal balance equation.

3. The method for optimizing power restoration in industrial parks considering ice storm icing prediction and multi-source coordination as described in claim 2, characterized in that, The physical loss term is composed of a weighted sum of the residuals of the heat conduction equation, the convection heat transfer equation, and the latent heat conservation equation for phase change. The weighting coefficients are determined through a grid search experiment on the validation set that jointly minimizes the prediction error and the physical violation.

4. The optimized method for restoring power supply to industrial parks considering ice storm icing prediction and multi-source synergy as described in claim 3, is characterized in that, The ice load degradation coefficient is defined as the ratio of the actual current carrying capacity to the rated current carrying capacity of the line under the current icing conditions. The value ranges from 0 to 1 and is calculated by substituting the time series prediction results of icing thickness into the line heat balance equation.

5. The optimized method for restoring power supply to industrial parks considering ice storm icing prediction and multi-source coordination as described in claim 4, characterized in that, The input of the causal dual-flow fault tracing and localization model includes a temporal flow and a spatial flow. The temporal flow is segmented and encoded using a patch embedding encoding method, while the spatial flow is encoded using a graph Transformer. The temporal flow and spatial flow are fused in the middle layer of the model through a cross-attention mechanism.

6. The optimized method for restoring power supply to industrial parks considering ice storm icing prediction and multi-source synergy as described in claim 5, is characterized in that, The causal dual-flow fault tracing and localization model is connected to a differentiable causal graph module based on a structural causal model after the cross-attention mechanism layer. The differentiable causal graph module maps the features extracted by the network to a causal variable chain consisting of icing thickness, mechanical stress, wire breakage and protection action. The intervention effect is calculated by do-calculus differentiable approximation to identify the root cause of the fault.

7. The method for optimizing power restoration in industrial parks considering ice storm icing prediction and multi-source synergy as described in claim 6, characterized in that, The causal dual-flow fault tracing and localization model is internally set with an iterative refinement loop. In each loop, the fault probability vector corresponding to the current output fault location set is scaled by a temperature coefficient to form a soft label mask, which is then superimposed on the feature vector of the input node and the inference is re-executed. The loop is terminated early when the confidence level meets the confidence level threshold.

8. The method for optimizing power restoration in industrial parks considering ice storm icing prediction and multi-source synergy as described in claim 7, characterized in that, The proposed minimum cost power restoration path algorithm after an ice storm maps the power restoration problem to a constrained Steiner tree problem. It uses the Dreyfus-Wagner dynamic programming algorithm to solve the problem accurately on the key subgraph. After using the minimum spanning tree approximation as the initial solution for the entire graph, it improves the solution through local search iteratively. The lower the ice load degradation coefficient, the higher the edge weight of the corresponding line.

9. The optimized method for restoring power supply to industrial parks considering ice storm icing prediction and multi-source coordination as described in claim 8, characterized in that, In the Benders decomposition method, the topology recovery binary decision master problem uses a graph neural network to generate a hot-start initial solution, and the graph neural network shares the spatial flow encoder parameters of the causal dual-flow fault tracing and localization model.

10. The method for optimizing power restoration in industrial parks considering ice storm icing prediction and multi-source synergy as described in claim 9, characterized in that, The power flow feasibility verification subproblem adopts the linearized DistFlow equation, with branch active power, reactive power and node voltage squares as state variables. After ignoring second-order minor quantities, a linear equation system is formed to complete the power flow verification.

Images