Time-domain numerical recursive calculation method for power system transient stability based on hhl algorithm

Abstract

The present application relates to the technical fields of quantum computing and power system, especially to a power system transient stability time domain numerical recursive calculation method based on HHL algorithm. The present application firstly establishes the mathematical model of power system components; adopts implicit trapezoidal integration method to discretize the state equation, simplifies the difference variable through elimination method, and reduces the network nodes of network topology; then, the difference equation and algebraic equation are solved simultaneously, and the quantum Newton-Raphson algorithm based on HHL algorithm is introduced for efficient solution. The present application introduces quantum algorithm into the transient stability time domain simulation combined with implicit trapezoidal integration method and Newton-Raphson method, utilizes the exponential acceleration characteristics of HHL algorithm in solving linear equations, effectively reduces the complexity of model solution, thereby greatly improves the solving efficiency of residual equation set and overall simulation, and effectively alleviates the "dimension disaster" problem in large-scale power system transient stability analysis.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of quantum computing and power system, and particularly relates to a time-domain numerical recursion calculation method for power system transient stability based on HHL algorithm. BACKGROUND

[0002] The time-domain calculation of power system transient stability usually needs to perform numerical recursion calculation on large-scale, strong nonlinear, differential algebraic equations in a short time. The traditional time-domain simulation of large-scale power system transient stability usually differentiates the state equation by using an implicit trapezoidal integration method with A-stability and strong adaptability to "rigid" equations to obtain a difference equation, and then solves the system difference equation and algebraic equation by using the Newton-Raphson method to avoid the "interface error" caused by the alternately solving of differential equations and algebraic equations and the interface error of components and network.

[0003] With the continuous expansion of the scale of the power system and the rapid growth of the proportion of new energy access, the number of system nodes and state variables increases exponentially. In the process of solving the equation set by using the Newton-Raphson method, the size of the Jacobian matrix increases rapidly, and the solving process of the residual equation set becomes the main calculation bottleneck, which causes the overall simulation efficiency to decrease greatly due to the "dimension disaster" problem.

[0004] Although the existing method introduces parallel computing, node reduction, and sparse matrix optimization to reduce the calculation pressure, it still cannot meet the high efficiency and high precision requirements of large-scale systems at the same time.

[0005] Under this background, quantum computing technology provides a new technical path for the analysis of power system transient stability. SUMMARY

[0006] In order to overcome the "dimension disaster" problem of the time-domain calculation of power system transient stability in the above-mentioned prior art, the present application provides a time-domain numerical recursion calculation method for power system transient stability based on HHL algorithm, which can reduce the complexity of the time-domain calculation of power system transient stability and significantly improve the calculation efficiency.

[0007] The time-domain numerical recursion calculation method for power system transient stability based on HHL algorithm provided by the present application comprises the following steps:

[0008] S1, constructing a mathematical model of each component in the transient stability calculation, the mathematical model comprising a state equation and an algebraic equation; the components comprising: a generator, a load, an excitation system, a prime mover and a governor system thereof;

[0009] S2, calculating the initial value of the transient stability simulation by using the power flow calculation result at the steady state of the power system and the original parameters of each component;

[0010] S3, combine the initial value calculation result with the mathematical model of each element, difference the state equation set and perform elimination and order reduction to obtain a reduced equation set; convert the algebraic equation set under the dq coordinate axis to the xy coordinate axis; and reduce the network nodes;

[0011] S4, form a to-be-solved equation set containing to-be-solved variables by comprehensively processing the equation after S3 and the mathematical model, and construct an iteration equation set of the quantum Newton-Raphson method based on the HHL algorithm as a to-be-solved quantum equation set for the to-be-solved equation set;

[0012] S5, solve the quantum circuit of the to-be-solved quantum equation set constructed by the quantum Newton-Raphson algorithm, output the solving result to a classical computer to correct errors, and obtain a correction result as a to-be-solved variable solving result of the current time step.

[0013] Preferably, step S5 is followed by the following steps:

[0014] S6, compare the correction result with a set stability condition, if unstable, end the simulation and output the result; if the stability condition is met, execute step S7;

[0015] S7, check whether the system state changes; if not, update the non-stationary variables, and then execute step S8; if changes occur, modify the network structure and calculate the algebraic jump variables, then update the non-stationary variables, and then execute step S8; the non-stationary variables are intermediate variables in the difference equation obtained by difference processing of the state equation and updated with recursion;

[0016] S8, update the time t to t+Δt, where Δt is the time step, and update the to-be-solved equation set and the to-be-solved quantum equation set in step S4 synchronously according to the updated non-stationary variables, and then return to step S5.

[0017] Preferably, the initial values of the transient stability simulation in step S2 include: the initial values of the variables to be solved in the quantum Newton-Raphson method, the non-stationary variables, and the variables constant in each stage of stability calculation.

[0018] Preferably, the variables constant in each stage of stability calculation include: the generator capacity conversion ratio K m , the load system capacity conversion ratio K H , and constant impedance loads.

[0019] Preferably, in step S3, the network nodes are reduced by: incorporating the constant impedance load part in the load into the node admittance matrix of the power system, then reducing the network nodes, eliminating the intermediate network nodes, and obtaining a reduced network structure that only retains the generator nodes and dynamic load nodes.

[0020] Preferably, in step S3, the operation mode of eliminating and reducing the order of the state equation is as follows: firstly, extracting the state variables appearing in the algebraic equation as target variables; using the implicit trapezoidal integral method to differentiate the state equation of the generator, the model of the excitation system, the model of the prime mover and its governor, and the state equation of the dynamic load; and performing variable elimination on the obtained differential equation set to obtain an elimination equation set only retaining the target variables and the strong nonlinear slip ratio s.

[0021] Preferably, in step S7, the calculation mode of the algebraic jump variable is as follows:

[0022] S7.1, correcting the node admittance matrix before reduction, then incorporating the corrected node admittance matrix into the original node network of the constant impedance load, and then performing network node reduction to form a new reduced network admittance matrix;

