CN103353921A - Method for determining power system parallel load flow based on heuristic automatic network partitioning - Google Patents

Abstract

The invention relates to a method for determining a power system parallel load flow based on heuristic automatic network partitioning, which comprises the following steps: partitioning distribution network; generating the admittance matrixes of the partitioning nodes, the power flow calculation constant-jacobian matrix and the voltage-power sensitivity matrixes of virtual generators in partitions, sending the voltage-power sensitivity matrixes to a coordination layer; sending initial voltage magnitude and phase angle of virtual engines at the split points by the coordination layer; carrying out power flow calculation on the partitions, sending output power of the virtual engines to a coordination server; judging whether the boundary coordination quantity precision is met or not according to the mismatch quantity of the sum of output power of two virtual engines at the split points by the coordination layer, if no, correcting the voltage of the partitioning points, sending the corrected result to sub-areas, skipping to the last step, and otherwise, skipping to the next step; outputting the data result. The determining method has good convergence, high calculation speed, less memory occupied rate and calculation convenience and flexibility.

Description

A kind ofly determine method based on heuristic automatic network subregion electric system parallel flow

Technical field

The present invention relates to electric power system tide and determine method, be specifically related to a kind of electric system parallel flow and determine method.

Background technology

The effect that electric power system tide calculates is to the running status of given service condition and network structure calculating whole system, belongs to and finds the solution large-scale nonlinear system of equations category.Classical method for solving comprises Newton method, PQ decomposition method, the front Dai Fa etc. that pushes back.The major requirement of trend computational algorithm has convenience, dirigibility of convergence, computing velocity, memory usage and calculating etc.

On the one hand, along with the characteristics of modern power systems large system strong nonlinearity and multicomponent become increasingly conspicuous, its calculated amount and computation complexity sharply increase, on the other hand, along with parallel machine and parallel computing development and maturation, so that parallel computation becomes the effective way that solves electric power system tide.Progressively ripe and use along with the trunking computer expanded of high performance price ratio and distributed parallel routine library can realize that the distributed parallel of electric power system tide finds the solution by computer cluster.

For the continuous expansion of electric system scale with to on-line analysis and improving constantly that control requires, adopt the relatively large electric power system tide of centralized algorithm analytical calculation and Optimal Power Flow problem often to exist calculator memory not enough, speed of convergence waits the dimension disaster problem slowly.

Summary of the invention

For the deficiencies in the prior art, the purpose of this invention is to provide and a kind ofly determine method based on heuristic automatic network subregion electric system parallel flow, the present invention is directed to the parallel actual demand of finding the solution of electric system problem, at first realize the auto-partition of electric power networks by a kind of heuristic rule, secondly realize zone Divided Parallel Calculation by interregional decoupling method, realize global convergence by the iteration convergence condition at last, be fit to trend electric system, that have stronger concurrency and the parallel method for solving of Optimal Power Flow thereby develop a cover.

The objective of the invention is to adopt following technical proposals to realize:

The invention provides and a kind ofly determine method based on heuristic automatic network subregion electric system parallel flow, its improvements are that described method comprises the steps:

A, distribution network is carried out subregion;

The voltage of B, the admittance matrix that generates the subregion node, normal Jacobian matrix that trend is calculated and the virtual generator of subregion-power sensitivity matrix, and voltage-power sensitivity matrix is sent to cooperation layer;

C, sent initial voltage amplitude and the phase angle of each virtual generator in split point place by cooperation layer;

D, each subregion is carried out trend calculate, and the output power of each virtual generator is sent to coordination server;

E, cooperation layer are according to the amount of mismatch of each split point place two virtual generated output power sum, judge that whether the border coordination accuracy of measurement satisfies, and does not then revise split point voltage if do not satisfy, and correction result is sent to all subregion, jump to step D, otherwise jump to step F;

F, parallel load flow algorithm finish, the output data result.

Wherein, in the described steps A, the power distribution network subregion is based on heuristic automatic network subregion, comprises the steps:

1) network is carried out topological analysis, with the network chord virtual for off-state so that network becomes radial pattern, Radial network is numbered and degree of depth traversal, and the numbering of each frontier node is stored in the array;

2) definition comprises frontier node in the D of all interior nodes value, makes D _iBe to comprise not total node number of subregion node of node i self and downstream thereof, each node D value of initialization is-1;

3) choose in the frontier node array arbitrary frontier node and recall as present node, setting frontier node D value is 1, simultaneously this frontier node of rejecting in the frontier node array;

4) recall from present node, judge each child node D of its father node _jWhether all definitely be worth: if, find the solution the D value of father node according to following expression (13), and with this father node as present node; If not, return step 3);

$D_k=1+\underset{j\∈k}{\mathrm{\Σ}}{D}_j---\left(13\right);$

5) the D value of judging frontier node whether satisfy D 〉=[N/m] and the present node child node whether unique: if, carry out step 6), if not, then select D value in present node and the child node thereof near the node of [N/m] as present node, and carry out step 6); Wherein N represents nodes, and m is the number of partitions;

The branch road that 6) will connect present node and father node thereof is as the cutting branch road, and search present node downstream D value is not all nodes of 0 or-1, is combined as new subregion with present node;

