CN106410811A - Trend calculating method of firstly iterating small impedance branch end point to change Jacobian matrix - Google Patents

Abstract

The invention discloses a trend calculating method of firstly iterating a small impedance branch end point to change a Jacobian matrix. During first iteration, a small impedance branch end point adopts ai and bi values calculated through given values Pis and Qis to calculate a Jacobian matrix element. A normal branch end point during the first iteration and all the nodes during subsequent iterations adopt a traditional method to calculate the Jacobian matrix element. In the invention, during a first iteration process, two end nodes of a small impedance branch adopt a Jacobian matrix calculating method which is different from a method used during later iteration processes so that a convergence problem of rectangular coordinate Newton method trend calculating during small impedance branch electrical power system analysis is solved. When conventional rectangular coordinate Newton method trend calculating is not converged, by using the method, reliable convergence can be performed and iteration frequencies are less than that in an existing patent technology. In the invention, trend calculating can be performed on a normal electric power system and there is no harmful effects.

Description

Iteration small impedance branches end points changes the tidal current computing method of Jacobian matrix first

Technical field

The present invention relates to the tidal current computing method containing small impedance branches power system, particularly a kind of right angle of power system 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 given according to power system Fixed service condition and network structure determine the running status of whole power system.Load flow calculation is also other analyses of power system Basis, such as safety analysis, transient stability analysis etc. will use Load flow calculation.Due to there is convergence reliability, calculating speed relatively Advantage fast and that memory requirements is moderate, Newton method becomes the main stream approach of current Load flow calculation.Newton method is divided into polar coordinate and straight Two kinds of forms of angular coordinate, the Newton Power Flow of two kinds of forms calculates and is all widely used in power system.

In rectangular coordinate Newton Power Flow calculates, the voltage of node i is expressed as using rectangular coordinate：

To normal electric power networks, rectangular coordinate Newton Power Flow calculates has good convergence, but runs into containing little resistance During the Ill-conditioned network of anti-branch road, rectangular coordinate Newton Power Flow calculates and may dissipate.Power system small impedance branches can be divided into Little impedance line and little impedance transformer branch road, on mathematical model, circuit is considered as no-load voltage ratio is 1:1 transformator, therefore under Analyze only during surface analysis taking little impedance transformer branch road as a example.Little impedance transformer model is shown in Fig. 1, the non-standard no-load voltage ratio of transformator K is located at node i side, and impedance is located at standard no-load voltage ratio side.Transformer impedance z_ij=r_ij+jx_ijVery little, admittance is

In formula, y_ij、g_ij、b_ijIt is respectively the admittance of small impedance branches, conductance and susceptance between node i and node j；r_ij、 x_ijIt is respectively the resistance of small impedance branches and reactance between node i and node j.

Due to small impedance branches l_i-jImpedance very little, the voltage drop also very little of branch road, the therefore electricity of transformator two end node Pressure should meet：

As shown in Fig. 2 existing rectangular coordinate Newton load flow calculation method, mainly include the following steps that：

A, input initial data and initialization voltage

According to the feature of power system node, Load flow calculation is divided into 3 classes power system node：Node active power and nothing 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.

Voltage initialization is started using flat, and that is, the voltage real part of PV node and balance nodes draws definite value, the electricity of PQ node Compacting portion takes 1.0；The imaginary part of all voltages all takes 0.0.Here unit adopts perunit value.

B, formation bus admittance matrix

If node i and the original self-conductance of node j be respectively G from susceptance_i0、B_i0、G_j0、B_j0, increase by one between them Self-admittance after bar small impedance branches and transadmittance are respectively：

In formula, Y_ii、Y_jjIt is respectively node i and the self-admittance of node j；Y_ijFor the transadmittance between node i and node j； r_ij、x_ijIt is respectively the resistance of small impedance branches and reactance between node i and node j；K is little impedance between node i and node j The no-load voltage ratio (if power transmission line branch road, no-load voltage ratio is 1) of branch road；

C, calculating power and voltage deviation

The power deviation computing formula of PQ node is：