[0023] S7.2, extracting the algebraic quantity at the moment of system mutation and the state quantity at the previous moment, substituting them into the equation set to be solved in step S4, then updating the equation set to be solved, and then executing step S8.

[0024] Preferably, the mathematical model of the generator adopts a three-order practical model, a five-order practical model, or a classical two-order model.

[0025] The mathematical model of the load includes a constant impedance load model and a dynamic load model, and the dynamic load model adopts a dynamic load model considering the mechanical transient process of an induction motor, a dynamic load model considering the electromechanical transient process of an induction motor, or a dynamic load model considering the electromagnetic transient process of an induction motor.

[0026] The application provides a power system transient stability time domain numerical recursive calculation system based on an HHL algorithm.

[0027] The application provides a storage medium storing a computer program, and the computer program is used to implement the power system transient stability time domain numerical recursive calculation method based on the HHL algorithm.

[0028] The application has the following advantages:

[0029] (1) The application provides a power system transient stability time domain numerical recursive calculation method based on an HHL algorithm, which comprises the following steps: firstly, a mathematical model of power system elements such as a generator, a load, an excitation system, a prime mover and a speed regulator system is established; an implicit trapezoidal integration method is used to discretize the state equation, the difference variable is simplified by the elimination method, and the network topology is reduced; then, the difference equation and the algebraic equation are combined, and a quantum Newton-Raphson algorithm based on the HHL algorithm is introduced to solve the equation.

[0030] (2) In this process, the calculation accuracy of the application is controlled by the preset allowable residual error, so that the same accuracy as the classical algorithm can be ensured while the acceleration is realized. Meanwhile, the application can be combined with various generator models, load models, excitation and speed regulation system models, and is suitable for power systems of different sizes and structures, and has good universality and adaptability. BRIEF DESCRIPTION OF DRAWINGS

[0031] Figure 1 A flowchart of the power system transient stability time domain numerical recursive calculation method based on the HHL algorithm is provided.

[0032] Figure 2 A schematic diagram of the IEEE9 node system topology and fault point in the embodiment is provided.

[0033] Figure 3 A relative swing angle between the generators G1 and G2 in the embodiment is provided.

[0034] Figure 4 A relative swing angle between the generators G1 and G3 in the embodiment is provided.

[0035] Figure 5 A relative swing angle between the generators G2 and G3 in the embodiment is provided. DETAILED DESCRIPTION

[0036] The technical solutions in the embodiments of the application will be described clearly and completely in combination with the drawings in the embodiments of the application. Obviously, the described embodiments are only part of the embodiments of the application, rather than all the embodiments of the application. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the application.

[0037] AsFigure 1 As shown, the power system transient stability time domain numerical recursive calculation method based on the HHL algorithm proposed in the embodiment comprises the following steps S1-S8.

[0038] S1, construct a mathematical model of each element in the transient stability calculation, including: generator, load, excitation system, prime mover and its governor system, and the corresponding mathematical model of each element is respectively denoted as synchronous generator model, load model, excitation system model, prime mover and its governor model.

[0039] The mathematical model of each element includes state equation (also known as differential equation) and algebraic equation.

[0040] The synchronous generator model can be selected from three-order practical model, five-order practical model and classical two-order model, etc.

[0041] The load model can be selected from constant impedance load model, dynamic load model considering mechanical transient process of induction motor, dynamic load model considering electromechanical transient process of induction motor, and dynamic load model considering electromagnetic transient process of induction motor, etc.

[0042] The excitation system model and the prime mover and its governor model also have various optional models. Although the model selection is various, the processing method is unified.

[0043] S2, calculate the initial value of transient stability simulation by using the power flow calculation results and the original parameters of each element in the steady state of the power system;

[0044] In this step S2, the power flow calculation results used include the node power and node voltage, and the nodes include generator nodes, load nodes and intermediate network nodes.

[0045] S i =P i +jQ i ;

[0046] U i =U xi +jU yi ;

[0047] Wherein, S i represents the power of node i, P i is the active power of node i, Q i is the reactive power of node i, j is an imaginary number; U i is the voltage of node i, U xi is the x-axis voltage component of node i, and U yi is the y-axis voltage component of node i.

[0048] The simulation starting stage is steady state, and the network node initial value needs to calculate the injection current of each node of the network, which can be obtained by the steady state network admittance matrix and the voltage of each node to form the initial injection current vector :

[0049] = x +j y =YU ;

[0050] wherein, x is the x-axis current vector, y is the y-axis current vector, Y is the admittance matrix, U is the voltage vector;

[0051] The initial value of the transient stability simulation mainly includes three parts:

[0052] 1. The initial value of the variable to be solved in the quantum Newton-Raphson method;

[0053] 2. The intermediate variable updated with recursion in the difference equation obtained by difference processing of the state equation of each element, referred to as the non-constant variable;

[0054] 3. The constant variable in each stage of stability calculation, including: the generator capacity conversion ratio K' m calculated according to the generator capacity and system capacity, the load system capacity conversion ratio K H calculated according to the dynamic load active power proportion coefficient K, and the constant impedance load.

[0055] In the stability calculation stage of the system, the constant impedance load is incorporated into the network admittance matrix for stability calculation.

[0056] According to the common knowledge in the field, the initial values of the above variables can be calculated by setting the state variables of the model at steady state to 0.

[0057] S3, process the mathematical model of each element in step S1, difference process and eliminate the state equation in the mathematical model, and convert the algebraic equation group in the mathematical model from dq coordinate axis to xy coordinate axis; and reduce the network nodes.

[0058] In this step S3, the specific operation mode of eliminating and reducing the state equation is: first, extract the state variables appearing in the algebraic equation as target variables; for the state equation of the generator in S1, the excitation system model, the prime mover and the governor model and the dynamic load T mThe state equation is discretized by using an implicit trapezoidal integration method to obtain a synchronous generator differential equation set, an excitation system differential equation set, a prime mover and a governor differential equation set, and a dynamic load differential equation set; the differential equation set is subjected to variable reduction to obtain a reduced equation set that only retains a target variable and a slip s with strong nonlinearity;

