Multi-section power transmission limit calculation method based on multi-target migratable reinforcement learning

Abstract

The invention relates to the technical field of electrical engineering, and particularly discloses a multi-section power transmission limit calculation method based on multi-target migratable reinforcement learning, and the method comprises the following steps: considering a transient power angle and voltage stability constraint, and building a reverse search model for power transmission limit calculation; decomposing a multi-section coupling quota calculation problem of the reverse search model into sub-problem models by adopting a decomposition strategy; different sub-problem models are converted into a Markov decision process, and solving is carried out in combination with a neighborhood-based parameter migration strategy and an improved DDPG algorithm; and the end-to-end output of the Pareto frontier of the multi-section coupling quota is realized by training the model. Compared with a traditional multi-target meta-heuristic algorithm, the method has the advantage of better convergence and diversity.

Description

Technical Field

[0001] The present invention relates to the field of electrical engineering technology, and in particular to a method for calculating multi-section power transmission limits based on multi-objective transferable reinforcement learning. Background Art

[0002] With the continuous increase in the penetration rate of new energy and the substantial growth in load, the uncertainty of the power system has become increasingly strong, and the system operation mode has become more complicated. It is increasingly difficult for the dispatching department to accurately identify extreme scenarios and accurately grasp the section limit. In order to reserve a part of the safety margin, the engineering often actively lowers the transmission limit, which further wastes the section transmission potential and hinders the absorption of new energy.

[0003] The calculation of section transmission quota belongs to the problem of solving high-dimensional non-convex nonlinear models. This means that the calculation of transmission quota in actual projects can only be reduced to traditional engineering thinking, relying on intensive manpower and operating experience to perform "trial and error" verification of limited extreme methods, and its accuracy is difficult to verify. At present, there are few methods and theories for calculating section transmission quota, mainly including expert experience method and computational geometry method. The rules formulated by the expert experience method may be too subjective, and it is impossible to find the best mechanism to plan the overall situation. In the process of matching the quota in the expert database, conservative thresholds and rough division of typical operating modes are considered, which reduces the section transmission capacity. The computational geometry method often loses some features related to the operating mode, and its quota calculation results are also conservative, and there are limitations when facing coupled sections with more than three dimensions.

[0004] The multi-section transmission limit setting problem evolves into a non-convex nonlinear multi-objective optimization problem with differential equations. The solution to this problem is not limited to a single optimal solution, but a Pareto front composed of multiple optimal solutions. This Pareto front can be considered as a stable convex boundary in the section flow cut set space. In the process of solving the non-convex nonlinear model, there may be multiple local optimal solutions, it is difficult to ensure convergence, and the computational complexity is too high. In addition, it is affected by multiple operating mode adjustments and operating state constraints. The actual Pareto frontier of this problem is often not a convex boundary in the usual sense, which further increases the difficulty of solving the model. Summary of the invention

[0005] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a multi-section transmission limit calculation method based on multi-objective transferable reinforcement learning.

[0006] The object of the present invention is achieved by the following technical solution: a multi-section transmission limit calculation method based on multi-objective transferable reinforcement learning, the method comprising the following steps:

[0007] Step S1, considering transient power angle and voltage stability constraints, establishing a reverse search model for transmission limit calculation;

[0008] Step S2, using a decomposition strategy to decompose the multi-section coupling limit calculation problem of the reverse search model into sub-problem models;

[0009] Step S3, transform different sub-problem models into Markov decision processes, and solve them by combining the neighborhood-based parameter migration strategy with the improved DDPG algorithm;

[0010] Step S4: Implement the Pareto frontier (stable convex boundary) of the end-to-end output multi-section coupling limit by training the sub-problem model.

[0011] Specifically, the reverse search model for transmission limit calculation takes minimizing the external power transmission of section instability as the objective function, and its expression is:

[0012]

[0013] Among them, P ij For contact line I ij The transmission power in the specified direction along the upper edge; Z is the set of tie lines;

[0014] Specifically, the system power flow equation of the reverse search model for transmission limit calculation is:

[0015]

[0016] Among them, i∈N, N is the set of all nodes; P Gi , Q Gi are the active power and reactive power output by the generator respectively; Q Ci Inject power into the reactive power compensation device; P Li , Q Li are the node active and reactive loads respectively; U i ,θ i is the node voltage amplitude and phase angle; Y ij , α ij is the magnitude and phase angle of the node admittance matrix;

[0017] The synchronous machine rotor motion equation of the reverse search model for the transmission limit calculation is:

[0018]

[0019] Where i∈S G , S G for all generators aggregate; and are the rotor speed at time τ and the rated speed of the rotor respectively; is the rotor angle; E′ di and E′ qiare the d-axis and q-axis components of the generator voltage respectively; M i is the inertia constant; T Mi is the output mechanical power; D i is the damping torque coefficient; I di and I qi are the d-axis and q-axis components of the generator current respectively; X di With X qi is the synchronous reactance; X′ di and X′ qi is the transient reactance; T′ d0i and T′ q0i are the open circuit time constants of the d-axis and q-axis respectively; E fdi Specifically, the constraints of the reverse search model for the power transmission quota calculation include steady-state operation constraints and transient instability constraints. The steady-state operation constraints are expressed as:

[0020]

[0021] Among them, P Gi , Q Gi are the active power and reactive power output by the generator respectively; U i is the node voltage amplitude; I ij For line l ij Current after a fault occurs; For line I ij Current upper limit; S G represents the set of generators; N represents the set of nodes; Z represents the set of tie lines.

[0022] Specifically, the transient instability constraint is:

[0023]

[0024] Among them, x(τ) is the state variable; y(τ) is the algebraic variable; u is the control variable; x(τ0) and y(τ0) represent the initial operation mode; τ is the simulation time step; (τ0,τ end ] indicates the transient process time period; Transient power angle instability constraints and transient voltage instability constraints;

[0025] Transient power angle stability evaluation index f TSI :

[0026]

[0027] Among them, f TSI <0, the system is transiently unstable, and f TSI The larger the system, the more stable it is.

[0028] Voltage stability evaluation index ηU,i :

[0029]