In formula, P_is、Q_isIt is respectively injection active power and the reactive power that node i gives, P_isFor power supply active power with The difference of load active power, Q_isDifference for power supply reactive power and reactive load power；a_i、b_iIt is respectively the calculating note of node i Enter real part and the imaginary part of electric current phasor, be

In formula, n is the nodes of power system.

During Load flow calculation convergence, Δ P in formula (6)_i、ΔQ_iAll level off to 0, therefore a_iAnd b_iIt is equal to by set-point P_isAnd Q_is The a calculating_isAnd b_is

The active power of PV node and voltage deviation computing formula are：

In formula, V_isThe voltage magnitude giving for node i.

Balance nodes are not involved in iterative calculation it is not necessary to calculate power deviation or voltage deviation.

Seek the value of maximum absolute value in each node power or voltage deviation, referred to as maximum amount of unbalance, if maximum uneven The absolute value weighed is less than given convergence precision, goes to step F, otherwise execution step D.

D, formation Jacobian matrix J

Element (during i ≠ j) computing formula of Jacobian matrix J is as follows：

Element (during i=j) computing formula of Jacobian matrix J is as follows：

PQ section knock type (16)-(19) calculate Jacobian matrix element；PV section knock type (16), (17), (20), (21) meter Calculate Jacobian matrix element；Balance nodes do not calculate Jacobian matrix element.

E, solution update equation and correction voltage real part e, imaginary part f

The fundamental equation (6) of Load flow calculation and (9) are Nonlinear System of Equations, are generally asked using successive Linearization Method iteration Solution.The equation that linearisation obtains is referred to as update equation, for seeking the correction of voltage real part and imaginary part.

Update equation is：

In formula, J is Jacobian matrix；Δ P and Δ Q is respectively active power and reactive power deviation column vector；ΔV²For Voltage magnitude deviation column vector；Δ e and Δ f is respectively real part and the imaginary part correction column vector of voltage phasor；For The local derviation matrix to voltage phasor real part column vector transposition for the active power departure function column vector, subscript T is transposition symbol.

Voltage correction formula is：

In formula, subscript t represents the t time iteration.

F, output node and branch data.

To normal electric power networks, Newton Power Flow calculates has good convergence, but runs into containing small impedance branches During Ill-conditioned network, Newton Power Flow calculates and may dissipate.And small impedance branches generally existing in power system, convergence is electricity The most important index of this kind of nonlinear problem of Force system Load flow calculation, calculating does not restrain and just cannot obtain non trivial solution.Therefore change Kind rectangular coordinate Newton Power Flow calculates and has very important significance for the convergence containing small impedance branches power system.

Chinese patent ZL201410299531.5 discloses a kind of modification conventional Cartesian coordinate Newton Power Flow that passes through and calculates The method of Jacobian matrix, improves the convergence of Load flow calculation.The method calculates during Jacobi's element using by set-point P_is And Q_isThe a calculating_iAnd b_iValue, efficiently solves the dissipating of small impedance branches electric power system tide calculating being 0 containing resistance and asks Topic.But when the resistance of small impedance branches is not 0, the method iterationses increase, and convergence is deteriorated, and does not even restrain.

Chinese patent ZL201410315785.1 proposes the rectangular coordinate Newton Power Flow that a kind of Jacobian matrix changes Computational methods, iteration and subsequently each iteration adopt different Jacobian matrix computational methods, iteration meter first to the method first Calculate during Jacobi's element using by set-point P_isAnd Q_isThe a calculating_iAnd b_iValue, during subsequently each iterative calculation Jacobi's element still Using traditional method, efficiently solve the divergence problem calculating containing the small impedance branches electric power system tide that resistance is not 0, but should When power system comprises the small impedance branches that a plurality of resistance is not 0, iterationses increase method, and convergence is deteriorated.

Chinese patent ZL201611094297.8 proposes a kind of cattle changing Jacobian matrix with iteration and node type Method tidal current computing method, the method first iteration when all PQ nodes with subsequently each iteration using different Jacobian matrixes Computational methods, during iteration, all PQ nodes calculate during Jacobi's element using by set-point P first_isAnd Q_isThe a calculating_iAnd b_i Value, first during iteration all PV nodes and during subsequently each iteration all nodes calculate and during Jacobi's element, still adopt tradition side Method, efficiently solves the divergence problem that power system comprises the Load flow calculation of the small impedance branches that a plurality of resistance is not 0, but iteration Number of times is still more, requires further improvement.