[0059] In the mathematical model, only the stator winding voltage equation of the synchronous generator is a dq coordinate axis equation, and thus, the algebraic equation set in the mathematical model is converted from the dq coordinate axis to the xy coordinate axis, which is actually converting the stator winding voltage equation of the synchronous generator from the dq coordinate axis to the xy synchronous rotating coordinate axis.

[0060] In this step S3, the network nodes are reduced in the following manner: the constant impedance load part (referred to as constant impedance load) in the load is incorporated into the node admittance matrix of the power system, the value of the constant impedance load is calculated from the initial value in step S2, then the network is reduced in nodes, the intermediate network nodes are eliminated, and only the generator nodes and the dynamic load nodes are retained; and the network node voltage equation is decomposed into two real linear algebraic equations in the xy coordinate system.

[0061] The node reduction can reduce the order of the equation set of the node admittance matrix, and the network node reduction needs to be performed once each time the network topology changes.

[0062] S4, the reduced equation set in step S3, the converted stator winding voltage equation of the synchronous generator, the real linear algebraic equation, and the algebraic equation set of the load model in the mathematical model are combined to form a to-be-solved equation set, and a quantum Newton-Raphson method iteration equation set based on the HHL algorithm is constructed for the to-be-solved equation set as a to-be-solved quantum equation set;

[0063] S5, a quantum circuit of the to-be-solved quantum equation set in S4 is constructed by using the quantum Newton-Raphson algorithm, and is solved, and the solution result is output to a classical computer to correct errors and obtain a correction result;

[0064] S6, the correction result is compared with a set stability condition, if instability occurs, the simulation is ended and the result is output; if the stability condition is met, subsequent steps are performed;

[0065] S7, whether a fault or an operation occurs in the system is checked, if not, the original reduced network admittance matrix is still used, and each non-steady variable is updated, and then step S8 is performed;

[0066] If a fault or an operation occurs, the network structure is modified and the algebraic jump variable is calculated, then each non-steady variable is updated, and then step S8 is performed;

[0067] S8, update the time t to t+Δt, Δt is the time step, update the equation group to be solved in step S4 and the quantum equation group to be solved according to the updated non-steady variable, and then return to step S5;

[0068] In this way, steps S5-S7 are repeated, and recursive calculation is performed when the system is stable to obtain the solution results at each time step at the specified time step.

[0069] In step S7, the specific operation of the algebraic jump variable calculation is as follows:

[0070] S7.1, according to the node admittance matrix before the system fault is corrected and reduced, then the corrected node admittance matrix is incorporated into the original node network of the constant impedance load, and then the network node is reduced to form a new reduced network admittance matrix;

[0071] S7.2, extract the algebraic quantity (algebraic quantity refers to the variable in the mathematical model that can have a mutation, and the state quantity that cannot have a mutation is called) at the mutation time of the system mutation and the state quantity at the previous moment, and substitute it into the equation group to be solved in step S4, then update the quantum equation group to be solved, and then execute step S8.

[0072] Hereinafter, the synchronous generator model adopts a five-order practical model, the load model adopts a dynamic load model considering the induction motor electromechanical transient process, and the excitation system model adopts a typical three-order model composed of a regulator, an exciter, and an excitation voltage soft feedback. The prime mover adopts a hydraulic turbine, i.e., the prime mover and its governor system adopt a dynamic hydraulic turbine mechanical governor system.

[0073] The practical five-order model of the generator includes state equations and algebraic equations. The state equations include rotor motion equations and rotor winding transient equations, and the state equation formula is as follows:

[0074] (1);

[0075] The algebraic equation includes a stator winding voltage equation, which is expressed as:

[0076] (2);

[0077] In formulas (1) and (2), ω is the angular velocity of the generator rotor, P m is the prime mover power, P e is the generator output power, D is the constant damping coefficient, Δω is the transition term, Δω=ω-1; ω B is the nominal value, indicating the rated speed of the system; T J is the inertia time constant of the generator unit; δ is the angular displacement of the q-axis leading synchronous rotating coordinate axis x-axis of the generator, which is simply referred to as the rotor angle of the generator; T' d0T is the open-circuit transient time constant for the d-axis of the generator, T' d0 T is the open-circuit sub-transient time constant for the d-axis of the generator, T' q0 E' is the open-circuit sub-transient EMF for the q-axis of the generator, E'' q E is the open-circuit transient EMF for the q-axis of the generator, E'' q E is the open-circuit super-transient EMF for the q-axis of the generator, E'' d E is the open-circuit super-transient EMF for the d-axis of the generator; E f X is the stator excitation EMF of the generator, X d X is the synchronous reactance for the d-axis of the generator, X' d X is the transient reactance for the d-axis of the generator, X'' d X is the super-transient reactance for the d-axis of the generator, X q X is the synchronous reactance for the q-axis of the generator, X'' q X is the super-transient reactance for the q-axis of the generator, U is the stator winding resistance of the generator, U d U and I d are the d-axis voltage and current of the generator, U q U and I q are the q-axis voltage and current of the generator;

[0078] Among the above variables, T J , t, T' d0 , T'' d0 and T'' q0 have the unit of s, i.e., second; ω B is determined according to actual conditions; and the rest of the variables are in per unit.

[0079] The dynamic load model is composed of a constant impedance part and an induction motor considering electromechanical transient process. The constant impedance part is incorporated into the network admittance matrix. The induction motor model can be equivalent to a synchronous motor model in a specific case, and the state equation of the induction motor model includes a rotor motion equation and a rotor winding transient equation, as follows:

[0080] (3) ;

[0081] The algebraic equation of the induction motor model includes two stator voltage equations:

[0082] (4) ;

[0083] The equivalent impedance Z L of the induction motor model is:

[0084] (5) ;