[0030] The voltage is stable when the bus voltage at the central point of the system is lower than 0.75pu and lasts for less than 1s, and critical voltage is stable when the voltage lasts for 0.8 to 1s; and voltage is unstable when the bus voltage at the central point is less than 0.75pu and lasts for more than 1s; where i∈N; U i,min The node voltage U is the transient process after the fault is cleared. i Minimum value; U i,cr is the voltage threshold of node i, which is 0.75 pu; T cr To allow V i,min <0.75pu duration, take 1s; T b The node voltage is lower than U i,cr The total duration of cr,j represents the jth U i,min =U i,cr The time step of U i,sta is the voltage when node i returns to normal; k τ is the critical voltage shift factor, which is taken as 0.75; when η U,i <0, it is determined that the busbar transient voltage is unstable;

[0031] The transient instability constraint satisfies the transient power angle instability constraint or the voltage instability constraint. The transient instability constraint is processed by the Big-M method. The model is easy to solve by rewriting equation (5) into the following equation:

[0032]

[0033] Among them, M1 and M2 are sufficiently small numbers; δ i is a 0-1 decision variable, when δ I =1, the i-th constraint is satisfied.

[0034] Specifically, when calculating the multi-section coupled transmission limit, multiple sections are optimized simultaneously, and the specific definition is as follows:

[0035]

[0036] in, Contains the state, algebra, and control variables of the power system; The objective function consists of m different independent section limit calculations; represent the system equality constraints and inequality constraints respectively, Including stable operation constraints and transient instability constraints.

[0037] Specifically, in step S2, the multi-section coupling limit calculation problem is converted into n scalar optimization sub-problems by weighted summation, and the converted scalar objective function is expressed as:

[0038]

[0039] in, represents the objective function of the ith subproblem, is the weight vector, and m is the number of cross sections to be solved. Specifically, the Markov decision process consists of five tuples:<S,A,r,P,γ> The composition is described as follows:

[0040] State space S:

[0041] S=[P G ,U G ,P L ,Q L ,P f ,f TSI ,η U ] (11)

[0042] Among them, P G is the active output of the generator; U G is the node voltage; P L is the load active power; Q L is the load reactive power; P f is the sum of the transmission flows of all tie lines in the selected section; f TSI is the transient power angle stability evaluation index, η U It is the transient voltage stability evaluation index;

[0043] The action space A includes the generator action space a gen and the total load action space a load , a gen is the generator active power adjustment ΔP g , its adjustment range is N represents the maximum active power of the generator and the number of training steps per round of RL; a load Load active power adjustment ΔP Load_all , its adjustment range is Indicates the maximum active power of total load fluctuation;

[0044] The reward function r is expressed as:

[0045] r con =-[f TSI +M1(1-δ1)]-[η U +M2(1-δ2)]+(δ1+δ2-1) (12)

[0046]

[0047] Among them, M1 and M2 are sufficiently small numbers; δ i is a 0-1 decision variable, when δ i =1, the i-th constraint is satisfied; r con is the reward value for the agent to explore the operation mode that meets the transient instability constraint, which is a continuous reward; g is the reward function calculation condition; r i- The total is the final reward value of subproblem i; r i-obj is the normalized reward value of the scalar objective function of subproblem i; ω s ,s=1,2 is the weight coefficient of the reward item; K is the penalty for power flow non-convergence / line overload / power flow reversal / generator over-limit;

[0048] The discount factor γ, γ∈[0,1] represents the discount on future rewards;

[0049] The state transfer matrix P is determined by the interaction environment.

[0050] Specifically, in step S3, the decomposition method of the multi-objective computing problem and the parameter migration strategy are combined to model the limit calculation sub-problem as a neural network through the MORL method. represents the network model parameters of the i-1th sub-problem, and defines [ω,b] as the network parameters that have not been optimized, [ω * ,b * ] is the optimized neural network parameters. Assuming that the i-1th limit calculation subproblem has been solved and its network parameters have reached the optimal or near-optimal state, when training the i-th subproblem, the empirical knowledge of subproblem i-1 can be used to quickly solve it, and the optimal network parameters of the latter can be As the starting parameters of the former neural network.

[0051] Specifically, the improved deep deterministic policy gradient method is used to train the MDP model of the sub-problem, allowing the agent to fully interact with the environment, increase the random exploration ability, and consider adding random noise to the action:

[0052] a t =μ(s t |θ μ )+ξ t (16)

[0053] Among them, μ(s t |θ μ ) is the Actor network, ξ t represents random noise;

[0054] Execute action a t After that, you will get an immediate reward. t+1 and the next state S t+1, forming a quaternion t ,a t ,r t+1 ,S t+1 >Stored in the experience replay pool R, randomly extract Batch_size quadruplets from R as the input of the current Actor network and Critic network, and calculate the true value Q i :

[0055]

[0056] Among them, γ is the discount factor.

[0057] Update the Critic network parameters θ by minimizing the loss function and calculating the policy gradient Q and Actor network parameters θ μ :

[0058]

[0059] Among them, L is the loss value when the Critic network is updated. By minimizing the loss function, the Q value of the Critic network is closer to the target Q value; N is the number of samples sampled from the experience replay pool for each training; The gradient of the Actor network is calculated based on the Q value of the Critic network to update the weight of the Actor network; It is the gradient of the Q value of the Critic network with respect to the action, reflecting the impact of action selection on the Q value.

[0060] Finally, the obtained (θ Q ,θ μ ) Update the target network parameters through a relatively stable soft update method:

[0061]

[0062] Among them, τ is the divergence factor. By repeatedly iterating equations (16) to (20), the network parameters are continuously updated and the algorithm eventually converges.

[0063] The present invention has the following advantages:

[0064] ​The present invention follows the idea of ​​solving multi-objective optimization problems, and proposes a method for calculating multi-section transmission limits based on multi-objective reinforcement learning, which realizes the end-to-end rapid output of the precise Pareto frontier of multi-section coupling limits. This frontier is the stable convex boundary that needs to be found in engineering. Compared with the traditional multi-objective meta-heuristic algorithm, the present invention saves 81.2% of the solution time while maintaining good convergence and diversity, and the obtained hypervolume index is improved by as little as 8.9%. In addition, the introduced neighborhood-based parameter migration strategy allows the current sub-problem to refer to the model parameters of the previous sub-problem, which greatly improves the training speed while improving the convergence of the algorithm. The multi-objective transferable reinforcement learning method adopted by the present invention provides a new automated solution for the multi-section transmission limit setting that is difficult to set in engineering. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 It is a schematic diagram of the calculation method flow of the present invention;