7) present node D value is modified to 0, and from then on node continues to recall as present node; Guarantee that all frontier nodes all recall, Radial network divides end of extent;

8) output subregion result.

Wherein, among the described step B, the power flow equation behind the subregion is as follows:

Wherein: formula (1) is the a-quadrant power flow equation, and first and second is respectively a-quadrant voltage status amount and each node injecting power in the variable, and third and fourth is respectively voltage and the injecting power amount of virtual generator that the a-quadrant connects; Formula (2) is B zone power flow equation, and first and second is respectively B zone voltage status amount and each node injecting power in the variable, and third and fourth is respectively voltage and the injecting power amount of B virtual generator that the zone connects; Formula (3), (4) are the equation relation of Coordination;

All there is following equation set in each node in the subregion:

$\left\{\begin{array}{cc}{P}_i=U_i\underset{j\∈C}{\mathrm{\Σ}}{U}_j(G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}})& i,j=1,...N\\ Q_i=U_i\underset{j\∈C}{{\mathrm{\Σ}}_j}(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})& i,j=1,...N\end{array}\right.---\left(5\right);$

Wherein: G, B are real part and the imaginary part of the bus admittance matrix of partition network; C is set of network nodes, P _i, Q _iBe respectively the meritorious and reactive power of injection at node i place; U _i, U _j, θ _IjBe respectively node i, the voltage magnitude of j and node i, the phase difference of voltage of j;

System of equations (5) is carried out the single order Taylor expansion at the approximate solution place, and is deformed into following increment equation form:

$\left[\begin{array}{cc}H& N\\ M& L\end{array}\right]\left[\begin{array}{c}\mathrm{U\Δ\θ}\\ \mathrm{\ΔU}\end{array}\right]=\left[\begin{array}{c}\mathrm{\ΔP}/U\\ \mathrm{\ΔQ}/U\end{array}\right]---\left(6\right);$

Matrix H, N, M, L is respectively

Wherein each element of matrix is:

$\left\{\begin{array}{c}{H}_{\mathrm{ij}}=B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}+w(i,j)\×Q_i/U_i^2\\ N_{\mathrm{ij}}=-G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-w(i,j)\×P_i/U_i^2\\ M_{\mathrm{ij}}=G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-w(i,j)\×P_i/U_i^2\\ L_{\mathrm{ij}}=B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-w(i,j)\×Q_i/U_i^2\end{array}\right.---\left(7\right);$

In the formula: w (i, j) expression: its value is 1 when i=j; Its value is 0 when i ≠ j; U, θ are respectively amplitude vector and the phase angle vector of Radial network node voltage; Δ U, Δ θ are respectively the first increment of node voltage amplitude and phase angle in the power flow equation; Δ P, Δ Q are respectively the first increment that each node injects meritorious and reactive power;

System of equations (6) is carried out piecemeal:

$\left[\begin{array}{c}\mathrm{\Δ}{S}_n\\ \mathrm{\Δ}{S}_g\end{array}\right]=\left[\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right]\×\left[\begin{array}{c}\mathrm{\Δ}{U}_n\\ \mathrm{\Δ}{U}_g\end{array}\right]---\left(8\right);$

Wherein: Δ S=[Δ P/U Δ Q/U] ^T, Δ U=[U Δ θ Δ U] ^TΔ S _n, Δ S _gBe respectively non-virtual generator node corresponding to Δ S and virtual generator node section, Δ U _n, Δ U _gBe respectively non-virtual generator node corresponding to Δ U matrix and virtual generator node section; J ₁₁, J ₁₂, J ₂₁, J ₂₂Be respectively the matrix H of distinguishing behind virtual generator node and the non-virtual generator node, N, J, the piecemeal combination of L;

Formula (8) equal sign left side is power increment, and the equal sign right side is voltage increment; Virtual generator node voltage U when subregion carries out trend calculating, θ is constant, and the Jacobi matrix of update equation is J ₁₁

By the proper Δ S of formula (8) _n=0 o'clock, Δ S _gWith Δ U _gRelation:

$\mathrm{\Δ}{S}_g=(J_{22}-J_{21}{J}_{11}^{-1}{J}_{12})\mathrm{\Δ}{U}_g---\left(9\right);$

When having the system balancing node, subregion removes matrix J ₂₂, J ₂₁, J ₁₂The ranks that middle balance node is corresponding; When subregion contains PV type node, from matrix, remove PV node Δ U _PVCorresponding row and Δ Q _PVCorresponding row.

Wherein, among the described step C, the areas of disconnection interconnection, with the impedance mean allocation on the interconnection in the partitioned area of its contact, the breakpoint place all adds virtual generator node, and all be set to U θ node (U and θ are known), described U θ node mark the one Initial Voltage Value is made as 1, and the voltage phase angle initial value is set to 0.

Wherein, in the described step e, formula (9) is characterizing each level of coordinating to the idle increment of exerting oneself of each virtual generated power after the virtual generator voltage correction; The amount of mismatch of through type (9) and current each virtual generator power equation in split vertexes place, will be so that amount of mismatch be more and more less to each split vertexes voltage correction, when equation (4) equal sign left side is worth less than permissible error ε, then equation (4) satisfies, iteration finishes, jump to step F, otherwise return step D.

