Power system Lyapunov function approximation and stability domain estimation method and device

By using a Lyapunov function approximation and stability domain estimation method based on Gaussian process regression for power systems, the problem of insufficient speed and accuracy of traditional methods in power grid stability analysis is solved. This method achieves fast and accurate stability domain estimation and is applicable to power systems with complex dynamic characteristics.

CN119298052APending Publication Date: 2025-01-10SOUTH CHINA UNIV OF TECH
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202411259324.7
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2024-09-09
Publication Date
2025-01-10

AI Technical Summary

Technical Problem

In the existing technology, traditional transient stability analysis methods are difficult to meet the speed requirements of online safety and stability assessment of modern power grids, and existing Lyapunov function construction methods are difficult to achieve a good balance in large-scale power systems with limited computing resources, and cannot be effectively applied to power systems with complex dynamic characteristics.

Method used

A method for approximating the Lyapunov function and estimating the stability region of a power system based on Gaussian process regression is adopted. By establishing a mathematical model of the power system, the bounded inverse Lyapunov function is learned using the Gaussian process regression model, and combined with a stochastic direction search algorithm, the boundary of the stability region is quickly estimated.

Benefits of technology

It achieves fast and accurate stability domain estimation in large-scale power systems, with high computation speed and accuracy, unaffected by the increase in the dimension of state variables, applicable to nonanalytic Lyapunov functions, and suitable for stability analysis of nonlinear dynamic systems.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN119298052A_ABST
    Figure CN119298052A_ABST
Patent Text Reader

Abstract

The invention discloses an electric power system Lyapunov function approximation and stability domain estimation method and device, and the method comprises the steps: building an electric power system mathematical model, and solving a stable equilibrium point; randomly generating a group of initial sample points with stable labels near the stable equilibrium point, and calculating a corresponding bounded inverse Lyapunov function value; constructing and obtaining an initial sample data set; taking the initial sample data set as input, learning a bounded inverse Lyapunov function of the power system by using a Gaussian process regression model, searching a state point closest to the currently estimated stability domain boundary, and calculating a corresponding bounded inverse Lyapunov function value; iteratively training the model; and when it is detected that an iteration stopping condition is reached, stopping the iteration process, and outputting an estimation result of the stability domain. According to the method, the time complexity is not increased along with the increase of the dimensionality of the state variable of the power system, the form of a function does not need to be limited, and the method has higher calculation speed and higher stability domain estimation precision.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] The present invention relates to the technical field of power system stability domain estimation, and in particular to a power system Lyapunov function approximation and stability domain estimation method and device. Background Art

[0002] With the widespread access to renewable energy and the continuous expansion of the scale of grid interconnection, the dynamic behavior of my country's power system is becoming increasingly complex. Traditional transient stability analysis methods (i.e., time-domain simulation methods) require a large number of simulations under pre-selected fault scenarios to obtain conclusions on whether the system can recover to a stable state after a large disturbance, which is difficult to meet the speed requirements of online safety and stability assessment of modern power grids. The stability domain method developed from Lyapunov's stability theory describes the area in the state space that can converge to the equilibrium point as a whole from the perspective of the domain. It can not only quickly and intuitively give a comprehensive evaluation of the current operating status of the system, but also provide key information such as stability margin and optimization control direction to power grid personnel, and has good application prospects.

[0003] For large-scale power systems with complex dynamic characteristics, how to construct Lyapunov functions that meet the theorem conditions has always been a difficult problem that hinders the further development of stability domain methods. In the existing literature, the methods for constructing Lyapunov functions for nonlinear dynamic systems mainly include the Sum of Square method and the neural network method. Among them, the former is based on the Positivstellensatz theorem and the S-procedure theory, converting the conditions of positive definiteness of the function and negative definiteness of the Lie derivative that the Lyapunov function needs to satisfy into the sum of square conditions, and then using the semidefinite programming solver to solve it; while the latter uses a neural network with a specific structure to learn the Lyapunov function of the dynamic system from the operating data obtained by simulation or historical measurement, and introduces the SMT falsifier in the learning process to ensure the validity of the obtained function. The above two methods are essentially searching for analytical Lyapunov functions that meet the conditions in a specific function space (generally limited to polynomial function space). Therefore, under the condition of limited computing resources, it is difficult to achieve a good balance between computing speed and stability domain estimation accuracy, and cannot be effectively applied to the transient stability analysis problem of large-scale power systems. Summary of the invention

[0004] In order to solve at least one of the technical problems existing in the prior art to a certain extent, the object of the present invention is to provide a method and device for approximating and estimating the stability domain of a power system Lyapunov function based on Gaussian process regression.

[0005] The first technical solution adopted by the present invention is:

[0006] A method for approximating a Lyapunov function and estimating a stability region of a power system comprises the following steps:

[0007] According to the topological structure of the power system after the fault is removed, a mathematical model of the power system is established, and the stable equilibrium point of the power system after the fault is removed is solved;

[0008] Based on the established mathematical model of the power system, a set of initial sample points with stable labels are randomly generated near the stable equilibrium point, and the corresponding bounded inverse Lyapunov function values ​​are calculated;

[0009] The initial sample data set is constructed based on the initial sample points and the calculated bounded inverse Lyapunov function values;

[0010] Taking the initial sample data set as input, the Gaussian process regression model is used to learn the bounded inverse Lyapunov function of the power system;

[0011] Combined with the random direction search algorithm, the state point closest to the currently estimated stability domain boundary is searched on the randomly generated direction vector, and the corresponding bounded inverse Lyapunov function value is calculated; the searched state point and the corresponding bounded inverse Lyapunov function value are added to the sample data set as new samples and participate in the next iteration of model training;

[0012] When it is detected that the iteration stop condition is reached, the iteration process is stopped and the estimation result of the stable domain is output.

[0013] Furthermore, the establishment of a mathematical model of the power system includes:

[0014] Based on the transient and sub-transient operating parameters of each dynamic device in the power system, a mathematical model of the power system suitable for transient stability analysis is established; the mathematical model of the power system suitable for transient stability analysis is a set of differential algebraic equations.

[0015] Further, when only the second-order dynamic characteristics of the synchronous generator are considered, the power system mathematical model includes the rotor motion equation of the synchronous generator, the stator voltage equation of the synchronous generator and the power flow balance equation of the network;

[0016] The rotor motion equation of the synchronous generator is as follows:

[0017]

