CN106356859B - A kind of rectangular coordinate Newton load flow calculation method based on Matlab - Google Patents

Abstract

The invention discloses a kind of rectangular coordinate Newton load flow calculation methods based on Matlab, using matrix operation and complex operation.And Jacobian matrix and calculate node power are formed according to the programming feature of Matlab, comprising the following steps: calculate Jacobi's initial calculation matrix；Calculate node complex power and Injection Current phasor；Calculate initial Jacobian matrix piecemeal submatrix；Jacobian matrix submatrix diagonal element is corrected to voltage local derviation with node Injection Current phasor and voltage deviation；Jacobian matrix is formed by Jacobian matrix submatrix and is adjusted.The present invention is realized in Matlab platform, uses the Matlab various tools provided and function test and analysis result convenient for scientific research personnel.Program is write, modify and is debugged it was verified that the present invention had both facilitated scientific research personnel, while calculating speed provides an outstanding analysis tool also substantially close to the speed realized on C language platform for scientific research personnel.

Description

A Cartesian Coordinate Newton Method Power Flow Calculation Method Based on Matlab

technical field

The invention relates to a power system Newton method power flow calculation method, in particular to a Cartesian coordinate Newton method power flow calculation method suitable for research purposes.

Background technique

Power system power flow calculation is a basic calculation for studying the steady-state operation of the power system. It determines the operating state of the entire network according to the given operating conditions and network structure. The power flow calculation is also the basis of other analysis of the power system, such as safety analysis, transient stability analysis, etc., all use power flow calculation. Cartesian coordinate Newton method power flow calculation method is the most commonly used power flow calculation method, and researchers often conduct further research on the basis of Cartesian coordinate Newton method power flow calculation. Practical commercial software employs advanced techniques such as sparse matrix techniques and node-optimized numbering. Although these technologies can greatly increase the speed of power flow calculation and reduce memory usage, programming is very cumbersome and difficult to modify and maintain, and it is not easy to add new functions, so it is not suitable for researchers to use for research purposes.

Matlab software uses matrix as the most basic data unit, which can easily handle various matrix and vector operations, and can also handle complex number types conveniently and naturally. Its instruction expressions are very close to those commonly used in mathematics, and there are a large number of common and practical The function brings great convenience to programming. Matlab software is easy to use, short in code, easy to operate, easy to program and debug, powerful in calculation, and also has very powerful visual graphics processing and interactive functions, providing an efficient programming tool for scientific research and engineering applications. It has become the basic tool and preferred platform in many scientific fields, and has been widely used in various scientific and engineering computing fields. In order to meet the needs of more and more researchers who need to carry out further research based on the Cartesian Newton method power flow calculation on the Matlab platform, there is an urgent need for a Cartesian coordinate Newton method power flow calculation based on Matlab software that is easy to program, modify and debug. Calculation method.

According to the characteristics of the power system nodes, the power flow calculation divides the power system nodes into three categories: nodes whose active power and reactive power are known, nodes whose voltage amplitude and voltage phase angle are unknown are called PQ nodes; nodes whose active power and voltage amplitude are unknown Nodes with known node reactive power and unknown voltage phase angle are called PV nodes; nodes with known node voltage amplitude and voltage phase angle but unknown node active power and reactive power are called balanced nodes.

The Newton method power flow calculation is divided into two categories: the node voltage in the Newton method power flow calculation uses polar coordinates The calculation method when expressed is called the polar coordinate Newton method power flow calculation method; the node voltage in the Newton method power flow calculation adopts rectangular coordinates The calculation method when expressed is called Cartesian coordinate Newton method power flow calculation method. Cartesian Newton method power flow calculation main equation is as follows:

The nodal admittance matrix is:

In the formula, Y _ik is the node admittance matrix element. When the subscript i≠k, it is the mutual admittance of node i and node k. When the subscript i=k, it is the self-admittance of node i; n is the number of nodes.

The node power equation is:

In the formula, P _i and Q _i are the node active power and reactive power of node i respectively; e _i and e _k are the node voltage phasor real parts of node i and node k respectively; f _i and f _k are node i and the node voltage phasor imaginary part of node k; G _ik , B _ik are the real part and imaginary part of node admittance matrix element Y _ik respectively.