Wherein, adopt following permanent Jacobian matrix to carry out the ectonexine iteration, comprising:

The voltage of each sub regions-power sensitivity matrix is sent to cooperation layer, obtains following matrix:

$\left[\begin{array}{c}\mathrm{\Δ}{S}_g^{\′}\\ \mathrm{\Δ}{S}_g^{\′\′}\end{array}\right]=\left[\begin{array}{c}{J}^{\′}\\ J^{\′\′}\end{array}\right]\mathrm{\Δ}{U}_g---\left(10\right);$

Each split point update equation is formula (11):

[ΔS' _g+ΔS″ _g]＝[J'+J'']×ΔU _g （11）；

When adopting Newton-Laphson method to carry out Load Flow Solution, all to again form Jacobi matrix in each iteration, again Jacobi matrix be formed factor table to be used for finding the solution of Incremental Equation when finding the solution Incremental Equation; Jacobi matrix is as follows in the abbreviation formula (7):

$\left[\begin{array}{cc}H& N\\ M& L\end{array}\right]\≈\left[\begin{array}{cc}B& -G\\ G& B\end{array}\right]---\left(12\right);$

The speed of convergence that adopts the method for permanent Jacobi matrix is single order.

Compared with the prior art, the beneficial effect that reaches of the present invention is:

1. this parallel flow determines that the method convergence is good, and speed of convergence is first-order linear;

2. the method can be deployed on many computing machines relatively easily, carries out the distributed parallel trend and calculates, and can guarantee faster speed;

3. only need to store the constant Jacobian matrix matrix at cooperation layer, memory usage is few, can flexible expansion.

Description of drawings

Fig. 1 is composition decomposition computing block diagram provided by the invention;

Fig. 2 is the network diagram in the A of being decomposed into provided by the invention, B zone;

Fig. 3 is the process flow diagram of determining method based on heuristic automatic network subregion electric system parallel flow provided by the invention;

Fig. 4 is provided by the invention based on didactic electric system auto-partition process flow diagram.

Embodiment

Below in conjunction with accompanying drawing the specific embodiment of the present invention is described in further detail.

One, parallel load flow algorithm comprises:

With network partition, each subregion becomes independently subnet and coordinates by interconnection and all the other subregions from the structure of distribution network in the present invention.The principle of subregion be each subregion scale should be quite and interconnection quantity can not be too much.Each subregion carries out trend calculating independently when calculating, and after finishing the parameter relevant with Coordination is sent to coordination server, and coordination server sends it to each subregion after calculating Coordination.The block diagram that composition decomposition calculates as shown in Figure 1.

Finish owing to all need each subsystem all to calculate during each calculating of coordinating level, so each subsystem scale can not have big difference, otherwise the mutual stand-by period that can occur growing is caused the wasting of resources.Coordinate in addition mainly handle data communication of level, its computation burden should be tried one's best little otherwise can be affected whole counting yield, and the calculated amount of coordinating level depends primarily on the number of Coordination, secondly is the selection of tuning algorithm.For better parallel effect can be arranged, partitioned mode fine or not most important.After network is cut apart each subnet keep cutting branch road half, namely cut branch road and split from mid point punishment, and with the voltage of split point as Coordination.This dividing method is node split in essence, but different from general node split, and the node of division belongs to dummy node.The below provides the computing method of composition decomposition take the network that is split into two districts as example.

As shown in Figure 2, with the generator of split vertexes place access of virtual injecting power and the node voltage with equivalent split point, behind the voltage of given virtual generator or power, just can carry out independently trend to each subregion and calculate.General distribution network carries out may occurring only having a zone that the situation of balance node is arranged behind the subregion, do not have the zone of balance node only to know that the injecting power of interconnection can't carry out trend and calculate, the voltage of split point then can not occurred carrying out trend as Coordination calculate this situation.

Provided by the inventionly determine that based on heuristic automatic network subregion electric system parallel flow method flow as shown in Figure 3, comprises the steps:

A, distribution network is carried out subregion, the process flow diagram of distribution network subregion comprises the steps: as shown in Figure 4

Carry out parallel load flow algorithm, if will obtain higher parallel speed-up ratio, should guarantee each subregion scale quite and all subregion interconnection sum should lack as far as possible.For meeting this requirement, the present invention has designed a kind of new electric system auto-partition algorithm based on heuristic rule.By this algorithm, whole electric system can be decomposed into a plurality of subregions that are applicable to parallel load flow algorithm fast and effectively.

The quality of network partition method is applied to different situations, and its evaluation criterion is also different.The subregion that here will do is mainly used in parallel load flow algorithm.At present, the quality of subregion there is no clear and definite standard and weighs, but different partitioned mode is when carrying out follow-up calculating, the computer resource that consumes, and working time is different.Analyze the distributed parallel calculation process after the subregion, partitioned mode should be so that institute divide each regional calculation scale should be close preferably as can be known, otherwise can waste computer resource owing to mutual wait; In addition, the flow process of an iterative computation was roughly " distributed parallel-coordinations calculating-distributed parallel " during Distributed Power Flow and optimal load flow calculated, and from the angle consideration of operation time, it is as far as possible little to make each time interative computation always calculate scale (time scale) as far as possible.Therefore when carrying out subregion for the distributed parallel trend of electric power networks and optimal load flow calculating, should mainly consider above-mentioned two aspects.For the quality of parallel efficiency calculation behind the reflection subregion that can be more definite, defined two parameters of the total complexity of resource utilization and parallel computation, and set up corresponding mathematical model and make it to quantize.

