CN106294993A

CN106294993A - A kind of transient energy function analysis method considering that inverter current is saturated

Info

Publication number: CN106294993A
Application number: CN201610647261.1A
Authority: CN
Inventors: 孙震宇; 马迪; 刘馨月; 张梦月; 王俊; 杨阳; 卜京; 姚娟; 谢云云; 殷明慧
Original assignee: Nanjing University of Science and Technology
Current assignee: Nanjing University of Science and Technology
Priority date: 2016-08-09
Filing date: 2016-08-09
Publication date: 2017-01-04
Anticipated expiration: 2036-08-09
Also published as: CN106294993B

Abstract

The invention provides a transient energy function analysis method considering inverter current saturation, based on the inverter virtual synchronous generator control strategy, by establishing the virtual synchronous generator rotor motion equation and considering the inverter current saturation factor , use the first integration method to construct the transient energy function, and use the linear path method to approximate the non-integrable items related to the path in the transient energy function, and then use the BCU method to obtain the critical energy value and the critical cut-off time. The transient energy function analysis method considering the virtual moment of inertia and inverter current saturation in the present invention can quantitatively evaluate the stability of the power system containing the virtual synchronous generator, and the constructed transient energy function takes into account the inverter Saturation factor, which is closer to reality, solves the problem of system instability caused by traditional online evaluation conservativeness, and is an important supplement to the time domain simulation method.

Description

A Transient Energy Function Analysis Method Considering Inverter Current Saturation

技术领域technical field

本发明设计一种电力系统稳定与控制技术，特别是一种考虑逆变器电流饱和的暂态能量函数分析方法。The invention designs a power system stability and control technology, in particular a transient energy function analysis method considering inverter current saturation.

背景技术Background technique

目前，随着分布式电源渗透率的不断增加以及将分布式电源与电网相连接电力电子器件采用数字电路控制的方式，暂态响应速度快，且几乎没有惯性，当DG占到一定比例时，电网小的扰动都可以造成安全稳定问题。在此背景下，国内外学者提出了虚拟同步发电机技术，通过控制逆变器模拟同步发电机的运行原理实现了分布式电源的友好接入。At present, with the increasing penetration rate of distributed power and the use of digital circuit control for power electronic devices connected to distributed power to the grid, the transient response speed is fast and there is almost no inertia. When DG accounts for a certain proportion, Small disturbances in the power grid can cause security and stability problems. In this context, scholars at home and abroad have proposed virtual synchronous generator technology, which realizes the friendly access of distributed power sources by controlling the inverter to simulate the operating principle of synchronous generators.

对虚拟同步发电机的研究国内外学者关注点在并网逆变器模拟同步发电机惯性和阻尼，虽然有些学者考虑虚拟转动惯量逆变器暂态能量函数分析，但是均忽略了带有虚拟惯量的下垂控制逆变器的电流饱和特性对暂态稳定性的影响，这样得到的临界切除时间太过于保守，当系统在线评估时会因为切除时间过长从而导致系统不稳定这一重大问题，直接影响到国民经济，因此，针对保守性导致的不稳定问题需要考虑逆变器电流饱和因素对暂态稳定的影响。Research on virtual synchronous generators domestic and foreign scholars focus on the grid-connected inverter to simulate the inertia and damping of synchronous generators. Although some scholars consider the transient energy function analysis of virtual moment of inertia inverters, they all ignore the virtual inertia The impact of the current saturation characteristics of the droop control inverter on the transient stability, the critical cut-off time obtained in this way is too conservative, and when the system is evaluated online, it will lead to the major problem of system instability due to the long cut-off time, directly It affects the national economy. Therefore, in view of the instability caused by conservatism, it is necessary to consider the influence of the inverter current saturation factor on the transient stability.

发明内容Contents of the invention

本发明的目的在于提供一种考虑逆变器电流饱和的暂态能量函数分析方法，解决了传统的考虑虚拟转动惯量逆变器暂态能量函数分析方法得到的临界切除时间太过于保守，当系统在线评估时会因为切除时间过长从而导致系统不稳定这一重大问题。The purpose of the present invention is to provide a transient energy function analysis method considering inverter current saturation, which solves the problem that the critical cut-off time obtained by the traditional transient energy function analysis method considering the virtual moment of inertia inverter is too conservative. There is a major problem of system instability due to the long resection time during online evaluation.

本方法包括以下步骤：This method comprises the following steps:

步骤1，对逆变器建立含有虚拟同步发电机的电力系统的详细数学模型，该模型包括类似于传统发电机的两阶模型转子运动方程以及非线性负载方程和网络方程；Step 1. Establish a detailed mathematical model of the power system with a virtual synchronous generator for the inverter, which includes a two-order model rotor motion equation similar to a traditional generator, as well as nonlinear load equations and network equations;