The power deviation and voltage deviation equations are:

In the formula, ΔP _i and _{ΔQ i} are the node active power deviation and reactive power deviation of node i respectively; _{ΔU i} ² is the node voltage amplitude deviation of node i; Active power and injected reactive power; U _iS is the node voltage amplitude given by node i; m is the number of PQ nodes.

The power deviation and voltage deviation equations can also be expressed as:

In the formula, a _i and b _i are the real part and imaginary part of the injected current phasor of node i respectively, which is

The corrected equations are:

In the formula, J is the Jacobian matrix, and H, N, M, L, R, K are the block sub-matrices of the Jacobian matrix. The calculation formula of each element of the Jacobian matrix is:

when j≠i

H _ij ＝-G _ij e _i -B _ij f _i (7)

N _ij ＝B _ij e _i -G _ij f _i (8)

M _ij ＝B _ij e _i -G _ij f _i (9)

L _ij =G _ij e _i +B _ij f _i (10)

R _ij =0 (11)

K _ij =0 (12)

when j=i

H _ii ＝-G _ii e _i -B _ii f _i -a _i (13)

N _ii =B _ii e _i -G _ii f _i -b _i (14)

M _ii ＝B _ii e _i -G _ii f _i +b _i (15)

L _ii =G _ii e _i +B _ii f _i -a _i (16)

R _ii =-2e _i (17)

K _ii =-2f _i (18)

As shown in Figure 1-2, the existing Cartesian coordinate Newton method power flow calculation method mainly includes the following steps:

A. Raw data input and voltage initialization;

The original data includes line and transformer branch data, node injection active power and reactive power, node voltage amplitude, node reactive power compensation data, as well as convergence accuracy and maximum number of iterations.

The voltage initialization adopts a flat start, that is, the real part of the node voltage of the PV node and the balance node takes a given value, and the real part of the node voltage of the PQ node takes 1.0; the imaginary part of all node voltages takes 0.0. The units here are per unit.

B. Form a node admittance matrix;

Form the node admittance matrix shown in formula (1) according to the input line and transformer branch data.

C. Form the Jacobian matrix;

Calculate each element of the Jacobian matrix according to formula (7) ~ formula (18) and formula (5).

D. Calculate node power, power deviation and voltage deviation;

Calculate the node power according to formula (2), and calculate the node power deviation and node voltage deviation according to formula (3).

E. Solve the equation and correct the real part e and imaginary part f of the node voltage;

Solve the correction equation group (6), and obtain the column vector Δe of the correction amount of the real part of the voltage and the column vector Δf of the correction amount of the imaginary part of the voltage.

The voltage correction formula is:

In the formula, the superscript (t) represents the value of the t-th iteration; Δe _i and Δf _i are the correction amount of the real part of the node voltage and the correction amount of the imaginary part of the node voltage of node i, respectively.

F. Determine whether the maximum unbalance |ΔP| _max , |ΔQ| _max and |ΔU ² | _max are all less than the convergence precision ε; if they are all less than the convergence precision ε, proceed to step G, otherwise return to step C for the next iteration;

G. Calculate the active power and reactive power of the balance node and the reactive power of the PV node, calculate the active power and reactive power of each branch, and end.

Cartesian coordinate Newton method power flow calculation software directly using the above principles has a slow calculation speed, and commercially used Cartesian coordinate Newton method power flow calculation software uses sparse matrix technology and node optimization numbering technology, which is relatively complicated and is not suitable for scientific researchers to use this as a basis for further research. Conduct scientific research. Therefore, Chinese patent ZL201010509556.5 proposes a Newton method power flow calculation method suitable for research purposes, and provides an easy-to-modify and maintain Newton method power flow calculation method for researchers who conduct further research based on polar coordinate Newton method power flow calculations , whose characteristics are as follows:

(1) Does not use sparse matrix technology and node optimization numbering, which greatly reduces the difficulty of implementing the algorithm;

(2) Avoid unnecessary calculations through simple logic judgments, and improve the calculation speed of power flow calculations.