[0018] Where n g represents the number of synchronous generators in the system; δ i and ω i represent the power angle and electrical angular velocity of the i-th synchronous generator respectively; f 0 Indicates the base frequency of the system; P m,i and Pe,i Denote the input mechanical power and output electromagnetic power of the i-th synchronous generator respectively; D i represents the damping coefficient of the i-th synchronous generator; T J,i represents the inertia time constant of the i-th synchronous generator;

[0019] The stator voltage equation of the synchronous generator is as follows:

[0020]

[0021] Where U x,i and U y,i I represents the x-axis and y-axis terminal voltage of the i-th synchronous generator respectively; x,i and I y,i They represent the x-axis and y-axis output currents of the i-th synchronous generator respectively; E′ i represents the transient internal potential of the i-th synchronous generator; x′ d,i represents the d-axis transient reactance of the i-th synchronous generator;

[0022] The power flow balance equation of the network is as follows:

[0023]

[0024] In the formula, and They represent the nth g Line n g The real and imaginary parts of the corresponding elements in the column.

[0025] Further, the method of solving the stable equilibrium point of the power system after the fault is removed includes:

[0026] The analytical expressions of the terminal voltage and output current are solved from the stator voltage equation and power flow balance equation of the synchronous generator. They are functions of the power angle δ of all synchronous generators. Substituting them into the rotor motion equation of the synchronous generator, the following simplified mathematical model of the power system is obtained:

[0027]

[0028] All differential terms in the simplified mathematical model of the power system are set to zero, and then the Newton method is used to solve the stable equilibrium point of the power system after the fault is removed.

[0029] Furthermore, the expression of the bounded inverse Lyapunov function is:

[0030]

[0031] Where V(x) represents the bounded inverse Lyapunov function of the power system; x represents the n-dimensional state variable in the power system, such as δ and ω in the above-mentioned rotor motion equation; x eq represents the stable equilibrium point of the dynamic system; α() represents a Γ-type function; φ(t,x) represents the system operation trajectory calculated based on the mathematical model of the dynamic system with x as the initial state point, which is a non-analytical function about time t; R A Represents the actual stability domain of the power system. Among them, the power system is a power system, and the Lyapunov function can be applied to power systems and other power systems. Therefore, the power system belongs to the power system and can be analyzed using the bounded Lyapunov function.

[0032] Furthermore, constructing an initial sample data set according to the initial sample points and the calculated bounded inverse Lyapunov function values ​​includes:

[0033] A1. At the stable equilibrium point x eq A sample point x is randomly generated nearby (i+1) , and calculate its running trajectory φ(t,x (i+1) );

[0034] A2. If the sample point x (i+1) The running trajectory of s The internal energy returns to the stable equilibrium point x eq If it is near, the corresponding approximate bounded inverse Lyapunov function value is calculated using the following trapezoidal integral formula to update the initial sample data set; otherwise, return to step A1;

[0035]

[0036] A3. When the number of samples in the initial sample data set reaches the set maximum value N init When , the calculation ends; otherwise, return to step A1.

[0037] Furthermore, the method of learning the bounded inverse Lyapunov function of the power system using the Gaussian process regression model includes:

[0038] Assuming a bounded inverse Lyapunov function, let V(x) be a Gaussian process, and use the trapezoidal integral formula to calculate the approximate value is the observation of the Gaussian process, and the difference between the two is regarded as observation noise ε, and ε is assumed to follow a univariate Gaussian distribution have to:

[0039]

[0040] In the formula, represents the Gaussian process regression model, σ is the standard deviation of the univariate Gaussian distribution, which determines the width and discreteness of the observed noise distribution;

[0041] Given the mean function, covariance function, and the distribution of the observation noise, Gaussian process regression can be used to find the optimal solution in the observation sample. Constructing a data-driven model of bounded inverse Lyapunov function in

[0042] The Gaussian process regression model is trained, and the prediction formula of the bounded inverse Lyapunov function value of any set of state points in the state space by the trained Gaussian process regression model is:

[0043]

[0044] Where V * Indicates that in the prediction sample X * The bounded inverse Lyapunov function value predicted in ; represents the value of the approximate bounded inverse Lyapunov function observed in the observed sample X; K N Represents the covariance matrix between the observed sample and the predicted sample; K * Represents the covariance matrix between prediction samples; I is the identity matrix; N is the number of samples.

[0045] Furthermore, the random direction search algorithm is combined to search for a state point closest to the currently estimated stability domain boundary on the randomly generated direction vector, and calculate the corresponding bounded inverse Lyapunov function value; the searched state point and the corresponding bounded inverse Lyapunov function value are added to the sample data set as new samples, and participate in the next iteration of model training, including:

[0046] In each iteration, based on the trained Gaussian process regression model, the vector d is randomly generated. (i+1) Search for the state point x that is closest to the currently estimated stability region boundary (i+1) , and calculate the state point x (i+1) The corresponding approximate bounded inverse Lyapunov function value After that, it is added as a new sample point into the sample data set to participate in the next iteration of model training.

[0047] Furthermore, the estimation formula of the stability domain is:

[0048]

[0049] Where E(x) and Var(x) represent the mean and variance of the Gaussian distribution that point x obeys, respectively; Φ -1 Represents a standard Gaussian distribution The inverse cumulative distribution function; ρ represents the confidence level, which is used to measure the estimated state point in the stable domain at the set t s It can return to the stable equilibrium point x within time eq The probability near ; c represents the level set parameter, which is used to measure the size of the estimated stable domain.

[0050] The second technical solution adopted by the present invention is:

[0051] An electronic device comprises a processor and a memory, wherein the memory stores at least one instruction, at least one program, a code set or an instruction set, and the at least one instruction, the at least one program, the code set or the instruction set is loaded and executed by the processor to implement a method for approximating a Lyapunov function and estimating a stability domain of a power system as described above.

[0052] The third technical solution adopted by the present invention is:

[0053] A computer-readable storage medium stores at least one instruction, at least one program, a code set or an instruction set, wherein the at least one instruction, the at least one program, the code set or the instruction set is loaded and executed by a processor to implement a method for approximating a Lyapunov function and estimating a stability domain of a power system as described above.

[0054] The fourth technical solution adopted by the present invention is:

[0055] A computer program product or a computer program includes computer instructions stored in a computer-readable storage medium. A processor of a computer device can read the computer instructions from the computer-readable storage medium, and the processor executes the computer instructions, so that the computer device executes the above method.