[0066] Figure 2 It is a schematic diagram of decomposing the sub-problems of the present invention;

[0067] Figure 3 It is a schematic diagram of the overall framework of the MORL algorithm calculation of the present invention;

[0068] Figure 4 A schematic diagram of the neighborhood-based parameter migration strategy process of the present invention;

[0069] Figure 5 It is a schematic diagram of the system structure and cross-section division of the present invention;

[0070] Figure 6 A schematic diagram of a reward value curve during the training process of the first sub-problem of the present invention;

[0071] Figure 7 It is a schematic diagram of the tidal current change of Section 1 and Section 2 after the reward value curve of the present invention is flattened;

[0072] Figure 8 A Pareto frontier diagram of the MORL of the present invention on the two-section transmission quota problem;

[0073] Fig. 9 It is a schematic diagram of stable and unstable working conditions in two cross-section tidal current spaces of the present invention;

[0074] Fig.10 It is a schematic diagram of the search trajectory of the cross-section power transmission limit of the present invention;

[0075] Fig.11 A schematic diagram showing the comparison of the Pareto fronts obtained by the MORL algorithm of the present invention and the NSGA-II algorithm;

[0076] Fig.12 A schematic diagram showing the comparison of the Pareto fronts obtained by the MORL and MOEA / D algorithms of the present invention;

[0077] Fig.13 The training performance comparison diagram of the parameter migration strategy adopted by the present invention and the original strategy (I);

[0078] Fig.14 This is a comparison chart of the training performance between the parameter migration strategy adopted by the present invention and the original strategy (II). DETAILED DESCRIPTION

[0079] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention, that is, the embodiments described are only part of the embodiments of the present invention, rather than all of the embodiments. The components of the embodiments of the present invention described and shown in the drawings herein can be arranged and designed in various different configurations.

[0080] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making creative work are within the scope of protection of the present invention.

[0081] It should be noted that relational terms such as "first" and "second" are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of further restrictions, the elements defined by the sentence "including a..." do not exclude the existence of other identical elements in the process, method, article or device including the elements.

[0082] The present invention is further described below in conjunction with the accompanying drawings, but the protection scope of the present invention is not limited to the following description.

[0083] like Figures 1 to 14 As shown, a multi-section transmission limit calculation method based on multi-objective transferable reinforcement learning includes the following steps:

[0084] Step S1, considering transient power angle and voltage stability constraints, establishing a reverse search model for transmission limit calculation;

[0085] Step S2, using a decomposition strategy to decompose the multi-section coupling limit calculation problem of the reverse search model into sub-problem models;

[0086] Step S3, transform different sub-problem models into Markov decision processes, and solve them by combining the neighborhood-based parameter migration strategy with the improved DDPG algorithm;

[0087] Step S4: Implement the Pareto frontier (stable convex boundary) of the end-to-end output multi-section coupling limit by training the sub-problem model.

[0088] The traditional method for calculating the transmission limit is to first determine the local operation boundary and then select a limited operation mode to perform a traversal calculation of the TTC, and take its upper (lower) limit to obtain the limit value. This method is greatly affected by the operation mode and has low calculation efficiency. Whether the calculation result is accurate also needs to be considered. In order to improve the calculation efficiency of the transmission limit as much as possible, the present invention proposes a reverse search model for the section transmission limit. The section transmission limit can be obtained by reverse searching from the instability space. The reverse search model for the transmission limit calculation takes minimizing the transmission power of the section instability as the objective function, and its expression is:

[0089]

[0090] Among them, P ij For contact line I ij The transmission power in the specified direction along the upper edge; Z is the set of tie lines;

[0091] Furthermore, the system power flow equation of the reverse search model for the transmission limit calculation is:

[0092]

[0093] Among them, i∈N, N is the set of all nodes; P Gi , Q Gi are the active power and reactive power output by the generator respectively; Q Ci Inject power into the reactive power compensation device; Li , Q Li are the node active and reactive loads respectively; U i ,θ i is the node voltage amplitude and phase angle; Y ij , α ij is the magnitude and phase angle of the node admittance matrix;

[0094] The synchronous machine rotor motion equation of the reverse search model for the transmission limit calculation is:

[0095]

[0096] Where i∈S G , S G for all generators aggregate; and are the rotor speed at time τ and the rated speed of the rotor respectively; is the rotor angle; E′ di and E′ qi are the d-axis and q-axis components of the generator voltage respectively; M i is the inertia constant; T Mi is the output mechanical power; D i is the damping torque coefficient; I di and I qi are the d-axis and q-axis components of the generator current respectively; X di With X qi is the synchronous reactance; X′ di and X′ qi is the transient reactance; T′ d0i and T′ q0i are the open circuit time constants of the d-axis and q-axis respectively; E fdi is the excitation output voltage. Further, the constraints of the reverse search model for the transmission limit calculation include steady-state operation constraints and transient instability constraints, and the steady-state operation constraints are expressed as:

[0097]

[0098] Among them, P Gi , Q Gi are the active power and reactive power output by the generator respectively; U i is the node voltage amplitude; I ij For line l ij Current after a fault occurs; For line I ij Current upper limit; S G represents the generator set; N represents the node set; Z represents the tie line set; Formula (4) represents the upper and lower limits of power source active output, the upper and lower limits of power source reactive power, the upper and lower limits of node voltage and the thermal stability constraints of the line respectively;

[0099] The transient instability constraint is:

[0100]

[0101] Among them, x(τ) is the state variable; y(τ) is the algebraic variable; u is the control variable; x(τ0) and y(τ0) represent the initial operation mode; τ is the simulation time step; (τ0,τ end ] indicates the transient process time period; Transient power angle instability constraints and transient voltage instability constraints;

[0102] Transient power angle stability evaluation index f TSI :

[0103]

[0104] Among them, f TSI <0, the system is transiently unstable, and f TSI The larger the system, the more stable it is;