步骤2，利用首次积分法构建暂态能量函数，将转子运动方程写成一阶方程组，利用首次积分法对两边同时进行积分，并将逆变器电流饱和因素考虑在内构造计及阻尼的能量型李雅普诺夫函数，对暂态能量函数中和路径有关的不可积项采用线性路径方法近似；Step 2. Use the first integral method to construct the transient energy function, write the rotor motion equation as a first-order equation system, use the first integral method to integrate both sides simultaneously, and take the inverter current saturation factor into account to construct the energy that takes damping into account type Lyapunov function, the non-integrable items related to the path in the transient energy function are approximated by the linear path method;

步骤3，采用BCU方法获得临界能量值：首先计算故障时轨线至故障清除时刻及至出口点，利用故障清除时刻的参数计算该时刻的能量，接着计算故障后系统的最小梯度点并根据求解功率平衡方程获取主导不稳定不平衡点，并将该点状态参数获得临界能量值；Step 3, using the BCU method to obtain the critical energy value: first calculate the trajectory from the fault time to the fault clearing time and to the exit point, use the parameters at the fault clearing time to calculate the energy at this time, then calculate the minimum gradient point of the system after the fault and solve it according to the power The balance equation obtains the dominant unstable and unbalanced point, and obtains the critical energy value of the state parameter at this point;

步骤4，当暂态能量等于主导不稳定平衡点处的能量时求解获得临界切除时间。Step 4, when the transient energy is equal to the energy at the dominant unstable equilibrium point, the critical cut-off time is obtained by solving.

本发明求得临界故障切除时间并与传统的忽略逆变器饱和得到的临界能量进行对比更贴近于实际。The invention obtains the critical fault removal time and compares it with the critical energy obtained by ignoring the saturation of the inverter, which is closer to reality.

下面结合说明书附图对本发明做进一步描述。The present invention will be further described below in conjunction with the accompanying drawings.

附图说明Description of drawings

图1是考虑虚拟转动惯量逆变器暂态能量函数分析方法流程图。Fig. 1 is a flow chart of the analysis method for the transient energy function of the inverter considering the virtual moment of inertia.

图2是三台虚拟同步发电机九节点系统仿真图。Figure 2 is a simulation diagram of a nine-node system with three virtual synchronous generators.

图3是BCU方法获取临界能量和CCT流程图。Fig. 3 is a flow chart of BCU method to obtain critical energy and CCT.

图4是系统总能量及动能和势能仿真结果图。Figure 4 is the simulation results of total energy, kinetic energy and potential energy of the system.

具体实施方式detailed description

本发明公开了一种考虑虚拟转动惯量和逆变器电流饱和的暂态能量函数分析方法，基于逆变器虚拟同步发电机控制策略基础上，通过建立虚拟同步发电机转子运动方程并考虑逆变器的电流饱和因素，利用首次积分法构建暂态能量函数，针对暂态能量函数中和路径有关的不可积项采用线性路径方法近似，接着采用BCU方法获得临界能量值，最后通过获取临界切除时间和时域仿真获取的临界时间进行了对比。本发明的考虑虚拟转动惯量和逆变器电流饱和的暂态能量函数分析方法能够定量的对含有虚拟同步发电机的电力系统的稳定性进行评估，所构造的暂态能量函数考虑了逆变器的饱和因素，更贴近实际，解决了传统的在线评估保守性导致系统不稳定的问题，是对时域仿真方法的一种重要补充。具体实施方式如下(如图1所示)：The invention discloses a transient energy function analysis method considering the virtual moment of inertia and inverter current saturation. Based on the inverter virtual synchronous generator control strategy, by establishing the virtual synchronous generator rotor motion equation and considering The current saturation factor of the transformer, using the first integral method to construct the transient energy function, using the linear path method to approximate the non-integrable items related to the path in the transient energy function, then using the BCU method to obtain the critical energy value, and finally by obtaining the critical cut-off time Compared with the critical time obtained by time domain simulation. The transient energy function analysis method of the present invention considering the virtual moment of inertia and inverter current saturation can quantitatively evaluate the stability of the power system containing the virtual synchronous generator, and the constructed transient energy function takes the inverter The saturation factor is closer to reality, and it solves the problem of system instability caused by traditional online evaluation conservativeness, and is an important supplement to the time domain simulation method. The specific implementation is as follows (as shown in Figure 1):

具体地specifically

第一步，建立含有电力系统的详细数学模型，其中包括：In the first step, a detailed mathematical model containing the power system is established, which includes:

①建立虚拟同步发电机的转子运动方程。逆变器采用虚拟同步发电机控制策略，建立类似于传统发电机的两阶模型转子运动方程如下；①Establish the rotor motion equation of the virtual synchronous generator. The inverter adopts the virtual synchronous generator control strategy, and establishes a two-order model rotor motion equation similar to the traditional generator as follows;