[0056] Compared with the prior art, the present invention has the following advantages and effects:

[0057] (1) Compared with the analytical Lyapunov function construction and stability domain estimation methods commonly used in existing literature based on sum of squares and neural networks, the time complexity of the method of the present invention does not increase with the increase of the dimension of the power system state variables, and there is no need to restrict the specific form of the constructed function (that is, a non-analytical Lyapunov function can be obtained). Therefore, the method of the present invention has a faster calculation speed and higher stability domain estimation accuracy.

[0058] 2. The bounded inverse Lyapunov function proposed in the present invention is defined in the entire state space, while the traditional inverse Lyapunov function takes positive infinity at state points outside the actual stable domain. Therefore, the former is more suitable for training and prediction of data-driven models than the latter.

[0059] 3. The present invention introduces an active learning strategy based on random direction search, which allows the Gaussian process regression model to complete the sampling of key sample points in the state space while learning the bounded inverse Lyapunov function of the power system. Compared with the traditional random sampling, the estimation results of the stable domain can be quickly corrected.

[0060] 4. The present invention is not only applicable to the field of transient stability analysis of power systems, but also to the field of stability analysis of other nonlinear dynamic systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the embodiments of the present invention or the drawings of related technical solutions in the prior art are introduced below. It should be understood that the drawings introduced below are only for the convenience of clearly describing some embodiments of the technical solutions of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative work.

[0062] Figure 1 It is a flowchart of a method for Lyapunov function approximation and stability region estimation of a power system in Embodiment 1 of the present invention;

[0063] Figure 2 is a flowchart of generating an initial sample data set in step 2 of embodiment 2 of the present invention;

[0064] Figure 3 It is a flowchart of the method for approximating the Lyapunov function and estimating the stability domain of a dynamic system based on Gaussian process regression in step 3 of embodiment 3 of the present invention;

[0065] Figure 4 The stability domain result estimated by the Lyapunov function approximation and stability domain estimation method of the power system based on Gaussian process regression proposed by the present invention is in δ 1 -δ 2 Schematic diagram of projection in state space;

[0066] Figure 5 The stability domain result is estimated by using the analytical Lyapunov function construction based on the sum of squares and the stability domain estimation method in δ 1 -δ 2 Schematic diagram of projection in state space;

[0067] Figure 6 The stability domain result is estimated by using the analytical Lyapunov function construction based on neural network and the stability domain estimation method in δ 1 -δ 2 Schematic diagram of the projection in state space. DETAILED DESCRIPTION

[0068] The embodiments of the present invention are described in detail below, and examples of the embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and are not to be construed as limitations of the present invention. For the step numbers in the following embodiments, they are only provided for the convenience of explanation, and the order between the steps is not limited in any way, and the execution order of each step in the embodiment can be adaptively adjusted according to the understanding of those skilled in the art.

[0069] In the description of the present invention, it should be understood that descriptions involving orientations, such as up, down, front, back, left, right, etc., and orientations or positional relationships indicated are based on the orientations or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be understood as a limitation on the present invention.

[0070] In the description of the present invention, "several" means one or more, "more" means more than two, "greater than", "less than", "exceed" etc. are understood as not including the number itself, and "above", "below", "within" etc. are understood as including the number itself. If there is a description of "first" or "second", it is only used for the purpose of distinguishing the technical features, and cannot be understood as indicating or implying the relative importance or implicitly indicating the number of the indicated technical features or implicitly indicating the order of the indicated technical features.

[0071] In the description of the present invention, unless otherwise clearly defined, terms such as setting, installing, connecting, etc. should be understood in a broad sense, and technicians in the relevant technical field can reasonably determine the specific meanings of the above terms in the present invention based on the specific content of the technical solution.

[0072] Aiming at the stability domain estimation problem of large-scale power systems with complex dynamic characteristics, the present invention proposes a power system Lyapunov function approximation and stability domain estimation scheme based on Gaussian process regression. Compared with the analytical Lyapunov function construction and stability domain estimation methods based on Sum of Square and neural network commonly used in existing literature, the time complexity of this scheme does not increase with the increase of the dimension of the power system state variables, and there is no need to restrict the specific form of the constructed function (that is, a non-analytical Lyapunov function can be obtained), so it has a faster calculation speed and higher stability domain estimation accuracy.

[0073] Example 1

[0074] like Figure 1As shown, this embodiment provides a method for approximating a Lyapunov function and estimating a stability region of a power system, comprising the following steps:

[0075] S1. According to the topological structure of the power system after the fault is removed, a mathematical model of the power system is established, and the stable equilibrium point of the power system after the fault is removed is solved.

[0076] Specifically, a mathematical model of the power system suitable for transient stability analysis is established based on the topological structure of the power system after fault removal and the transient and sub-transient operating parameters of each dynamic device in the system, and the stable equilibrium point of the power system after fault removal is solved by setting the differential terms in the model to zero.

[0077] In some embodiments, the power system mathematical model suitable for transient stability analysis is essentially a set of differential algebraic equations. Taking the second-order dynamic characteristics of the synchronous generator as an example, it specifically includes the rotor motion equation of the synchronous generator, the stator voltage equation of the synchronous generator, and the power flow balance equation of the network.

[0078] The rotor motion equation of the synchronous generator is as follows:

[0079]

[0080] Where n g represents the number of synchronous generators in the system; δ i and ω i represent the power angle and electrical angular velocity of the i-th synchronous generator respectively; f 0 Indicates the base frequency of the system; P m,i and P e,i Denote the input mechanical power and output electromagnetic power of the i-th synchronous generator respectively; D i represents the damping coefficient of the i-th synchronous generator; T J,i represents the inertia time constant of the i-th synchronous generator.

[0081] The stator voltage equation of the synchronous generator is as follows:

[0082]

[0083] Where U x,i and U y,i I represents the x-axis and y-axis terminal voltage of the i-th synchronous generator respectively; x,i and I y,i They represent the x-axis and y-axis output currents of the i-th synchronous generator respectively; E′ i represents the transient internal potential of the i-th synchronous generator; x′ d,i represents the d-axis transient reactance of the i-th synchronous generator.

[0084] The power balance equation of the network is as follows:

[0085]

[0086] In the formula, and They represent the nth g Line n g The real and imaginary parts of the corresponding elements in the column.