[0105] Voltage stability evaluation index η U,i :

[0106]

[0107] The voltage is stable when the bus voltage at the central point of the system is lower than 0.75pu and lasts for less than 1s, and critical voltage is stable when the voltage lasts for 0.8 to 1s; and voltage is unstable when the bus voltage at the central point is less than 0.75pu and lasts for more than 1s; where i∈N; U i,min The node voltage U is the transient process after the fault is cleared. i Minimum value; U i,cr is the voltage threshold of node i, which is 0.75 pu; T cr To allow V i,min <0.75pu duration, take 1s; T b The node voltage is lower than U i,cr The total duration of cr,j represents the jth U i,min =U i,cr The time step of U i,sta is the voltage when node i returns to normal; k τ is the critical voltage shift factor, which is taken as 0.75; when η U,i <0, the busbar transient voltage is judged to be unstable; hereinafter, η U Represents the voltage stability margin of the system, η U is η U,i The minimum value of .

[0108] Furthermore, the transient instability constraint adopted by the section transmission limit calculation model is an "or" form constraint, that is, the reverse search model of the transmission limit calculation satisfies the transient power angle instability constraint or the voltage instability constraint. The Big-M method is used to process the transient instability constraint. By rewriting equation (5) into the following equation, the model is easy to solve:

[0109]

[0110] Among them, M1 and M2 are sufficiently small numbers; δi is a 0-1 decision variable, when δ I = 1, the i-th constraint is satisfied, that is, This method can be used to design RL reward functions;

[0111] Furthermore, complex power systems often need to consider the impact of multiple sections coupled with each other. Therefore, it is necessary to calculate the multi-section coupled transmission limit. This problem belongs to a multi-objective optimization problem. When calculating the multi-section coupled transmission limit, multiple sections are optimized at the same time. The specific definition is as follows:

[0112]

[0113] in, Contains the state, algebra, and control variables of the power system; The objective function consists of m different independent section limit calculations; represent the system equality constraints and inequality constraints respectively, Including stable operation constraints and transient instability constraints. The goal of this function is to find a set of optimal solutions that weigh different objectives, that is, Pareto optimal solutions. Furthermore, the present invention intends to introduce a decomposition strategy to model the multi-section coupled limit calculation problem, decomposing it into a set of scalar optimization sub-problems, and solving the sub-problems by means of coupled section collaborative optimization to obtain the optimal solutions corresponding to different sub-problems. These optimal solutions constitute the Pareto frontier of the original problem. The present invention uses the well-known weighted sum (Weighted sum, WS) method to decompose the multi-section limit calculation problem. In step S2, the multi-section coupled limit calculation problem is converted into n scalar optimization sub-problems by means of weighted sum. The converted scalar objective function is expressed as:

[0114]

[0115] in, represents the objective function of the ith subproblem, is the weight vector, and m is the number of sections to be solved. The Pareto frontier of the multi-section transmission quota problem can be described by the n optimal solutions of the above n sub-problems. The larger the value of n, the more thorough the decomposition and the more complete the description. The process can be described by Figure 2 Intuitive explanation.

[0116] Furthermore, the multi-section quota calculation problem can be decomposed into multiple sub-problems with different weight vectors, and finally the optimal solutions of all sub-problems constitute the Pareto frontier of the original problem. Therefore, the solution of the sub-problems is crucial, and the objective function of the model is unique. First, the section quota calculation model is reconstructed into the form of a Markov decision process (MDP), which consists of a five-tuple<S,A,r,P,γ> The composition is described as follows:

[0117] State space S:

[0118]

[0119] Among them, P G is the active output of the generator; U G is the node voltage; P L is the load active power; Q L is the load reactive power; P f is the sum of the transmission flows of all tie lines in the selected section; f TSI is the transient power angle stability evaluation index, η U It is the transient voltage stability evaluation index; S t represents the state vector at time t;

[0120] The action space A includes the generator action space a gen and the total load action space a load , a gen is the generator active power adjustment ΔP g , its adjustment range is N represents the maximum active power of the generator and the number of training steps per round of RL; a load Load active power adjustment ΔP Load_all , its adjustment range is Indicates the maximum active power of the total load fluctuation; in order to improve the stability of the training, the output of the generators other than the balancing machine is constrained to The load total active power is constrained within In the following, a i =[a gen ,a load] represents the action vector reward function r at time t. The reward function r is a hard indicator for evaluating the quality of the agent's actions and plays a decisive role in the solution process of MDP. The reward function is designed by combining discrete rewards and continuous rewards. First, the constraints are processed. When the power flow does not converge, the line is overloaded, the power flow is reversed, etc., a larger penalty is given to increase the training speed; for transient instability constraints, the transient power angle stability index and the voltage stability margin are selected as key goals, and the continuous reward value is designed in combination with the Big-M processing method. Then, the scalar objective function of the sub-problem is normalized and added to the total reward function. For sub-problem i, the reward function can be expressed by the following formula:

[0121] r con =-[f TSI +M1(1-δ1)]-[η U +M2(1-δ2)]+(δ1+δ2-1) (12)

[0122]

[0123] Among them, M1 and M2 are sufficiently small numbers; δ i is a 0-1 decision variable, when δ i =1, the i-th constraint is satisfied; r con is the reward value for the agent to explore the operation mode that meets the transient instability constraint, which is a continuous reward; g is the reward function calculation condition; r i- The total is the final reward value of subproblem i; r i-obj is the normalized reward value of the scalar objective function of subproblem i; ω s ,s=1,2 is the weight coefficient of the reward item; K is the penalty for power flow non-convergence / line overload / power flow reversal / generator over-limit;

[0124] The discount factor γ, γ∈[0,1] represents the discount on future rewards;

[0125] The state transfer matrix P is determined by the interaction environment.

[0126] The MORL computational framework constructed by coordinating the optimization of the MDP model and multiple sub-problems is as follows: Figure 3 shown.

[0127] When there are many subproblems, solving the above multiple MDP model takes a lot of time. This is because for different subproblems, the agent has to re-explore the optimal strategy, and does not make good use of RL's inherent self-learning and knowledge transfer capabilities. According to formula (10), two adjacent subproblems have similar weight vectors, so the optimal solutions of the two are very likely to be similar. Based on this, a subproblem can be solved by drawing on the knowledge of its adjacent subproblems.