What resource utilization characterized is the whole utilization ratio of server, and n is the number of partitions of network.P _iIt is the interstitial content of i subregion.Depend on its node scale because the subregion calculation scale is approximate, use P here _i ²Characterize the calculation scale of subregion i.After calculating separately and all finish, each subregion need to utilize each child partition server computational data to carry out the calculating of Coordination by coordination server.If each child partition calculation scale has big difference, can cause each subregion to be waited for mutually and cause the waste of resource.η is more approaching the closer to 1 each subregion calculation scale of expression, and the level of resources utilization is also higher.

The total complexity of parallel computation: O=max (P _i ²)+M * l ²＜2 〉;

The estimation equation of above-mentioned parallel computation complexity is obtained by preamble parallel computation process.The step of parallel computation can be sketched and be following three steps: at first each child servers is carried out the separately trend calculating in zone, and Coordination relevant data result sent to coordination server, then the data that send by each child servers of coordination server calculate the modified value of each Coordination, the data that to coordinate at last value send to each child servers, the power flow equation of each each child partition of child servers parallel computation.Calculating in the general electric system need to be carried out repeatedly iteration, only need to repeat above-mentioned steps again until satisfy needed precision.The parallel computation complexity of first corresponding first step in the above-mentioned O formula; The computation complexity of second corresponding coordination server.The node scale is P _iThe computation complexity of child partition i be P _i ², l is the Coordination number, l ²For coordinating the calculation scale of level server.Coordinating level server Main Resources is used for being responsible for and the communicating by letter of each child partition, calculation scale is compared each child partition calculation server should be little many, second of total complexity O formula be multiply by coefficient M(2～5 here) be higher than the calculation cost of each child servers to characterize the cooperation layer server, the parallel computation total scale of total employed partition method of the larger representative of complexity O is larger.Because the contact amount number C in each subregion and the external world _iMuch smaller than P _i, therefore total the 3rd of complexity O formula is less with respect to first, causes Coordination too much so that does not satisfy l=min (P if subregion is extremely unreasonable ₁, P ₂KP _N-1, P _n) time, total the 3rd of complexity O formula will account for larger specific gravity, and the calculated amount than subregion not is large on the contrary to make scoring area.Comprehensively above-mentioned, the scale of each child partition should approach, and this conclusion is consistent about the conclusion of resource utilization η with preamble.

As can be known described by preamble, with after the network division, all subregion node scale of division should be as far as possible close, so just can make scoring area after, index η is large and complexity O is less.For the electric power networks of a N node, if it will be divided into m subregion, then the interstitial content of each subregion should be about N/m, and for guaranteeing to round, this paper chooses the desirable nodes of each subregion and is [N/m], and this value is near the integer of N/m.In order to realize preferably above-mentioned purpose, the heuristic search that adopts node to recall here carries out subregion.Its basic ideas are to recall forward from the network frontier node, thereby in only place the branch road disconnection are formed a new subregion, and new number of partitions should be suitable with average [N/m].Its primitive rule is when dating back to a certain node, if subregion node sum is suitable with average [N/m] for this node and downstream thereof, its next-door neighbour's upstream branch road is disconnected; Otherwise continue to recall toward the upstream.

As shown in Figure 2, now be described below based on the detailed process of the network auto-partition algorithm of heuristic rule:

1) at first will become the radiativity network for network after the off-state with the chord of network is virtual.The Radial network topology is the simplest connected network, it is divided also relatively easily realize.

2) secondly the radiativity network is numbered, the method for son is carried out degree of depth traversal behind the employing late father, is father node for two node serial number smallers on the branch road like this.Record the numbering of each frontier node, define the D value of each node, D _iBe to comprise not total node number of subregion node of node i self and downstream thereof.Its meaning is that if will connect the Branch cutting of i and father node thereof, the new subregion interstitial content that comprises the i node that splits off so is D _iIndividual.

3) by after the above-mentioned degree of depth traversal numbering all frontier nodes being stored in an array, select any one frontier node to begin to recall, in the frontier node set, reject this node simultaneously.

4) D that sets again frontier node is 1, to any one node k, has:

$D_k=1+\underset{j\∈k}{\mathrm{\Σ}}{D}_j---\left(13\right);$

J is the downstream node that directly links to each other with node k.D appears when recalling the k node _kWhen 〉=[N/m] and k node child node are unique, select to connect the branch road of k node and father node thereof as cutting branch road splitting network, produce a new subregion, so do can guarantee generation new subregion interstitial content near average [N/m], then with the D value D of k node _kBe set to 0, continue to recall toward the upstream of k node;

5) if when this moment, k node child node was not unique, choose D value in k node and its child node near the node of [N/m], will this node of connection and the branch cutting of its father node;

6) this node D value is set to 0, and from then on node continues to recall.

7) when the D that dates back to certain its child node of node does not also all obtain, select another frontier node of also recalling to recall, can guarantee that all frontier nodes have all been recalled after, the D value of all nodes of network all can be obtained.