[0085] where s is the slip, and T eT is the electromagnetic torque of the motor; m T is the mechanical torque of the motor, E x E is the transient EMF of the motor in x-axis, E y X is the synchronous reactance of the motor, X is the transient reactance of the motor, T is the open-circuit transient time constant of the motor stator, K H K is the capacity conversion ratio of the motor; U Lx and I Lx are the voltage and current components in the x-axis of the motor, respectively, U Ly and I Ly are the voltage and current components in the y-axis of the motor, respectively; r s and X s are the resistance and leakage reactance of the stator winding, respectively, X m is the mutual inductance between the stator and rotor, r r and X r are the equivalent resistance and leakage reactance of the rotor winding, respectively.

[0086] The model of the excitation system is:

[0087] (6) ;

[0088] wherein, U R is the output voltage of the voltage regulator of the excitation system, T A and K A are the time constant and amplification factor of the inertia amplification link of the excitation system, respectively, U ref is the reference voltage of the excitation system, U t is the amplitude of the terminal voltage of the generator, U F is the negative feedback voltage of the excitation system;

[0089] T L and K L are the time constant and amplification factor of the exciter, respectively, E f is the excitation EMF of the generator stator, S E is the saturation coefficient of the exciter, and S E is related to E f and can be processed by piecewise linearization S E E f = K1E f -K2, K1 and K2 are experimental constants related to the saturation characteristics of the exciter; T F and K F are the time constant and amplification factor of the negative feedback link of the excitation system, respectively.

[0090] The model of the hydraulic turbine and its governor (i.e., the prime mover and its governor system) is:

[0091] (7) ;

[0092] where µ is the servomotor displacement, K δ is the magnification of the centrifugal pendulum speed measurement component, T s is the servomotor time constant, ω r is the reference speed, ω is the generator rotor angular speed, δ i is the static droop coefficient δ i = K i / K δ , K i is the hard feedback magnification, P1 is the transition term; ε represents the dead zone range, ± is determined by the direction of the dead zone, and when the dead zone is not considered, it is set to 0;

[0093] β is the soft feedback coefficient, β = K β / K δ , T i and K β are the soft feedback time constant and magnification, respectively; T w is the water flow time constant, K' m is the prime mover capacity conversion coefficient; P m is the prime mover power.

[0094] S2, the initial value calculation of the transient stability simulation is carried out by using the steady-state power flow calculation results and the original parameters of each element, and the power flow calculation results used include the power of each node and the node voltage .

[0095] The simulation starts in a steady state, and the initial value of each node of the network needs to be calculated. The initial injection current of each node can be calculated by the steady-state network admittance matrix and the voltage of each node:

[0096] = x +j y =YU (8) ;

[0097] wherein, x is the x-axis current vector, y is the y-axis current vector, Y is the admittance matrix, U is the voltage vector;

[0098] The initial value calculation of each element mainly includes three parts:

[0099] 1. The initial value of the variable to be solved in the quantum Newton-Raphson method, including the q-axis super transient electromotive force E'' q0 of the generator part, the d-axis super transient electromotive force E'' d0Rotor angle δ0, voltage and current U along the x and y axes x0 U y0 I x0 and I y0 The initial values ​​of the seven variables; and the slip s0, transient electromotive force, voltage, and current E' of the dynamic load component. x0 、E' y0 U Lx0 U Ly0 I Lx0 and I Ly0 Seven initial values ​​for variables, etc.; a total of 14 initial values ​​for variables in this part can be obtained by simulating after setting the state variables to 0 when the system is in steady state.

[0100] 2. Unstationary variables; see the variable set in the subsequent difference equations for details. .

[0101] 3. Constant variables in each stage of the stability calculation, including the generator capacity reduction ratio calculated based on generator capacity and system capacity. According to the dynamic load active power ratio coefficient Calculate the load system capacity conversion ratio and constant impedance loads.

[0102] The initial values ​​of the above variables can be calculated by setting the state variables of the model to 0 in steady state.

[0103] S3. Process the differential (state) algebraic equations, mainly including differentiating the state equation system and performing elimination to reduce its order, transforming the algebraic equation system from the dq coordinate axis to the xy coordinate axis, and reducing the number of network nodes. The specific operations are as follows:

[0104] S3.1. The generator state equation (1), excitation system model equation (6), turbine and governor model equation (7), and dynamic load state equation (3) in S1 are differentially processed using the implicit trapezoidal integral method. Then, variable elimination is performed. Based on the input and output variables between each component, further elimination is performed, retaining only the state variables appearing in the algebraic equations and the slip rate with strong nonlinearity. After processing, the following system of difference equations is obtained:

[0105] Synchronous generator differential equations:

[0106] (9);

[0107] Dynamic load difference equations:

[0108] (10);

[0109] In equations (9) and (10), P e,n+1 The generator output power P at time step n+1 e U x,n+1 and I x,n+1 These represent the generator x-axis voltage and current at time step n+1, respectively. y,n+1 and I y,n+1 These represent the generator's y-axis voltage and current at time step n+1, respectively, δ n+1 Let δ be the rotor angle of the generator at time step n+1, and T be the rotor angle of the generator. m,n+1 The mechanical torque T of the motor at time step n+1 m T e,n+1 The electromagnetic torque T of the motor at time step n+1 e s n+1 Let s be the slip rate at time step n+1, and E' be the slip rate at time step n+1. y,n+1 The transient electromotive force E' of the motor along the y-axis at time step n+1 y E' x,n+1 Let E' be the transient electromotive force E' of the motor along the x-axis at time step n+1. x I Ly,n+1 Let I be the current component on the y-axis of the motor at time step n+1. Lx,n+1 The current component on the x-axis of the motor at time step n+1;

[0110] There are two sets of coefficients in equations (9) and (10). and b ? These are all parameter variables obtained after simplification following differentiation, among which the variables... When the time step is fixed, b is a constant, and the variable is b. ? These are non-stationary variables that are updated as the process continues.