The method proposed in Chinese patent ZL201010509556.5 provides an easy-to-modify and maintain polar coordinate Newton method power flow calculation method for researchers who conduct further research on the basis of polar coordinate Newton method power flow calculation. This method is fast when implemented in a compiled programming language such as C language, but the calculation speed is very slow when implemented in an interpreted programming language such as Matlab. At the same time, this patent does not make full use of the characteristics of Matlab that is good at matrix operations and complex operations. Therefore, a Cartesian coordinate Newton method power flow calculation method that makes full use of the characteristics of Matlab and calculates quickly is needed for researchers who conduct scientific research on the Matlab platform.

Contents of the invention

In order to solve the above-mentioned problems existing in the prior art, the present invention proposes a Matlab-based Cartesian Newton method power flow calculation method, which can make full use of Matlab's unique characteristics of being good at matrix operations and complex number operations, and has the advantage of faster calculation speed Cartesian coordinate Newton method power flow calculation method.

In order to achieve the above object, the technical solution of the present invention is as follows: a Matlab-based Cartesian Newton method power flow calculation method, using matrix operations and complex operations. Include the following steps:

A. Raw data input and voltage initialization;

B. Form a node admittance matrix;

C. Form Jacobian matrix and calculate node power;

Matlab is good at matrix operations and complex number operations, so using Matlab programming, the Jacobian matrix calculation method based on matrix operations and complex number operations should be derived.

The elements of the Jacobian matrix are related to the node type. When forming the Jacobian matrix, the node type must be judged, and which nodes need to form the Jacobian matrix elements are determined according to the node type. Such processing is easy to handle for an algorithm realized by loops, but it is not suitable for a method of forming a Jacobian matrix by an overall matrix operation. Therefore, when the present invention forms the Jacobian matrix, it does not judge the node type, but forms Jacobian matrix elements for all nodes, and then removes unnecessary rows and columns.

The formulas for calculating Jacobian matrix elements and node powers using matrix operations are derived below.

Analysis of formula (7) ~ formula (10), we can get

M _ij = N _ij (21)

L _ij ＝-H _ij (22)

Therefore, H _ij and N _ij are calculated first, and after H _ij and N _ij are calculated, M _ij and L _ij can be obtained naturally.

Before deriving the formula for computing the Jacobian matrix by matrix operations, let's see how each element of the Jacobian matrix is represented by a complex number or a phasor.

From formula (7) and formula (8), it can be seen that each element of the Jacobian matrix should be composed of Y _ij and the conjugate of the product of get. Therefore, H _ij , N _ij , M _ij , and L _ij can be generated by the following formula:

where the superscript (^) represents the conjugate of a complex number.

Equation (23) can be regarded as the element in row i and column j of matrix J ⁰ obtained by multiplying the node admittance matrix with the corresponding position element of the following formula:

Equation (24) can be realized by a matrix filling function repmat of Matlab, so the block submatrix of Jacobian matrix can be obtained by the following formula:

In the formula, J ⁰ is the Jacobian initial calculation matrix; repmat is the matrix filling function of Matlab, and 1 in the function means filling time The row is repeated once, and n in the function indicates when filling The column of is repeated n times; .* indicates that the elements of the corresponding rows and columns of the two matrices are multiplied.

The block sub-matrix of the initial Jacobian matrix obtained from J ⁰ is:

H ⁰ ＝-Re(J ⁰ ) (26)

N ⁰ =Im(J ⁰ ) (27)

M ⁰ =Im(J ⁰ ) (28)

L ⁰ =Re(J ⁰ ) (29)

In the formula, H ⁰ , N ⁰ , M ⁰ , and L ⁰ are block sub-matrixes of the initial Jacobian matrix; Re means to take the real part of the matrix elements; Im means to take the imaginary parts of the matrix elements.

The off-diagonal elements of the block sub-matrix of the initial Jacobian matrix obtained from formulas (26) to (29) are already Jacobian matrix elements, and the diagonal elements need to be corrected.

The first and second items on the right side of the equations in formula (13) ~ formula (16) are the diagonal elements of H ⁰ , N ⁰ , M ⁰ , L ⁰ , so only the obtained H ⁰ , N ⁰ , M ⁰ and L ⁰ are corrected by a _i and b _i , and the diagonal elements of the block sub-matrix of the Jacobian matrix can be obtained.

Replace the control variable summed in the nodal power equation (2) by j, and rearrange as:

Equation (23) expands to:

Observing formula (30) and formula (31), the complex power can be obtained:

In the formula, is a complex power column vector.

The injected current phasor is:

In the formula, is the column vector of the injection current phasor; is the node voltage column vector.

Since a _i and b _i are the real part and imaginary part of the injected current phasor of node i respectively, so

In the formula, is the injected current phasor of node i.

The diagonal elements of the block submatrix of the Jacobian matrix are corrected with the injected current phasor as follows:

The off-diagonal elements of the voltage deviation to the voltage partial conduction are 0, and the diagonal elements are:

R _ii =-2e _i i = m+1,...,n-1 (40)

K _ii =-2f _i i = m+1,...,n-1 (41)

The elements from the m+1 to n-1 rows of R and K can be placed in the corresponding rows of M and L to form M' and L':

M _i ′ _j =M _ij i=1,...,mj=1,...,n-1 (42)

L _i ′ _j =L _ij i=1,...,mj=1,...,n-1 (43)

M _i ′ _j =R _ij i=m+1,...,n-1 j=1,...,n-1 (44)

L _i ′ _j =K _ij i=m+1,...,n-1 j=1,...,n-1 (45)

Forming the Jacobian matrix and calculating node power includes the following steps:

C1. Calculating the Jacobian initial calculation matrix J ⁰ ;

C2. Compute node complex power;

C3, calculate the initial Jacobian matrix block sub-matrix by J ⁰ ;

C4. Calculate the injection current phasor;

C5, using the injected current phasor to modify the diagonal elements of the Jacobian matrix block submatrix;

C6, calculate R, K, form M', L';

C7, the Jacobian matrix is formed by the Jacobian matrix block sub-matrix;

C8. Adjust the Jacobian matrix, delete the rows corresponding to the active power deviation and the reactive power deviation of the balance node; delete the columns corresponding to the correction amount of the real part of the voltage at the balance node and the correction amount of the imaginary part of the voltage, and end.

D. Calculate node power deviation and node voltage deviation;

Equation (4) calculates the node power deviation and voltage deviation equation and writes the matrix operation as follows:

In the formula, ΔP and _ΔQ are the column vectors of node active power deviation and reactive power deviation respectively; _ΔU ² is the column vector of node voltage amplitude deviation; Inject reactive power column vector; U _S is the column vector of node voltage given value; is the node voltage column vector.

The calculated node power deviation vectors ΔP and ΔQ remove the PV node reactive power deviation and the balance node active power deviation and reactive power deviation; the node voltage amplitude deviation column vector removes the node voltage amplitude deviation of the PQ node and the balance node.

E. Solve the equation and correct the real part e and imaginary part f of the voltage;

Directly invoke the algorithm of solving linear equations of Matlab software to solve the correction equations (6), and obtain the column vector Δe of the correction amount of the real part of the node voltage and the column vector Δf of the correction amount of the imaginary part of the node voltage.

The formulas (19) and (20) for voltage correction are rewritten into matrix form as:

e ^(t+1) = e ^(t) -Δe ^(t) (48)

f ^(t+1) = f ^(t) -Δf ^(t) (49)

In the formula, the superscript (t) represents the value of the t-th iteration; Δe and Δf are the column vectors of the correction amount of the real part of the node voltage and the column vector of the correction amount of the imaginary part of the node voltage, respectively.

F. Determine whether the maximum unbalance |ΔP| _max , |ΔQ| _max and |ΔU ² | _max are all less than the convergence precision ε; if they are all less than the convergence precision ε, proceed to step G, otherwise return to step C for the next iteration.

G. Calculate the active power and reactive power of the balance node and the reactive power of the PV node, calculate the active power and reactive power of each branch, and end.

Compared with the prior art, the present invention has the following beneficial effects:

1, the method that the present invention proposes realizes in Matlab platform, is convenient to scientific research personnel to use various tools and the function that Matlab provides to test and analyze calculation result.

2. The method proposed by the present invention adopts matrix operation and complex number operation, reduces the program code, simplifies programming, makes the program clearer, and is convenient for scientific researchers to modify the program, debug and improve the program, and add new functions.