[0128] In step S3, the decomposition method of the multi-objective computing problem and the parameter migration strategy are combined to model the limit calculation sub-problem as a neural network through the MORL method. represents the network model parameters of the i-1th sub-problem, and defines [ω,b] as the network parameters that have not been optimized, [ω * ,b * ] is the optimized neural network parameters. Assuming that the i-1th limit calculation subproblem has been solved and its network parameters have reached the optimal or near-optimal state, when training the i-th subproblem, the empirical knowledge of subproblem i-1 can be used to quickly solve it, and the optimal network parameters of the latter can be As the starting parameters of the former neural network. In this way, by transferring the model parameters of the previous adjacent subproblem to the current subproblem as the starting parameters, the training time cost of the model is effectively reduced. In simple terms, the network model parameters are deployed to the adjacent subproblems in an orderly manner along the direction of weight reduction (increase). This pattern is as follows Figure 4 The strategy used makes it possible for MORL to quickly solve the multi-section limit problem. Combining the decomposition method of multi-objective calculation problems and the parameter migration strategy, the multi-section transmission limit calculation process based on MORL is shown in Table 1:

[0129] Table 1 Pseudo code of the overall calculation process of the algorithm

[0130]

[0131]

[0132] Furthermore, the improved deep deterministic policy gradient method is used to train the MDP model of the sub-problem. DDPG draws on the idea of ​​DQN replaying the experience library to break the correlation between data, improves the stability of training, and realizes the rapid convergence of the algorithm through the Actor-Critic network. The original DDPG contains 4 neural networks, namely, Actor network μ(s|θ μ ), Critic network Q(s,a|θ Q ) and the copied target network μ tar (s|θ μ ) and Q tar (s,a|θ Q ), let the agent fully interact with the environment, increase the ability of random exploration, and consider adding random noise to the action:

[0133] a t =μ(s t |θ μ )+ξ t (16)

[0134] Among them, μ(s t |θ μ ) is the Actor network, ξ t represents random noise;

[0135] Execute action a t After that, you will get an immediate reward. t+1 and the next state S t+1 , forming a quaternion t ,a t ,r t+1 ,S t+1 >Stored in the experience replay pool R, randomly extract Batch_size quadruplets from R as the input of the current Actor network and Critic network, and calculate the true value Q i :

[0136]

[0137] Among them, γ is the discount factor.

[0138] Update the Critic network parameters θ by minimizing the loss function and calculating the policy gradient Q and Actor network parameters θ μ :

[0139]

[0140] Among them, L is the loss value when the Critic network is updated. By minimizing the loss function, the Q value of the Critic network is closer to the target Q value; N is the number of samples sampled from the experience replay pool for each training; The gradient of the Actor network is calculated based on the Q value of the Critic network to update the weight of the Actor network; It is the gradient of the Q value of the Critic network with respect to the action, reflecting the impact of action selection on the Q value.

[0141] Finally, the obtained (θ Q ,θ μ ) Update the target network parameters through a relatively stable soft update method:

[0142]

[0143] Among them, τ is the divergence factor. By repeatedly iterating equations (16) to (20), the network parameters are continuously updated and the algorithm eventually converges.

[0144] ​The improved method of DDPG adopted by the present invention uses the experience replay method for training in the original DDPG algorithm. This method fully considers the temporal correlation between samples and reuses samples. However, this method ignores the importance of different samples. Randomly extracting batch samples is likely to cause a large number of low reward value samples to be trained, which in turn causes DDPG to converge slowly. The present invention adopts a classified experience replay method based on instant rewards for DDPG, and puts samples with instant rewards less than the average reward value into experience buffer pool 1, and samples with instant rewards greater than the average reward value into experience buffer pool 2. When training batch samples in proportion, more samples are extracted from buffer pool 2 each time, and a small number of samples are extracted from buffer pool 1, in order to maintain sample diversity while improving convergence performance.

[0145] The above is the conventional training process of the improved DDPG algorithm for single-objective problems. The training process of the sub-problem model after decomposition of the multi-section transmission limit calculation problem is shown in Table 2:

[0146] Table 2 Pseudocode of sub-problem training algorithm

[0147]

[0148]

[0149] Example system:

[0150] The present invention uses the IEEE 39-bus system as a test example. The coupling sections investigated are the transmission section 1 composed of the AC tie lines 1-39, 2-3 and 18-3 and the transmission section 2 composed of the AC tie line 15-16. The system structure and section division are as follows: Figure 5 As shown. Considering the thermal stability constraints, transient power angle stability constraints and transient voltage stability constraints, the key fault of the section is set as a three-phase short circuit fault of the tie line at 0.1s, and the fault line is removed at 0.15s; the transient stability simulation time is 5s. The system benchmark power is 100MW, the generator adopts the classic model, considering the generator excitation control, and the load adopts the constant impedance model.

[0151] The continuous action space of the generator and the total load trained by the MORL algorithm is set to: a_bound = [-0.6, 0.6], p L_all _bound=[-6,6]. The system uses a three-layer BP neural network to construct the Actor network and the Critic network. During the training process, the playback mechanism uses a classified buffer pool, in which the Batch_size extraction ratio of the two types of buffer pools is 0.7, and random noise is added to the action selection. The other network parameter settings are shown in Table 3:

[0152] Table 3 Parameter settings of network model

[0153]

[0154]

[0155] Validity analysis of multi-section coupling limit calculation method:

[0156] The simulation starts with a randomly selected system instability mode. The initial sub-problem is trained first, and the model weights of the subsequent sub-problems are transferred using the field parameter migration strategy. The number of training steps for the first sub-problem is 60,000 steps, and data statistics are performed with 20 steps as one round. The data during the training process is normalized. The reward value curve during the training process of the first sub-problem is as follows: Figure 6 shown.