[0111] variable Including intermediate variables ' δ , δ , ω , ' pm , 1. 2. 3. p1 , pm , µ , q1 , ' q , '' q , EF E EF UR E UF q2 d d s Lx

[0112]

[0113] where ω B is the nominal value; h is the time step, T J is the inertia time constant of the generator unit, D is the constant damping coefficient of the generator, β is the soft feedback coefficient, T i is the soft feedback time constant, T w is the water flow time constant, K' m is the capacity conversion coefficient of the prime mover, K δ is the magnification of the centrifugal flyweight speed measuring part, T s is the servomotor time constant, δ i is the static regulation difference coefficient, T' d0 is the d-axis open circuit transient time constant, T'' d0 is the d-axis open circuit sub-transient time constant of the generator, T A and K A ​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​T and K are respectively the time constant and the amplification factor of the inertial amplification link L and K L are respectively the time constant and the amplification factor of the exciter, K1 is an experimental constant related to the saturation characteristic of the exciter, T F and K F are respectively the time constant and the amplification factor of the exciter negative feedback link, X d is the synchronous reactance of the generator d-axis, X' d is the transient reactance of the generator d-axis, X'' d is the super-transient reactance of the generator d-axis, X q is the synchronous reactance of the generator q-axis, X'' q is the super-transient reactance of the generator q-axis, T'' q0 is the open-circuit super-transient time constant of the generator q-axis, T' is the open-circuit transient time constant of the motor stator;

[0114] non-steady variable b ? includes the intermediate variables b' δ , b' pm , b p1 , b µ , b pm , b ω , b δ , b' Ef , b UR , b UF , b Ef , b q , b' q , b'' q , b d , b s , b Lx , b Ly ;

[0115] , , , , , , ;

[0116] , , , ;

[0117] , ,

[0118] ;

[0119] ; ;

[0120] ;

[0121] ;

[0122] where ω B is the nominal value; Δω r is the intermediate term, Δω r = ω r - 1, ω r is the reference speed; ε represents the range of the prime mover dead zone, U ref is the excitation system reference voltage, K2 is an experimental constant related to the saturation characteristic of the exciter, X is the synchronous reactance of the motor, X' is the transient reactance of the motor, K H is the motor capacity conversion ratio, P 1,n is the transition term P1 at time step n, µ n is the prime mover actuator displacement µ at time step n, Δω n is the transition term Δω at time step n, P m,n is the prime mover power P m at time step n, P e,n is the generator output power P e at time step n, δ n is the rotor angle δ of the generator at time step n, U t,n is the generator terminal voltage amplitude U t at time step n, U F,n is the excitation negative feedback voltage U F at time step n, U R,n is the voltage regulator output voltage U R at time step n, E f,n is the generator stator excitation electromotive force E f at time step n, E' q,n is the generator q-axis transient electromotive force E' q at time step n, I d,n is the generator d-axis current I d at time step n, E'' q,n is the generator q-axis super transient electromotive force E'' q at time step n, I q,n is the generator q-axis current I q at time step n, E'' d,n is the generator d-axis super transient electromotive force E'' d at time step n, T m,n is the motor mechanical torque T m at time step n, T e,n is the motor electromagnetic torque T e at time step n, s nLet s be the slip rate at time step n, and E' be the slip rate at time step n. y,n The transient electromotive force E' of the motor on the y-axis at time step n y E' x,n The transient electromotive force E' of the motor on the x-axis at time step n x I Ly,n I represents the current component on the y-axis of the motor at time step n. Lx,n The current component on the x-axis of the motor at time step n; the obtained parameters and variable values ​​at time step n are used to obtain the solution at the next time step n+1; when time step n is time t, time step n+1 is time t+Δt.

[0123] S3.2, Transform the generator algebraic equation, i.e., the stator winding voltage equation (2), into the transformation matrix. Transform from the dq coordinate axis to the xy synchronously rotated coordinate axis for network-based solution:

[0124] (11);

[0125] Among them, E'' d E'' is the d-axis ultratransient electromotive force of the generator. q X'' is the q-axis ultratransient electromotive force of the generator. d X'' is the d-axis transient reactance of the generator. q For the generator's q-axis ultratransient reactance, U x U is the generator x-axis voltage. y I is the generator's y-axis voltage. x I is the x-axis current of the generator. y This represents the generator's y-axis current. δ is the generator stator winding resistance, and δ is the generator rotor angle.

[0126] S3.3. Incorporate the constant impedance load portion of the load into the node admittance matrix. The value of the constant impedance load is calculated from the initial value. Then, perform node reduction on the network, eliminating intermediate network nodes and retaining only generator nodes and dynamic load nodes to reduce the order of the equation system. Network node reduction needs to be performed once every time the network topology changes. Then, decompose the network node voltage equations into two real linear algebraic equations in the xy coordinate system:

[0127] (12);

[0128] Where i and j represent node numbers, I xi and I yi U represents the x-axis and y-axis currents injected at node i, respectively. xi and U yi The x-axis and y-axis voltages of the i-th node are respectively. ij=G ij +jB ij , Y ij represents admittance between nodes i and j, G ij represents conductance between nodes i and j, B ij represents susceptance between nodes i and j; when i=j, Y ij , G ij and B ij represent self-admittance, conductance and susceptance of node i respectively; j∈i represents that the parameters of node j are summed to the parameters of node i.

[0129] S4, construct the iteration equation set of quantum Newton-Raphson method based on HHL algorithm. Seven equations corresponding to formula (9), formula (11) and formula (12) form the equation set to be solved for a single generator node, and seven equations corresponding to formula (4), formula (10) and formula (12) form the equation set to be solved for a single dynamic load node, and let the number of generator nodes and the number of dynamic load nodes of the reduced node network be N G and N L respectively, then the total number of linear equation sets to be solved simultaneously is N=7N G +7N L , and the transient stability residual equation set and the Jacobian matrix are solved by using Newton-Raphson method for the N equation sets, which are converted into the quantum equation set to be solved shown in formula (13), also known as the residual equation set:

[0130] (13);

