CN106356859A - Matlab-based rectangular coordinate newton power flow calculation method - Google Patents

Abstract

The invention discloses a Matlab-based rectangular coordinate newton power flow calculation method. Matrix operation and complex operation are adopted, a Jacobi matrix is formed according to the programming characteristic of Matlab, and the power of a node is calculated. The method comprises the following steps: calculating an initial Jacobi calculation matrix; calculating the complex power and the injection current phasor of the node; calculating partitioned sub-matrixes of the initial Jacobi matrix; correcting diagonal elements of the sub-matrixes of the Jacobi matrix by virtue of the injection current phasor of the node and a partial derivative of a voltage deviation to a voltage; forming the Jacobi matrix by virtue of the sub-matrixes of the Jacobi matrix, and performing regulation. The method is implemented on a Matlab platform, so that a researcher can conveniently test and analyze a calculation result by virtue of various tools and functions provided by Matlab. Practices show that convenience is brought to program compiling, modification and debugging of the researcher, the calculation speed is substantially approximate to the speed achieved on a C language platform and an excellent analysis tool is provided for the researcher.

Description

A kind of rectangular coordinate Newton load flow calculation method based on matlab

Technical field

The present invention relates to a kind of power system Newton load flow calculation method, particularly a kind of suitable research purpose use Rectangular coordinate Newton load flow calculation method.

Background technology