8) according to the method described above, divided m-1 rear minute end of extent, output subregion result.

Radial network has been divided into m zone after carrying out above-mentioned steps, and disconnected branches becomes interregional interconnection.Owing to also disconnected before this chord of network, also to restore after subregion is finished in addition, if two nodes of chord at different subregions, then this chord also is interregional interconnection; If two nodes of chord are at same subregion, then this chord also belongs to this subregion.The subregion of network is all finished by this.

Adopt said method to carry out to ensure higher resource utilization behind the subregion, and the total complexity of parallel computation along with having one more, the change of number of partitions reduce first the large trend of rear change.This is that this can make cooperation layer bear larger communication pressure and calculating pressure because subregion too much can cause the interconnection number also to become many, too much interconnections and mean too much Coordination.For the number of partitions situation of in advance appointment not, can be since two subregions, the enhancing number of partitions.The total complexity of the previously described parallel computation of employing was come subregion is assessed after subregion was finished, increase number of partitions, until stop increasing number of partitions when the total complexity parameter of parallel computation begins to descend during a certain number of partitions, the number of partitions of this moment be the number of partitions of the best in theory.

The voltage of B, the admittance matrix that generates the subregion node, normal Jacobian matrix that trend is calculated and the virtual generator of subregion-power sensitivity matrix, and voltage-power sensitivity matrix is sent to cooperation layer:

Power flow equation behind the subregion is as follows:

Formula (1) is the a-quadrant power flow equation, and first and second is respectively a-quadrant voltage status amount and each node injecting power in the variable, and third and fourth is respectively voltage and the injecting power amount of virtual generator that the a-quadrant connects; Formula (2) is B zone power flow equation, and first and second is respectively B zone voltage status amount and each node injecting power in the variable, and third and fourth is respectively voltage and the injecting power amount of B virtual generator that the zone connects; Formula (3), (4) are the equation relation of Coordination.

Carry out assignment for before each subregion calculates beginning voltage and the phase angle of the virtual generator of each subregion interconnection punishment knick point, the reference angle of phase angle is total system balance node angle.Each subregion can begin to carry out separately at this moment trend and calculate, and obtain the injecting power of virtual generator, since each split point voltage and virtual voltage there are differences every regional interconnection on two virtual generator injecting power summations connecing non-vanishing, be that formula (4) does not satisfy, amount of mismatch by calculating formula (4) is to the voltage correction of each virtual generator in interconnection place, constantly iteration until amount of mismatch with enough precision near 0.Therefore in whole computation process Chinese style (1), set up all the time (2), (3), and formula (4) is then constantly being set up with certain precision in the iteration.

The key of finding the solution is how cooperation layer revises voltage and the angle values of the virtual generator in split vertexes place.Since need to be to the correction of split point voltage so that two generated powers in each split point place, reactive power sum are zero, therefore need to establish the virtual generator voltage variable quantity in split point place and concern with equation between its injecting power variable quantity.

The node type in the diffServ network not, all there is following equation in each node:

$\left\{\begin{array}{cc}{P}_i=U_i\mathrm{\Σ}{U}_j(G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}})& i,=1,...N\\ Q_i=U_i\mathrm{\Σ}{U}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})& i,=1,...N\end{array}\right.---\left(5\right);$

G, B are real part and the imaginary part of the bus admittance matrix of network.System of equations (5) is carried out the single order Taylor expansion near certain approximate solution, and is deformed into following increment equation form:

$\left[\begin{array}{cc}H& N\\ M& L\end{array}\right]\left[\begin{array}{c}\mathrm{V\Δ\θ}\\ \mathrm{\ΔV}\end{array}\right]=\left[\begin{array}{c}\mathrm{\ΔP}/V\\ \mathrm{\ΔQ}/V\end{array}\right]---\left(6\right);$

Wherein each element of matrix is:

$\left\{\begin{array}{c}{H}_{\mathrm{ij}}=B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}+w(i,j)\×Q_i/{V_i}^2\\ N_{\mathrm{ij}}=-G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-w(i,j)\×P_i/{V_i}^2\\ M_{\mathrm{ij}}=G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-w(i,j)\×P_i/{V_i}^2\\ L_{\mathrm{ij}}=B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-w(i,j)\×Q_i/{V_i}^2\end{array}\right.---\left(7\right);$

W (i, j) expression in the formula: its value is 1 when i=j; Its value is 0 when i ≠ j; U, θ are respectively amplitude vector and the phase angle vector of Radial network node voltage; Δ U, Δ θ are respectively the first increment of node voltage amplitude and phase angle in the power flow equation; Δ P, Δ Q are respectively the first increment that each node injects meritorious and reactive power.

Be the relation between analyzing virtual generator voltage increment and the injecting power increment, system of equations (6) carried out piecemeal:

$\left[\begin{array}{c}\mathrm{\Δ}{S}_n\\ \mathrm{\Δ}{S}_g\end{array}\right]=\left[\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right]\×\left[\begin{array}{c}\mathrm{\Δ}{U}_n\\ \mathrm{\Δ}{U}_g\end{array}\right]---\left(8\right);$