[0157] Depend on Figure 6 It can be seen that the reward value curve fluctuates greatly before about 500 rounds, and basically reaches convergence after 1000 rounds of training. The reason for the oscillation is that the decision-making subject explores the optimal strategy through continuous "trial and error", but this does not affect the convergence process. Before 500 rounds, the reward value curve oscillates violently. At this time, the system is often in a stable state or the power flow does not converge, which cannot reflect the process of the intelligent agent searching for the transmission limit boundary from the unstable operation mode (the power flow when the system is stable is often relatively small). After the curve tends to stabilize, according to the definition of the reward function, the operation mode of the system at this time is basically unstable, which meets the constraint conditions of the limit calculation. Figure 7 The flow changes of Section 1 and Section 2 after the reward value curve is flattened are given. The section flow corresponding to each round is the cumulative average of 20 update steps. It can be seen that both target values ​​have undergone huge changes within 1000 to 2000 training rounds. The flow of Section 1 dropped sharply from 800MW and finally converged to around 600MW; the flow of Section 2 dropped from 500MW to 350MW after repeated rises and falls and reached stability. The above shows that the decision result of the MORL method is reasonable. For the sub-problem under a certain weight, the section transmission limit can be well searched. Figure 7 The last converged section power flow value is the coupled section transmission limit value of the first sub-problem. However, in actual scheduling, the section transmission limit value under a certain weight is not of concern, and it is even more unrealistic to impose preferences on sections. The present invention hopes to find the convex boundary domain of the mutually restrained section transmission limits to guide the optimization control of multiple sections.

[0158] For the first sub-problem model trained, the parameter migration strategy based on the neighborhood is used to transfer model parameters to the neighboring sub-problems. A small number of training rounds are performed with 20 update steps to generate 50 sets of weight vectors, i.e., 50 scalar optimization sub-problems. The Pareto frontier of the multi-objective transferable reinforcement learning algorithm on the two-section transmission limit problem is obtained. The visualization results are shown in the figure below. Figure 8 As shown. It can be seen that the Pareto frontier of the section transmission limit problem is different from the Pareto frontier of the conventional two-objective optimization (the conventional PF is a "slide" type), and its shape is similar to an "S" curve. This is because the two sections selected by the present invention are mutually clamped, and are not a simple relationship of increasing and decreasing together. This is determined by the non-convex nonlinear model of the section limit calculation. This Pareto frontier is the safe convex boundary of the two section limits. Through rough division, it can be determined that the system in the lower left corner flow space is roughly stable, and the upper right corner flow space is roughly unstable.

[0159] In order to verify the effectiveness of the above MORL method in searching the Pareto frontier, the present invention verifies the operation mode in the two cross-section flow space by time domain simulation to determine the stability boundary. The generator output and total load are set to fluctuate randomly within their upper and lower limits, and the terminal voltage fluctuates within the value range. Latin hypercube sampling is used to generate 12,000 operating conditions, and a three-phase short-circuit fault is set as the expected accident set, and a transient stability check is performed on it. Fig. 9 A schematic diagram of stable and unstable operating conditions in the tidal space of two sections is given. In the figure, circles represent stable modes and crosses represent unstable modes. Since stable and unstable operating modes are confused, the unstable mode is overlaid on the stable mode to better distinguish the safety boundary. Figure 8 and Fig. 9 It can be seen that the convex boundary of the cross-section limit searched by MORL is basically consistent with the stable boundary obtained by sampling, which verifies the effectiveness of the method proposed in this invention. In addition, the boundary searched by MORL is more conservative than the sampling boundary, which may be caused by the insufficient number of samples, further illustrating that the search accuracy of the method of this invention is high.

[0160] To intuitively describe the relationship between training data and results, Fig.10 The search trajectory of the intelligent agent for the transmission limit operation point of the section is given. As can be seen from the figure, different transmission limit boundary points will be obtained by using different training weights under the initial operating condition of point A. This is because the selection of different weight vectors will change the reward function in MORL, which will lead to certain differences in the searched transmission limits. In addition, when the intelligent agent searches for a transmission limit similar to that obtained at the initial point A under the initial operating condition of point B, the weight vectors used by the two are basically the same, which proves that the proposed method has good adaptability to the changing initial operating mode.

[0161] Comparison results with meta-heuristic algorithms:

[0162] In order to better verify the superiority of the MORL algorithm, this paper uses the most widely used multi-objective metaheuristic algorithm as the comparison algorithm, and selects the NSGA-II and MOEA / D algorithms that perform best in this field for testing. To ensure the accuracy of the test results, the example is repeated 20 times each time, and the best result is taken as the comparison data. Fig.11 and 12 The Pareto fronts of section transmission limit obtained by different iteration numbers of NSGA-II algorithm and MOEA / D algorithm are compared with the search results of MORL.

[0163] Depend on Fig.11 It can be seen that for NSGA-II that executes the normal number of iterations (100 times), the Pareto front obtained is quite different from the actual result, and it is easy to fall into the local optimum. As the number of iterations gradually increases, the convergence level of NSGA-II gradually improves. When it increases to 8 times the initial number, the Pareto front of the section transmission limit searched by NSGA-II is close to the MORL algorithm, but it is slightly lower than the MORL algorithm in terms of convergence and diversity. Fig.12 In the figure, when the number of iterations increases to 800, the local Pareto front searched by the MOEA / D algorithm is very similar to the MORL search results and even performs better locally. However, the MOEA / D algorithm is significantly inferior to MORL and NSGA-II in terms of diversity, and all the non-dominated solutions obtained are concentrated in a small area, which is unacceptable for constructing a complete transmission limit stability convex boundary. From the above, it can be seen that the MORL algorithm proposed in the present invention is superior to the NSGA-II and MOEA / D algorithms in terms of convergence and diversity methods, which verifies that the MORL algorithm has a stronger global optimal search performance.

[0164] In order to further intuitively demonstrate the specific performance of each algorithm, Table 4 gives the hypervolume (HV) index and computational time of different algorithms (including different iterations of the same algorithm). It can be seen that MORL has the largest HV index among all algorithms, that is, MORL has the best convergence and diversity. In terms of computational time, the average running time of MORL is 117.93s, which is much lower than the 800 iterations of NSGA-II and MOEA / D algorithms with better convergence. Although the NSGA-II algorithm takes less time to solve the problem with 100 iterations, its convergence is extremely poor, while MORL only takes about 40s more to search for the Pareto frontier with the best performance in all aspects, showing the high efficiency of the MORL algorithm. In general, the NSGA-II algorithm has better diversity but poor convergence; the MOEA / D algorithm has better convergence but poor diversity, and it takes a lot of time to solve. The MORL algorithm improves the HV index by 8.9% and saves 81.2% of the solution time compared with the above optimal multi-objective metaheuristic algorithm.