[0131] where, Δf is the residual, J (k) is the transient Jacobian matrix of the kth iteration; x (k) is the variable value of the kth iteration of Newton-Raphson method, x (k) is the calculation result of the previous time step; Δf(x (k) ) is the transient residual of the kth iteration, Δx (k) is the transient correction amount Δx of the kth iteration, and the variable x includes 14 kinds of variables of E'' q , E'' d , δ, U x , U y , I x , I y of synchronous generators and s, E' x , E' y , U Lx , U Ly , I Lx , I Ly of induction motors.

[0132] S5. Constructing the quantum circuit of the transient stability residual equation of quantum Newton-Raphson method to solve efficiently, outputting the result to the classical computer to correct the error, iteratively solving the result of the next time step, and repeating the process for recursive calculation, as follows:

[0133] S5.1, check whether the transient stability residual equation formula (13) meets the calculation requirements of HHL algorithm, including:

[0134] 1. The transient stability Jacobian matrix J meets the Hermitian condition;

[0135] 2. The equation order is an integer power of 2;

[0136] 3. The transient stability residual Δf(x (k) ) is a unit vector.

[0137] When the conditions are not met, the following pretreatment is needed for the residual equation set:

[0138] 1. If the transient stability Jacobian matrix J does not meet the Hermitian condition, the following formula is used to perform Hermitian processing on formula (13);

[0139] (14);

[0140] Where J H is the conjugate transpose matrix of J; 01 is an N-order 0 vector added by Hermitian, that is, a vector composed of N 0 elements; J (k) is the transient Jacobian matrix of the kth iteration, J (k)H is the conjugate transpose matrix of J (k) .

[0141] 2. If the equation order of formula (13) does not meet the integer power of 2, the following formula (15) is used to perform order upgrading on formula (13):

[0142] (15);

[0143] Where, is an identity matrix whose order is , represents the integer part of the number; represents the eigenvalue of the kth iteration, which is an integer power of 2 within the range of the absolute value of the eigenvalue ; 01 is an N-order 0 vector added by Hermitian; 02 is a 0 vector with the same order as I; where |λ (k) | min represents the minimum value of the eigenvalue |λ(k) | max denotes maximum value of eigenvalue.

[0144] 3. If the transient stability residual error Δf(x (k) ) is not a unit vector;

[0145] Specifically, formula (13) can be normalized by formula (16):

[0146] (16);

[0147] Formula (15) can be normalized by formula (16.1):

[0148] (16.1);

[0149] wherein, denotes the modulus value.

[0150] In this way, through the above processing, it can be ensured that formula (13) satisfies conditions 1-3 in step S5.1 after processing. If formula (13) satisfies conditions 1-3, a quantum circuit of formula (13) is constructed, and formula (13) is solved; otherwise, formula (13) is processed to ensure that the processed formula satisfies conditions 1-3, and a quantum circuit is constructed for the processed formula, so as to solve formula (13).

[0151] It is assumed that formula (13) satisfies conditions 1-3 in step S5.1, and the solving process of formula (13) is described.

[0152] S5.2, constructing a HHL algorithm quantum circuit for fast solving of the residual error equation set shown in formula (13), specifically including the following sub-steps S5.2.1-S5.2.5.

[0153] S5.2.1, setting three groups of quantum bits q a , q b , q c needed by the HHL algorithm:

[0154] q b is a quantum register for storing phase information of eigenvalues, and the number of quantum bits n b depends on the range of absolute values of eigenvalues of the transient stability Jacobian matrix and the storage accuracy; for every n pre bit of decimal accuracy, n quantum bits need to be added, and represents rounding up; in specific implementation, n b =n0+ , This represents the floor function, where n0 is the number of qubits corresponding to the base precision.

[0155] q a As an auxiliary qubit with controlled rotation, its function is to convert q b The eigenvalue information stored in the q is transferred to the amplitude. a The number of qubits n a Fixed at 1;

[0156] q c For storing transient residuals Δf(x) (k) ), and is used to output the transient correction quantity Δx, whose number of qubits n c =log2N.

[0157] S5.2.2, Regarding the transient residual -Δf(x) (k) Amplitude encoding is performed to obtain |-Δf(x) (k) The amplitude is encoded as |-Δf(x)>. (k) Input quantum register q c In the middle, simultaneously on the quantum register q b Applying a Hadamard gate to each qubit constructs a quantum superposition state. Then, the transient stability Jacobian matrix is ​​used. Preparation of unitary matrix j is an imaginary number; t0 is a set value used to ensure the eigenvalue λ i (k) phase The value of λ is between 0 and 1; i (k) J is the Jacobian matrix (k) eigenvalues; with q b To control the bit in q c A series of power gates of U are applied, followed by a quantum inverse Fourier transform, to transform the eigenvalues ​​λ. i (k) The phase information is stored in the quantum register q. b middle.

[0158] S5.2.3, using quantum register q b To control the bit, in the auxiliary bit q a Apply a series of genus gates R y (θ i ), and the eigenvalue λ i (k) The phase information is transferred to the amplitude.

[0159] Specifically: ;

[0160] Where, θi For transition, θ i =2arcsin , C is a normalization constant to ensure the sum of amplitude square is 1.

[0161] S5.2.4, apply inverse quantum phase estimation (IQPE) to remove q b 's entanglement with q c , to get the quantum state shown in equation (17).

[0162] This step is specifically implemented by applying the adjoint operators of the gates used in quantum phase estimation part in reverse order, at this time the quantum state is obtained:

[0163] (17);

[0164] Where, β i (k) is the coefficient of amplitude encoding |-Δf(x (k) )> on the corresponding eigenvector, |u i (k) is a set of eigenvectors of the temporary Jacobian matrix J (k) of the kth iteration.

[0165] S5.2.5, apply quantum measurement to q a and q c at the same time, when the measurement result of q a is , the output of q c is:

[0166] (18);

[0167] Where, is the approximate solution of the temporary correction Δx (k) of the kth iteration, is the quantum state of , and after inverse encoding and inverse normalization processing of , we can get .

