CN115693754A - Broadband harmonic instability analysis method for distributed power station - Google Patents

Abstract

The invention provides a broadband harmonic instability analysis method for a distributed power station, which comprises the following steps: firstly, establishing a distributed power station equivalent frequency domain impedance model according to a grid-connected inverter control strategy and power station element parameters, and deducing a transfer function matrix of current to a point of common coupling when an inverter side and a power grid side respectively act; secondly, indirectly judging whether a right half-plane pole exists in the transfer function matrix based on a generalized Nyquist stability criterion; finally, if the transfer function matrix has a right half-plane pole, the system is unstable. And if the system is stable, determining the cut-off frequency and the phase margin of the system through the intersection point of the Nyquist curve and the unit circle, and determining the potential harmonic amplification point of the system. The method can accurately judge the stability of the system and accurately position the potential broadband harmonic instability point of the system under the consideration of the dq axis coupling and line distribution parameter characteristics, and can provide important theoretical guidance for harmonic management and control.

Description

Broadband harmonic instability analysis method for distributed power station

Technical Field

The invention relates to the technical field of distributed power stations, in particular to a broadband harmonic instability analysis method for a distributed power station.

Background

Under the drive of a double-carbon target, the power generation of a distributed power station represented by photovoltaic and wind power is rapidly developed. However, while each unit of distributed generation feeds clean energy into the grid through the converter, the unit also injects wide-frequency-domain and high-frequency subharmonics into the grid. Due to the distribution parameter characteristics of the lines, broadband amplification is easily caused in the transmission process of harmonic waves through the lines, and the phenomenon is particularly prominent when the distributed power station is merged into a weak power grid through a long line. If the harmonic wave is continuously amplified, the voltage and current waveform is seriously deteriorated, and the harmonic wave is induced to be unstable. Harmonic instability can cause equipment damage and endanger the stable operation of the power grid. Therefore, the stability of the distributed power station and the potential broadband harmonic amplification point are accurately analyzed, and important theoretical guidance can be provided for harmonic management and control. However, the existing harmonic instability research is mainly focused on low frequency bands, and the line is usually simplified into a lumped parameter model, namely series impedance or concentrated pi type in the analysis process. Due to the fact that the lumped parameter model cannot describe the characteristics of the circuit outside the port in a wide frequency range, analysis of high-frequency harmonic amplification points is inaccurate, and even the analysis cannot reach the high-frequency harmonic amplification points. Furthermore, ignoring the coupling between the dq axes may not accurately determine system stability. Therefore, it is necessary to consider dq axis coupling and line distribution parametric models for distributed power station broadband harmonic instability analysis.

Disclosure of Invention

The invention solves the problem that the stability of a system cannot be accurately judged by neglecting the coupling between dq axes in the prior art, provides a distributed power station broadband harmonic instability analysis method, can accurately judge the stability of the system by considering the characteristics of the dq axis coupling and line distribution parameters, accurately positions potential broadband harmonic instability points of the system, and provides important theoretical guidance for harmonic management and control.

In order to realize the purpose, the following technical scheme is provided:

a broadband harmonic instability analysis method for a distributed power station comprises the following steps:

s1, acquiring power grid parameters, wherein the power grid parameters comprise a grid-connected inverter control strategy and power station element parameters;

s2, establishing a distributed power station equivalent frequency domain impedance model according to a grid-connected inverter control strategy and power station element parameters, and calculating a transfer function matrix of current to a point of common coupling when an inverter side and a power grid side act respectively;

s3, indirectly judging whether a right half-plane pole exists in the transfer function matrix based on a generalized Nyquist stability criterion; if yes, judging that the system is unstable; if not, entering S4;

and S4, determining that the system is stable, determining the cut-off frequency and the phase margin of the system through the intersection point of the Nyquist curve and the unit circle, and determining the potential harmonic amplification point of the system.

The method can accurately judge the stability of the system under the consideration of the dq axis coupling and the line distribution parameter characteristics, accurately position the potential broadband harmonic instability point of the system, and provide important theoretical guidance for harmonic management and control and management.

Preferably, the S2 specifically includes the following steps:

in a distributed power station in which n grid-connected inverters are connected to the grid through a single line, the current expression of a common coupling point is as follows:

will I _pcc The second term on the right is noted