[0165] Table 4. Super volume index and computational time of the algorithm for solving the transmission quota problem

[0166]

[0167] It is worth mentioning that the MORL algorithm of the present invention took 13.69 hours to train the transmission limit stability boundary. This is because the intelligent agent needs to constantly change the generator output action to perform power flow calculation and transient stability analysis and solve multiple sub-problems at the same time. The trained model can be directly used to output the Pareto frontier.

[0168] Effectiveness analysis of parameter migration strategy:

[0169] The present invention first compares the model training conditions of the original strategy and the migration strategy. The original strategy trains each sub-problem for 3000 rounds with 20 update steps per round. Fig.13 The Pareto frontiers of the outputs of the two training strategies are given. It can be found that when the same number of rounds of training is used, the model convergence and diversity of the original strategy are very poor. This is because when training a sub-problem from scratch each time, the agent cannot use the knowledge learned in the past and loses the key information of the neighborhood.

[0170] In addition, in order to eliminate the result error caused by the number of training times, the present invention further expands the number of training rounds under the original strategy to 12,000, updates 20 steps per round, and displays the results after training on Fig.14As can be seen from the figure, if the parameters of the neighborhood are not transferred, even if the number of training times increases by 3 times, the performance of the original strategy is still not good enough. Therefore, the application of the neighborhood-based parameter transfer strategy to the MORL algorithm is very practical and efficient, which greatly reduces the training time while improving the convergence.

[0171] The present invention follows the idea of ​​solving multi-objective optimization problems, and proposes a method for calculating multi-section transmission limits based on multi-objective reinforcement learning, which realizes the end-to-end rapid output of the precise Pareto frontier of multi-section coupling limits. This frontier is the stable convex boundary that needs to be found in engineering. Compared with the traditional multi-objective metaheuristic algorithm, the proposed method saves 81.2% of the solution time while maintaining good convergence and diversity, and the obtained hypervolume index is improved by as little as 8.9%. In addition, the introduced neighborhood-based parameter migration strategy allows the current sub-problem to refer to the model parameters of the previous sub-problem, which greatly improves the training speed while improving the convergence of the algorithm. The multi-objective transferable reinforcement learning method adopted by the present invention provides a new automated solution for the multi-section transmission limit setting that is difficult to set in engineering.

[0172] The above is only a preferred embodiment of the present invention and does not limit the present invention in any form. Any technician familiar with the art can make many possible changes and modifications to the technical solution of the present invention by using the above-mentioned technical content without departing from the scope of the technical solution of the present invention, or modify it into an equivalent embodiment of equivalent changes. Therefore, any changes, modifications, equivalent changes and modifications made to the above embodiments based on the technology of the present invention without departing from the content of the technical solution of the present invention belong to the protection scope of the present technical solution.

Claims

1. A multi-section transmission limit calculation method based on multi-objective transferable reinforcement learning, characterized by: The method comprises the following steps: Step S1, considering transient power angle and voltage stability constraints, establishing a reverse search model for transmission limit calculation; Step S2, using a decomposition strategy to decompose the multi-section coupling limit calculation problem of the reverse search model into sub-problem models; Step S3, transform different sub-problem models into Markov decision processes, and solve them by combining the neighborhood-based parameter migration strategy with the improved DDPG algorithm; Step S4: Implement the Pareto frontier of the end-to-end output multi-section coupling limit by training the sub-problem model.

2. The multi-section transmission limit calculation method based on multi-objective transferable reinforcement learning according to claim 1 is characterized in that: The reverse search model for transmission quota calculation takes minimizing the external power transmission of section instability as the objective function, and its expression is: Among them, P ij For contact line I ij The transmission power in the specified direction; Z is the collection of tie lines.

3. The multi-section transmission limit calculation method based on multi-objective transferable reinforcement learning according to claim 1 is characterized in that: The system power flow equation of the reverse search model for transmission limit calculation is: Among them, i∈N, N is the set of all nodes; P Gi , Q Gi are the active power and reactive power output by the generator respectively; Q Ci Inject power into the reactive power compensation device; P Li , Q Li are the node active and reactive loads respectively; U i ,θ i is the node voltage amplitude and phase angle; Y ij , α ij is the magnitude and phase angle of the node admittance matrix; The synchronous machine rotor motion equation of the reverse search model for the transmission limit calculation is: Where i∈S G , S G for all generators aggregate; and are the rotor speed at time τ and the rated speed of the rotor respectively; is the rotor angle; E′ di and E′ qi are the d-axis and q-axis components of the generator voltage respectively; M i is the inertia constant; T Mi is the output mechanical power; D i is the damping torque coefficient; I di and I qi are the d-axis and q-axis components of the generator current respectively; X di With X qi is the synchronous reactance; X′ di and X′ qi is the transient reactance; T′ d0i and T′ q0i are the open circuit time constants of the d-axis and q-axis respectively; E fdi is the excitation output voltage.

4. The method for calculating multi-section power transmission quota based on multi-objective transferable reinforcement learning according to claim 3 is characterized in that: The constraints of the reverse search model for transmission quota calculation include steady-state operation constraints and transient instability constraints. The steady-state operation constraints are expressed as: Among them, P Gi , Q Gi are the active power and reactive power output by the generator respectively; U i is the node voltage amplitude; I ij For line l ij Current after a fault occurs; For line I ij Current upper limit; S G represents the generator set; N represents the node set; Z represents the tie line set; The transient instability constraint is: Among them, x(τ) is the state variable; y(τ) is the algebraic variable; u is the control variable; x(τ0) and y(τ0) represent the initial operation mode; τ is the simulation time step; (τ0,τ end ] indicates the transient process time period; Transient power angle instability constraints and transient voltage instability constraints; Transient power angle stability evaluation index f TSI : Among them, f TSI <0, the system is transiently unstable, and f TSI The larger the system, the more stable it is; Voltage stability evaluation index η U,i : The voltage is stable when the bus voltage at the central point of the system is lower than 0.75pu and lasts for less than 1s, and critical voltage is stable when the voltage lasts for 0.8 to 1s; and voltage is unstable when the bus voltage at the central point is less than 0.75pu and lasts for more than 1s; where i∈N; U i,min The node voltage U is the transient process after the fault is cleared. i Minimum value; U i,cr is the voltage threshold of node i, which is 0.75 pu; T cr To allow V i,min <0.75pu duration, take 1s; T b The node voltage is lower than U i,cr The total duration of cr,j represents the jth U i,min =U i,cr The time step of U i,sta is the voltage when node i returns to normal; k τ is the critical voltage shift factor, which is taken as 0.75; when η U,i <0, the busbar transient voltage is determined to be unstable.