[0087] When only the second-order dynamic characteristics of the synchronous generator are considered, the analytical expressions of the terminal voltage and output current can be further solved from the stator voltage equation and the power flow balance equation. They are functions of the power angle δ of all synchronous generators. Substituting them into the rotor motion equation, the following simplified mathematical model of the power system can be obtained:

[0088]

[0089] All differential terms in the simplified mathematical model of the power system are set to zero, and then the Newton method is used to solve the stable equilibrium point of the power system after the fault is removed.

[0090] S2. Based on the mathematical model of the power system after fault removal obtained in step S1, a set of initial sample points with stable labels are randomly generated near the stable equilibrium point, and the corresponding bounded inverse Lyapunov function value is calculated.

[0091] The bounded inverse Lyapunov function is obtained by performing nonlinear transformation on the basis of the traditional inverse Lyapunov function. Its specific expression is:

[0092]

[0093] Where V(x) represents the bounded inverse Lyapunov function of the dynamical system; represents the n-dimensional state variables in the power system, such as δ and ω in the above rotor motion equation; x eq represents the stable equilibrium point of the dynamic system; α(·):[0,∞)→[0,∞) represents the Γ-type function; φ(t,x) represents the system operation trajectory calculated based on the mathematical model of the dynamic system with x as the initial state point, which is a non-analytical function about time t; R A Represents the actual stability region of the dynamical system.

[0094] Among them, the Γ-type function α(·) refers to any continuous function that satisfies the following three conditions:

[0095]

[0096] The actual stability region R of the dynamic system AThe specific definition is:

[0097]

[0098] The bounded inverse Lyapunov function V(x) can satisfy the three conditions of the Lyapunov function:

[0099] V(x eq )=0

[0100]

[0101] For the state points within the actual stable domain, the value of the bounded inverse Lyapunov function is less than 1, while for the state points outside the actual stable domain, the value of the bounded inverse Lyapunov function is equal to 1. Therefore, based on the bounded inverse Lyapunov function V(x), the stable domain of the dynamic system can be estimated:

[0102]

[0103] S3. Construct an initial sample data set based on the initial sample points and the calculated bounded inverse Lyapunov function value.

[0104] As an optional implementation, the specific process of generating the initial sample data set is as follows:

[0105] S31, at the stable equilibrium point x eq A random sample point is generated nearby Based on the mathematical model of the power system, the running trajectory φ(t,x (i+1) ).

[0106] S32, if the sample point x (i+1) The running trajectory of s The internal energy returns to the stable equilibrium point x eq Nearby, that is, ||φ(t s ,x (i+1) )-x eq ||<ξ, then use the following trapezoidal integral formula to calculate the corresponding approximate bounded inverse Lyapunov function value and update the initial sample data set; otherwise, return to step S31.

[0107]

[0108] S33, when the number of samples in the initial sample data set reaches the set maximum value N init , end the calculation; otherwise, return to step S31.

[0109] S4. Taking the initial sample data set as input, the Gaussian process regression model is used to learn the bounded inverse Lyapunov function of the power system.

[0110] Gaussian process regression is essentially a data-driven model that uses Gaussian process to solve regression problems. The Gaussian process is specifically defined as any finite set of points {x 1 ,x 2 ,…,x N} have random variables {f(x 1 ),f(x 2 ),…,f(x N )} is a random process that obeys an N-dimensional Gaussian distribution. In the function space perspective, it can be represented by a mean function m(x) and a covariance function k(x,x′):

[0111]

[0112] Using the Gaussian process regression model to learn the bounded inverse Lyapunov function of the dynamical system, we first need to assume that the bounded inverse Lyapunov function is a Gaussian process, and the approximate value calculated using the trapezoidal integral formula is is the observation of the Gaussian process, and the difference between the two is regarded as observation noise ε, and ε is assumed to follow a univariate Gaussian distribution have to:

[0113]

[0114] Given the mean function (generally fixed to zero), the covariance function, and the distribution of the observation noise, Gaussian process regression can be used to find the optimal value for the observed sample. A data-driven model of the bounded inverse Lyapunov function is constructed in , and further, the value of the bounded inverse Lyapunov function at any point in the state space is predicted.

[0115] The training of Gaussian process regression model essentially refers to learning the mean function, covariance function and hyperparameters in the observation noise in the data, which generally needs to be solved by maximum likelihood estimation and conjugate gradient method. When the mean function is fixed to zero, the calculation formula of the log-likelihood function and its gradient to each hyperparameter is as follows:

[0116]

[0117] in, Indicates that in the observed sample X = [x 1 ,x 2 ,…,x N ] is the approximate bounded inverse Lyapunov function value observed in ; K represents the covariance matrix between the observed samples; θ represents the hyperparameter to be learned in the covariance function.

[0118] After obtaining the trained Gaussian process regression model, the prediction formula for the bounded inverse Lyapunov function value of any set of state points in the state space is:

[0119]

[0120] Among them, V * =[V *1 ,V *2 ,…,V *M ] T Indicates that in the prediction sample X * =[x *1 ,x *2 ,…,x *M ] is the value of the bounded inverse Lyapunov function predicted in N Represents the covariance matrix between the observed samples and the predicted samples; K * Represents the covariance matrix between prediction samples.

[0121] Furthermore, the estimation formula of the stability domain is:

[0122]

[0123] Where E(x) and Var(x) represent the mean and variance of the Gaussian distribution that point x obeys respectively; Φ -1 (·) represents standard Gaussian distribution The inverse cumulative distribution function; ρ represents the confidence level, which is used to measure the estimated state point in the stable domain at the set t s It can return to the stable equilibrium point x within time eq The probability near ; c represents the level set parameter, which is used to measure the size of the estimated stable domain.

[0124] S5. Combined with the random direction search algorithm, search for the state point closest to the currently estimated stability domain boundary on the randomly generated direction vector, and calculate the corresponding bounded inverse Lyapunov function value; add the searched state point and the corresponding bounded inverse Lyapunov function value as new samples to the sample data set and participate in the next iteration of model training.

[0125] The time complexity of training the Gaussian process regression model is closely related to the number of training samples. In the case of random sampling, the training efficiency and generalization ability are relatively poor, so an active learning strategy needs to be introduced. Specifically, in each iteration, based on the trained Gaussian process regression model, the random direction vector d is generated. (i+1) Search for the state point x that is closest to the current estimated stability region boundary (i+1) , and calculate the corresponding approximate bounded inverse Lyapunov function value After that, it is added as a new sample point into the sample data set to participate in the next iteration of model training.

[0126]