Content of the invention

For solving the problems referred to above that prior art exists, the present invention will propose one kind, and iteration small impedance branches end points changes first Become the tidal current computing method of Jacobian matrix, the method can improve the small impedance branches electric power that its analysis is not 0 containing resistance The convergence rate of system.

The method of traditional calculating Jacobian matrix is to derive from the ultimate principle of Newton method, normal impedance branch road End points using traditional method calculate Jacobian matrix element be suitable, but the end points of small impedance branches adopt traditional method meter Calculating Jacobian matrix element then can lead to Load flow calculation to dissipate.First during iteration, voltage is at the beginning of the voltage of flat startup method setting Value, the branch power that normal impedance branch road calculates is more or less the same with actual value, the injecting power of the node being connected with these branch roads Value of calculation is close with set-point, therefore first the end points of iteration normal impedance branch road calculated using injecting power value of calculation refined can Than matrix element also relatively rationally.Small impedance branches due to its impedance very little, the inconsistent band of both end voltage initial value and actual value The voltage difference of the very little come will calculate very big branch power, the injecting power value of calculation of the node being connected with this branch road Very big, Load flow calculation can be led to dissipate, therefore the end points of iteration small impedance branches should not use injecting power value of calculation first Calculate Jacobian matrix element, and Jacobian matrix element should be calculated using injecting power set-point or initial value.

To achieve these goals, the present invention proposes a kind of rectangular coordinate Newton load flow calculation method to improve trend Computational convergence.During the iteration first of the present invention, small impedance branches end points is using by set-point P_isAnd Q_isThe a calculating_iAnd b_iValue Calculate Jacobian matrix element, first during iteration normal leg endpoint and during subsequently each iteration all nodes then adopt tradition Method calculates Jacobian matrix element.

Technical scheme is as follows：Iteration small impedance branches end points changes the Load flow calculation side of Jacobian matrix first Method, comprises the following steps：

A, input initial data and initialization voltage；

B, determine connected branch type T of two end nodes according to the size of branch resistance and reactance

Form comprising the following steps that of node institute's chord road type array：

B1, reading branch data, arrange small resistor threshold value r_minWith low reactance threshold value x_min；

B2, node institute's chord road type array T reset；

B3, make m=1；

B4, first and last node number i and j taking branch road m, resistance r, reactance x；

B5, judge whether to meet r≤r_minAnd x≤x_minCondition, if be unsatisfactory for, go to step B7；

B6, make T_i=1, T_j=1；

B7, make m=m+1；

B8, judge that whether m is more than circuitry number l, if m is not more than l and goes to step B4；Otherwise go to step C；

C, formation bus admittance matrix；

D, setting iteration count t=0；

E, calculating power and voltage deviation, seek maximum amount of unbalance Δ W_max；

F, the maximum amount of unbalance absolute value of judgement | Δ W_max| whether it is less than convergence precision ε；If less than convergence precision ε, hold Row step J；Otherwise, execution step G；

G, formation Jacobian matrix；

In addition to iteration first, Jacobian matrix computational methods still adopt traditional method.The Jacobian matrix meter of iteration first Calculation method adopts distinct methods according to the type that node connects branch road.For the end points of small impedance branches, because adopting tradition side Method calculates Jacobian matrix and Load flow calculation can be led to dissipate, so calculate calculating injection using formula (8) during Jacobian matrix element The real part of electric current phasor and imaginary part effect are preferable；For the end points of normal impedance branch road, still traditionally calculate Jacobi Matrix element, i.e. the real part a of the node i injection current phasor in Jacobian matrix computing formula_iWith imaginary part b_iCalculate by formula (7).

Form comprising the following steps that of Jacobian matrix element：

G1, by formula (10)-(15) calculate i ≠ j when Jacobian matrix element；