Wherein:

a 2 × 1 matrix representing dq axes, and other bold italic variables all represent 2 × 2 matrices; I.C. A _ref And U _grid Respectively representing the control reference current and the power grid voltage of a single inverter; i represents a 2 × 2 identity matrix; represents G _eq 、Y _eq And Z _{eq_g} Respectively express inverter side equivalent current source coefficient, output admittance and electric wire netting equivalent impedance, and the expression respectively is:

wherein the content of the first and second substances,

G _inv and Y _inv Respectively representing the output equivalent current source coefficient and the equivalent admittance of a single inverter; z _c γ and l represent line characteristic impedance, propagation coefficient and length, respectively; r ₀ 、L ₀ And C ₀ Respectively arranging a resistor, an inductor and a capacitor under the unit length of the line; s represents the laplacian operator; w is a ₁ Representing the fundamental angular velocity; z is a linear or branched member _DT Representing the equivalent impedance of the inverter side step-up transformer; l is _MT 、L _g Respectively representing box transformer and grid equivalent inductance.

Preferably, the S3 specifically includes the following steps:

system stability is represented by a transfer function matrix G _eq 、Y _eq Co-determined with TF; transfer function matrix G _eq 、Y _eq TF indirectly judges determinant | G through generalized Nyquist stability criterion _eq |、|Y _eq And if the | TF | and the | TF | have a right half-plane pole, judging whether a system stability condition is met, if not, judging that the system is unstable, and if so, performing S4.

Preferably, the system stabilizing conditions are:

the first condition is as follows: determinant | G _eq I and Y _eq All I have no right half-plane pole;

and a second condition: the determinant | TF | has no right half-plane pole;

preferably, the transfer function matrix G _eq 、Y _eq And TF indirectly judges the determinant | G through generalized Nyquist stability criterion _eq |、|Y _eq Whether the | and | TF | have a right half-plane pole specifically includes the following steps:

step S301: calculate | G _inv |、|Y _inv I and I + Y _inv Z _DT If there is a right half-plane pole, distributing zero poles of | and if there is a right half-plane pole;

step S302: determination | G _eq I and Y _eq Whether a right half-plane pole exists;

when | Y _inv If there is no right half-plane pole, G _eq The rate matrix of (c) can be written as:

when | I + Y _inv Z _DT If there is no right half-plane pole, G _eq The rate matrix of (c) can be written as:

from the rate matrix, we derive:

1)|I+Y _inv Z _DT the right half-plane pole of | and | Y _inv The sum of the zero points of the right half plane of the L is the determinant L | G _eq The number of poles of the right half-plane;

2)|I+Y _inv Z _DT zero point of right half-plane of | and | Y _inv The sum of the poles of the right half plane of the matrix L' G _eq The number of right half plane poles;

if and only if L | G _eq Or L' | G _eq The characteristic function Nyquist curve of (1, j 0) is wound counterclockwise on the s plane ₊ -2N _- ) And determinant L | G _eq Or L' | G _eq If the number of poles of the right half-plane is the same, the condition one is satisfied, otherwise, the condition one is not satisfied; n is a radical of hydrogen ₊ And N _- Respectively representing the positive and negative crossing times of a Nyquist curve;

when w increases, if the Nyquist curve starts from-1 and passes downwards or upwards with a left negative real axis, the Nyquist curve is called 0.5 times of positive or negative passes when (2N) ₊ -2N _- ) The number of the right half-plane poles is the same as that of the return rate matrix determinant | L |, and the right half-plane pole does not exist in the determinant | I/(I + L) |;

step S303: it is determined whether | TF | has a right half-plane pole.

Preferably, the step S303 specifically includes the following steps:

calculating the rate matrix of the TF as follows:

if the first condition is satisfied, it can be determined that the L | TF has no right half plane pole, and if the nyquist curve of the rate matrix L | TF does not surround the (-1, j 0) point or clockwise surrounding is the same as the number of turns of counterclockwise surrounding, the L | TF transfer function matrix satisfies the GNC criterion, that is, the second stable condition is satisfied.

The invention has the beneficial effects that: the method can accurately judge the stability of the system under the consideration of the dq axis coupling and the line distribution parameter characteristics, accurately position the potential broadband harmonic instability point of the system, and provide important theoretical guidance for harmonic management and control and management.