3. Since Matlab optimizes matrix operations, using matrix operations is much faster than programming by matrix element loop operations. At the same time, directly calling Matlab's equation solving algorithm also greatly improves the calculation speed. Practice has proved that the method of the present invention is convenient for scientific research personnel to write, modify and debug programs, and at the same time, the calculation speed is basically close to the speed realized on the C language platform, providing an excellent analysis tool for scientific research personnel .

Description of drawings

The present invention has 4 accompanying drawings. in:

Fig. 1 is a flow chart of the existing Cartesian coordinate Newton method power flow calculation.

Fig. 2 is a flow chart of forming a Jacobian matrix by the existing rectangular coordinate Newton method.

Fig. 3 is a flow chart of the flow calculation of the Cartesian Newton method of the present invention.

Fig. 4 is a flow chart of forming Jacobian matrix and calculating node power in the present invention.

Detailed ways

The present invention will be further described below in conjunction with the accompanying drawings, and a modified 445-node actual system calculation example is calculated according to the flow shown in Figures 3-4.

445 nodes The actual large-scale power system has 445 nodes and 544 branches, including a large number of small impedance branches. In order to compare various methods, these small impedance branches are changed to normal impedance branches to meet the requirements of various methods.

The calculation example of the 445-node actual system is calculated by using the present invention and several comparison methods, and the convergence precision is 0.00001 during the calculation. Several power flow calculation algorithms are:

Method 1: The method of Chinese patent ZL201010509556.5 is changed to the form of Cartesian coordinates and implemented by Matlab language.

Method 2: The method of the Chinese patent ZL201010509556.5 is changed to the form of Cartesian coordinates and implemented in Matlab language, but the equation solving algorithm of Matlab is directly called.

Method 3: The method of Chinese patent ZL201010509556.5 is changed to Cartesian coordinates and implemented in Matlab language, but the equation solving algorithm directly invokes Matlab's equation solving algorithm, and at the same time pre-defines the dimension of the matrix variable.

Method 4: the method of the present invention.

The calculation time of several methods is shown in Table 1, and the calculation time does not include the time for data input and output and branch power calculation.

Table 1 Comparison of calculation time of several Cartesian Newton methods for power flow calculation

Power flow calculation method Calculation time (s) method 1 20.167 Method 2 4.730 Method 3 1.027 Method 4 0.562

As can be seen from Table 1, when the method of Chinese patent ZL201010509556.5 is changed to Cartesian coordinates and directly adopts Matlab to realize, the calculation time is very long; method, the calculation speed can be greatly improved. Matlab's equation-solving method is relatively mature and stable, which is conducive to the stability of the algorithm. The call of Matlab's equation-solving method is also very simple, which simplifies programming and makes the program clearer; Defining the dimension and avoiding the continuous expansion of the size of the matrix during the execution of the program can greatly improve the calculation speed. The calculation result of the present invention shows that the calculation speed can be further greatly improved by adopting the matrix operation technology when forming the Jacobian matrix, and the program is further simplified.

The method of the present invention can be implemented in any version of the MATLAB programming language, but it is recommended to use a newer version of the MATLAB language.

The present invention is not limited to this embodiment, and any equivalent ideas or changes within the technical scope disclosed in the present invention are listed in the protection scope of the present invention.

Claims (1)