G2, make i=1；

G3, judge whether to meet t=0 and T simultaneously_i=1 condition, if be unsatisfactory for this condition to go to step G4；If full Foot, then press the real part a of the injection current phasor of formula (8) calculate node i_iWith imaginary part b_i, then go to step G5；

G4, by formula (7) calculate node i injection current phasor real part a_iWith imaginary part b_i；

G5, by formula (16)-(21) calculate i=j when Jacobian matrix element；

G6, make i=i+1；

G7, judge that whether i is more than nodes n, if i is not more than n and goes to step G3；Otherwise go to step H；

H, solution update equation and correction voltage real part e, imaginary part f；

I, make t=t+1, return to step E carries out next iteration；

J, output node and branch data.

Compared with prior art, the invention has the advantages that：

1st, the present invention is adopted and later each iterative process by small impedance branches two ends node in iterative process first Different Jacobian matrix computational methods, solve rectangular coordinate Newton Power Flow and calculate and contain small impedance branches electric power in analysis Convergence problem during system.Using conventional Cartesian coordinate Newton Power Flow calculate do not restrain when, this method can reliable conveyance, And it is fewer than existing patented technology iterationses.

2nd, because the present invention can not only prop up containing little impedance in the calculating analysis of effectively solving conventional Cartesian coordinate Newton Power Flow The convergence problem of road power system, also can carry out Load flow calculation to normal power system simultaneously, not have harmful effect.

Brief description

The present invention has 6, accompanying drawing.Wherein：

Fig. 1 is power system little impedance transformer model schematic.

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

Fig. 3 is the flow chart that patented method 1 rectangular coordinate Newton Power Flow calculates.

Fig. 4 is the flow chart that patented method 2 rectangular coordinate Newton Power Flow calculates.

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

Fig. 6 is the flow chart that the present invention forms node institute's chord road type array.

Specific embodiment

Below in conjunction with the accompanying drawings the present invention is described further.Little impedance transformer model according to Fig. 1, adopts The flow chart being calculated with the rectangular coordinate Newton Power Flow shown in Fig. 5-6, has carried out Load flow calculation to an actual large-scale power grid. This actual large-scale power grid has 445 nodes, containing substantial amounts of small impedance branches.Wherein, the small impedance branches of x≤0.001 have 49 Bar, the small impedance branches of x≤0.0001 have 41, and the small impedance branches of x≤0.00001 have 22.Wherein resistance value is minimum It is the small impedance branches l between node 118 and node 125_118-125For x=0.00000001, no-load voltage ratio k=0.9565, k is located at section Point 118 sides.The convergence precision of Load flow calculation is 0.00001.In order to verify that present invention calculating is not 0 small impedance branches containing resistance The convergence of power system, small impedance branches l_118-125、l_60-122And l_287-310Resistance be changed to r=0.0001.

As a comparison, using following 3 kinds of control methods, Load flow calculation has been carried out to this actual large-scale power grid simultaneously：

Conventional method：Conventional rectangular coordinate Newton Power Flow method；

Patented method 1：The patented method of Patent No. ZL201410315785.1；

Patented method 2：The patented method of Application No. ZL201611094297.8.

Iterationses the results are shown in Table 1.

The iteration result of the different trend method of table 1

Method	Conventional method	Patented method 1	Patented method 2	The inventive method
					Iteration result	Do not restrain	7 convergences	6 convergences	5 convergences

From table 1, for amended 445 node practical power systems examples, conventional Cartesian coordinate Newton Power Flow Method does not restrain, and the inventive method and existing patented method can restrain, but the iterationses of the inventive method are more special than existing Few 2 times of sharp method 1, fewer than existing patented method 21 time.

Different each iteration maximum amount of unbalances of tidal current computing method are shown in Table 2.Unit is perunit value.

The different each iteration maximum amount of unbalance of trend method of table 2

Iteration sequence number	Conventional method	Patented method 1	Patented method 2	The inventive method
					0	-4754.570367135	-4754.570367135	-4754.570367135	-4754.570367135