Drawings

FIG. 1 is a schematic flow diagram of the present invention;

FIG. 2 is a single inverter topology and control strategy;

FIG. 3 is a topological diagram of a distributed power station structure and system parameters of n grid-connected inverters connected through a single line;

FIG. 4 shows the example of the method of the present invention _eq And L' | G _eq A nyquist plot of the characteristic function of;

FIG. 5 is a Nyquist plot of the L | TF signature function obtained using the method of the present invention in an exemplary embodiment;

FIG. 6 is a schematic diagram of a current waveform of a PCC point obtained by simulation in the embodiment of the present invention;

fig. 7 is a frequency spectrum of PCC point current a phase obtained by simulation in the embodiment of the present invention.

Detailed Description

The embodiment is as follows:

the embodiment provides a distributed power station harmonic instability analysis method considering line distribution parameters, and with reference to fig. 1, the method includes the following steps:

step S1: acquiring power grid parameters, wherein the power grid parameters comprise grid-connected inverter control strategies and power station element parameters; reference is made to the grid-connected inverter control strategy shown in fig. 2 and the distributed power station structure topology and system parameters shown in fig. 3.

Step S2: and establishing a distributed power station equivalent frequency domain impedance model according to a grid-connected inverter control strategy and power station element parameters, and deducing a transfer function matrix of PCC current to a point of common coupling when an inverter side and a power grid side respectively act. In a distributed power station with n grid-connected inverters connected through a single line, the PCC point current expression is as follows:

will I _pcc The second term on the right side is

For the sake of simplicity of expression, the voltage-to-current variables of all bold italics are represented as a 2 x 1 matrix of dq axes, e.g.,

the other bold italic variables all represent a 2 x 2 matrix. In the above formula, I _ref And U _grid Respectively representing the control reference current of a single inverter and the voltage of a power grid; i denotes a 2 × 2 identity matrix. Represents G _eq 、Y _eq And Z _{eq_g} The method comprises the following steps of respectively representing an inverter side equivalent current source coefficient, an output admittance and a power grid equivalent impedance, wherein the expressions are respectively as follows:

wherein, the first and the second end of the pipe are connected with each other,

in the above formula, G _inv And Y _inv Respectively representing the output equivalent current source coefficient and the equivalent admittance of a single inverter, and deriving the output equivalent current source coefficient and the equivalent admittance according to the control strategy and the control parameters of the inverter; z is a linear or branched member _c γ and l represent line characteristic impedance, propagation coefficient and length, respectively; r ₀ 、L ₀ And C ₀ Respectively arranging a resistor, an inductor and a capacitor under the unit length of the line; s represents the laplacian operator; w is a ₁ Representing the fundamental angular velocity; z _DT Representing the equivalent impedance of the inverter side step-up transformer; l is a radical of an alcohol _MT 、L _g Respectively representing box transformer and grid equivalent inductance.

And step S3: indirectly judging whether a right half-plane pole exists in a transfer function matrix based on a generalized Nyquist stability criterion; if yes, judging that the system is unstable; g _eq I _ref And Y _eq U _grid Can be respectively regarded as the harmonic emission source of the inverter side and the network side, and TF can be regarded as the mutual common action term of harmonic waves of both sides. If not, entering S4;

the step S3 specifically includes the following steps:

system stability is represented by a transfer function matrix G _eq 、Y _eq And TF, the system stability conditions are as follows:

the first condition is as follows: determinant | G _eq I and Y _eq All I have no right half-plane pole;

and (2) carrying out a second condition: the determinant | TF | has no right half-plane pole.

The transfer function matrix G obtained in step S2 _eq 、Y _eq And TF indirectly determines determinant | G through generalized Nyquist stability criterion (GNC) _eq |、|Y _eq If there is a right half-plane pole for | and | TF |.

Determining determinant | G _eq |、|Y _eq Whether or not | and | TF | exist for the right half-plane pole includes the steps of:

step S301: calculate | G _inv |、|Y _inv I and I + Y _inv Z _DT If there is a right half-plane pole, the system is unstable.

Step S302: determination | G _eq I and Y _eq If there is a right half-plane pole. Condition-requirement determinant | G for system stability _eq I and Y _eq All I have no right half plane pole, and the observation formula G _eq And Y _eq As can be seen, | G _eq I and Y _eq The | pole equations are consistent, so that only one of the two stability determinations needs to be determined.