It is the basic calculating that research power system mesomeric state runs that electric power system tide calculates, and it is according to given operation Condition and network structure determine the running status of whole network.Load flow calculation is also the basis of other analyses of power system, such as pacifies Complete analysis, transient stability analysis etc. will use Load flow calculation.Rectangular coordinate Newton load flow calculation method is a kind of the most frequently used Tidal current computing method, scientific research personnel often by rectangular coordinate Newton Power Flow calculating based on further studied.Real Business software adopts the advanced techniques such as sparse matrix technology and node optimizing code.Although these technology can increase substantially The speed of Load flow calculation, reduction EMS memory occupation amount, but programming bothers and is difficult to change and safeguards very much, is difficult to increase new work( Can, thus be not suitable for scientific research personnel and use for research purposes.

Matlab software, with matrix for most basic data unit, can easily process various matrixes and vector operation, Easily can also naturally process complex data type, its instruction expression formula and conventional form in mathematics are very close to also in a large number Common practical function, brings convenience to programming.Matlab software is easy to use, code is short and small easy to operate it is easy to program And debugging, computing function is powerful, also has very powerful visualized graphs simultaneously and processes and Interactive function, is scientific research And engineer applied provides a kind of efficient programming tool, have become as basic tool and the first-selection of many scientific domains at present Platform, is widely used in various science and engineering calculation field.In order to adapt to increasing scientific research personnel's needs The demand further studied based on the calculating of rectangular coordinate Newton Power Flow on matlab platform, in the urgent need to one Plant the rectangular coordinate Newton load flow calculation method being easily programmed, change and debugging based on matlab software.

According to the feature of power system node, Load flow calculation is divided into 3 classes: node active power and nothing power system node Work(power is known, node voltage amplitude and the unknown node of voltage phase angle are referred to as pq node；Node active power and voltage magnitude Known, node reactive power and the unknown node of voltage phase angle are referred to as pv node；Node voltage amplitude and voltage phase angle are it is known that save Point active power and the unknown node of reactive power are referred to as balance nodes.

Newton Power Flow calculates and is divided into two classes: Newton Power Flow calculates interior joint voltage and adopts polar coordinate Computational methods during expression, referred to as polar coordinate Newton load flow calculation method；Newton Power Flow calculates interior joint voltage using straight Angular coordinateComputational methods during expression, referred to as rectangular coordinate Newton load flow calculation method.Rectangular coordinate newton The main equation of method Load flow calculation is as follows:

Bus admittance matrix is:

In formula, y_ikFor bus admittance matrix element, as subscript i ≠ k, it is node i and the transadmittance of node k, when subscript i During=k, it is the self-admittance of node i；N is nodes.

Node power equation is:

$\left\{\begin{array}{cc}{p}_i=e_i\underset{k=1}{\overset{n}{\mathrm{\σ}}}(g_{ik}{e}_k-b_{ik}{f}_k)+f_i\underset{k=1}{\overset{n}{\mathrm{\σ}}}(g_{ik}{f}_k+b_{ik}{e}_k)& i=1,...,n\\ q_i=f_i\underset{k=1}{\overset{n}{\mathrm{\σ}}}(g_{ik}{e}_k-b_{ik}{f}_k)-e_i\underset{k=1}{\overset{n}{\mathrm{\σ}}}(g_{ik}{f}_k+b_{ik}{e}_k)& i=1,...,n\end{array}\right.---\left(2\right)$

In formula, p_i、q_iIt is respectively node active power and the reactive power of node i；e_i、e_kIt is respectively node i and node k Node voltage phasor real part；f_iAnd f_kIt is respectively node i and the node voltage phasor imaginary part of node k；g_ik、b_ikIt is respectively node Admittance matrix element y_ikReal part and imaginary part.

Power deviation and voltage deviation equation are:

$\left\{\begin{array}{cc}{\mathrm{\δp}}_i=p_{is}-p_i=p_{is}-e_i\underset{k=1}{\overset{n}{\mathrm{\σ}}}(g_{ik}{e}_k-b_{ik}{f}_k)-f_i\underset{k=1}{\overset{n}{\mathrm{\σ}}}(g_{ik}{f}_k+b_{ik}{e}_k)& i=1,...,n-1\\ {\mathrm{\δq}}_i=q_{is}-q_i=q_{is}-f_i\underset{k=1}{\overset{n}{\mathrm{\σ}}}(g_{ik}{e}_k-b_{ik}{f}_k)+e_i\underset{k=1}{\overset{n}{\mathrm{\σ}}}(g_{ik}{f}_k+b_{ik}{e}_k)& i=1,...,m\\ {\mathrm{\δu}}_i^2=u_{is}^2-(e_i^2+f_i^2)& i=m+1,...,n-1\end{array}\right.---\left(3\right)$

In formula, δ p_i、δq_iIt is respectively node active power deviation and the reactive power deviation of node i；δu_i ²For node i Node voltage amplitude deviation；p_is、q_isIt is respectively node injection active power and the injection reactive power that node i gives；u_isFor The node voltage amplitude that node i gives；M is pq nodes.

Power deviation and voltage deviation equation can also be expressed as:

$\left\{\begin{array}{c}{\mathrm{\δp}}_i=p_{is}-e_i{a}_i-f_i{b}_i\\ {\mathrm{\δq}}_i=q_{is}-f_i{a}_i+e_i{b}_i\\ {\mathrm{\δu}}_i^2=u_{is}^2-(e_i^2+f_i^2)\end{array}\right.---\left(4\right)$

In formula, a_i、b_iIt is respectively the real part of injection current phasor and the imaginary part of node i, be

$\left\{\begin{array}{c}{a}_i=\underset{k=1}{\overset{n}{\σ}}(g_{ik}{e}_k-b_{ik}{f}_k)\\ b_i=\underset{k=1}{\overset{n}{\σ}}(g_{ik}{f}_k+b_{ik}{e}_k)\end{array}\right.---\left(5\right)$

Update equation group is:

$\left[\begin{array}{c}\mathrm{\δ}p\\ \mathrm{\δ}q\\ \mathrm{\δ}{u}^2\end{array}\right]=j\left[\begin{array}{c}\mathrm{\δ}e\\ \mathrm{\δ}f\end{array}\right]=\left[\begin{array}{cc}h& n\\ m& l\\ r& k\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}e\\ \mathrm{\δ}f\end{array}\right]---\left(6\right)$

In formula, j is Jacobian matrix, and h, n, m, l, r, k are the piecemeal submatrix of Jacobian matrix.The each unit of Jacobian matrix Plain computing formula is:

As j ≠ i

h_ij=-g_ije_i-b_ijf_i(7)

n_ij=b_ije_i-g_ijf_i(8)

m_ij=b_ije_i-g_ijf_i(9)

l_ij=g_ije_i+b_ijf_i(10)

r_ij=0 (11)

k_ij=0 (12)

As j=i

h_ii=-g_iie_i-b_iif_i-a_i(13)

n_ii=b_iie_i-g_iif_i-b_i(14)

m_ii=b_iie_i-g_iif_i+b_i(15)

l_ii=g_iie_i+b_iif_i-a_i(16)

r_ii=-2e_i(17)

k_ii=-2f_i(18)

As shown in Figure 1-2, existing rectangular coordinate Newton load flow calculation method, mainly comprises the steps that

A, initial data input and voltage initialization；

Raw data packets vinculum road and transformer branch data, node injection active power and reactive power, node voltage Amplitude, node reactive-load compensation data, and convergence precision, maximum iteration time.

Voltage initialization is started using flat, and that is, the node voltage real part of pv node and balance nodes draws definite value, pq node Node voltage real part take 1.0；The imaginary part of all node voltages all takes 0.0.Here unit adopts perunit value.

B, formation bus admittance matrix；

Circuit according to input and transformer branch data form the bus admittance matrix as shown in formula (1).

C, formation Jacobian matrix；

Calculate each element of Jacobian matrix by formula (7)～formula (18) and formula (5).

D, calculate node power and power deviation and voltage deviation；

By formula (2) calculate node power, by formula (3) calculate node power deviation and node voltage deviation.

E, solve equation and revise node voltage real part e and imaginary part f；

Solution update equation group (6), obtains voltage real part correction column vector δ e and voltage imaginary part correction column vector δ f.

Voltage correction formula is:

$\begin{array}{cc}{e}_i^{(t+1)}=e_i^{\left(t\right)}-{\mathrm{\δe}}_i^{\left(t\right)}& i=1,...,n-1\end{array}---\left(19\right)$

$\begin{array}{cc}{f}_i^{(t+1)}=f_i^{\left(t\right)}-{\mathrm{\δf}}_i^{\left(t\right)}& i=1,...,n-1\end{array}---\left(20\right)$

In formula, subscript (t) represents the value of the t time iteration；δe_iWith δ f_iIt is respectively the node voltage real part correction of node i Amount and node voltage imaginary part correction.

F, maximum amount of unbalance | the δ p | of judgement_max、|δq|_maxWith | δ u²|_maxWhether it is both less than convergence precision ε；If all Less than convergence precision ε, carry out step g, otherwise return to step c carries out next iteration；

The reactive power of g, the active power of calculated equilibrium node and reactive power and pv node, calculates each branch road wattful power Rate and reactive power, terminate.

Directly adopt the rectangular coordinate Newton Power Flow software for calculation calculating speed that above-mentioned principle is realized slower, commercially use Rectangular coordinate Newton Power Flow software for calculation adopt sparse matrix technology and node optimizing code technology, more complicated, uncomfortable Close scientific research personnel and carry out scientific research based on this further.Therefore, Chinese patent zl201010509556.5 proposes one Plant the Newton load flow calculation method that suitable research purpose uses, be to be carried out into one based on the calculating of polar coordinate Newton Power Flow The scientific research personnel of step research provides a Newton load flow calculation method being easy to change and safeguard, its feature is as follows:

(1) sparse matrix technology and node optimizing code are not adopted, greatly reduce algorithm realizes difficulty；

(2) judge to avoid unnecessary computing by simple logic, improve the calculating speed of Load flow calculation.

Chinese patent zl201010509556.5 proposed method is to be carried out based on the calculating of polar coordinate Newton Power Flow The scientific research personnel studying further provides a polar coordinate Newton load flow calculation method being easy to change and safeguard.The method When being realized using the compiled programming language such as c language speed quickly, but when being realized using matlab this kind of explanation type programming language Calculating speed is then very slow, and this patent does not make full use of the feature that matlab is good at matrix operationss and complex operation yet simultaneously.Cause This needs a feature making full use of matlab and the quick rectangular coordinate Newton load flow calculation method of calculating supplies The scientific research personnel carrying out scientific research on matlab platform uses.

Content of the invention

For solving the problems referred to above that prior art exists, the present invention will propose a kind of rectangular coordinate newton based on matlab Method tidal current computing method, can make full use of the distinctive feature being good at matrix operationss and complex operation of matlab, have simultaneously again The rectangular coordinate Newton load flow calculation method of very fast calculating speed.

To achieve these goals, technical scheme is as follows: a kind of rectangular coordinate Newton method based on matlab Tidal current computing method, using matrix operationss and complex operation.Comprise the following steps:

A, initial data input and voltage initialization；

B, formation bus admittance matrix；

C, formation Jacobian matrix and calculate node power；

Matlab is good at matrix operationss and complex operation, therefore using matlab programming it should derive based on matrix fortune Calculate the Jacobian matrix computational methods with complex operation.

Jacobian matrix element is relevant with node type, wants decision node type during conventional formation Jacobian matrix, according to Node type determines which node needs to form Jacobian matrix element.So process for the algorithm to be realized by circulation, It is easily processed, but do not fit through the method that matrix integral operation forms Jacobian matrix.Therefore, the present invention forms Jacobean matrix During battle array, not decision node type, Jacobian matrix element is all formed to all nodes, removes unwanted row and column afterwards again.

Derive below and calculate the formula of Jacobian matrix element and node power using matrix operationss.

Analysis to formula (7)～formula (10), can obtain

m_ij=n_ij(21)

l_ij=-h_ij(22)

Therefore, first seek h_ijAnd n_ij, obtain h_ijAnd n_ijAfterwards, m can naturally just be obtained_ijAnd l_ij.

Derive before calculating the formula of Jacobian matrix by matrix operationss, first look at Jacobian matrix each element how with again Number or phasor representation.

From formula (7) and formula (8), Jacobian matrix each element should be by y_ijWithConjugationProduct, that is, Obtain.Therefore h_ij、n_ij、m_ij、l_ijCan be generated by following formula:

$j_{ij}^0=y_{ij}{\hat{u}}_i---\left(23\right)$

In formula, subscript (^) represents the conjugation of plural number.

Formula (23) can regard the matrix j being obtained with following formula correspondence position element multiplication by bus admittance matrix as⁰The i-th row The element of jth row:

Formula (24) then can use a matrix fill-in function repmat of matlab to realize, the therefore sub- square of Jacobian matrix piecemeal Battle array can be obtained by following formula:

$j^0=repmat(\hat{u},1,n).*y---\left(25\right)$

In formula, j⁰For Jacobi's initial calculation matrix；Repmat is the matrix fill-in function of matlab, 1 expression in function During fillingRow be repeated 1 times, when in function, n represents fillingRow repeat n time；.* represent the element of the corresponding ranks of two matrixes It is multiplied.

By j⁰Obtaining initial Jacobian matrix piecemeal submatrix is:

h⁰=-re (j⁰) (26)

n⁰=im (j⁰) (27)

m⁰=im (j⁰) (28)

l⁰=re (j⁰) (29)

In formula, h⁰、n⁰、m⁰、l⁰Piecemeal submatrix for initial Jacobian matrix；Re represents the real part taking matrix element；im Represent the imaginary part taking matrix element.

The nondiagonal element of the initial Jacobian matrix piecemeal submatrix being obtained by formula (26)～formula (29) has been Jacobi Matrix element, diagonal element also needs to revise.

The the 1st and the 2nd on the right side of equation in formula (13)～formula (16) is exactly h⁰、n⁰、m⁰、l⁰Diagonal element, it is right therefore only to need The h obtaining⁰、n⁰、m⁰、l⁰Use a_iAnd b_iRevise, just can obtain Jacobian matrix piecemeal submatrix diagonal element.

The control variable of summation in node power equation (2) is changed into j, and is rearranged as:

$\left\{\begin{array}{cc}{p}_i=\underset{j=1}{\overset{n}{\σ}}(e_i{g}_{ij}+f_i{b}_{ij})e_j-\underset{j=1}{\overset{n}{\σ}}(e_i{b}_{ij}-f_i{g}_{ij})f_j& i=1,...,n\\ q_i=-\underset{j=1}{\overset{n}{\σ}}(e_i{b}_{ij}-f_i{g}_{ij})e_j-\underset{j=1}{\overset{n}{\σ}}(e_i{g}_{ij}+f_i{b}_{ij})f_j& i=1,...,n\end{array}\right.---\left(30\right)$

Formula (23) expands into:

$j_{ij}^0=(e_i{g}_{ij}+f_i{b}_{ij})+j(e_i{b}_{ij}-f_i{g}_{ij})---\left(31\right)$

Observation type (30) and formula (31), can obtain complex power:

$\tilde{s}=p+jq={\hat{j}}^0\hat{u}---\left(32\right)$

In formula,For complex power column vector.

Injection current phasor is:

$\stackrel{\·}{i}=y\stackrel{\·}{u}---\left(33\right)$

In formula,For injection current phasor column vector；For node voltage column vector.

Due to a_i、b_iIt is respectively the real part of injection current phasor and the imaginary part of node i, therefore

$\begin{array}{cc}{a}_i=\mathrm{re}\left({\stackrel{\·}{i}}_i\right)& i=1,...,n\end{array}---\left(34\right)$

$\begin{array}{cc}{b}_i=\mathrm{im}\left({\stackrel{\·}{i}}_i\right)& i=1,...,n\end{array}---\left(35\right)$

In formula,Injection current phasor for node i.

With injection current phasor, Jacobian matrix piecemeal submatrix diagonal element is modified as follows:

$\begin{array}{cc}{h}_{ii}=h_{ii}^0-\mathrm{re}\left(i_i\right)& i=1,...,n-1\end{array}---\left(36\right)$

$\begin{array}{cc}{n}_{ii}=n_{ii}^0-\mathrm{im}\left(i_i\right)& i=1,...,n-1\end{array}---\left(37\right)$

$\begin{array}{cc}{m}_{ii}=m_{ii}^0+\mathrm{im}\left(i_i\right)& i=1,...,m\end{array}---\left(38\right)$

$\begin{array}{cc}{l}_{ii}=l_{ii}^0-\mathrm{re}\left(i_i\right)& i=1,...,m\end{array}---\left(39\right)$

Voltage deviation is 0 to the nondiagonal element of voltage local derviation, and diagonal element is:

r_ii=-2e_iI=m+1 ..., n-1 (40)

k_ii=-2f_iI=m+1 ..., n-1 (41)

The element of r, k m+1 row to n-1 row can be put into the corresponding row of m, l, formation m ', l ':

m_i′_j=m_ijI=1 ..., m j=1 ..., n-1 (42)

l_i′_j=l_ijI=1 ..., m j=1 ..., n-1 (43)

m_i′_j=r_ijI=m+1 ..., n-1 j=1 ..., n-1 (44)

l_i′_j=k_ijI=m+1 ..., n-1 j=1 ..., n-1 (45)

Form Jacobian matrix and calculate node power, comprise the following steps:

C1, calculating Jacobi initial calculation matrix j⁰；

C2, calculate node complex power；

C3, by j⁰Calculate initial Jacobian matrix piecemeal submatrix；

C4, calculating injection current phasor；

C5, with injection current phasor, Jacobian matrix piecemeal submatrix diagonal element is modified；

C6, calculating r, k, form m ', l '；

C7, Jacobian matrix is formed by Jacobian matrix piecemeal submatrix；

$j=\left[\begin{array}{cc}h& n\\ m^{\′}& l^{\′}\end{array}\right]---\left(46\right)$

C8, Jacobian matrix is adjusted, removes balance nodes active power deviation and reactive power deviation is corresponding OK；Remove balance nodes voltage real part correction and the corresponding row of voltage imaginary part correction, terminate.

D, calculate node power deviation and node voltage deviation；

Formula (4) calculate node power deviation and voltage deviation equation are write as being formed as of matrix operationss:

$\left\{\begin{array}{c}\mathrm{\δ}p=p_s-\mathrm{re}\left(\tilde{s}\right)\\ \mathrm{\δ}q=q_s\mathrm{im}\left(\tilde{s}\right)\\ {\mathrm{\δu}}^2=u_s^2-|\stackrel{\·}{u}{|}^2\end{array}\right.---\left(47\right)$

In formula, δ p, δ q are respectively node active power deviation column vector and reactive power deviation column vector；δu²For section Point voltage amplitude deviation column vector；p_s、q_sBe respectively node give injection active power column vector and injection reactive power arrange to Amount；u_sFor node voltage set-point column vector；For node voltage column vector.

Remove pv node reactive power deviation and balance nodes in calculated node power bias vector δ p and δ q Active power deviation and reactive power deviation；Node voltage amplitude deviation column vector removes the node electricity of pq node and balance nodes Pressure amplitude value deviation.

E, solve equation and revise voltage real part e and imaginary part f；

Directly invoke solution system of linear equations algorithm solution update equation group (6) of matlab software, obtain node voltage real part Correction column vector δ e and node voltage imaginary part correction column vector δ f.

The formula (19) that voltage is modified and formula (20) are rewritten into matrix form and are:

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

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

In formula, subscript (t) represents the value of the t time iteration；δ e and δ f is respectively node voltage real part correction column vector With node voltage imaginary part correction column vector.

F, maximum amount of unbalance | the δ p | of judgement_max、|δq|_maxWith | δ u²|_maxWhether it is both less than convergence precision ε；If all Less than convergence precision ε, carry out step g, otherwise return to step c carries out next iteration.

The reactive power of g, the active power of calculated equilibrium node and reactive power and pv node, calculates each branch road wattful power Rate and reactive power, terminate.

Compared with prior art, the method have the advantages that

1st, method proposed by the present invention is realized in matlab platform, is easy to the various works that scientific research personnel uses matlab to provide Tool and function pair result of calculation are tested and are analyzed.

2nd, method proposed by the present invention adopts matrix operationss and complex operation, decreases program code, simplifies programming, makes Calling program becomes apparent from, and is easy to scientific research personnel's modification program, program is debugged and is improved, adds New function.

3rd, optimize because matlab carries out to matrix operationss, using matrix operationss than by the programming of matrix element loop computation Faster, directly invoke the equation solution algorithm of matlab simultaneously, also substantially increase calculating speed.It was verified that the present invention Method both facilitated scientific research personnel and program write, changed and is debugged, calculating speed is also substantially close in c language simultaneously The speed realized on speech platform, is that the research work of scientific research personnel provides an outstanding analytical tool.

Brief description

The present invention has 4, accompanying drawing.Wherein:

Fig. 1 is the flow chart that existing rectangular coordinate Newton Power Flow calculates.

Fig. 2 is the flow chart that existing rectangular coordinate Newton method forms Jacobian matrix.

Fig. 3 is the flow chart that rectangular coordinate Newton Power Flow of the present invention calculates.

Fig. 4 is the flow chart that the present invention forms Jacobian matrix and calculate node power.

Specific embodiment

Below in conjunction with the accompanying drawings the present invention is described further, amended to one according to flow process shown in Fig. 3-4 445 node real system examples are calculated.

The actual large-scale power system of 445 nodes has 445 nodes, and 544 branch roads, containing substantial amounts of small impedance branches.For Various methods are compared, these small impedance branches are changed to normal impedance branch road to meet the requirement of various methods.

Using the present invention and several control methods, 445 node real system examples are calculated, during calculating, restrained essence Spend for 0.00001.Several Load flow calculation algorithms are respectively as follows:

Method 1: Chinese patent zl201010509556.5 method is changed to Cartesian form, real using matlab language Existing.

Method 2: Chinese patent zl201010509556.5 method is changed to Cartesian form, real using matlab language Existing, but solve equation the equation solving algorithm directly invoking matlab.

Method 3: Chinese patent zl201010509556.5 method is changed to Cartesian form, real using matlab language Existing, but solve equation the equation solving algorithm directly invoking matlab, dimension is predefined to matrix variables simultaneously.

Method 4: the inventive method.

The calculating time of several method is shown in Table 1, the calculating time do not include data read in and output and branch power calculate Time.

Several rectangular coordinate Newton Power Flow of table 1 calculating calculating time compares

Tidal current computing method	Calculating time (s)
		Method 1	20.167
Method 2	4.730
		Method 3	1.027
Method 4	0.562

As seen from Table 1, Chinese patent zl201010509556.5 method is changed to directly adopt during Cartesian form Matlab realizes, and calculates the time very long；Patent zl201010509556.5 method is changed to Cartesian form and directly adopts When matlab realizes, if directly invoking the method that solves equation of matlab, calculating speed just can greatly improve, the solution side of matlab Cheng Fangfa comparative maturity is stable, is conducive to algorithm stability, matlab solves equation calling also very simply of method, simplifies Programming, makes program become apparent from；If dimension is predefined to matrix variables, it is to avoid program constantly expands in the process of implementation simultaneously The size of matrix, can greatly improve calculating speed.The result of calculation of the present invention shows to adopt matrix when forming Jacobian matrix Computing can greatly improve calculating speed further, and program also simplifies further.

The inventive method can be realized in the matlab programming language of any version, but suggestion is using more recent version Matlab language.

The present invention is not limited to the present embodiment, any equivalent concepts in the technical scope of present disclosure or change Become, be all classified as protection scope of the present invention.

Claims (1)

1. a kind of rectangular coordinate Newton load flow calculation method based on matlab, comprises the following steps:

A, initial data input and voltage initialization；

Raw data packets vinculum road and transformer branch data, node injection active power and reactive power, node voltage amplitude, Node reactive-load compensation data, and convergence precision and maximum iteration time；

Voltage initialization is started using flat, and that is, the node voltage real part of pv node and balance nodes draws definite value, the section of pq node Point voltage real part takes 1.0；The imaginary part of all node voltages all takes 0.0；Here unit adopts perunit value；

Described pq node is that active power and reactive power be known, voltage magnitude and the unknown node of voltage phase angle, described Pv node is that active power and voltage magnitude be known, reactive power and the unknown node of voltage phase angle, and described balance nodes are Voltage magnitude and voltage phase angle are it is known that the active power of node and the unknown node of reactive power；

B, formation bus admittance matrix；

Circuit according to input and transformer branch data form the bus admittance matrix as shown in formula (1)；

In formula, y_ikFor bus admittance matrix element, as subscript i ≠ k, it is node i and the transadmittance of node k, when subscript i=k When, it is the self-admittance of node i；N is nodes；

It is characterized in that:

C, formation Jacobian matrix and calculate node power；

C1, calculating Jacobi initial calculation matrix j⁰；

$j^0=repmat(\hat{u},1,n).*y---\left(2\right)$

In formula, j⁰For Jacobi's initial calculation matrix；Column vector for node voltage conjugation；Repmat is the matrix of matlab Stuffing function, in function during 1 expression fillingRow be repeated 1 times, when in function, n represents fillingRow repeat n time；.* represent The element multiplication of the corresponding ranks of two matrixes；

C2, calculate node complex power；

The node complex power column vector being made up of each node complex power is:

$\tilde{s}={\hat{j}}^0\hat{u}---\left(3\right)$

In formula,For node complex power column vector；Subscript (^) represents the conjugation of plural number；

C3, initial Jacobian matrix piecemeal submatrix is calculated by j0；

By j⁰Obtaining initial Jacobian matrix piecemeal submatrix is:

h⁰=-re (j⁰) (4)

n⁰=im (j⁰) (5)

m⁰=im (j⁰) (6)

l⁰=re (j⁰) (7)

In formula, h⁰、n⁰、m⁰、l⁰Piecemeal submatrix for initial Jacobian matrix；Re represents the real part taking matrix element；Im represents Take the imaginary part of matrix element；

C4, calculating injection current phasor；

$\stackrel{\·}{i}=y\stackrel{\·}{u}---\left(8\right)$

In formula,For injection current phasor column vector；

C5, with injection current phasor, Jacobian matrix piecemeal submatrix diagonal element is modified as follows:

$h_{ii}=h_{ii}^0-\mathrm{re}\left(i_i\right),i=1,...,n-1---\left(9\right)$

$n_{ii}=n_{ii}^0-\mathrm{im}\left(i_i\right),i=1,...,n-1---\left(10\right)$

$m_{ii}=m_{ii}^0+\mathrm{im}\left(i_i\right),i=1,...,m---\left(11\right)$

$l_{ii}=l_{ii}^0-\mathrm{re}\left(i_i\right),i=1,...,m---\left(12\right)$

In formula,Injection current phasor for node i；M is pq nodes；

C6, calculating r, k, form m ', l '；

Voltage deviation is 0 to the nondiagonal element of voltage local derviation, and diagonal element is:

r_ii=-2e_iI=m+1 ..., n-1 (13)

k_ii=-2f_iI=m+1 ..., n-1 (14)

Replace the corresponding row of m with the element of r m+1 row to n-1 row, form m '；Replace l with the element of k m+1 row to n-1 row Corresponding row, formed l '；

C7, Jacobian matrix is formed by Jacobian matrix piecemeal submatrix；

$j=\left[\begin{array}{cc}h& n\\ m^{\′}& l^{\′}\end{array}\right]---\left(15\right)$

In formula, j is Jacobian matrix, and h, n, m ', l ' are the piecemeal submatrix of Jacobian matrix；

C8, Jacobian matrix is adjusted, removes balance nodes active power deviation and the corresponding row of reactive power deviation；Go Fall balance nodes voltage real part correction and the corresponding row of voltage imaginary part correction, terminate；

D, calculate node power deviation and node voltage deviation；

It is calculated as follows node power deviation and node voltage deviation:

$\left\{\begin{array}{c}\mathrm{\δ}p=p_s-\mathrm{re}\left(\tilde{s}\right)\\ \mathrm{\δ}q=q_s-\mathrm{im}\left(\tilde{s}\right)\\ {\mathrm{\δu}}^2=u_s^2-|\stackrel{\·}{u}{|}^2\end{array}\right.---\left(16\right)$

In formula, δ p, δ q are respectively node active power deviation column vector and reactive power deviation column vector；δu²For node electricity Pressure amplitude value deviation column vector；p_s、q_sIt is respectively injection active power column vector and the injection reactive power column vector that node gives； u_sFor node voltage set-point column vector；For node voltage column vector；

Remove pv node reactive power deviation in calculated node power bias vector δ p and δ q and balance nodes are active Power deviation and reactive power deviation；Node voltage amplitude deviation column vector removes the node voltage width of pq node and balance nodes Value deviation；

E, solve equation and revise voltage real part e and imaginary part f；

Jacobian matrix j is obtained by step c and step d obtains node power bias vector δ p and δ q and node voltage amplitude is inclined Difference vector δ u², the update equation being configured to Load flow calculation is as follows:

$\left[\begin{array}{c}\mathrm{\δ}p\\ \mathrm{\δ}q\\ \mathrm{\δ}{u}^2\end{array}\right]=j\left[\begin{array}{c}\mathrm{\δ}e\\ \mathrm{\δ}f\end{array}\right]---\left(17\right)$

Directly invoke solution system of linear equations algorithm solution update equation group (17) of matlab software, obtain voltage real part correction to Amount δ e and voltage imaginary part correction vector δ f；

As the following formula node voltage real part and imaginary part are modified:

e^(t+1)=e^(t)-δe^(t)(18)

f^(t+1)=f^(t)-δf^(t)(19)

In formula, subscript (t) represents the value of the t time iteration；

F, maximum amount of unbalance | the δ p | of judgement_max、|δq|_maxWith | δ u²|_maxWhether it is both less than convergence precision ε；If both less than Convergence precision ε, carries out step g, and otherwise return to step c carries out next iteration；

The reactive power of g, the active power of calculated equilibrium node and reactive power and pv node, calculate each branch road active power and Reactive power, terminates.