[0127] x (i+1) =β·d (i+1)

[0128] S6. When it is detected that the iteration stop condition is reached, the iteration process is stopped and the estimation result of the stable domain is output.

[0129] As an optional implementation, the estimation result of the stable domain is obtained by the following steps:

[0130] S61. Train a Gaussian process regression model based on the current sample data set.

[0131] S62, using the trained Gaussian process regression model to determine whether the state points in the validation set are within the estimated stable domain. If the model can completely identify all state points with unstable labels in the validation set, the iterative process ends and the estimated result of the stable domain is output.

[0132] S63. If the number of samples in the sample data set reaches the set maximum value N max , the iteration process ends and the estimation result of the stable domain is output.

[0133] S64. Randomly generate a direction vector d in the state space (i+1) , and then use the trained Gaussian process regression model to search for the state point x in this direction that is closest to the currently estimated stability domain boundary (i+1) , and calculate its running trajectory φ(t,x (i+1) ) and the approximate bounded inverse Lyapunov function value

[0134] S65: The new state point x obtained in step S64 is (i+1) and the corresponding approximately bounded inverse Lyapunov function value Add it to the sample data set and return to step S61.

[0135] Example 2

[0136] This embodiment provides a method for approximating the Lyapunov function and estimating the stability region of a power system based on Gaussian process regression, which specifically includes the following steps:

[0137] Step 1: Based on the topological structure of the power system after fault removal and the transient and sub-transient operating parameters of each dynamic device in the system, a mathematical model of the power system suitable for transient stability analysis is established, and the stable equilibrium point of the power system after fault removal is solved by setting the differential terms in the model to zero.

[0138] The mathematical model of the power system suitable for transient stability analysis is essentially a set of differential algebraic equations. Taking the second-order dynamic characteristics of the synchronous generator as an example, it specifically includes the rotor motion equation of the synchronous generator, the stator voltage equation of the synchronous generator, and the power flow balance equation of the network.

[0139] The rotor motion equation of the synchronous generator is as follows:

[0140]

[0141] Among them, n g represents the number of synchronous generators in the system; δ i and ω i represent the power angle and electrical angular velocity of the i-th synchronous generator respectively; f 0 Indicates the base frequency of the system; P m,i and P e,i Denote the input mechanical power and output electromagnetic power of the i-th synchronous generator respectively; D i represents the damping coefficient of the i-th synchronous generator; T J,i represents the inertia time constant of the i-th synchronous generator.

[0142] The stator voltage equation of the synchronous generator is as follows:

[0143]

[0144] Among them, U x,i and U y,i I represents the x-axis and y-axis terminal voltage of the i-th synchronous generator respectively; x,i and I y,i They represent the x-axis and y-axis output currents of the i-th synchronous generator respectively; E i ′ represents the transient internal potential of the i-th synchronous generator; x′ d,i represents the d-axis transient reactance of the i-th synchronous generator.

[0145] The power balance equation of the network is as follows:

[0146]

[0147] Among them, G ij and B ij They represent the real and imaginary parts of the element corresponding to the i-th row and j-th column in the node admittance matrix.

[0148] When only the second-order dynamic characteristics of the synchronous generator are considered, the analytical expressions of the terminal voltage and output current can be further solved from the stator voltage equation and the power flow balance equation. They are functions of the power angle δ of all synchronous generators. Substituting them into the rotor motion equation, the following simplified mathematical model of the power system can be obtained:

[0149]

[0150] In the three-machine system used in this embodiment, the specific expression is:

[0151]

[0152] All differential terms in the above model are set to zero, and then the Newton method is used to solve the stable equilibrium point of the power system after the fault is removed. In the three-machine system used in this embodiment, it is specifically:

[0153] x eq =(δ 1 ,ω 1 ,δ 2 ,ω 2 )=(0.02,0,0.06,0)

[0154] Step 2: Based on the mathematical model of the power system after fault removal obtained in step 1, a set of initial sample points with stable labels are randomly generated near the stable equilibrium point, and the corresponding bounded inverse Lyapunov function values ​​are calculated. Together, they constitute the initial sample data set required for subsequent learning of the bounded inverse Lyapunov function of the power system.

[0155] (1) The bounded inverse Lyapunov function is obtained by performing nonlinear transformation on the traditional inverse Lyapunov function. Its specific expression is:

[0156]

[0157] Where V(x) represents the bounded inverse Lyapunov function of the dynamical system; represents the n-dimensional state variables in the power system, such as δ and ω in the above rotor motion equation; x eq represents the stable equilibrium point of the dynamic system; α(·):[0,∞)→[0,∞) represents the Γ-type function; φ(t,x) represents the system operation trajectory calculated based on the mathematical model of the dynamic system with x as the initial state point, which is a non-analytical function about time t; R A Represents the actual stability region of the dynamical system.

[0158] (2) A Γ-type function α(·) is any continuous function that satisfies the following three conditions:

[0159]

[0160] In this embodiment, the specific expression is:

[0161]

[0162] (3) Actual stability region R of the dynamic system A The specific definition is:

[0163]

[0164] (4) The bounded inverse Lyapunov function V(x) can satisfy the three conditions of the Lyapunov function:

[0165] V(x eq )=0

[0166]

[0167] (5) For the state points within the actual stable domain, the value of the bounded inverse Lyapunov function is less than 1, while for the state points outside the actual stable domain, the value of the bounded inverse Lyapunov function is equal to 1. Therefore, based on the bounded inverse Lyapunov function V(x), the stable domain of the dynamic system can be estimated:

[0168]

[0169] As an optional implementation, see Figure 2 , the steps to generate the initial sample data set specifically include:

[0170] Step 2.1: At the stable equilibrium point x eq A random sample point is generated nearby Based on the mathematical model of the power system, the running trajectory φ(t,x (i+1) ). In this embodiment, the specific sampling space is x∈[-5,5], and the initial sampling space is one tenth of the sampling space, that is, x∈[-0.5,0.5];

[0171] Step 2.2: If its running trajectory is at the set time t s The internal energy returns to the stable equilibrium point x eq Nearby, that is, ||φ(t s ,x (i+1) )-x eq ||<ξ, then use the following trapezoidal integral formula to calculate the corresponding approximate bounded inverse Lyapunov function value, otherwise return to step 2.1. In this embodiment, specifically take t s =100s,ξ=10 -4 ;

[0172]