Δ S=[Δ P/U Δ Q/U wherein] ^T, Δ U=[U Δ θ Δ U] ^TΔ S _n, Δ S _gBe respectively non-virtual generator node corresponding to Δ S and virtual generator node section, Δ U _n, Δ U _gBe respectively non-virtual generator node corresponding to Δ U matrix and virtual generator node section; J ₁₁, J ₁₂, J ₂₁, J ₂₂Be respectively the matrix H of distinguishing behind virtual generator node and the non-virtual generator node, N, J, the piecemeal combination of L;

Formula (8) left side is power increment, and the right side is voltage increment, system of equations the first behavior PQ node group, V θ node group in the second behavior network.V θ node voltage was constant when subregion carried out trend calculating, and the Jacobi matrix of update equation is J ₁₁And coordinate need to know when level is calculated when origin node all subregion in and load when constant its idle variation of exerting oneself of gaining merit when virtual generator voltage variation is only arranged.

Can be derived as Δ S by formula (8) _n=0 o'clock, Δ S _gWith Δ U _gRelation:

$\mathrm{\Δ}{S}_g=(J_{22}-J_{21}{J}_{11}^{-1}{J}_{12})\mathrm{\Δ}{U}_g---\left(9\right);$

When having the system balancing node, subregion removes matrix J ₂₂, J ₂₁, J ₁₂The ranks that middle balance node is corresponding; When subregion contains PV type node, from matrix, remove PV node Δ V _PVCorresponding row and Δ Q _PVCorresponding row.

Formula (9) is characterizing each level of coordinating to the idle increment of exerting oneself of each virtual generated power after the virtual generator voltage correction.

C, sent initial voltage amplitude and the phase angle of each virtual generator in split point place by cooperation layer: the areas of disconnection interconnection, with the impedance mean allocation on the interconnection in the partitioned area of its contact, the breakpoint place all adds virtual generator node, and all be set to V θ node, described V θ node mark the one Initial Voltage Value is made as 1, and the voltage phase angle initial value is set to 0.

D, each subregion is carried out trend calculate, and the output power of each virtual generator is sent to coordination server;

E, cooperation layer are according to the amount of mismatch of each split point place two virtual generated output power sum, judge that whether the border coordination accuracy of measurement satisfies, and does not then revise split point voltage if do not satisfy, and correction result is sent to all subregion, jump to step D, otherwise jump to step F:

Formula (9) is characterizing each level of coordinating to the idle increment of exerting oneself of each virtual generated power after the virtual generator voltage correction; The amount of mismatch of through type (9) and current each virtual generator power equation in split vertexes place, will be so that amount of mismatch be more and more less to each split vertexes voltage correction, when equation (4) equal sign left side is worth less than permissible error ε, then equation (4) satisfies, iteration finishes, jump to step F, otherwise return step D.Trend calculating in the subregion is called the internal layer iteration here, coordinates a level computing and be called external iteration.

The sensitivity matrix of each sub regions is sent to the coordination level, can obtain following matrix in the coordination level:

$\left[\begin{array}{c}\mathrm{\Δ}{S}_g^{\′}\\ \mathrm{\Δ}{S}_g^{\′\′}\end{array}\right]=\left[\begin{array}{c}{J}^{\′}\\ J^{\′\′}\end{array}\right]\mathrm{\Δ}{U}_g---\left(10\right);$

Because two groups of sensitivity matrix J ' that corresponding two the virtual generators of each split vertexes therefore can corresponding dimension equate and J " because what need to revise is each split point place power sum, also namely to makes and revise rear Δ S' _g+ Δ S " _g+ S' _g+ S " _g=0.Therefore each split point update equation is formula (11):

[ΔS' _g+ΔS″ _g]＝[J'+J'']×ΔU _g （11）；

When adopting traditional Newton-Laphson method to carry out Load Flow Solution, all to again form Jacobi matrix in each iteration, also need when finding the solution Incremental Equation again Jacobi matrix to be formed factor table to be used for finding the solution of Incremental Equation.Jacobi matrix kind element in the investigation formula (7) can carry out approximate processing to it.Generally speaking, each node phase angle difference is very little in the electric power networks, and θ is arranged _i-θ _j≈ 0, approximate thinks cos θ _Ij=1, sin θ _Ij=0; Have the 3rd during i=j in the Jacobi matrix element in the formula (7) in addition, the 3rd real imaginary part that approximates the load equivalent admittance since in the distribution network load equivalent admittance is former it can be omitted less than line admittance.Therefore Jacobi matrix can be done following simplification in the formula (7):

$\left[\begin{array}{cc}H& N\\ M& L\end{array}\right]\≈\left[\begin{array}{cc}B& -G\\ G& B\end{array}\right]---\left(12\right);$

Behind the permanent Jacobi matrix above adopting, the all subregion sensitivity matrix is also done corresponding correction in the formula (9), because each iteration all subregion Jacobi matrix is fixed, coordinate the matrix of coefficients of the update equation of level and also fix, having saved so the each external iteration that adopts accurate Jacobi matrix to cause all needs to recomputate and send separately the sensitivity matrix to computing time and the call duration time of cooperation layer.Because the update equation group matrix of coefficients of ectonexine iteration is all constant in whole computation process behind the permanent Jacobi matrix of employing, so just do not need repeatedly to generate the factor table of ectonexine iteration update equation group, greatly reduced the time of each iteration.Different is to adopt the method for accurate Jacobi matrix to have the characteristic of second order convergence, is single order and adopt the speed of convergence of the method for permanent Jacobi matrix.But facts have proved, be about ε=10 ^-6Under the precision of MVA, the method for change Jacobi matrix is lacked 1～2 time than the iterations of permanent jacobi method, and generally speaking, permanent jacobi method has higher counting yield.Among Fig. 3: e1 is the precision of subregion trend; E2 is the precision that boundary number is coordinated.