[0168] S5.3 set the convergence condition |Δx (k) |≤ε1, judge whether it is less than or equal to the convergence threshold ε1;

[0169] If not, modify the variable x (k+1) =x (k) + , update k to k+1, and repeat S5.2;

[0170] If it is satisfied, output the variable x.

[0171] Specifically, the stability index can be extracted from the variable x according to the set stability condition to determine the system stability. The variable x contains 14 variables of the synchronous generator, E q , E d , δ, U x , U y , I x , I y and the s, E' x , E' y , U Lx , U Ly , I Lx , I Ly of the induction motor.

[0172] S6, comparing the stability index in the solution result with the set stability condition;

[0173] If the system is unstable, i.e., the solution result does not satisfy the stability condition, the simulation is ended and the solution result is outputted;

[0174] If the solution result satisfies the stability condition, step S7 is performed.

[0175] S7, checking whether the system has a fault or operation, if not changed, the original reduced network admittance matrix is still used, and the non-steady variable b ? is updated; if the fault or operation occurs, the network structure is modified and the algebraic jump variable calculation is performed, and the specific operation is as follows:

[0176] S7.1, when the system change is three-phase symmetry, such as three-phase symmetrical metallic fault or three-phase line cut-off, only the admittance matrix of the positive sequence network of the system needs to be modified; if the fault is an asymmetrical fault, the negative sequence and zero sequence components are expressed by the positive sequence quantities according to the boundary conditions of the fault point, the positive sequence fault node admittance relationship is constructed, the effects of the negative sequence and zero sequence networks are equivalent by forming corresponding comprehensive impedance matrices, and the positive sequence admittance matrix is modified by the comprehensive impedance matrix. The modified node admittance matrix is incorporated into the original node network of the constant impedance load, and then the network node reduction is performed to form a new reduced network admittance matrix.

[0177] S7.2, at the time t n when the system topology changes, the state quantity of the system remains unchanged, i.e., the state quantity at the time t n+ is directly taken as the value at the time t n- , the time t n+ is adjacent to the time t n and located after the time t n , and the time t n- is adjacent to the time t n and located after the time t nprevious time step; while the population size is mutated, at which time t n- the state variables at time t are substituted into equations (4), (11), and (12), and then simultaneous solutions are obtained, so that the population sizes at time t are obtained n+ These population sizes are used as the predicted values at the next time together with the state variables calculated at the present time, and the intermediate variables (i.e., the unsteady variables required for calculation) at the present time are updated, and the unsteady variables at the next time are calculated.

[0178] S8, the current time t is updated to t+Δt, and the steps S5-S7 are repeated until the simulation time ends and the results are output, and Δt is the time step length.

[0179] In this embodiment, the classic method and the method of the present application are respectively used for transient stability time domain simulation on the IEEE9-node standard system shown in FIG. 1. Figure 2

[0180] The nodes of the IEEE9-node system are renumbered, nodes 1-3 are generator nodes, nodes 4-6 are load nodes, and nodes 7-9 are intermediate network nodes. The generator uses a five-order practical model, the load uses a dynamic load model considering the electromechanical transient process of an induction motor, the excitation system uses a typical three-order model, the prime mover and speed regulating system use a hydro-turbine speed regulating system, and the network uses a quasi-steady state model (the admittance value is referenced to the value of the IEEE9-node standard system). The parameters of each element are shown in Tables 1-4.

[0181] Table 1: Generator-related parameters

[0182] ;

[0183] Table 2: Excitation machine-related parameters

[0184] ;

[0185] Table 3: Hydro-turbine speed regulating system parameters

[0186] ;

[0187] Table 4: Dynamic load parameters

[0188] ;

[0189] K is a dynamic load active power proportionality coefficient.

[0190] The transient stability time domain simulation parameters are set as follows:

[0191] ​​The simulation duration is 3s, the simulation step h=0.02s, the first 0~0.2s is set as the pre-fault steady state, a three-phase short circuit occurs at the line 5-7 close to node 7 at 0.2s, and the fault is cleared after one time step; the convergence accuracy of the quantum Newton-Raphson method is set to , the solving variable is the change of the relative swing angle of the generator, and the stability condition is that the relative swing angle is less than 180 degrees.

[0192] The relative swing angle between the generators G1 and G2 is shown in Figure 3 , the relative swing angle between the generators G1 and G3 is shown in Figure 4 , and the relative swing angle between the generators G2 and G3 is shown in Figure 5 .

[0193] In this embodiment, the stability condition is that the peak value of the relative swing angle is not more than 180 degrees; from the simulation results, the peak value of each relative swing angle shows a decreasing trend, which is shown in Figures 3-5 , the second swing is less than the first swing, and it can be judged that the system is transiently stable. The first swing is the first peak value of the relative swing angle, and the second swing is the second peak value of the relative swing angle.

[0194] At the same time, the method of the present application highly coincides with the curve of the classical algorithm, proving the feasibility of the present application.

[0195] Of course, for those skilled in the art, the present application is not limited to the details of the above exemplary embodiments, but also includes the same or similar structures that can be realized in other specific forms without departing from the spirit or essential characteristics of the present application. Therefore, from any point of view, the embodiments should be regarded as exemplary and non-limiting, and the scope of the present application is defined by the appended claims rather than the above description, and therefore all changes falling within the meaning and scope of the equivalent elements of the claims are intended to be included in the present application. Any reference signs in the claims should not be regarded as limiting the claims involved.

[0196] In addition, it should be understood that although the present specification is described in terms of embodiments, not every embodiment contains only one independent technical solution, and the description manner of the specification is only for the sake of clarity, and those skilled in the art should regard the specification as a whole, and the technical solutions in each embodiment can also be properly combined to form other embodiments that those skilled in the art can understand.

[0197] The technologies, shapes, and structural parts not described in detail in the present application are well-known technologies.

Claims