[0173] Step 2.3: When the number of samples in the initial sample data set reaches the set maximum value N init The process ends when N is reached, otherwise, the process returns to step 2.1. init =10, the initial samples obtained are shown in Table 1.

[0174] Table 1 Initial sample

[0175]

[0176] Step 3: Take the initial sample data set obtained in step 2 as input. In each iteration, use the Gaussian process regression model to learn the bounded inverse Lyapunov function of the power system. Then, combine the random direction search algorithm to search for the state point closest to the currently estimated stability domain boundary on the randomly generated direction vector, and calculate the corresponding bounded inverse Lyapunov function value. This state point and the corresponding bounded inverse Lyapunov function value are added to the sample data set as new samples to participate in the model training of the next iteration. When the stability domain estimated by the Gaussian process regression model can fully identify all sample points with unstable labels in the validation set or the number of training samples reaches the set maximum value, stop the above iteration process and output the estimation result of the stability domain.

[0177] (1) Gaussian process regression is essentially a data-driven model that uses Gaussian processes to solve regression problems. The Gaussian process is specifically defined as any finite set of points {x 1 ,x 2 ,…,x N} have random variables {f(x 1 ),f(x 2 ),…,f(x N )} is a random process that obeys an N-dimensional Gaussian distribution. In the function space perspective, it can be represented by a mean function m(x) and a covariance function k(x,x′):

[0178]

[0179] (2) Using the Gaussian process regression model to learn the bounded inverse Lyapunov function of the dynamical system, we first need to assume that V(x) is a Gaussian process, and the approximate value calculated using the trapezoidal integral formula is is the observation of the Gaussian process, and the difference between the two is regarded as observation noise ε and is assumed to follow a univariate Gaussian distribution

[0180]

[0181] (3) Given a mean function (generally fixed to zero), a covariance function, and the distribution of observation noise, Gaussian process regression can be used to find the optimal solution for the observation sample. A data-driven model of the bounded inverse Lyapunov function is constructed, and further, the value of the bounded inverse Lyapunov function at any point in the state space is predicted. In this embodiment, the mean function m(x)≡0 is taken, and the covariance function k(x,x′) is the radial basis covariance function:

[0182]

[0183] Where s represents the signal variance; L represents the length scale matrix, and except for the diagonal elements, the remaining elements are zero. The diagonal elements in s and L together constitute the hyperparameters that need to be trained in the covariance function.

[0184] (4) The training of the Gaussian process regression model essentially refers to learning the mean function, covariance function, and hyperparameters in the observation noise in the data, which generally needs to be solved by maximum likelihood estimation and conjugate gradient method. When the mean function is fixed to zero, the calculation formula of the log-likelihood function and its gradient to each hyperparameter is as follows:

[0185]

[0186] in, Indicates that in the observed sample X = [x 1 ,x 2 ,…,x N ] is the approximate bounded inverse Lyapunov function value observed in ; K represents the covariance matrix between the observed samples; θ represents the hyperparameter to be learned in the covariance function.

[0187] (5) After obtaining the trained Gaussian process regression model, the prediction formula for the bounded inverse Lyapunov function value of any set of state points in the state space is:

[0188]

[0189] Among them, V * =[V *1 ,V *2 ,…,V *M ] T Indicates that in the prediction sample X * =[x *1 ,x *2 ,…,x *M ] is the value of the bounded inverse Lyapunov function predicted in N Represents the covariance matrix between the observed sample and the predicted sample; K* Represents the covariance matrix between prediction samples.

[0190] (6) Furthermore, the estimation formula of the stability domain is:

[0191]

[0192] Where E(x) and Var(x) represent the mean and variance of the Gaussian distribution that point x obeys respectively; Φ -1 (·) represents standard Gaussian distribution The inverse cumulative distribution function; ρ represents the confidence level, which is used to measure the estimated state point in the stable domain at the set t s It can return to the stable equilibrium point x within time eq The probability of being near; c represents the level set parameter, which is used to measure the size of the estimated stability domain. In this embodiment, ρ = 99.99%, c = 0.95. At this time, the standard Gaussian distribution The inverse cumulative distribution function of is:

[0193]

[0194] (7) The time complexity of training the Gaussian process regression model is closely related to the number of training samples. In the case of random sampling, the training efficiency and generalization ability are relatively poor, so it is necessary to introduce an active learning strategy, which is as follows: in each iteration, based on the trained Gaussian process regression model, the randomly generated direction vector d (i+1) Search for the state point x that is closest to the currently estimated stability region boundary (i+1) , and calculate the corresponding approximate bounded inverse Lyapunov function value After that, it is added as a new sample point into the sample data set to participate in the next iteration of model training.

[0195]

[0196] x (i+1) =β·d (i+1)

[0197] As an optional implementation, participate Figure 3 , the steps of the Lyapunov function approximation and stability region estimation method for dynamic systems based on Gaussian process regression are as follows:

[0198] Step 3.1: Train the Gaussian process regression model based on the current sample data set;

[0199] Step 3.2: Use the trained Gaussian process regression model to determine whether the state points in the validation set are within the estimated stable domain. If the model can fully identify all state points with unstable labels in the validation set, the iterative process ends and the estimation result of the stable domain is output. In this embodiment, the validation set is randomly generated in the sampling space x∈[-5,5], with a total of 10,000 samples, including 1,836 stable samples and 8,164 unstable samples;

[0200] Step 3.3: If the number of samples in the sample data set reaches the set maximum value N max , the iteration process ends and the estimation result of the stable domain is output. max =1000.

[0201] Step 3.4: Randomly generate a direction vector d in the state space (i+1) , and then use the trained Gaussian process regression model to search for the state point x in this direction that is closest to the currently estimated stability domain boundary (i+1) , and calculate its running trajectory φ(t,x (i+1) ) and the approximate bounded inverse Lyapunov function value

[0202] Step 3.5: The new state point x obtained in step 3.4 (i+1) and the corresponding approximately bounded inverse Lyapunov function value Add to the sample dataset and return to step 3.1.

[0203] In this embodiment, when the number of samples N reaches 300, the above iterative process ends. At this time, the Gaussian process regression model obtained can completely identify all state points with unstable labels in the validation set. The estimated results of the stable domain are projected onto δ 1 -δ 2 In the state space, Figure 4 As shown. Among them, the actual stable region is in δ 1 -δ 2 The projection in state space is obtained by scanning the grid points.