1	-3451593.823720038	-11.138394991	-3.264368583	23.913925681
					2	-886651.468310079	-6.163450054	-0.715148045	3.050019341
3	-222023.112200678	-1.441071252	-0.076847277	0.102604201
					4	-55754.415245002	-0.106199006	-0.002294590	-0.000454516
5	-13972.568194423	-0.006353455	-0.000017499	-0.000000012
					6	-6386.835620506	-0.000141863	-0.000000001
7	-6585.38761914	-0.000000062
					8	-378994.776907351
9	-98508.025841226
					10	-37917.863557986

As shown in Table 2,4 kinds of methods first before iteration maximum amount of unbalance identical and very big.First after iteration, existing patent Method and this patent method maximum amount of unbalance significantly reduce, existing 7 convergences of patented method 1 iteration；Existing patented method 2 changes 6 convergences of generation；This patent method maximum amount of unbalance reduces speed faster, 5 convergences of iteration；And conventional method is maximum uneven Weigh and then become big, finally dissipate.

The end points of small impedance branches is during each node power reactive power input value and the Load flow calculation convergence of PV node Value of calculation and initial calculation value are shown in Table 3.Unit is perunit value.

Value of calculation and initial calculation value when the power supply reactive power input value of table 3PV node and convergence

Node	Input value	Convergence value of calculation	Initial calculation value
				22	1.80000	1.36829	0.15000
400	0.10000	0.69586	204.87496
				439	0.80000	0.48861	525.04000
440	0.80000	0.48861	525.04000

From table 3, it is each node of PV node, power supply reactive power input value and trend for small impedance branches end points Calculating value of calculation during convergence has larger difference, but the difference of value of calculation when initial calculation value is restrained with Load flow calculation is more Greatly.Therefore first iteration when, the end points input value of small impedance branches calculates the real part of injection current phasor and imaginary part will more be closed Reason.

The present invention can be realized using any programming language and programmed environment, such as C language, C++, FORTRAN, Delphi etc..Development environment can adopt Visual C++, Borland C++Builder, Visual FORTRAN etc..

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. first iteration small impedance branches end points change Jacobian matrix tidal current computing method it is characterised in that：Including following Step：

A, input initial data and initialization voltage