When | Y _inv G if there is no right half-plane pole |, then _eq The return rate matrix of (c) can be written as:

when | I + Y _inv Z _DT G if there is no right half-plane pole |, then _eq The return rate matrix of (c) can be written as:

the observation of the two return rate matrixes shows that:

1)|I+Y _inv Z _DT the right half-plane pole point of | and | Y _inv The sum of the zero points of the right half plane of the L is the determinant L | G _eq The number of right half-plane poles.

2)|I+Y _inv Z _DT Zero point of right half-plane of | and | Y _inv The sum of the poles of the right half plane of the L is the determinant L' G _eq The number of right half-plane poles.

If and only if L | G _eq (or L' | G) _eq ) The total number of turns of the characteristic function Nyquist curve around the counterclockwise direction (-1, j 0) in the s plane or the positive and negative crossings (2N) of the characteristic function Nyquist curve of the positive frequency response Nyquist curve ₊ -2N _- ) And determinant L | G _eq (or L' | G) _eq ) And if the number of poles of the right half-plane is the same, the condition one is satisfied, otherwise, the condition one is not satisfied. The two return rate matrixes are in a reciprocal relation, which is actually the relation between the generalized Nyquist criterion and the inverse generalized Nyquist stability criterion, and the two judgment results are consistent.

N ₊ And N _- Respectively representing the positive and negative crossing times of the Nyquist curve. When w increases (from 0 to + ∞, positive frequency response), then the downward (upward) traversal of the nyquist curve starting with the negative real axis on the left from-1 is referred to as 0.5 positive traversal (negative traversal). When (2N) ₊ -2N _- ) The right half-plane pole is determined to be absent in the determinant | I/(I + L) |, i.e. the system is stable, otherwise the system is unstable.

Step S303: it is determined whether | TF | has a right half-plane pole. The second condition requires that the determinant | TF | has no right half-plane pole, that is, as long as the return rate matrix of TF satisfies the generalized nyquist stability criterion. The rate of return matrix for TF is:

observation formula G _eq 、Y _eq And L | TF, to know G _eq 、Y _eq And the L | TF pole equation is consistent, therefore, if the condition is satisfied, the L | TF can be judged to have no right half-plane pole.

On the premise that the first condition is satisfied, if none of nyquist curves of the rate matrix L | TF surrounds a point (-1, j 0) (or clockwise surrounding is the same as counterclockwise surrounding turns), the L | TF transfer function matrix satisfies the GNC criterion, that is, the second stable condition is satisfied. Finally, the cut-off frequency and the phase margin of the system can be given according to the Nyquist curve of L | TF.

And step S4: and when the system is judged to be stable and the harmonics on the two sides act together, the harmonic amplification characteristic at the PCC point depends on the cut-off frequency and the phase margin of the intersection point of the Nyquist curve of the characteristic function of the TF loop rate matrix and the unit circle. Harmonics at the cut-off frequency present a potential harmonic amplification risk, and the smaller the phase margin, the higher the harmonic amplification risk. In short, the cut-off frequency point corresponds to the system potential harmonic amplification frequency point.

The implementation case is as follows: the harmonic instability analysis method provided by the invention is verified in the embodiment, firstly, a system impedance model is constructed according to system parameters and is deduced to a current transfer function of a PCC point, and | G | is respectively calculated _inv |、|Y _inv I and I + Y _inv Z _DT And | distributing zero poles, and judging whether a right half plane zero pole exists or not, wherein the calculation result is shown in table 1.

TABLE 1

	Whether or not there is a right half-plane pole	Whether there is a right half-plane zero
			|G inv |	×	×
|Y inv |	×	V. presents one (6212.08, j0)
			|I+Y inv Z DT |	×	×

From table 1, it can be seen that: 1) Determinant | G _inv |、|Y _inv The right half plane pole does not exist in the I, and the inverter can stably operate to meet the design requirement; 2) Determinant L | G _eq There is one right half-plane pole; determinant L' | G _eq There is no right half-plane pole.