5. The method for calculating multi-section power transmission quota based on multi-objective transferable reinforcement learning according to claim 4 is characterized in that: The transient instability constraint satisfies the transient power angle instability constraint or the voltage instability constraint. The transient instability constraint is processed by the Big-M method. The model is easy to solve by rewriting equation (5) into the following equation: Among them, M1 and M2 are real numbers; δ i is a 0-1 decision variable, when δ i =1, the i-th constraint is satisfied.

6. The method for calculating multi-section power transmission quota based on multi-objective transferable reinforcement learning according to claim 5 is characterized in that: When calculating the multi-section coupled transmission limit, multiple sections are optimized simultaneously, and the specific definitions are as follows: in, Contains the state, algebra, and control variables of the power system; The objective function consists of m different independent section limit calculations; represent the system equality constraints and inequality constraints respectively, Including stable operation constraints and transient instability constraints.

7. The method for calculating multi-section power transmission limits based on multi-objective transferable reinforcement learning according to claim 6 is characterized in that: In step S2, the multi-section coupling limit calculation problem is converted into n scalar optimization sub-problems by weighted summation. The converted scalar objective function is expressed as: in, represents the objective function of the ith subproblem, is the weight vector, and m is the number of cross sections to be solved.

8. The method for calculating multi-section power transmission quota based on multi-objective transferable reinforcement learning according to claim 7 is characterized in that: The Markov decision process consists of five tuples<S,A,r,P,γ> The specific description is as follows: State space S: S=[P G ,U G ,P L ,Q L ,P f ,f TSI ,η U ] (11) Among them, P G is the active output of the generator; U G is the node voltage; P L is the load active power; Q L is the load reactive power; P f is the sum of the transmission flows of all tie lines in the selected section; f TSI is the transient power angle stability evaluation index, η U It is the transient voltage stability evaluation index; The action space A includes the generator action space a gen and the total load action space a load , a gen is the generator active power adjustment ΔP g , its adjustment range is N represents the maximum active power of the generator and the number of training steps per round of RL; a load Load active power adjustment ΔP Load_all , its adjustment range is Indicates the maximum active power of total load fluctuation; The reward function r is expressed as: r con =-[f TSI +M1(1-δ1)]-[η U +M2(1-δ2)]+(δ1+δ2-1) (12) Among them, r con is the reward value for the agent to explore the operation mode that satisfies the transient instability constraint, which is a continuous reward; g is the reward function calculation condition, M1 and M2 are real numbers; δ i is a 0-1 decision variable, when δ i = 1, the i-th constraint is satisfied; r i-总 is the final reward value of subproblem i; r i-obj is the normalized reward value of the scalar objective function of subproblem i; ω s ,s=1,2 is the weight coefficient of the reward item; K is the penalty for power flow non-convergence / line overload / power flow reversal / generator over-limit; The discount factor γ, γ∈[0,1] represents the discount on future rewards; The state transfer matrix P is determined by the interaction environment.

9. The method for calculating multi-section power transmission quota based on multi-objective transferable reinforcement learning according to claim 8 is characterized in that: In step S3, the decomposition method of the multi-objective computing problem and the parameter migration strategy are combined to model the limit calculation sub-problem as a neural network through the MORL method. represents the network model parameters of the i-1th sub-problem, and defines [ω,b] as the network parameters that have not been optimized, [ω * , b * ] is the optimized neural network parameters. Assuming that the i-1th limit calculation subproblem has been solved and its network parameters have reached the optimal or near-optimal state, when training the i-th subproblem, the empirical knowledge of subproblem i-1 can be used to quickly solve it, and the optimal network parameters of the latter can be As the starting parameters of the former neural network.

10. The method for calculating multi-section power transmission limits based on multi-objective transferable reinforcement learning according to claim 9 is characterized in that: The improved deep deterministic policy gradient method is used to train the MDP model of the sub-problem. The original DDPG contains four neural networks, namely Actor network μ(s|θ μ ), Critic network Q(s,a|θ Q ) and the copied target network μ tar (s|θ μ ) and Q tar (s,a|θ Q ); In order to allow the agent to fully interact with the environment and increase random exploration capabilities, consider adding random noise to the actions: a t =μ(s t |θ μ )+ξ t (16) Among them, μ(s t |θ μ ) is the Actor network, ξ t represents random noise; Execute action a t After that, you will get an immediate reward. t+1 and the next state S t+1 , forming a quaternion t ,a t ,r t+1 ,S t+1 >Stored in the experience replay pool R, randomly extract Batch_size quadruplets from R as the input of the current Actor network and Critic network, and calculate the true value Q i :​ Among them, γ is the discount factor; Update the Critic network parameters θ by minimizing the loss function and calculating the policy gradient Q and Actor network parameters θ μ : Among them, L is the loss value when the Critic network is updated. By minimizing the loss function, the Q value of the Critic network is closer to the target Q value; N is the number of samples sampled from the experience replay pool for each training; The gradient of the Actor network is calculated based on the Q value of the Critic network to update the weight of the Actor network; It is the gradient of the Q value of the Critic network with respect to the action, reflecting the impact of action selection on the Q value; Finally, the obtained (θ Q ,θ μ ) Update the target network parameters through a relatively stable soft update method: Among them, τ is the divergence factor. By repeatedly iterating equations (16) to (20), the network parameters are continuously updated and the algorithm eventually converges.

Images