$\begin{matrix} \{\begin{matrix} {\overset{\cdot \cdot}{θ θ}}_{i i} = = {\overset{~ ~}{ω ω}}_{i i} \\ {M m}_{i i} {\overset{\cdot \cdot}{\overset{~ ~}{ω ω}}}_{i i} = = {P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I} \end{matrix} \\ {P P}_{e e m m i i} = = \frac{{E E.}_{i i} {E E.}_{j j}}{\sqrt{{E E.}_{i i}^{22} - - 22 {E E.}_{i i} {E E.}_{j j} cos cos (({θ θ}_{i i} - - {θ θ}_{j j})) + + {E E.}_{j j}^{22}}} \\ {I I}_{i i j j} = = \{\begin{matrix} {I I}_{i i j j},, {I I}_{i i j j} \leq \leq {I I}_{max max} \\ {I I}_{max max},, {I I}_{i i j j} > > {I I}_{max max} \end{matrix} \end{matrix} - - - - - - ((11))$

式中，E代表节点电压，θ_i、分别为第i台发电机相对于惯性中心的转子角度以及角频率，单位分别为rad、rad·s-¹、M_i为第i台发电机的惯性常数，单位为s²·rad-¹，P_Mi代表第i台虚拟同步发电机的机械功率，P_emi代表第i台虚拟同步发电机的电磁功率，逆变器最大电流为I_max且θ_i＝δ_i-δ₀， In the formula, E represents the node voltage, θ _i , are the rotor angle and angular frequency of the i-th generator relative to the center of inertia, in units of rad, rad·s- ¹ , and Mi is the inertia constant of the _i -th generator, in units of s ² ·rad- ¹ , P _Mi represents the mechanical power of the i-th virtual synchronous generator, P _emi represents the electromagnetic power of the i-th virtual synchronous generator, the maximum current of the inverter is I _max and θ _i =δ _i -δ ₀ ,

式中，δ_i为第i台发电机的转子角速度，单位为rad，ω_i为角频率，单位为rad·s^-1,P_COI为惯性中心的加速功率，m为虚拟同步发电机的数量，δ₀、ω₀分别为惯性中心的转子角速度和角频率。In the formula, δ _i is the angular velocity of the rotor of the ith generator in rad, ω _i is the angular frequency in rad s ^-1 , P _COI is the acceleration power at the center of inertia, and m is the number of virtual synchronous generators , δ ₀ , ω ₀ are the angular velocity and angular frequency of the rotor at the center of inertia, respectively.

②非线性负载模型② Nonlinear load model

$\begin{matrix} {P P}_{l l i i} = = {f f}_{p p i i} (({E E.}_{i i})) \\ {Q Q}_{l l i i} = = {f f}_{q q i i} (({E E.}_{i i})) \end{matrix} - - - - - - ((22))$

式中，P_li，Q_li分别为负荷吸收的无功和有功功率。In the formula, P _li , Q _li are the reactive power and active power absorbed by the load, respectively.

③网络方程③Network equation

从节点i注入网络中的有功功率及无功功率为：The active power and reactive power injected into the network from node i are:

$\begin{matrix} 00 = = {P P}_{i i} + + {Σ Σ}_{j j = = 11}^{n no} {B B}_{i i j j} {E E.}_{i i} {E E.}_{j j} sin sin (({θ θ}_{i i} - - {θ θ}_{j j})) \\ 00 = = {Q Q}_{i i} + + {Σ Σ}_{j j = = 11}^{n no} {B B}_{i i j j} {E E.}_{i i} {E E.}_{j j} cos cos (({θ θ}_{i i} - - {θ θ}_{j j})) \\ i i = = m m + + 11,, ... ...,, n no \end{matrix} - - - - - - ((33))$

式中，P_i，Q_i分别为节点i注入网络中的有功功率及无功功率。In the formula, P _i and Q _i are the active power and reactive power injected into the network by node i respectively.

第二步，得到暂态能量函数。The second step is to obtain the transient energy function.

利用首次积分法构建暂态能量函数。将转子运动方程写成一阶方程组，利用首次积分法对两边同时进行积分，并将逆变器电流饱和因素考虑在内构造计及阻尼的能量型李雅普诺夫函数如式(4)～(7)，对暂态能量函数中和路径有关的不可积项采用线性路径方法近似，具体过程如下：The transient energy function is constructed using the first integration method. Write the rotor motion equation as a first-order equation group, use the first integral method to integrate both sides at the same time, and take the inverter current saturation factor into account to construct the energy-type Lyapunov function including damping, such as formulas (4)-(7 ), the non-integrable items related to the path in the transient energy function are approximated by the linear path method, and the specific process is as follows:

$V V (({\overset{~ ~}{ω ω}}_{g g i i},, θ θ,, E E.)) = = {V V}_{K K} + + {V V}_{P P} - - - - - - ((44))$

V_P＝V_P1+V_P2+V_P3+V_P4+V_damping (5)V _P ＝V _P1 +V _P2 +V _P3 +V _P4 +V _damping (5)

V_P4＝V_P41+V_P42+V_P43 (6)V _P4 ＝V _P41 +V _P42 +V _P43 (6)

$\{\begin{matrix} {V V}_{K K} = = \frac{11}{22} {Σ Σ}_{i i = = 11}^{m m} {M m}_{i i} {\overset{~ ~}{ω ω}}_{g g i i}^{22} \\ {V V}_{P P 11} = = - - {Σ Σ}_{i i = = 11}^{m m} {P P}_{{M m}_{i i}} (({θ θ}_{i i} - - {θ θ}_{i i}^{s the s})) \\ {V V}_{P P 22} = = {Σ Σ}_{i i = = 11}^{n no + + m m} {P P}_{i i} (({θ θ}_{i i} - - {θ θ}_{i i}^{s the s})) \\ {V V}_{P P 33} = = {Σ Σ}_{i i = = 11}^{n no + + m m} \frac{Q Q i i}{a a {(({E E.}_{i i}^{s the s}))}^{a a}} (({(({E E.}_{i i}))}^{a a} - - {(({E E.}_{i i}^{s the s}))}^{a a})) \\ {V V}_{P P 44} = = {Σ Σ}_{i i = = 11}^{m m} {Σ Σ}_{j j = = i i + + 11}^{m m} \frac{{E E.}_{i i} {E E.}_{j j}}{\sqrt{{E E.}_{i i}^{22} - - 22 {E E.}_{i i} {E E.}_{j j} cos cos (({θ θ}_{i i} - - {θ θ}_{j j})) + + {E E.}_{j j}^{22}}} {I I}_{i i j j} \\ - - {Σ Σ}_{i i = = 11}^{m m} {Σ Σ}_{j j = = i i + + 11}^{m m} \frac{{E E.}_{i i}^{s the s} {E E.}_{j j}^{s the s}}{\sqrt{{(({E E.}_{i i}^{s the s}))}^{22} - - 22 {E E.}_{i i}^{s the s} {E E.}_{j j}^{s the s} cos cos (({θ θ}_{i i}^{s the s} - - {θ θ}_{j j}^{s the s})) + + {(({E E.}_{j j}^{s the s}))}^{22}}} {I I}_{i i j j} \\ - - {Σ Σ}_{i i = = m m}^{n no + + m m - - 11} {Σ Σ}_{j j = = i i + + 11}^{n no + + m m} {B B}_{i i j j} (({E E.}_{i i} {E E.}_{j j} cos cos (({θ θ}_{i i} - - {θ θ}_{j j})) - - {E E.}_{i i}^{s the s} {E E.}_{j j}^{s the s} cos cos (({θ θ}_{i i}^{s the s} - - {θ θ}_{j j}^{s the s})))) \\ - - \frac{11}{22} {Σ Σ}_{i i = = 11}^{n no + + m m} {B B}_{i i i i} (({E E.}_{i i}^{22} - - {(({E E.}_{i i}^{s the s}))}^{22})) \\ {V V}_{d d a a m m p p i i n no g g} = = {Σ Σ}_{i i = = 11}^{m m} {&Integral; &Integral;}_{{t t}_{s the s}}^{t t} {D D.}_{i i} {ω ω}_{i i} \frac{{dθ dθ}_{i i}}{d d t t} = = {Σ Σ}_{i i = = 11}^{m m} {&Integral; &Integral;}_{{θ θ}_{i i,, s the s}}^{{θ θ}_{i i}} {D D.}_{i i} {ω ω}_{i i} {dθ dθ}_{i i} \end{matrix} - - - - - - ((77))$

其中，V为所述的能量函数，V_K为虚拟同步发电机的动能，V_P为系统总的势能，V_P1为全部虚拟同步发电机机械功率输入引起的转子势能，V_P2为全部有功负载引起的势能，V_P3为全部无功负载引起的势能，V_P4为储藏与网络中的势能，s代表稳定平衡点，i、j为节点的索引值，n为节点个数，D_i为第i台发电机的阻尼，B_ij、B_ii分别为节点i、j间的互导纳和节点i的自导纳，E_i、E_j分别是节点i、j处的电压，为第i台发电机的角频率，θ_i为第i台发电机相对于惯性中心的转子角度，其中a为常数取值一般为2。Among them, V is the energy function mentioned above, V _K is the kinetic energy of the virtual synchronous generator, V _P is the total potential energy of the system, V _P1 is the rotor potential energy caused by the mechanical power input of all virtual synchronous generators, and V _P2 is the total active load V _P3 is the potential energy caused by all reactive loads, V _P4 is the potential energy in storage and network, s represents the stable equilibrium point, i, j are the index values of nodes, n is the number of nodes, D _i is the The damping of generator i, B _ij and B _ii are the mutual admittance between nodes i and j and the self-admittance of node i respectively, E _i and E _j are the voltages at nodes i and j respectively, is the angular frequency of the i-th generator, θ _i is the rotor angle of the i-th generator relative to the center of inertia, where a is a constant and the value is generally 2.

第三步，得到临界能量值。The third step is to obtain the critical energy value.

结合图3，采用BCU方法获得临界能量值步骤如下Combined with Figure 3, the steps to obtain the critical energy value using the BCU method are as follows

步骤3-1，针对公式(1)电力系统模型，写出其收缩系统如式(8)Step 3-1, for the power system model of formula (1), write out its contraction system as formula (8)

$\{\begin{matrix} {\overset{\cdot &Center Dot;}{θ θ}}_{i i} = = {P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I} = = {f f}_{i i} ((θ θ)) \\ {θ θ}_{m m} = = \frac{{Σ Σ}_{i i = = 11}^{m m - - 11} {M m}_{i i} {θ θ}_{i i}}{{M m}_{m m}} \end{matrix} - - - - - - ((88))$

接着运用故障时的轨线求取出口点θ^EP,它是投影轨线存在收缩系统稳定边界上的一点，由公式(1)可以得到其故障时轨线。检测出口点θ^EP是由投影轨线到达第一个局部势能最大值。也就是由公式(1)得到的θ值带入故障后功率偏差量方程：Then use the trajectory at the time of failure Obtain the exit point θ ^EP , which is a point where the projected trajectory exists on the stable boundary of the shrinkage system, and its fault-time trajectory can be obtained from formula (1). The detection exit point θ ^EP is reached by the projected trajectory to the first local potential energy maximum. That is, the θ value obtained by formula (1) is brought into the power deviation equation after the fault:

${P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I} = = {f f}_{i i} ((θ θ)) - - - - - - ((99))$

当满足条件f_i*dθ＝0，也即后获得θ^EP。为了能够正确获取出口点，检测精度选取10^-5。When the condition f _i *dθ=0 is satisfied, that is Then get θ ^EP . In order to obtain the exit point correctly, the detection accuracy is selected as 10 ^-5 .

步骤3-2，以出口点θEP为初始点，对式(8)所得的收缩系统进行积分，沿着积分的曲线去寻找式(10)所示的第一个最小值:Step 3-2, taking the exit point θEP as the initial point, integrate the shrinkage system obtained by formula (8), and find the first minimum value shown in formula (10) along the integral curve:

$F f ((θ θ)) = = {Σ Σ}_{i i = = 11}^{m m} {f f}_{i i}^{22} ((θ θ)) = = {Σ Σ}_{i i = = 11}^{m m} {[[{P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I}]]}^{22} - - - - - - ((1010))$

得到的第一个最小值称为最小梯度点θ^MGP。The first minimum obtained is called the minimum gradient point θ ^MGP .

步骤3-3，以最小梯度点θ^EP为初值，用Newton-Raphson方法迭代求解(m-1)个故障后功率偏差量方程:Step 3-3, with the minimum gradient point θ ^EP as the initial value, use the Newton-Raphson method to iteratively solve (m-1) post-fault power deviation equations:

$\{\begin{matrix} {f f}_{i i} ((θ θ)) = = {P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I} = = 00 \\ {θ θ}_{m m} = = \frac{{Σ Σ}_{i i = = 11}^{m m - - 11} {M m}_{i i} {θ θ}_{i i}}{{M m}_{m m}} \end{matrix} - - - - - - ((1111))$

得到收缩系统的主导不稳定平衡点CUEP。将CUEP的状态参量代入式(7)暂态能量函数计算临界能量值V_cr。Get the dominant unstable equilibrium point CUEP of the contraction system. Substitute the state parameters of CUEP into the transient energy function of formula (7) to calculate the critical energy value V _cr .

第四步，得到临界切除时间并对暂态稳定性进行评估。The fourth step is to obtain the critical cut-off time and evaluate the transient stability.

令式(7)所得的暂态能量值等于第三步中所获得的临界能量值V_cr，从而求得临界故障切除时间并进行稳定性评估。Make the transient energy value obtained by formula (7) equal to the critical energy value V _cr obtained in the third step, so as to obtain the critical fault removal time and perform stability evaluation.

通过仿真来分析本发明提出的考虑虚拟转动惯量逆变器暂态能量函数分析方法的有效性。基于传统的IEEE 3机9节点系统的发电机换成虚拟同步发电机并含有电流限幅装置上，分别在不同的线路设置三相短路故障进行分析，得到不同的临界切除时间和系统能量包括系统的动能和势能见附图4，并与时域仿真得到的CCT结果进行对比。仿真系统图见附图2，虚拟同步发电机的参数如表1所示，得到的临界切除时间如表2所示。The effectiveness of the method for analyzing the transient energy function of the inverter considering the virtual moment of inertia proposed by the present invention is analyzed through simulation. The generator based on the traditional IEEE 3-machine 9-node system is replaced by a virtual synchronous generator with a current limiting device, and three-phase short-circuit faults are set up on different lines for analysis, and different critical cut-off times and system energy including system The kinetic and potential energies are shown in Figure 4, and compared with the CCT results obtained by time domain simulation. The simulation system diagram is shown in Figure 2, the parameters of the virtual synchronous generator are shown in Table 1, and the obtained critical cut-off time is shown in Table 2.

表1不同虚拟同步发电机参数Table 1 Parameters of different virtual synchronous generators

参数parameter 大小size 滤波电容filter capacitor 35uF35uF 寄生电阻Parasitic resistance 0.5′Ω0.5'Ω 负载load 50′Ω50'Ω 直流侧电压DC side voltage 650V650V 电网电压power voltage 380V380V 线路阻抗line impedance 0.01+j0.377′Ω0.01+j0.377'Ω PR比例常数kpPR proportionality constant kp 150150 PR积分时间常数kiPR integral time constant ki 44 电网频率f₀ Grid frequency f ₀ 50Hz50Hz 无功调节系数k_q Reactive adjustment coefficient k _q 1×10^-4 1×10 ^-4 电压调节系数k_v Voltage adjustment coefficient k _v 3.5×10^-2 3.5×10 ^-2 额定容量S₁ Rated capacity S ₁ 247.5MVA247.5MVA 额定容量S₂ Rated capacity S ₂ 192.0MVA192.0MVA 额定容量S₃ Rated capacity S ₃ 128.0MVA128.0MVA

表2考虑逆变器电流饱和不同的线路三相短路故障得到CCTTable 2 Considering the three-phase short-circuit fault of the line with different inverter current saturation, the CCT is obtained

故障节点failed node 故障线路fault line CCT暂态能量函数CCT transient energy function CCT时域仿真CCT Time Domain Simulation 44 4-54-5 0.1600.160 0.1620.162 44 4-64-6 0.1550.155 0.1600.160 55 5-75-7 0.1650.165 0.1680.168 66 6-96-9 0.1740.174 0.1790.179 77 7-87-8 0.1750.175 0.180.18 88 8-98-9 0.1300.130 0.1350.135

表3不考虑逆变器电流饱和不同的线路三相短路故障得到CCTTable 3 does not consider the three-phase short-circuit fault of the different lines of the inverter current saturation to get the CCT

故障节点failed node 故障线路fault line CCT暂态能量函数CCT transient energy function CCT时域仿真CCT Time Domain Simulation 44 4-54-5 0.1640.164 0.1680.168 44 4-64-6 0.1590.159 0.1640.164 55 5-75-7 0.1690.169 0.1730.173 66 6-96-9 0.1780.178 0.1820.182 77 7-87-8 0.1800.180 0.1840.184 88 8-98-9 0.1360.136 0.1400.140

根据表2和表3对比可知，传统的不考虑逆变器电流饱和因素的暂态能量函数分析方法得到的临界切除时间要比考虑逆变器电流饱和因素的时间长，如果故障切除时间长就会造成系统的不稳定，并且利用暂态能量函数法得到的临界切除时间(CCT)和时域仿真得到的误差很小，也验证本发明考虑虚拟转动惯量和逆变器电流饱和的暂态能量函数分析方法的有效性。因此，本发明方法更贴近实际，解决了传统的在线评估保守性导致系统不稳定的问题，是对时域仿真方法的一种重要补充。According to the comparison of Table 2 and Table 3, it can be seen that the critical cut-off time obtained by the traditional transient energy function analysis method without considering the inverter current saturation factor is longer than the time considering the inverter current saturation factor. If the fault cut-off time is long Can cause the instability of the system, and the critical cut-off time (CCT) that utilizes transient energy function method to obtain and the error that time domain simulation obtains is very little, also verifies that the present invention considers the transient energy of virtual moment of inertia and inverter current saturation Effectiveness of functional analysis methods. Therefore, the method of the invention is closer to reality, solves the problem of system instability caused by traditional online evaluation conservatism, and is an important supplement to the time domain simulation method.

Claims

1. A transient energy function analysis method considering inverter current saturation, comprising the following steps:

Step 1. Establish a detailed mathematical model of the power system with a virtual synchronous generator for the inverter, which includes a two-order model rotor motion equation similar to a traditional generator, as well as nonlinear load equations and network equations;

Step 2. Use the first integral method to construct the transient energy function, write the rotor motion equation as a first-order equation system, use the first integral method to integrate both sides simultaneously, and take the inverter current saturation factor into consideration to construct the energy that takes damping into account type Lyapunov function, the non-integrable items related to the path in the transient energy function are approximated by the linear path method;

Step 3, using the BCU method to obtain the critical energy value: first calculate the trajectory from the fault time to the fault clearing time and to the exit point, use the parameters at the fault clearing time to calculate the energy at this time, then calculate the minimum gradient point of the system after the fault and solve it according to the power The balance equation obtains the dominant unstable and unbalanced point, and obtains the critical energy value of the state parameter at this point;

Step 4, when the transient energy is equal to the energy at the dominant unstable equilibrium point, the critical cut-off time is obtained by solving.

2. method according to claim 1, is characterized in that, the mathematical model in step 1 is

\begin{matrix} \{\begin{matrix} {\overset{\cdot \cdot}{θ θ}}_{i i} = = {\overset{~ ~}{ω ω}}_{i i} \\ {M m}_{i i} {\overset{\cdot &Center Dot;}{\overset{~ ~}{ω ω}}}_{i i} = = {P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I} \end{matrix} \\ {P P}_{e e m m i i} = = \frac{{E E.}_{i i} {E E.}_{j j}}{\sqrt{{E E.}_{i i}^{22} - - 22 {E E.}_{i i} {E E.}_{j j} cos cos (({θ θ}_{i i} - - {θ θ}_{j j})) + + {E E.}_{j j}^{22}}} \\ {I I}_{i i j j} = = \{\begin{matrix} {I I}_{i i j j},, {I I}_{i i j j} \leq \leq {I I}_{max max} \\ {I I}_{max max},, {I I}_{i i j j} > > {I I}_{max max} \end{matrix} \end{matrix} - - - - - - ((11))

In the formula, E represents the node voltage, θ _i , are the rotor angle and angular frequency of the _i -th generator relative to the center of inertia, Mi is the inertia constant of the i-th generator, P _Mi represents the mechanical power of the i-th virtual synchronous generator, P _emi represents the i-th generator The electromagnetic power of the virtual synchronous generator, I _ij is the line current between node i and node j, and the maximum current of the inverter is I _max ;

θ _i =δ _i -δ ₀ , In the formula, δ _i is the rotor angular velocity of the i-th generator, ω _i is the angular frequency, P _COI is the acceleration power at the center of inertia, m is the number of virtual synchronous generators, δ ₀ and ω ₀ are the rotors at the center of inertia respectively Angular velocity and angular frequency.

3. The method according to claim 2, characterized in that, in step 2, the energy-type Lyapunov function considering the damping of the inverter current saturation factor into consideration is:

V V (({\overset{~ ~}{ω ω}}_{g g i i},, θ θ,, E E.)) = = {V V}_{K K} + + {V V}_{P P} - - - - - - ((22))

V _P ＝V _P1 +V _P2 +V _P3 +V _P4 +V _damping (3)

\{\begin{matrix} {V V}_{K K} = = \frac{11}{22} {Σ Σ}_{i i = = 11}^{m m} {M m}_{i i} {\overset{~ ~}{ω ω}}_{g g i i}^{22} \\ {V V}_{P P 11} = = - - {Σ Σ}_{i i = = 11}^{m m} {P P}_{{M m}_{i i}} (({θ θ}_{i i} - - {θ θ}_{i i}^{s the s})) \\ {V V}_{P P 22} = = {Σ Σ}_{i i = = 11}^{n no + + m m} {P P}_{i i} (({θ θ}_{i i} - - {θ θ}_{i i}^{s the s})) \\ {V V}_{P P 33} = = {Σ Σ}_{i i = = 11}^{n no + + m m} \frac{Q Q i i}{a a {(({E E.}_{i i}^{s the s}))}^{a a}} (({(({E E.}_{i i}))}^{a a} - - {(({E E.}_{i i}^{s the s}))}^{a a})) \\ {V V}_{P P 44} = = {Σ Σ}_{i i = = 11}^{m m} {Σ Σ}_{j j = = i i + + 11}^{m m} \frac{{E E.}_{i i} {E E.}_{j j}}{\sqrt{{E E.}_{i i}^{22} - - 22 {E E.}_{i i} {E E.}_{j j} cos cos (({θ θ}_{i i} - - {θ θ}_{j j})) + + {E E.}_{j j}^{22}}} {I I}_{i i j j} \\ - - {Σ Σ}_{i i = = 11}^{m m} {Σ Σ}_{j j = = i i + + 11}^{m m} \frac{{E E.}_{i i}^{s the s} {E E.}_{j j}^{s the s}}{\sqrt{{(({E E.}_{i i}^{s the s}))}^{22} - - 22 {E E.}_{i i}^{s the s} {E E.}_{j j}^{s the s} cos cos (({θ θ}_{i i}^{s the s} - - {θ θ}_{j j}^{s the s})) + + {(({E E.}_{j j}^{s the s}))}^{22}}} {I I}_{i i j j} \\ - - {Σ Σ}_{i i = = 11}^{n no + + m m - - 11} {Σ Σ}_{j j = = i i + + 11}^{n no + + m m} {B B}_{i i j j} (({E E.}_{i i} {E E.}_{j j} cos cos (({θ θ}_{i i} - - {θ θ}_{j j})) - - {E E.}_{i i}^{s the s} {E E.}_{j j}^{s the s} cos cos (({θ θ}_{i i}^{s the s} - - {θ θ}_{j j}^{s the s})))) \\ - - \frac{11}{22} {Σ Σ}_{i i = = 11}^{n no + + m m} {B B}_{i i i i} (({E E.}_{i i}^{22} - - {(({E E.}_{i i}^{s the s}))}^{22})) \\ {V V}_{d d a a m m p p i i n no g g} = = {Σ Σ}_{i i = = 11}^{m m} {&Integral; &Integral;}_{{t t}_{s the s}}^{t t} {D D.}_{i i} {ω ω}_{i i} \frac{{dθ dθ}_{i i}}{d d t t} = = {Σ Σ}_{i i = = 11}^{m m} {&Integral; &Integral;}_{{θ θ}_{i i,, s the s}}^{{θ θ}_{i i}} {D D.}_{i i} {ω ω}_{i i} {dθ dθ}_{i i} \end{matrix} - - - - - - ((44))

Among them, V is the energy function mentioned above, V _K is the kinetic energy of the virtual synchronous generator, V _P is the total potential energy of the system, V _P1 is the rotor potential energy caused by the mechanical power input of all virtual synchronous generators, and V _P2 is the total active load V _P3 is the potential energy caused by all reactive loads, V _P4 is the potential energy in storage and network, s represents the stable equilibrium point, i, j are the index values of nodes, n is the number of nodes, D _i is the The damping of generator i, B _ij and B _ii are the mutual admittance between nodes i and j and the self-admittance of node i respectively, E _i and E _j are the voltages at nodes i and j respectively, is the angular frequency of the i-th generator, θ _i is the rotor angle of the i-th generator relative to the center of inertia, where a is a constant and the value is generally 2.

4. method according to claim 3, is characterized in that, the concrete process of described step 3 is:

Step 3.1, obtain the contraction system of formula (1) as follows

\{\begin{matrix} {\overset{\cdot \cdot}{θ θ}}_{i i} = = {P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I} = = {f f}_{i i} ((θ θ)) \\ {θ θ}_{m m} = = \frac{{Σ Σ}_{i i = = 11}^{m m - - 11} {M m}_{i i} {θ θ}_{i i}}{{M m}_{m m}} \end{matrix} - - - - - - ((55))

Step 3.2, use the trajectory at the time of failure Find the exit point θ ^EP , θ ^EP is a point on the stable boundary of the contraction system where the projected trajectory exists, specifically as

The trajectory of θ ^EP at fault is obtained from the detailed mathematical model of the power system, and the value of θ obtained by formula (1) is brought into the power deviation equation after the fault. The detection exit point θ ^EP is obtained from the projection trajectory to the first local potential value; where the deviation equation is

{P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I} = = {f f}_{i i} ((θ θ)) - - - - - - ((66))

θ ^EP is obtained when the condition f _i *dθ=0 is satisfied;

Step 3.3, taking the exit point θ ^EP as the initial point, integrate the shrinkage system obtained by formula (6), and find the first minimum value shown in formula (7) along the integral curve

F f ((θ θ)) = = {Σ Σ}_{i i = = 11}^{m m} {f f}_{i i}^{22} ((θ θ)) = = {Σ Σ}_{i i = = 11}^{m m} {[[{P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I}]]}^{22} - - - - - - ((77))

The first minimum obtained is the minimum gradient point θ ^MGP .

Step 3.4, with the minimum gradient point θ ^MGP as the initial value, use the Newton-Raphson method to iteratively solve (m-1) post-fault power deviation equations

\{\begin{matrix} {f f}_{i i} ((θ θ)) = = {P P}_{{M m}_{i i}} - - {P P}_{e e m m i i} - - \frac{{M m}_{i i}}{{M m}_{T T}} {P P}_{C C O o I I} = = 00 \\ {θ θ}_{m m} = = \frac{{Σ Σ}_{i i = = 11}^{m m - - 11} {M m}_{i i} {θ θ}_{i i}}{{M m}_{m m}} \end{matrix} - - - - - - ((88))

Get the dominant unstable equilibrium point CUEP of the contraction system;

In step 3.4, the state parameters θ and ω of CUEP are substituted into formula (2) to obtain the critical energy value V _cr .

5. The method according to claim 4, characterized in that, in step 4, set the critical energy value _Vcr obtained in step 3.4 to be equal to the transient energy value of formula (2) gained, obtain corresponding θ, ω, when In domain simulation, it is found that the corresponding time values of θ and ω obtained are critical cutting time.