F, parallel load flow algorithm finish, the output data result.

A kind of electric system parallel flow based on heuristic automatic network subregion provided by the invention is determined method, for the parallel actual demand of finding the solution of electric system problem, at first realize the auto-partition of electric power networks by a kind of heuristic rule, secondly realize zone Divided Parallel Calculation by interregional decoupling method, realize global convergence by the iteration convergence condition at last, be fit to trend electric system, that have stronger concurrency and the parallel method for solving of Optimal Power Flow thereby develop a cover.

Should be noted that at last: above embodiment is only in order to illustrate that technical scheme of the present invention is not intended to limit, although with reference to above-described embodiment the present invention is had been described in detail, those of ordinary skill in the field are to be understood that: still can make amendment or be equal to replacement the specific embodiment of the present invention, and do not break away from any modification of spirit and scope of the invention or be equal to replacement, it all should be encompassed in the middle of the claim scope of the present invention.

Claims (6)

1. determine method based on heuristic automatic network subregion electric system parallel flow for one kind, it is characterized in that described method comprises the steps:

A, distribution network is carried out subregion;

The voltage of B, the admittance matrix that generates the subregion node, normal Jacobian matrix that trend is calculated and the virtual generator of subregion-power sensitivity matrix, and voltage-power sensitivity matrix is sent to cooperation layer;

C, sent initial voltage amplitude and the phase angle of each virtual generator in split point place by cooperation layer;

D, each subregion is carried out trend calculate, and the output power of each virtual generator is sent to coordination server;

E, cooperation layer are according to the amount of mismatch of each split point place two virtual generated output power sum, judge that whether the border coordination accuracy of measurement satisfies, and does not then revise split point voltage if do not satisfy, and correction result is sent to all subregion, jump to step D, otherwise jump to step F;

F, parallel load flow algorithm finish, the output data result.

2. electric system parallel flow as claimed in claim 1 is determined method, it is characterized in that, in the described steps A, the power distribution network subregion is based on heuristic automatic network subregion, comprises the steps:

1) network is carried out topological analysis, with the network chord virtual for off-state so that network becomes radial pattern, Radial network is numbered and degree of depth traversal, and the numbering of each frontier node is stored in the array;

2) definition comprises frontier node in the D of all interior nodes value, makes D _iBe to comprise not total node number of subregion node of node i self and downstream thereof, each node D value of initialization is-1;

3) choose in the frontier node array arbitrary frontier node and recall as present node, setting frontier node D value is 1, simultaneously this frontier node of rejecting in the frontier node array;

4) recall from present node, judge each child node D of its father node _jWhether all definitely be worth: if, find the solution the D value of father node according to following expression (13), and with this father node as present node; If not, return step 3);

$D_k=1+\underset{j\∈k}{\mathrm{\Σ}}{D}_j---\left(13\right);$

5) the D value of judging frontier node whether satisfy D 〉=[N/m] and the present node child node whether unique: if, carry out step 6), if not, then select D value in present node and the child node thereof near the node of [N/m] as present node, and carry out step 6); Wherein N represents nodes, and m is the number of partitions;

The branch road that 6) will connect present node and father node thereof is as the cutting branch road, and search present node downstream D value is not all nodes of 0 or-1, is combined as new subregion with present node;

7) present node D value is modified to 0, and from then on node continues to recall as present node; Guarantee that all frontier nodes all recall, Radial network divides end of extent;

8) output subregion result.

3. electric system parallel flow as claimed in claim 1 is determined method, it is characterized in that, among the described step B, the power flow equation behind the subregion is as follows:

Wherein: formula (1) is the a-quadrant power flow equation, and first and second is respectively a-quadrant voltage status amount and each node injecting power in the variable, and third and fourth is respectively voltage and the injecting power amount of virtual generator that the a-quadrant connects; Formula (2) is B zone power flow equation, and first and second is respectively B zone voltage status amount and each node injecting power in the variable, and third and fourth is respectively voltage and the injecting power amount of B virtual generator that the zone connects; Formula (3), (4) are the equation relation of Coordination;

All there is following equation set in each node in the subregion:

$\left\{\begin{array}{cc}{P}_i=U_i\underset{j\∈C}{\mathrm{\Σ}}{U}_j(G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}})& i,j=1,...N\\ Q_i=U_i\underset{j\∈C}{{\mathrm{\Σ}}_j}(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})& i,j=1,...N\end{array}\right.---\left(5\right);$

Wherein: G, B are real part and the imaginary part of the bus admittance matrix of partition network; C is set of network nodes, P _i, Q _iBe respectively the meritorious and reactive power of injection at node i place; U _i, U _j, θ _IjBe respectively node i, the voltage magnitude of j and node i, the phase difference of voltage of j;