According to the feature of power system node, Load flow calculation is divided into 3 classes power system node：Node active power and idle work( Rate is known, node voltage amplitude and the unknown node of voltage phase angle are referred to as PQ node；Known to node active power and voltage magnitude, 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 node has Work(power and the unknown node of reactive power are referred to as balance nodes；

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

B, determine connected branch type T of two end nodes according to the size of branch resistance and reactance

Form comprising the following steps that of node institute's chord road type array：

B1, reading branch data, arrange small resistor threshold value r_minWith low reactance threshold value x_min；

B2, node institute's chord road type array T reset；

B3, make m=1；

B4, first and last node number i and j taking branch road m, resistance r, reactance x；

B5, judge whether to meet r≤r_minAnd x≤x_minCondition, if be unsatisfactory for, go to step B7；

B6, make T_i=1, T_j=1；

B7, make m=m+1；

B8, judge that whether m is more than circuitry number l, if m is not more than l and goes to step B4；Otherwise go to step C；

C, formation bus admittance matrix

If node i and the original self-conductance of node j be respectively G from susceptance_i0、B_i0、G_j0、B_j0, increase by is little between them Self-admittance after impedance branch and transadmittance are respectively：

$Y_{ii}=\left(G_{i0}+\frac{r_{ij}}{k^2(r_{ij}^2+x_{ij}^2)}\right)+j\left(B_{i0}-\frac{x_{ij}}{k^2(r_{ij}^2+x_{ij}^2)}\right)---\left(1\right)$

$Y_{jj}=\left(G_{j0}+\frac{r_{ij}}{(r_{ij}^2+x_{ij}^2)}\right)+j\left(B_{j0}-\frac{x_{ij}}{(r_{ij}^2+x_{ij}^2)}\right)---\left(2\right)$

$Y_{ij}=-\frac{r_{ij}}{k(r_{ij}^2+x_{ij}^2)}+j\frac{x_{ij}}{k(r_{ij}^2+x_{ij}^2)}---\left(3\right)$

In formula, Y_ii、Y_jjIt is respectively node i and the self-admittance of node j；Y_ijFor the transadmittance between node i and node j；r_ij、x_ij It is respectively the resistance of small impedance branches and reactance between node i and node j；K is small impedance branches between node i and node j No-load voltage ratio, if power transmission line branch road, no-load voltage ratio k is 1；

D, setting iteration count t=0；

E, calculating power and voltage deviation, seek maximum amount of unbalance Δ W_max；

The power deviation computing formula of PQ node is：

$\left\{\begin{array}{c}{\mathrm{\ΔP}}_i=P_{is}-P_i=P_{is}-e_i{a}_i-f_i{b}_i\\ {\mathrm{\ΔQ}}_i=Q_{is}-Q_i=Q_{is}-f_i{a}_i+e_i{b}_i\end{array}\right.---\left(4\right)$

In formula, P_is、Q_isIt is respectively injection active power and the reactive power that node i gives, P_isFor power supply active power and load The difference of active power, Q_isDifference for power supply reactive power and reactive load power；e_iAnd f_iIt is respectively the voltage phasor of node i Real part and imaginary part；a_i、b_iIt is respectively the real part of calculating injection current phasor and the imaginary part of node i, be

$\left\{\begin{array}{c}{a}_i=\underset{m=1}{\overset{n}{\Σ}}(G_{im}{e}_m-B_{im}{f}_m)\\ b_i=\underset{m=1}{\overset{n}{\Σ}}(G_{im}{f}_m+B_{im}{e}_m)\end{array}\right.---\left(5\right)$

In formula, n is the nodes of power system；

During Load flow calculation convergence, Δ P in formula (4)_i、ΔQ_iAll level off to 0, therefore a_iAnd b_iIt is equal to by set-point P_isAnd Q_isCalculate The a going out_isAnd b_is

$\left\{\begin{array}{c}{a}_i\≈a_{is}=\frac{e_i{P}_{is}+f_i{Q}_{is}}{e_i^2+f_i^2}\\ b_i\≈b_{is}=\frac{f_i{P}_{is}-e_i{Q}_{is}}{e_i^2+f_i^2}\end{array}\right.---\left(6\right)$

The active power of PV node and voltage deviation computing formula are：

$\left\{\begin{array}{c}{\mathrm{\ΔP}}_i=P_{is}-P_i=P_{is}-e_i{a}_i-f_i{b}_i\\ {\mathrm{\ΔV}}_i^2=V_{is}^2-(e_i^2+f_i^2)\end{array}\right.---\left(7\right)$

In formula, V_isThe voltage magnitude giving for node i；

Balance nodes are not involved in iterative calculation it is not necessary to calculate power deviation or voltage deviation；

Seek the value of maximum absolute value in each node power or voltage deviation, referred to as maximum amount of unbalance Δ W_max；

F, the maximum amount of unbalance absolute value of judgement | Δ W_max| whether it is less than convergence precision ε；If less than convergence precision ε, execute step Rapid J；Otherwise, execution step G；

G, formation Jacobian matrix J

In addition to iteration first, Jacobian matrix computational methods still adopt traditional method；The Jacobian matrix calculating side of iteration first Method adopts distinct methods according to the type that node connects branch road；For the end points of small impedance branches, because adopting traditional method meter Calculating Jacobian matrix can lead to Load flow calculation to dissipate, so calculate calculating injection current using formula (6) during Jacobian matrix element The real part of phasor and imaginary part effect are preferable；For the end points of normal impedance branch road, still traditionally calculate Jacobian matrix Element, i.e. the real part a of the node i injection current phasor in Jacobian matrix computing formula_iWith imaginary part b_iCalculate by formula (5)；

Form comprising the following steps that of Jacobian matrix element：

Jacobian matrix element when G1, calculating i ≠ j；

As i ≠ j, the element computing formula of Jacobian matrix J is as follows：

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂e_j}=-G_{ij}{e}_i-B_{ij}{f}_i---\left(8\right)$

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂f_j}=B_{ij}{e}_i-G_{ij}{f}_i---\left(9\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂e_j}=B_{ij}{e}_i-G_{ij}{f}_i---\left(10\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂f_j}=G_{ij}{e}_i+B_{ij}{f}_i---\left(11\right)$

$\frac{\∂{\mathrm{\ΔV}}_i^2}{\∂e_j}=0---\left(12\right)$

$\frac{\∂{\mathrm{\ΔV}}_i^2}{\∂f_j}=0---\left(13\right)$

G2, make i=1；

G3, judge whether to meet t=0 and T simultaneously_i=1 condition, if be unsatisfactory for going to step G4；If it is satisfied, then pressing formula (6) The real part a of the injection current phasor of calculate node i_iWith imaginary part b_i, then go to step G5；

G4, by formula (5) calculate node i injection current phasor real part a_iWith imaginary part b_i；

Jacobian matrix element when G5, calculating i=j；

As i=j, the element computing formula of Jacobian matrix J is as follows：

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂e_i}=-a_i-G_{ii}{e}_i-B_{ii}{f}_i---\left(14\right)$

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂f_i}=-b_i+B_{ii}{e}_i-G_{ii}{f}_i---\left(15\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂e_i}=b_i+B_{ii}{e}_i-G_{ii}{f}_i---\left(16\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂f_i}=-a_i+G_{ii}{e}_i+B_{ii}{f}_i---\left(17\right)$

$\frac{\∂{\mathrm{\ΔV}}_i^2}{\∂e_i}=-2e_i---\left(18\right)$

$\frac{\∂{\mathrm{\ΔV}}_i^2}{\∂f_i}=-2f_i---\left(19\right)$

PQ section knock type (14)-(17) calculate Jacobian matrix element；PV section knock type (14), (15), (18), (19) calculating are refined Than matrix element；Balance nodes do not calculate Jacobian matrix element；

G6, make i=i+1；

G7, judge that whether i is more than nodes n, if i is not more than n and goes to step G3；Otherwise go to step H；

H, solution update equation and correction voltage real part e, imaginary part f

The fundamental equation (4) of Load flow calculation and (7) are Nonlinear System of Equations, using successive Linearization Method iterative；Linearly Change the equation obtaining and be referred to as update equation, for seeking the correction of voltage real part and imaginary part；

Update equation is：

$\left[\begin{array}{c}\mathrm{\Δ}P\\ \mathrm{\Δ}Q\\ {\mathrm{\ΔV}}^2\end{array}\right]=J\left[\begin{array}{c}\mathrm{\Δ}e\\ \mathrm{\Δ}f\end{array}\right]=\left[\begin{array}{cc}\frac{\∂\mathrm{\Δ}P}{\∂e^T}& \frac{\∂\mathrm{\Δ}P}{\∂f^T}\\ \frac{\∂\mathrm{\Δ}Q}{\∂e^T}& \frac{\∂\mathrm{\Δ}Q}{\∂f^T}\\ \frac{\∂{\mathrm{\ΔV}}^2}{\∂e^T}& \frac{\∂{\mathrm{\ΔV}}^2}{\∂f^T}\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}e\\ \mathrm{\Δ}f\end{array}\right]---\left(20\right)$

In formula, J is Jacobian matrix；Δ P and Δ Q is respectively active power and reactive power deviation column vector；ΔV²For voltage amplitude Value deviation column vector；Δ e and Δ f is respectively real part and the imaginary part correction column vector of voltage phasor；For wattful power The local derviation matrix to voltage phasor real part column vector transposition for the rate departure function column vector, subscript T is transposition symbol；

Voltage correction formula is：

$\left\{\begin{array}{c}{e}_i^{(t+1)}=e_i^{\left(t\right)}-{\mathrm{\Δe}}_i^{\left(t\right)}\\ f_i^{(t+1)}=f_i^{\left(t\right)}-{\mathrm{\Δf}}_i^{\left(t\right)}\end{array}\right.---\left(21\right)$

In formula, subscript t represents the t time iteration；

I, make t=t+1, return to step E carries out next iteration；

J, output node and branch data.