Due to L | G _eq The Nyquist curve of the method is complex, so that only the positive frequency response Nyquist curve is drawn; and L' | G _eq Its full frequency response nyquist curve is plotted. Obtaining L | G Using the method of the invention _eq And L' | G _eq The nyquist curve of (a) is shown in fig. 4. If and only if, L | G _eq Positive frequency nyquist curve of (N) positive crossing ₊ ) Negative crossing (N) _- ) The times are 0.5 times more; or, if and only if, L' | G _eq The full frequency nyquist curve of (1) is not enclosed (-1, j 0) or the clockwise enclosing number of turns is equal to the counterclockwise enclosing number of turns, then the condition one is true, otherwise it is not true. From FIG. 4, L | G _eq The number of positive crossings of the positive frequency Nyquist curve is more than that of negative crossings by 0.5 times; l' | G _eq OfNone of the stettes surrounds the point (-1, j 0), so condition one holds.

On the premise that the condition one is satisfied, it can be determined that L | TF has no right half plane pole, the broadband harmonic instability characteristic of the system depends on the condition that the nyquist curve of L | TF encloses (-1, j 0), if the nyquist curve encloses (-1, j 0), the system is unstable, and if the nyquist curve does not enclose, the system is stable. The nyquist curves of the L | TF characteristic functions obtained with the method of the present invention at different Short Circuit Ratios (SCR) are shown in fig. 5. As can be seen from FIG. 5, the SCR decreases, and the L | TF Nyquist curve gradually moves towards the enclosing (-1, j 0) direction, i.e., the smaller the SCR, the more easily the system is unstable; when SCR =4.71, the L | TF nyquist curve encloses (-1, j 0), and the system is unstable.

In the case of stable system, the intersection cut-off frequency of the L | TF nyquist curve obtained by the method of the present invention and the unit circle is shown in table 2, taking SCR =39.19, SCR =6.66, and SCR =5.52 as examples.

TABLE 2

Secondly, matlab software is used for building a distributed power station model of 8 grid-connected inverters which are connected in a grid mode through a single line and shown in the figure 3, and parameters are shown in a table 3; after the construction is finished, when the power grid is set at 0.2s, the SCR of the system is switched from 39.19 to 6.66 at 0.7s, is switched to 5.52 at 1.4s, and is switched to 4.71 at 1.4 s. The PCC point time domain simulation waveform is shown in fig. 6; the PCC point a phase spectrum for different SCRs is shown in fig. 7.

TABLE 3

As can be seen from fig. 6, when the SCRs are 39.19, 6.66 and 5.52, respectively, the current waveform is good, and when 1.4s is shifted to SCR =4.71, the current enters a constant amplitude oscillation link after diverging, the waveform quality is seriously deteriorated, and the system is unstable, which is consistent with the system instability when SCR =4.71 analyzed by the method of the present invention. Secondly, as can be seen from table 2 and fig. 7, under the condition of stable system, the harmonic current spectrum peaks of different SCRs are consistent with the theoretical analysis cut-off frequency result, which indicates that the potential harmonic method risk exists at the cut-off frequency, and also verifies the accuracy of the system potential harmonic amplification point obtained by the method of the present invention.

The above are only typical examples of the present invention, and besides, the present invention may have other embodiments, and all technical solutions formed by equivalent substitutions or equivalent transformations fall within the scope of the present invention as claimed.

Claims (6)

1. A broadband harmonic instability analysis method for a distributed power station is characterized by comprising the following steps:

s1, acquiring power grid parameters, wherein the power grid parameters comprise a grid-connected inverter control strategy and power station element parameters;

s2, establishing a distributed power station equivalent frequency domain impedance model according to a grid-connected inverter control strategy and power station element parameters, and calculating a transfer function matrix of current to a point of common coupling when an inverter side and a power grid side act respectively;

s3, indirectly judging whether a right half-plane pole exists in the transfer function matrix based on a generalized Nyquist stability criterion; if yes, judging that the system is unstable; if not, entering S4;

and S4, judging that the system is stable, determining the cut-off frequency and the phase margin of the system through the intersection point of the Nyquist curve and the unit circle, and determining the potential harmonic amplification point of the system.

2. The method according to claim 1, wherein the S2 specifically comprises the following steps:

in a distributed power station in which n grid-connected inverters are connected to the grid through a single line, the current expression of a common coupling point is as follows:

will I _pcc The second term on the right side is