[0204] The proposed method is compared with the analytical Lyapunov function construction and stability domain estimation method based on Sum of Square and neural network. The stability domain estimation results are 1 -δ 2 The projection in state space is Figure 5 and Figure 6 shown.

[0205] Furthermore, 10,000 test samples are randomly generated in the sampling space x∈[-5,5], of which 1,860 are stable samples and 8,140 are unstable samples. The training time, test time, and stability judgment accuracy of different methods are shown in Table 2.

[0206] Table 2 Comparison of different methods

[0207] method Training time / s Test time / s Accuracy / % Recall / % Based on Gaussian process regression 100.2122 0.1094 86.23 100 Based on Sum of Square 12025.5355 0.0042 82.61 100 Based on neural network 6231.0711 0.0052 82.66 100

[0208] contrast Figure 4 , Figure 5 and Figure 6 As can be seen from the data in Table 2, the power system Lyapunov function approximation and stability domain estimation method based on Gaussian process regression proposed in the present invention has faster calculation speed and higher stability domain estimation accuracy.

[0209] Example 3

[0210] An embodiment of the present invention further provides an electronic device, the electronic device comprising a processor and a memory, the memory storing at least one instruction, at least one program, a code set or an instruction set, the at least one instruction, the at least one program, the code set or the instruction set being loaded and executed by the processor to implement the following Figure 1 A method for approximating the Lyapunov function and estimating the stability region of a power system is shown.

[0211] It is understood that the memory may include a random access memory (RAM) or a read-only memory (ROM). Optionally, the memory includes a non-transitory computer-readable storage medium. The memory may be used to store instructions, programs, codes, code sets, or instruction sets. The memory may include a program storage area and a data storage area, wherein the program storage area may store instructions for implementing an operating system, instructions for at least one function, instructions for implementing the above-mentioned various method embodiments, etc.; the data storage area may store data created according to the use of the server, etc.

[0212] The processor may include one or more processing cores. The processor uses various interfaces and lines to connect the various parts of the entire server, and executes various functions of the server and processes data by running or executing instructions, programs, code sets or instruction sets stored in the memory, and calling data stored in the memory. Optionally, the processor can be implemented in at least one hardware form of digital signal processing (DSP), field programmable gate array (FPGA), and programmable logic array (PLA). The processor can integrate one or a combination of a central processing unit (CPU) and a modem. Among them, the CPU mainly processes the operating system and application programs; the modem is used to process wireless communications. It can be understood that the above-mentioned modem may not be integrated into the processor, but implemented separately through a chip.

[0213] Since the electronic device is an electronic device corresponding to the Lyapunov function approximation and stability domain estimation method of an electric power system in an embodiment of the present invention, and the principle of solving the problem by the electronic device is similar to that of the method, the implementation of the electronic device can refer to the implementation process of the above-mentioned method embodiment, and the repeated parts will not be repeated.

[0214] Example 4

[0215] The embodiment of the present invention further provides a computer-readable storage medium, wherein the storage medium stores at least one instruction, at least one program, code set or instruction set, and the at least one instruction, the at least one program, the code set or instruction set is loaded and executed by a processor to implement the following Figure 1 A method for approximating the Lyapunov function and estimating the stability region of a power system is shown.

[0216] Those skilled in the art can understand that all or part of the steps in the various methods of the above embodiments can be completed by instructing related hardware through a program, and the program can be stored in a computer-readable storage medium, and the storage medium includes a read-only memory (ROM), a random access memory (RAM), a programmable read-only memory (PROM), an erasable programmable read-only memory (EPROM), a one-time programmable read-only memory (OTPROM), an electronically erasable rewritable read-only memory (EEPROM), a compact disc (CD-ROM) or other optical disc storage, magnetic disk storage, magnetic tape storage, or any other computer-readable medium that can be used to carry or store data.

[0217] Since the storage medium is a storage medium corresponding to the Lyapunov function approximation and stability domain estimation method of an electric power system in an embodiment of the present invention, and the principle of solving the problem by the storage medium is similar to that of the method, the implementation of the storage medium can refer to the implementation process of the above-mentioned method embodiment, and the repeated parts will not be repeated.

[0218] Example 5

[0219] In some possible implementations, various aspects of the method of the embodiment of the present invention may also be implemented in the form of a program product, which includes a program code. When the program product is run on a computer device, the program code is used to enable the computer device to execute the steps of the method for approximating the Lyapunov function and estimating the stability domain of a power system according to various exemplary embodiments of the present application described above in this specification. Among them, the executable computer program code or "code" for executing each embodiment can be written in a high-level programming language such as C, C++, C#, Smalltalk, Java, JavaScript, Visual Basic, structured query language (e.g., Transact-SQL), Perl, or in various other programming languages.

[0220] It should be understood that the various parts of the present invention can be implemented by hardware, software, firmware or a combination thereof. In the above-mentioned embodiments, a plurality of steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented by hardware, as in another embodiment, it can be implemented by any one of the following technologies known in the art or their combination: a discrete logic circuit having a logic gate circuit for implementing a logic function for a data signal, a dedicated integrated circuit having a suitable combination of logic gate circuits, a programmable gate array (PGA), a field programmable gate array (FPGA), etc.

[0221] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of the different embodiments or examples, without contradiction.

[0222] The above embodiments are only for illustrating the technical concept and features of the present invention, and their purpose is to enable ordinary technicians in the field to understand the content of the present invention and implement it accordingly, and they cannot be used to limit the protection scope of the present invention. Any equivalent changes or modifications made based on the essence of the content of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for approximating Lyapunov functions and estimating stability regions of power systems, characterized in that: The following steps are involved: According to the topological structure of the power system after the fault is removed, a mathematical model of the power system is established, and the stable equilibrium point of the power system after the fault is removed is solved; Based on the established mathematical model of the power system, a set of initial sample points with stable labels are randomly generated near the stable equilibrium point, and the corresponding bounded inverse Lyapunov function values ​​are calculated; The initial sample data set is constructed based on the initial sample points and the calculated bounded inverse Lyapunov function values; Taking the initial sample data set as input, the Gaussian process regression model is used to learn the bounded inverse Lyapunov function of the power system; Combined with the random direction search algorithm, the state point closest to the currently estimated stability domain boundary is searched on the randomly generated direction vector, and the corresponding bounded inverse Lyapunov function value is calculated; the searched state point and the corresponding bounded inverse Lyapunov function value are added to the sample data set as new samples and participate in the next iteration of model training; When it is detected that the iteration stop condition is reached, the iteration process is stopped and the estimation result of the stable domain is output.