1. A time-domain numerical recursive calculation method for the transient stability of a power system based on the HHL algorithm, characterized in that, Includes the following steps: S1. Construct mathematical models for each component in transient stability calculation. The mathematical models include state equations and algebraic equations. The components include: generator, load, excitation system, prime mover and its governor system. S2. Calculate the initial values ​​for transient stability simulation using the power flow calculation results and the original parameters of each component in steady state of the power system; S3. Based on the initial value calculation results, process the mathematical model of each component, differentiate the state equation system and perform elimination to reduce the order to obtain the elimination equation system; transform the algebraic equation system under the dq coordinate axis to the xy coordinate axis; and reduce the network nodes. S4. The equations processed by S3 and the mathematical model are combined to form a set of equations to be solved, which includes variables to be solved. For the set of equations to be solved, an iterative set of quantum Newton-Raphson method based on the HHL algorithm is constructed as a set of quantum equations to be solved. S5. The quantum circuit of the quantum equation system to be solved is constructed using the quantum Newton-Raphson algorithm. The solution is output to the classical computer to correct the error. The corrected result is used as the solution result of the variable to be solved at the current time step. The set of quantum equations to be solved in step S4 is as follows: Where Δf is the residual, J (k) Let x be the transient Jacobian matrix of the k-th iteration; (k) Let be the variable value in the k-th iteration of the Newton-Raphson method. x at time (k) This is the calculation result from the previous time step; Δf(x) (k) Let Δx be the transient residual of the k-th iteration. (k) Let Δx be the transient correction value for the k-th iteration; In step S5, the quantum equations to be solved are first made to meet the computational requirements of the HHL algorithm, and then the variable x is iterated until the convergence condition |Δx| is met. (k) |≤ε1, where ε1 is the convergence threshold; the iteration method for variable x is as follows: construct the HHL algorithm quantum circuit to quickly solve the quantum equation system to be solved, and obtain Δx. (k) Approximate solution Correcting variable x (k+1) =x (k) + .

2. The time-domain numerical recursive calculation method for transient stability of power systems based on the HHL algorithm as described in claim 1, characterized in that, Step S5 is followed by the following steps: S6. Compare the correction result with the set stability condition. If instability occurs, the simulation ends and the result is output. If the stability condition is met, proceed to step S7. S7. Check if the system state has changed; if not, update each unsteady variable and then execute step S8; if it has changed, modify the network structure and perform algebraic jump variable calculation, then update each unsteady variable and then execute step S8; unsteady variables are intermediate variables that are updated as the state equation is recursively processed in the difference equation. S8. Update time t to t+Δt, where Δt is the time step. Synchronously update the system of equations to be solved and the system of quantum equations to be solved in step S4 according to the updated unsteady variables, and then return to step S5.

3. The time-domain numerical recursive calculation method for transient stability of power systems based on the HHL algorithm as described in claim 1, characterized in that, The initial values ​​for transient stability simulation in step S2 include: the initial values ​​of the variables to be solved in the quantum Newton-Raphson method, unsteady variables, and constant variables in each stage of stability calculation; the unsteady variables are the intermediate variables that are updated as the recursion proceeds in the difference equation after the state equation has been processed by the difference equation.

4. The time-domain numerical recursive calculation method for transient stability of power systems based on the HHL algorithm as described in claim 3, characterized in that, The constant variables in each stage of the stability calculation include: generator capacity conversion ratio K' m Load system capacity conversion ratio K H and constant impedance loads.

5. The time-domain numerical recursive calculation method for transient stability of power systems based on the HHL algorithm as described in claim 1, characterized in that, In step S3, the method for reducing network nodes is as follows: the constant impedance load portion of the load is incorporated into the node admittance matrix of the power system, and then the network nodes are reduced to eliminate intermediate network nodes, resulting in a reduced network structure that retains only generator nodes and dynamic load nodes.

6. The time-domain numerical recursive calculation method for transient stability of power systems based on the HHL algorithm as described in claim 1, characterized in that, In step S3, the operation method for eliminating variables and reducing the order of the state equation is as follows: First, extract the state variables that appear in the algebraic equation as the target variables; then, perform differential processing on the generator state equation, excitation system model, prime mover and its governor model, and dynamic load state equation using the implicit trapezoidal integral method. The resulting difference equations are reduced by variables to obtain a reduced equation set that retains only the objective variable and the highly nonlinear slip rate s.

7. The time-domain numerical recursive calculation method for transient stability of power systems based on the HHL algorithm as described in claim 1, characterized in that, The algebraic jump variable is calculated in step S7 as follows: S7.

1. Based on the system fault correction, the node admittance matrix before reduction is corrected. Then, the corrected node admittance matrix is ​​incorporated into the original node network of the constant impedance load. Then, the network nodes are reduced to form a new reduced network admittance matrix. S7.2 Extract the number of generations at which the system mutation occurs and the state variables of the previous moment, substitute them into the system of equations to be solved in step S4, then update the system of quantum equations to be solved, and then execute step S8.

8. The method for time-domain numerical recursive calculation of power system transient stability based on the HHL algorithm as described in any one of claims 1-7, characterized in that, The mathematical model of the generator adopts a third-order practical model, a fifth-order practical model, or a classical second-order model; The mathematical models of the load include constant impedance load models and dynamic load models. The dynamic load models adopt dynamic load models that consider the mechanical transient process of the induction motor, dynamic load models that consider the electromechanical transient process of the induction motor, or dynamic load models that consider the electromagnetic transient process of the induction motor.

9. A time-domain numerical recursive calculation system for transient stability of a power system based on the HHL algorithm, characterized in that, It includes a memory and a processor. The memory stores a computer program, and the processor is connected to the memory. The processor is used to execute the computer program to implement the time-domain numerical recursive calculation method for transient stability of power systems based on the HHL algorithm as described in any one of claims 1-8.

10. A storage medium, characterized in that, The system contains a computer program that, when executed, implements the time-domain numerical recursive calculation method for transient stability of power systems based on the HHL algorithm as described in any one of claims 1-8.

Images