Wherein:

a 2 × 1 matrix representing dq axes, and other bold italic variables all represent 2 × 2 matrices; i is _ref And U _grid Respectively representing the control reference current and the power grid voltage of a single inverter; i represents a 2 × 2 identity matrix; represents G _eq 、Y _eq And Z _{eq_g} Respectively express inverter side equivalent current source coefficient, output admittance and electric wire netting equivalent impedance, and the expression respectively is:

wherein the content of the first and second substances,

G _inv and Y _inv Respectively representing the output equivalent current source coefficient and the equivalent admittance of a single inverter; z _c γ and l represent line characteristic impedance, propagation coefficient and length, respectively; r ₀ 、L ₀ And C ₀ Respectively line unitResistance, inductance and capacitance under length; s represents the laplacian operator; w is a ₁ Representing the fundamental angular velocity; z _DT Representing the equivalent impedance of the inverter side step-up transformer; l is _MT 、L _g Respectively representing box transformer and grid equivalent inductance.

3. The method according to claim 2, wherein the S3 specifically comprises the following steps:

system stability is represented by a transfer function matrix G _eq 、Y _eq Co-determined with TF; transfer function matrix G _eq 、Y _eq TF indirectly judges determinant | G through generalized Nyquist stability criterion _eq |、|Y _eq And if the | TF | and the | TF | have a right half-plane pole, judging whether a system stability condition is met, if not, judging that the system is unstable, and if so, performing S4.

4. The method of claim 3, wherein the system stability conditions are:

the first condition is as follows: determinant | G _eq I and Y _eq All | have no right half-plane pole;

and (2) carrying out a second condition: the determinant | TF | has no right half-plane pole.

5. The method as claimed in claim 4, wherein the transfer function matrix G is a matrix of transfer functions _eq 、Y _eq TF indirectly judges determinant | G through generalized Nyquist stability criterion _eq |、|Y _eq Whether the | and | TF | have a right half-plane pole specifically includes the following steps:

step S301: calculate | G _inv |、|Y _inv I and I + Y _inv Z _DT Distributing zero poles of | if a right semiplane pole exists;

step S302: determination | G _eq I and Y _eq Whether a right half-plane pole exists;

when | Y _inv I do notIn the presence of the right half-plane pole, G _eq The return rate matrix of (c) can be written as:

when | I + Y _inv Z _DT G if there is no right half-plane pole |, then _eq The return rate matrix of (c) can be written as:

from the rate matrix:

1)|I+Y _inv Z _DT the right half-plane pole point of | and | Y _inv The sum of zero points of the right half plane of the L is the determinant L | G _eq The number of poles of the right half-plane;

2)|I+Y _inv Z _DT zero point of right half-plane of | and | Y _inv The sum of the poles of the right half plane of the L is the determinant L' G _eq The number of poles of the right half-plane;

if and only if L | G _eq Or L' | G _eq The total number of turns of the characteristic function Nyquist curve around the counterclockwise direction (-1, j 0) in the s plane or the positive and negative crossings (2N) of the characteristic function Nyquist curve of the positive frequency response Nyquist curve ₊ -2N _- ) And determinant L | G _eq Or L' | G _eq If the number of poles of the right half-plane is the same, the condition one is satisfied, otherwise, the condition one is not satisfied; n is a radical of ₊ And N _- Respectively representing the positive and negative crossing times of the Nyquist curve;

when w increases, if the Nyquist curve starts from-1 and passes downwards or upwards with a left negative real axis, the Nyquist curve is called 0.5 times of positive or negative passes when (2N) ₊ -2N _- ) The number of the right half-plane poles is the same as that of the return rate matrix determinant | L |, and the right half-plane poles do not exist in the determinant | I/(I + L) |;

step S303: it is determined whether | TF | has a right half-plane pole.

6. The method as claimed in claim 5, wherein said S303 includes the following steps:

calculating the rate matrix of the TF as follows:

if the first condition is satisfied, it can be determined that the L | TF has no right half plane pole, and if the nyquist curve of the rate matrix L | TF does not surround the (-1, j 0) point or clockwise surrounding is the same as the number of turns of counterclockwise surrounding, the L | TF transfer function matrix satisfies the GNC criterion, that is, the second stable condition is satisfied.

Images