1. A Matlab-based Cartesian Newton method power flow calculation method, comprising the following steps: A. Raw data input and voltage initialization; Raw data includes line and transformer branch data, node injected active power and reactive power, node voltage amplitude, node reactive power compensation data, as well as convergence accuracy and maximum number of iterations; The voltage initialization adopts a flat start, that is, the real part of the node voltage of the PV node and the balance node takes a given value, and the real part of the node voltage of the PQ node takes 1.0; the imaginary part of all node voltages takes 0.0; the unit here is the per unit value; The PQ node is a node with known active power and reactive power, unknown voltage amplitude and voltage phase angle, and the PV node is known with active power and voltage amplitude, and unknown reactive power and voltage phase angle node, the balanced node is a node whose voltage amplitude and voltage phase angle are known, and whose active power and reactive power are unknown; B. Form a node admittance matrix; Form the node admittance matrix shown in formula (1) according to the input line and transformer branch data; In the formula, Y _ik is the node admittance matrix element, when the subscript i≠k, it is the mutual admittance of node i and node k, when the subscript i=k, it is the self-admittance of node i; n is the number of nodes; It is characterized in that it also includes the following steps: C. Form Jacobian matrix and calculate node power; C1. Calculating the Jacobian initial calculation matrix J ⁰ ; In the formula, J ⁰ is the Jacobian initial calculation matrix; is the column vector of the node voltage conjugate; repmat is the matrix filling function of Matlab, and 1 in the function means when filling The row is repeated once, and n in the function indicates when filling The column of is repeated n times; .* indicates that the elements of the corresponding rows and columns of the two matrices are multiplied; C2. Compute node complex power; The node complex power column vector composed of the complex power of each node is: In the formula, is the node complex power column vector; the superscript ^ represents the conjugate of the complex number; C3, calculate the initial Jacobian matrix block sub-matrix by J0; The block sub-matrix of the initial Jacobian matrix obtained from J ⁰ is: H ⁰ ＝-Re(J ⁰ ) (4) N ⁰ =Im(J ⁰ ) (5) M ⁰ =Im(J ⁰ ) (6) L ⁰ =Re(J ⁰ ) (7) In the formula, H ⁰ , N ⁰ , M ⁰ , L ⁰ are block sub-matrixes of the initial Jacobian matrix; Re means to take the real part of the matrix element; Im means to take the imaginary part of the matrix element; C4. Calculate the injected current phasor; In the formula, is the column vector of the injection current phasor; C5. Use the injected current phasor to modify the diagonal elements of the block submatrix of the initial Jacobian matrix as follows: In the formula, is the injection current phasor of node i; m is the number of PQ nodes; C6, calculate R, K, form M', L'; The off-diagonal elements of the voltage deviation to the voltage partial conduction are 0, and the diagonal elements are: R _ii =-2e _i i = m+1,...,n-1 (13) K _ii =-2f _i i = m+1,...,n-1 (14) Replace the corresponding row of M with the elements from row m+1 to row n-1 of R to form M'; replace the corresponding row of L with elements from row m+1 to row n-1 of K to form L'; C7, the Jacobian matrix is formed by the Jacobian matrix block sub-matrix after modification; In the formula, J is the Jacobian matrix, H, N, M', L' are the block sub-matrixes of the modified Jacobian matrix; C8, adjust the Jacobian matrix, remove the row corresponding to the active power deviation and reactive power deviation of the balance node; remove the column corresponding to the correction amount of the real part of the voltage at the balance node and the correction amount of the imaginary part of the voltage, and end; D. Calculate node power deviation and node voltage deviation; Calculate the node power deviation and node voltage deviation as follows: In the formula, ΔP and _ΔQ are the column vectors of node active power deviation and reactive power deviation respectively; _ΔU ² is the column vector of node voltage amplitude deviation; Inject reactive power column vector; U _S is the column vector of node voltage given value; is the node voltage column vector; The PV node reactive power deviation and the balance node active power deviation and reactive power deviation are removed from the calculated node power deviation vectors ΔP and ΔQ; the node voltage amplitude deviation column vector removes the node voltage amplitude deviation of the PQ node and the balance node; E. Solve the equation and correct the real part e and imaginary part f of the voltage; The Jacobian matrix J obtained in step C and the node power deviation vectors ΔP and ΔQ and the node voltage amplitude deviation vector ΔU ² are obtained in step D, and the correction equation for power flow calculation is constructed as follows: Directly invoke the linear equations algorithm of Matlab software to solve the correction equations (17), and obtain the voltage real part correction vector Δe and the voltage imaginary part correction vector Δf; The real and imaginary parts of the node voltage are corrected as follows: e ^(t+1) = e ^(t) -Δe ^(t) (18) f ^(t+1) = f ^(t) -Δf ^(t) (19) In the formula, the superscript (t) represents the value of the tth iteration; F. Judging whether the maximum unbalance |ΔP| _max , |ΔQ| _max and |ΔU ² | _max are all less than the convergence precision ε; if they are all less than the convergence precision ε, proceed to step G, otherwise return to step C for the next iteration; G. Calculate the active power and reactive power of the balance node and the reactive power of the PV node, calculate the active power and reactive power of each branch, and end.