2. A method for approximating the Lyapunov function and estimating the stability region of a power system according to claim 1, characterized in that: The establishing of the power system mathematical model comprises: Based on the transient and sub-transient operating parameters of each dynamic device in the power system, a mathematical model of the power system suitable for transient stability analysis is established; the mathematical model of the power system suitable for transient stability analysis is a set of differential algebraic equations.

3. A method for approximating the Lyapunov function and estimating the stability region of a power system according to claim 2, characterized in that: When only the second-order dynamic characteristics of the synchronous generator are considered, the power system mathematical model includes the rotor motion equation of the synchronous generator, the stator voltage equation of the synchronous generator and the power flow balance equation of the network; The rotor motion equation of the synchronous generator is as follows: Where n g represents the number of synchronous generators in the system; δ i and ω i represent the power angle and electrical angular velocity of the ith synchronous generator respectively; f0 represents the reference frequency of the system; P m,i and P e,i Denote the input mechanical power and output electromagnetic power of the i-th synchronous generator respectively; D i represents the damping coefficient of the i-th synchronous generator; T J,i represents the inertia time constant of the i-th synchronous generator; The stator voltage equation of the synchronous generator is as follows: Where U x,i and U y,i I represents the x-axis and y-axis terminal voltage of the i-th synchronous generator respectively; x,i and I y,i They represent the x-axis and y-axis output currents of the i-th synchronous generator respectively; E′ i represents the transient internal potential of the i-th synchronous generator; x′ d,i represents the d-axis transient reactance of the i-th synchronous generator; The power flow balance equation of the network is as follows: In the formula, and They represent the nth g Line n g The real and imaginary parts of the corresponding elements in the column.

4. A method for approximating the Lyapunov function and estimating the stability region of a power system according to claim 3, characterized in that: The method of solving the stable equilibrium point of the power system after the fault is removed includes: The analytical expressions of the terminal voltage and output current are solved from the stator voltage equation and power flow balance equation of the synchronous generator. They are functions of the power angle δ of all synchronous generators. Substituting them into the rotor motion equation of the synchronous generator, the following simplified mathematical model of the power system is obtained: All differential terms in the simplified mathematical model of the power system are set to zero, and then the Newton method is used to solve the stable equilibrium point of the power system after the fault is removed.

5. The method for Lyapunov function approximation and stability region estimation of a power system according to claim 1, characterized in that: The expression of the bounded inverse Lyapunov function is: Where V(x) represents the bounded inverse Lyapunov function of the dynamical system; x represents the n-dimensional state variable in the dynamical system; x eq represents the stable equilibrium point of the dynamic system; α() represents a Γ-type function; φ(t,x) represents the system operation trajectory calculated based on the mathematical model of the dynamic system with x as the initial state point, which is a non-analytical function of time r; R A Represents the actual stability region of the dynamical system.

6. A method for approximating the Lyapunov function and estimating the stability region of a power system according to claim 5, characterized in that: The initial sample data set is constructed based on the initial sample points and the calculated bounded inverse Lyapunov function values, including: A1. At the stable equilibrium point x eq A sample point x is randomly generated nearby (i+1) , and calculate its running trajectory φ(t,x (i+1) ); A2. If the sample point x (i+1) The running trajectory of s The internal energy returns to the stable equilibrium point x eq If it is near, the corresponding approximate bounded inverse Lyapunov function value is calculated using the following trapezoidal integral formula to update the initial sample data set; otherwise, return to step A1; A3. When the number of samples in the initial sample data set reaches the set maximum value N init When , the calculation ends; otherwise, return to step A1.

7. A method for approximating the Lyapunov function and estimating the stability region of a power system according to claim 1, characterized in that: The method of learning the bounded inverse Lyapunov function of the power system by using the Gaussian process regression model includes: Assuming a bounded inverse Lyapunov function, let V(x) be a Gaussian process, and use the trapezoidal integral formula to calculate the approximate value is the observation of the Gaussian process, and the difference between the two is regarded as observation noise ε, and ε is assumed to follow a univariate Gaussian distribution have to: In the formula, represents the Gaussian process regression model, σ is the standard deviation of the univariate Gaussian distribution; Given the mean function, covariance function, and the distribution of the observation noise, Gaussian process regression can be used to find the optimal solution in the observation sample. Constructing a data-driven model of bounded inverse Lyapunov function in The Gaussian process regression model is trained, and the prediction formula of the bounded inverse Lyapunov function value of any set of state points in the state space by the trained Gaussian process regression model is: Where V * Indicates that in the prediction sample X * The bounded inverse Lyapunov function value predicted in ; represents the value of the approximate bounded inverse Lyapunov function observed in the observed sample X; K N Represents the covariance matrix between the observed sample and the predicted sample; K * Represents the covariance matrix between prediction samples; I is the identity matrix; N is the number of samples.

8. A method for approximating the Lyapunov function and estimating the stability region of a power system according to claim 1, characterized in that: The combined random direction search algorithm searches for the state point closest to the currently estimated stability region boundary on the randomly generated direction vector, and calculates the corresponding bounded inverse Lyapunov function value; The searched state points and the corresponding bounded inverse Lyapunov function values ​​are added to the sample data set as new samples and participate in the next iteration of model training, including: In each iteration, based on the trained Gaussian process regression model, the vector d is randomly generated. (i+1) Search for the state point x that is closest to the currently estimated stability region boundary (i+1) , and calculate the state point x (i+1) The corresponding approximate bounded inverse Lyapunov function value After that, it is added as a new sample point into the sample data set to participate in the next iteration of model training.

9. A method for approximating the Lyapunov function and estimating the stability region of a power system according to claim 1, characterized in that: The estimation formula of the stability region is: Where E(x) and Var(x) represent the mean and variance of the Gaussian distribution that point x obeys, respectively; Φ -1 Represents a standard Gaussian distribution The inverse cumulative distribution function; ρ represents the confidence level, which is used to measure the estimated state point in the stable domain at the set t s It can return to the stable equilibrium point x within time eq The probability near ; c represents the level set parameter, which is used to measure the size of the estimated stable domain.

10. An electronic device, characterized in that: The electronic device includes a processor and a memory, wherein the memory stores at least one instruction, at least one program, a code set or an instruction set, and the at least one instruction, the at least one program, the code set or the instruction set is loaded and executed by the processor to implement the method described in any one of claims 1 to 9.