System of equations (5) is carried out the single order Taylor expansion at the approximate solution place, and is deformed into following increment equation form:

$\left[\begin{array}{cc}H& N\\ M& L\end{array}\right]\left[\begin{array}{c}\mathrm{U\Δ\θ}\\ \mathrm{\ΔU}\end{array}\right]=\left[\begin{array}{c}\mathrm{\ΔP}/U\\ \mathrm{\ΔQ}/U\end{array}\right]---\left(6\right);$

Matrix H, N, M, L is respectively

Wherein each element of matrix is:

$\left\{\begin{array}{c}{H}_{\mathrm{ij}}=B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}+w(i,j)\×Q_i/U_i^2\\ N_{\mathrm{ij}}=-G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-w(i,j)\×P_i/U_i^2\\ M_{\mathrm{ij}}=G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-w(i,j)\×P_i/U_i^2\\ L_{\mathrm{ij}}=B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-w(i,j)\×Q_i/U_i^2\end{array}\right.---\left(7\right);$

In the formula: w (i, j) expression: its value is 1 when i=j; Its value is 0 when i ≠ j; U, θ are respectively amplitude vector and the phase angle vector of Radial network node voltage; Δ U, Δ θ are respectively the first increment of node voltage amplitude and phase angle in the power flow equation; Δ P, Δ Q are respectively the first increment that each node injects meritorious and reactive power;

System of equations (6) is carried out piecemeal:

$\left[\begin{array}{c}\mathrm{\Δ}{S}_n\\ \mathrm{\Δ}{S}_g\end{array}\right]=\left[\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right]\×\left[\begin{array}{c}\mathrm{\Δ}{U}_n\\ \mathrm{\Δ}{U}_g\end{array}\right]---\left(8\right);$

Wherein: Δ S=[Δ P/U Δ Q/U] ^T, Δ U=[U Δ θ Δ U] ^TΔ S _n, Δ S _gBe respectively non-virtual generator node corresponding to Δ S and virtual generator node section, Δ U _n, Δ U _gBe respectively non-virtual generator node corresponding to Δ U matrix and virtual generator node section; J ₁₁, J ₁₂, J ₂₁, J ₂₂Be respectively the matrix H of distinguishing behind virtual generator node and the non-virtual generator node, N, J, the piecemeal combination of L;

Formula (8) equal sign left side is power increment, and the equal sign right side is voltage increment; Virtual generator node voltage U when subregion carries out trend calculating, θ is constant, and the Jacobi matrix of update equation is J ₁₁

By the proper Δ S of formula (8) _n=0 o'clock, Δ S _gWith Δ U _gRelation:

$\mathrm{\Δ}{S}_g=(J_{22}-J_{21}{J}_{11}^{-1}{J}_{12})\mathrm{\Δ}{U}_g---\left(9\right);$

When having the system balancing node, subregion removes matrix J ₂₂, J ₂₁, J ₁₂The ranks that middle balance node is corresponding; When subregion contains PV type node, from matrix, remove PV node Δ U _PVCorresponding row and Δ Q _PVCorresponding row.

4. electric system parallel flow as claimed in claim 1 is determined method, it is characterized in that, among the described step C, the areas of disconnection interconnection, in the partitioned area of its contact, the breakpoint place all adds virtual generator node with the impedance mean allocation on the interconnection, and all is set to U θ node, described U θ node mark the one Initial Voltage Value is made as 1, and the voltage phase angle initial value is set to 0.

5. electric system parallel flow as claimed in claim 1 is determined method, it is characterized in that, in the described step e, formula (9) is characterizing each level of coordinating to the idle increment of exerting oneself of each virtual generated power after the virtual generator voltage correction; The amount of mismatch of through type (9) and current each virtual generator power equation in split vertexes place, will be so that amount of mismatch be more and more less to each split vertexes voltage correction, when equation (4) equal sign left side is worth less than permissible error ε, then equation (4) satisfies, iteration finishes, jump to step F, otherwise return step D.

6. electric system parallel flow as claimed in claim 5 is determined method, it is characterized in that, adopts following permanent Jacobian matrix to carry out the ectonexine iteration, comprising:

The voltage of each sub regions-power sensitivity matrix is sent to cooperation layer, obtains following matrix:

$\left[\begin{array}{c}\mathrm{\Δ}{S}_g^{\′}\\ \mathrm{\Δ}{S}_g^{\′\′}\end{array}\right]=\left[\begin{array}{c}{J}^{\′}\\ J^{\′\′}\end{array}\right]\mathrm{\Δ}{U}_g---\left(10\right);$

Each split point update equation is formula (11):

[ΔS' _g+ΔS″ _g]＝[J'+J'']×ΔU _g （11）；

When adopting Newton-Laphson method to carry out Load Flow Solution, all to again form Jacobi matrix in each iteration, again Jacobi matrix be formed factor table to be used for finding the solution of Incremental Equation when finding the solution Incremental Equation; Jacobi matrix is as follows in the abbreviation formula (7):

$\left[\begin{array}{cc}H& N\\ M& L\end{array}\right]\≈\left[\begin{array}{cc}B& -G\\ G& B\end{array}\right]---\left(12\right);$

The speed of convergence that adopts the method for permanent Jacobi matrix is single order.

Images