CN104578157A - Load flow calculation method of distributed power supply connection power grid - Google Patents

Abstract

A load flow calculation method of a distributed power supply connection power grid comprises the steps that S1. initial data of an electric power system are read; S2. sampling frequency N and the dimensions s of input random variables are determined; S3. an s * N order sampling matrix is generated; S4. sampling frequency is initialized, namely n is equal to 1; S5. whether n is larger than the sampling frequency N is judged, and yes, the probability statistics results of the variables are directly output; otherwise, S6 is carried out; S6. a wind power and photovoltaic power generation output model is determined, and a load random model is determined; S7. a load flow calculation model is determined; S8. an optimized economic model is determined; S9. load flow calculation is carried out; S10. data such as voltage, branch power and power generation cost of a <n>th node group are determined; and S11. a next round of load flow calculation is carried out, t is equal to t + 1, and S5 is carried out. The probability distribution of the output random variables can be estimated well, the uncertainty problem in an electricity market can be well solved, debugging manpower and material resources are saved, and production cost is lowered.

Description

A kind of tidal current computing method of distributed power source access electrical network

Technical field

The present invention relates to the application of system for distribution network of power, particularly relate to the tidal current computing method of a kind of distributed power source access electrical network.

Background technology

Wind energy, solar energy are green clean energy resourcies, greatly develop wind-powered electricity generation, photovoltaic is conducive to reducing Fossil fuel consumption, reducing carbon emission level.But because it has feature that is intermittent and randomness, control to propose requirements at the higher level to power system operation.In recent years, the development of China's wind-powered electricity generation, photovoltaic generation rapidly, installed capacity increases sharply, dissolve and encounter difficulties, the characteristic of wind energy turbine set, photovoltaic plant is more fully considered in planning and design of power system and operation control, grasp its fluctuation pattern, significant to the fail safe and economy improving system.

Wind-powered electricity generation, photovoltaic are exerted oneself very large by the impact of natural weather condition, when having extensive wind energy and photovoltaic access in system, its fluctuation of exerting oneself can thermoelectricity relatively in the past, water power scheduling of exerting oneself different.How under the condition meeting the system power equilibrium of supply and demand, priority scheduling new forms of energy, when considering new forms of energy fluctuation, allow fired power generating unit bear base load, and Different periods change of exerting oneself is less; Allow water power regulate peak load, between Different periods, larger fluctuation can be had; Consider meritorious generating optimization and idle generating optimization simultaneously and system is always generated electricity network minimal, these are all had higher requirement to the modeling of optimal load flow.Electrical network optimal load flow has very high practical value, its first time, by economy and fail safe, the combining of meritorious and idle work optimization almost Perfect, meets that Iarge-scale system interconnects, the requirement of programming and planning personnel, traffic control personnel after electrical network popularization.

Exert oneself due to new forms of energy and have uncertainty, probability optimal load flow calculates and also becomes more complicated.At present, existing pertinent literature is studied the optimal load flow containing new forms of energy.The uncertainty that blower fan is exerted oneself considered by document " considering the Stochastic Optimal Power Flow Approach of injecting power distribution ", establishes the optimal load flow model of chance constraint.Document " the wind power integration capability analysis based on probability optimal load flow " uses the particle swarm optimization algorithm probability optimal load flow model of stochastic technique, assesses the feasibility of wind power integration ability and validity.But above-mentioned study general only considers wind energy turbine set, and seldom study wind energy turbine set and photovoltaic plant simultaneously connecting system on the impact of optimal load flow.Wind energy turbine set and photovoltaic plant are exerted oneself and are all had randomness and power producing characteristics is different, add the uncertain factor in electricity market.

At present, consider that random sex optimal load flow computational methods mainly comprise Monte Carlo method, the Cumulant Method Using, point estimations, ant group algorithm etc.Monte Carlo method can well study the impact of random factor on system optimal trend, but the method needs thousands of the different running statuses of analogue system just can obtain rational result, and computing time is long, committed memory is large.Input stochastic variable separate or meet linear relationship prerequisite under, the Cumulant Method Using Gram-Charlier launches progression, Cornish-Fisher launches progression etc. and carries out matching, thus obtain the probability density function exporting stochastic variable, improve computational efficiency.Although point estimations has computational speed faster, its High Order Moment error exporting stochastic variable is larger.Ant group algorithm operand directly affects computational speed very greatly.

Summary of the invention

In order to solve the problem, the invention provides the tidal current computing method of a kind of distributed power source access electrical network, comprising the following steps:

S1: the primary data reading electric power system;

S2: the dimension s determining sampling number N and input stochastic variable;

S3: according to following 3 steps, generates s × N rank sampling matrixs, to form in point range the individual point (j=1 ..., s; N=1 ...) step as follows:

S3-1: N-1 integer 2 system numbers are represented, i.e. formula (1)

N-1＝a _R-1a _R-2…a ₂a ₁(1)

Wherein a _n∈ Z _b, Z _b=0,1 ..., b-1}, R are for meeting b ^rthe maximum of the r of≤N;

S3-2: to N-1=a _r-1a _r-2a ₂a ₁sort, obtain the sequence [d after sorting ₁d ₂d _nd _r] ^tfor formula (2)

${[d_1{d}_2\&CenterDot;\&CenterDot;\&CenterDot;d_n\&CenterDot;\&CenterDot;\&CenterDot;d_R]}^T=C_{N\&times;N}^i{[a_1{a}_2\&CenterDot;\&CenterDot;\&CenterDot;a_n\&CenterDot;\&CenterDot;\&CenterDot;a_{R-1}]}^T---\left(2\right)$

Wherein, for generator matrix, 0≤d _n≤ b-1; Introduce generator matrix to reset a ₁a ₂a _na _r-1in the position of each numeral; The position of numeral is after over-reset, and the Digital size of other dimension of each peacekeeping is identical, but the difference that puts in order, thus ensure that the uniformity of result;

S3-3: through the calculating of S3-2 step, 2 binary form of formula (3) can be expressed as:

$x_n^j=0.d_1{d}_2\&CenterDot;\&CenterDot;\&CenterDot;d_R---\left(3\right)$

Finally, 2 systems are represented 10 system numbers are converted into according to formula (2);

S4: by sampling number initialization: make n=1;

S5: the size judging n and sampling number N, if n > is N, the probability statistics result of direct output variable; If n≤N, turn S6;

S6: determine that wind-powered electricity generation and photovoltaic generation are exerted oneself model, determine Stochastic Load Model

S6-1: wind speed obeys Wei Buer distribution, active power of wind power field P _wprobability density function can be expressed as formula (4):

$f\left(P_w\right)=\frac{k}{k_{1}c}\left(\frac{P_w-k_2}{k_{1}c}\right)\exp [-{\left(\frac{P_w-k_2}{k_{1}c}\right)}^k]---\left(4\right)$

In formula: k, c are respectively form parameter and the scale parameter of Wei Buer distribution, p _rfor blower fan rated power, v _r, v _cibe respectively rated wind speed and incision wind speed;

Wind-powered electricity generation is treated to PQ node, makes power of fan factor in Load flow calculation invariable, then reactive power presses following formula (5) calculating:

In formula: for power-factor angle, for grid-connected blower fan, generally be positioned at fourth quadrant, for negative value.

S6-2: photovoltaic is exerted oneself stochastic model

In certain hour section, Intensity of the sunlight can think that obeying beta distributes, then photovoltaic plant power output P _pvprobability density function be expressed as formula (6):

$f\left(P_{\mathrm{pv}}\right)=\frac{\mathrm{\&Gamma;}(\mathrm{\&alpha;}+\mathrm{\&beta;})}{\mathrm{\&Gamma;}\left(\mathrm{\&alpha;}\right)+\mathrm{\&Gamma;}\left(\mathrm{\&beta;}\right)}{\left(\frac{P_{\mathrm{pv}}}{R_{\mathrm{pv}}}\right)}^{\mathrm{\&alpha;}-1}{\left(\frac{P_{\mathrm{pv}}}{R_{\mathrm{pv}}}\right)}^{\mathrm{\&beta;}-1}---\left(6\right)$

In formula: R _pv=A η γ _maxfor emulation peak power output, A is the solar cell emulation gross area, and η emulates total photoelectric conversion efficiency, γ _maxfor the maximum intensity of illumination in a period of time, Γ is Gamma function, and α, β are the form parameter of beta distribution;

Identical with wind-powered electricity generation, in Load flow calculation using photovoltaic plant also as PQ node;

S6-3: Stochastic Load Model

Load has time variation, and the method that a lot of related documents is proposed region load is predicted obtains its probability distribution; And as medium-term and long-term load prediction results, the probability distribution rule of load accords with normal distribution substantially; Its average and variance all can be obtained by a large amount of historical statistical datas; Like this, the probability density function that is meritorious and reactive power of load is respectively formula (7) and (8):

$f\left(P\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&delta;}}_P}\exp (-\frac{{(P-{\mathrm{\&mu;}}_P)}^2}{{{2\mathrm{\&delta;}}_P}^2})---\left(7\right)$

$f\left(P\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&delta;}}_Q}\exp (-\frac{{(Q-{\mathrm{\&mu;}}_Q)}^2}{{{2\mathrm{\&delta;}}_Q}^2})---\left(8\right)$

In formula: μ _pfor the average of active power, δ _p ²for the variance of active power, μ _qfor the average of reactive power, δ _q ²for the variance of reactive power;

S7: determine power flow algorithm

The present invention sets up all kinds of energy and gains merit and the minimum optimal load flow model for target function of reactive power generation total cost, adjustment generator output and reactive source exert oneself to meet load needs and system cloud gray model constraint as far as possible, guarantee current loads need and meet each node voltage bound and transmission line transmission power limit situation under, find the minimum generator output arrangement of feasible and total generating expense and electric network swim distribution;

S7-1: target function

The generation optimization model that the present invention builds is as follows:

$\min F=\underset{i\&Element;\mathrm{Ng}}{\mathrm{\&Sigma;}}[C_{\mathrm{Gpi}}\left(P_{\mathrm{Gi}}\left(t\right)\right)+C_{\mathrm{Gqi}}\left(Q_{\mathrm{Gi}}\left(t\right)\right)+\underset{j\&Element;\mathrm{Nq}}{\mathrm{\&Sigma;}}{C}_{\mathrm{gqj}}\left(Q_{\mathrm{gj}}\left(t\right)\right)---\left(9\right)$

C in formula (9) _gpi, C _gqifor the meritorious of unit i and reactive power generation cost function, C _gqj, C _gqjfor the reactive power generation cost function of reactive power compensator j, P _gi(t), Q _git () to be i-th generating set exert oneself at period t meritorious and idlely exert oneself, Q _gjt () is exerted oneself at the idle of period t for jth platform reactive power compensator; N _g, N _qfor generator node number and reactive-load compensation equipment number; Target function is that the generating expense of each period system is minimum;

S7-2: equality constraint

Equality constraint is the node trend Constraints of Equilibrium of each period:

$P_{\mathrm{Gi}}-P_{\mathrm{Di}}-V_i\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}{\cos \mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})=0---\left(10\right)$

$Q_{\mathrm{Gi}}+Q_{\mathrm{gi}}-{Q_{\mathrm{Di}}-V}_i\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}{\sin \mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}})=0---\left(11\right)$

In formula (10), formula (11): V _i, θ _ifor node voltage and phase angle, θ _ij=θ _i-θ _j; P _di, Q _difor burden with power and load or burden without work; G _ij, B _ijfor conductance and the susceptance of node admittance matrix;

S7-3: inequality constraints formula (12)

$\left\{\begin{array}{c}\underset{\&OverBar;}{P_{\mathrm{Gi}}}\&le;P_{\mathrm{Gi}}\&le;\stackrel{\&OverBar;}{P_{\mathrm{Gi}}}\\ \underset{\&OverBar;}{Q_{\mathrm{Gi}}}\&le;Q_{\mathrm{Gi}}\&le;\stackrel{\&OverBar;}{Q_{\mathrm{Gi}}}\\ \underset{\&OverBar;}{Q_{\mathrm{gi}}}\&le;Q_{\mathrm{gi}}\&le;\stackrel{\&OverBar;}{Q_{\mathrm{gi}}}\\ \underset{\&OverBar;}{V_i}\&le;V_i\&le;\stackrel{\&OverBar;}{V_i}\\ S_i\&le;\stackrel{\&OverBar;}{S_i}\end{array}\right.,\left\{\begin{array}{c}i\&Element;N_g\\ i\&Element;N_g\\ i\&Element;N_q\\ i\&Element;N\\ i\&Element;N_b\end{array}\right.---\left(12\right)$

In formula, for the meritorious bound of exerting oneself of generator i; for the idle bound of exerting oneself of generator i; q _gifor the idle bound of exerting oneself of reactive-load compensation equipment i; for node voltage amplitude bound; for circuit ij continues the transmission capacity limit (MVA); N, N _bfor set of node, set of fingers;

P _GT,i(t+1)-P _GT,i≤R _i,up(13)

P _GT,i(t)-P _GT,i(t+1)≤R _i,down(14)

In formula (13), R _{i, up}it is the upwards creep speed of i-th fired power generating unit; In formula (14), R _{i, down}it is the downward creep speed of i-th fired power generating unit;

S8: determine optimal economic model

S8-1: the generating expense of thermal power plant

Meritorious the exerting oneself with coal consumption amount of Coal-fired group is that standard carries out charging, the meritorious cost function C that exerts oneself of unit i _gpicalculate with formula (15).A in formula _i, b _i, c _iit is the coal consumption cost coefficient of i-th fired power generating unit;

$C_{\mathrm{Gpi}}=a_i{P}_{\mathrm{Gi}}^2+b_i{P}_{\mathrm{Gi}}+c_i---\left(15\right)$

The Reactive Power Price of Generation Side is divided into two parts: reactive capability electricity price and capacity of idle power electricity price.Idle opportunity cost and the active loss expense of what capacity of idle power electricity price related generally to is generator, the present invention is using the total reactive power generation expense of idle opportunity cost as generator side;

Idle opportunity cost is the profit corresponding to active power generate output that this generator loses because of output reactive power; If ignore the Power generation limits of prime mover, and suppose this is idle opportunity cost C _op(Q _gi) can represent such as formula (16);

$C_{\mathrm{Gqi}}\left(Q_{\mathrm{Gi}}\right)=C_{\mathrm{op}}\left(Q_{\mathrm{Gi}}\right)=k[C_{\mathrm{Gpi}}\left(S_{\mathrm{Gi},\max }\right)-C_{\mathrm{Gpi}}\left(\sqrt{S_{\mathrm{Gi},\max }^2-Q_{\mathrm{Gi}}^2}\right)]---\left(16\right)$

Formula (15) is updated in formula (16), and carries out Taylor expansion, remain into , also arrange after ignoring high-order term and obtain formula (17);

$C_{\mathrm{Gqi}}\left(Q_{\mathrm{Gi}}\right)={\mathrm{dQ}}_{\mathrm{Gi}}^2+e---\left(17\right)$

C _gqi(Q _gi) be the idle cost function of exerting oneself of generating set i, S _{gi, max}for the specified apparent power of generating set i, Q _gifor generating set i idle go out force value, k is the profit margin of power plant production active power, is generally 5%-10%;

S8-2: the generating expense of hydroelectric plant

Current China water power operating cost is generally 4 ~ 9 points/kilowatt hour, and China's thermoelectricity operating cost is about 0.09-0.19 unit/kilowatt hour, the present invention adopts the form of water power Active Generation cost formula (15) to carry out charging, and the value of its design parameter is approximate and a in thermoelectricity _i, b _i, c _idoubly, m is the ratio of thermoelectricity operating cost and water power operating cost electricity price to difference m, a _i, b _i, c _ivalue has fine setting change, to distinguish the generating expense in similar power station; The reactive power generation expense of hydroelectric plant also adopts the idle expense of exerting oneself of similar thermal power plant, and according to the charging way of formula (16), C wherein _gpiget the Active Generation cost function of corresponding hydroelectric plant;

S8-3: photovoltaic plant and wind energy turbine set generating expense

The rate for incorporation into the power network of current photovoltaic plant and wind energy turbine set is still higher than traditional energy, but along with the reduction of photovoltaic apparatus and wind power equipment cost, and country is for the reinforcement of generation of electricity by new energy subsidy policy, the further reduction of the rate for incorporation into the power network of photovoltaic generation and wind power generation it is expected to.Preferentially to call new forms of energy to greatest extent for criterion in the present invention, after order subsidy, photovoltaic generation expense and wind power generation expense are lower than the online generating price of thermoelectricity and water power, and it is identical that the choosing of its meritorious cost function chooses mode with the meritorious expense of water power;

S8-4: the generating expense of reactive-load compensation equipment

With capacitor, reactor, synchronous compensator, SVC is idle expense for fixed cost expression formula (18):

$G_{\mathrm{gqj}}\left(Q_{\mathrm{gj}}\right)=\frac{C_f}{Y\&times;365\&times;24\&times;p}{Q}_{\mathrm{gi}}=f_q{Q}_{\mathrm{gj}}---\left(18\right)$

Wherein Y is the useful life of shunt capacitor, usually gets 15 years; P is average service rate, is similar to and is taken as 2/3, C _ffor the fixed cost of capacitor unit capacity, on average can be taken as 62500 yuan/MVar, calculate f with these data _q=1.97;

S9: Load flow calculation

Utilize Lagrangian method to process the equality constraint in optimization problem, thus the optimization problem with equality constraint is converted into unconfined optimization problem; Utilize the penalty function method process inequality constraints of logrithmic barrier function method, finally solve unconstrained optimization problem optimal solution by Newton method;

The following mathematical formulae of nonlinear problem is represented:

obj min.f(x)

s.t.h(x)＝0 (19)

$\underset{\&OverBar;}{g}\&le;g\left(x\right)\&le;\stackrel{\&OverBar;}{g}$

Wherein: min.f (x) is target function, be a nonlinear function; H (x)=[h ₁(x) ..., h _m(x)] ^tfor Nonlinear Equality Constrained condition, g (x)=[g ₁(x) ..., g _r(x)] ^tfor nonlinear complementary problem.Suppose total k variable in above model, m equality constraint, r inequality constraints.During with interior point method Solve problems (19), first inequality constraints is converted into equality constraint, constructs barrier function simultaneously; First introduce slack variable l > 0, u > 0, l ∈ R for this reason ^r, u ∈ R ^r, the inequality constraints of formula (19) is converted into equality constraint, and target function is transformed into barrier function, following optimization problem A can be obtained:

$\begin{array}{cc}\mathrm{obj}& \min .f\left(x\right)-\mathrm{\&mu;}(\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln k_j+\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln u_j)\end{array}$

s.t.h(x)＝0 (11)

$g\left(x\right)+u-\stackrel{\&OverBar;}{g}=0$

$g\left(x\right)-l-\underset{\&OverBar;}{g}=0$

Wherein Discontinuous Factors u > 0; Work as l _ior u _iwhen border, be tending towards infinitely great with superior function, the minimal solution therefore meeting above obstacle target function can not find on border, just can only may obtain optimal solution when meeting l > 0, u>0; Like this, just by the conversion of target function, the optimization problem limited containing inequality is become only containing the optimization problem A of equality constraint restriction, therefore can solve by direct method of Lagrange multipliers.

The Lagrangian of Optimized model A is:

$L=f\left(x\right)-y^{T}h\left(x\right)-z^T[g\left(x\right)-1-\underset{\&OverBar;}{g}]-w^T[g\left(x\right)+u-\stackrel{\&OverBar;}{g}]-\mathrm{\&mu;}(\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln l_j+\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln u_j)---\left(12\right)$

In formula: y=[y ₁..., y _m] ^t, z=[z ₁..., z _r] ^t, w=[w ₁..., w _r] ^tbe Lagrange multiplier; The necessary condition that this problem minimum exists is the partial derivative of Lagrangian to all variablees and multiplier is 0, thus constrained optimization is converted into unconstrained optimization, next can use Newton Algorithm of the prior art;

S10: record the data such as the n-th group node voltage, branch power and cost of electricity-generating;

S11: carry out next round Load flow calculation, t=t+1, turn S5;

The present invention compared with the existing technology, has the following advantages and beneficial effect:

1. computational speed of the present invention is fast, accuracy is high, and the probability statistics information obtained can react the operation conditions of electricity market all sidedly, effectively can process the uncertain problem in electricity market, have good engineering practical value;

2. emulated by embodiment and verify, the present invention can promote the safety and economic operation of distributed power source access electrical network, reduces network loss simultaneously, improves node voltage level, the economy effectively raised and practicality;

3. the wind energy turbine set that adopts of the present invention and the node electricity price of photovoltaic plant hybrid system, network loss and branch power fluctuation situation less than independent windfarm system, the photovoltaic capacity of connecting system is larger, node electricity price is lower, the effect of Load flow calculation and control can be given full play to more comprehensively, effectively, quickly, present invention saves the man power and material of debugging, the production cost reduced, has certain economic benefit.

Accompanying drawing explanation

Fig. 1 is flow chart of steps of the present invention;

When Fig. 2 is embodiment photovoltaic access electrical network, the desired value of node electricity price;

When Fig. 3 is embodiment photovoltaic access electrical network, the standard deviation of node electricity price.

Embodiment

A tidal current computing method for distributed power source access electrical network, comprises the following steps:

S1: the primary data reading electric power system;

S2: the dimension s determining sampling number N and input stochastic variable;

S3: according to following 3 steps, generates s × N rank sampling matrixs, to form in point range the individual point (j=1, s; N=1) step as follows:

S3-1: N-1 integer 2 system numbers are represented, i.e. formula (1)

N-1＝a _R-1a _R-2···a ₂a ₁(1)

Wherein a _n∈ Z _b, Z _b=0,1, b-1}, R are for meeting b ^rthe maximum of the r of≤N;

S3-2: to N-1=a _r-1a _r-2a ₂a ₁sort, obtain the sequence [d after sorting ₁d ₂d _nd _r] ^tfor formula (2)

${[d_1{d}_2\&CenterDot;\&CenterDot;\&CenterDot;d_n\&CenterDot;\&CenterDot;\&CenterDot;d_R]}^T=C_{N\&times;N}^i{[a_1{a}_2\&CenterDot;\&CenterDot;\&CenterDot;a_n\&CenterDot;\&CenterDot;\&CenterDot;a_{R-1}]}^T---\left(2\right)$

Wherein, for generator matrix, 0≤d _n≤ b-1; Introduce generator matrix to reset a ₁a ₂a _na _r-1in the position of each numeral; The position of numeral is after over-reset, and the Digital size of other dimension of each peacekeeping is identical, but the difference that puts in order, thus ensure that the uniformity of result;

S3-3: through the calculating of S3-2 step, 2 binary form of formula (3) can be expressed as:

$x_n^j=0.d_1{d}_2\&CenterDot;\&CenterDot;\&CenterDot;d_R---\left(3\right)$

Finally, 2 systems are represented 10 system numbers are converted into according to formula (2);

S4: by sampling number initialization: make n=1;

S5: the size judging n and sampling number N, if n > is N, the probability statistics result of direct output variable; If n≤N, turn S6;

S6: determine that wind-powered electricity generation and photovoltaic generation are exerted oneself model, determine Stochastic Load Model

S6-1: wind speed obeys Wei Buer distribution, active power of wind power field P _wprobability density function can be expressed as formula (4):

$f\left(P_w\right)=\frac{k}{k_{1}c}\left(\frac{P_w-k_2}{k_{1}c}\right)\exp [-{\left(\frac{P_w-k_2}{k_{1}c}\right)}^k]---\left(4\right)$

In formula: k, c are respectively form parameter and the scale parameter of Wei Buer distribution, p _rfor blower fan rated power, v _r, v _cibe respectively rated wind speed and incision wind speed;

Wind-powered electricity generation is treated to PQ node, makes power of fan factor in Load flow calculation invariable, then reactive power presses following formula (5) calculating:

In formula: for power-factor angle, for grid-connected blower fan, generally be positioned at fourth quadrant, for negative value;

S6-2: photovoltaic is exerted oneself stochastic model

In certain hour section, Intensity of the sunlight can think that obeying beta distributes, then photovoltaic plant power output P _pvprobability density function be expressed as formula (6):

$f\left(P_{\mathrm{pv}}\right)=\frac{\mathrm{\&Gamma;}(\mathrm{\&alpha;}+\mathrm{\&beta;})}{\mathrm{\&Gamma;}\left(\mathrm{\&alpha;}\right)+\mathrm{\&Gamma;}\left(\mathrm{\&beta;}\right)}{\left(\frac{P_{\mathrm{pv}}}{R_{\mathrm{pv}}}\right)}^{\mathrm{\&alpha;}-1}{(1-\frac{P_{\mathrm{pv}}}{R_{\mathrm{pv}}})}^{\mathrm{\&beta;}-1}---\left(6\right)$

In formula: R _pv=A η γ _maxfor emulation peak power output, A is the solar cell emulation gross area, and η emulates total photoelectric conversion efficiency, γ _maxfor the maximum intensity of illumination in a period of time, Γ is Gamma function, and α, β are the form parameter of beta distribution;

Identical with wind-powered electricity generation, in Load flow calculation using photovoltaic plant also as PQ node;

S6-3: Stochastic Load Model

Load has time variation, and the method that a lot of related documents is proposed region load is predicted obtains its probability distribution; And as medium-term and long-term load prediction results, the probability distribution rule of load accords with normal distribution substantially.Its average and variance all can be obtained by a large amount of historical statistical datas; Like this, the probability density function that is meritorious and reactive power of load is respectively formula (7) and (8):

$f\left(P\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&delta;}}_P}\exp (-\frac{{(P-{\mathrm{\&mu;}}_P)}^2}{{{2\mathrm{\&delta;}}_P}^2})---\left(7\right)$

$f\left(P\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&delta;}}_Q}\exp (-\frac{{(Q-{\mathrm{\&mu;}}_Q)}^2}{{{2\mathrm{\&delta;}}_Q}^2})---\left(8\right)$

In formula: μ _pfor the average of active power, δ _p ²for the variance of active power, μ _qfor the average of reactive power, δ _q ²for the variance of reactive power;

S7: determine power flow algorithm

The present invention sets up all kinds of energy and gains merit and the minimum optimal load flow model for target function of reactive power generation total cost, adjustment generator output and reactive source exert oneself to meet load needs and system cloud gray model constraint as far as possible, guarantee current loads need and meet each node voltage bound and transmission line transmission power limit situation under, find the minimum generator output arrangement of feasible and total generating expense and electric network swim distribution;

S7-1: target function

The generation optimization model that the present invention builds is as follows:

$\min F=\underset{i\&Element;\mathrm{Ng}}{\mathrm{\&Sigma;}}[C_{\mathrm{Gpi}}\left(P_{\mathrm{Gi}}\left(t\right)\right)+C_{\mathrm{Gqi}}\left(Q_{\mathrm{Gi}}\left(t\right)\right)+\underset{j\&Element;\mathrm{Nq}}{\mathrm{\&Sigma;}}{C}_{\mathrm{gqj}}\left(Q_{\mathrm{gj}}\left(t\right)\right)---\left(9\right)$

C in formula (9) _gpi, C _gqifor the meritorious of unit i and reactive power generation cost function, C _gqj, C _gqjfor the reactive power generation cost function of reactive power compensator j, P _gi(t), Q _git () to be i-th generating set exert oneself at period t meritorious and idlely exert oneself, Q _gjt () is exerted oneself at the idle of period t for jth platform reactive power compensator; N _g, N _qfor generator node number and reactive-load compensation equipment number; Target function is that the generating expense of each period system is minimum;

S7-2: equality constraint

Equality constraint is the node trend Constraints of Equilibrium of each period:

$P_{\mathrm{Gi}}-P_{\mathrm{Di}}-V_i\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}{\cos \mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})=0---\left(10\right)$

$Q_{\mathrm{Gi}}+Q_{\mathrm{gi}}-{Q_{\mathrm{Di}}-V}_i\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}{\sin \mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}})=0---\left(11\right)$

In formula (10), formula (11): V _i, θ _ifor node voltage and phase angle, θ _ij=θ _i-θ _j; P _di, Q _difor burden with power and load or burden without work; G _ij, B _ijfor conductance and the susceptance of node admittance matrix;

S7-3: inequality constraints formula (12)

$\left\{\begin{array}{c}\underset{\&OverBar;}{P_{\mathrm{Gi}}}\&le;P_{\mathrm{Gi}}\&le;\stackrel{\&OverBar;}{P_{\mathrm{Gi}}}\\ \underset{\&OverBar;}{Q_{\mathrm{Gi}}}\&le;Q_{\mathrm{Gi}}\&le;\stackrel{\&OverBar;}{Q_{\mathrm{Gi}}}\\ \underset{\&OverBar;}{Q_{\mathrm{gi}}}\&le;Q_{\mathrm{gi}}\&le;\stackrel{\&OverBar;}{Q_{\mathrm{gi}}}\\ \underset{\&OverBar;}{V_i}\&le;V_i\&le;\stackrel{\&OverBar;}{V_i}\\ S_i\&le;\stackrel{\&OverBar;}{S_i}\end{array}\right.,\left\{\begin{array}{c}i\&Element;N_g\\ i\&Element;N_g\\ i\&Element;N_q\\ i\&Element;N\\ i\&Element;N_b\end{array}\right.---\left(12\right)$

In formula, for the meritorious bound of exerting oneself of generator i; for the idle bound of exerting oneself of generator i; q _gifor the idle bound of exerting oneself of reactive-load compensation equipment i; for node voltage amplitude bound; for circuit ij continues the transmission capacity limit (MVA); N, N _bfor set of node, set of fingers;

P _GT,i(t+1)-P _GT,i≤R _i,up(13)

P _GT,i(t)-P _GT,i(t+1)≤R _i,down(14)

In formula (13), R _{i, up}it is the upwards creep speed of i-th fired power generating unit; In formula (14), R _{i, down}it is the downward creep speed of i-th fired power generating unit;

S8: determine optimal economic model

S8-1: the generating expense of thermal power plant

Meritorious the exerting oneself with coal consumption amount of Coal-fired group is that standard carries out charging, the meritorious cost function C that exerts oneself of unit i _gpicalculate with formula (15).A in formula _i, b _i, c _iit is the coal consumption cost coefficient of i-th fired power generating unit;

$C_{\mathrm{Gpi}}=a_i{P}_{\mathrm{Gi}}^2+b_i{P}_{\mathrm{Gi}}+c_i---\left(15\right)$

The Reactive Power Price of Generation Side is divided into two parts: reactive capability electricity price and capacity of idle power electricity price.Idle opportunity cost and the active loss expense of what capacity of idle power electricity price related generally to is generator, the present invention is using the total reactive power generation expense of idle opportunity cost as generator side;

Idle opportunity cost is the profit corresponding to active power generate output that this generator loses because of output reactive power; If ignore the Power generation limits of prime mover, and suppose this is idle opportunity cost C _op(Q _gi) can represent such as formula (16);

$C_{\mathrm{Gqi}}\left(Q_{\mathrm{Gi}}\right)=C_{\mathrm{op}}\left(Q_{\mathrm{Gi}}\right)=k[C_{\mathrm{Gpi}}\left(S_{\mathrm{Gi},\max }\right)-C_{\mathrm{Gpi}}\left(\sqrt{S_{\mathrm{Gi},\max }^2-Q_{\mathrm{Gi}}^2}\right)]---\left(16\right)$

Formula (15) is updated in formula (16), and carries out Taylor expansion, remain into , also arrange after ignoring high-order term and obtain formula (17);

$C_{\mathrm{Gqi}}\left(Q_{\mathrm{Gi}}\right)={\mathrm{dQ}}_{\mathrm{Gi}}^2+e---\left(17\right)$

C _gqi(Q _gi) be the idle cost function of exerting oneself of generating set i, S _{gi, max}for the specified apparent power of generating set i, Q _gifor generating set i idle go out force value, k is the profit margin of power plant production active power, is generally 5%-10%;

S8-2: the generating expense of hydroelectric plant

Current China water power operating cost is generally 4 ~ 9 points/kilowatt hour, and China's thermoelectricity operating cost is about 0.09-0.19 unit/kilowatt hour, the present invention adopts the form of water power Active Generation cost formula (15) to carry out charging, and the value of its design parameter is approximate and a in thermoelectricity _i, b _i, c _idoubly, m is the ratio of thermoelectricity operating cost and water power operating cost electricity price to difference m, a _i, b _i, c _ivalue has fine setting change, to distinguish the generating expense in similar power station.The reactive power generation expense of hydroelectric plant also adopts the idle expense of exerting oneself of similar thermal power plant, and according to the charging way of formula (16), C wherein _gpiget the Active Generation cost function of corresponding hydroelectric plant;

S8-3: photovoltaic plant and wind energy turbine set generating expense

The rate for incorporation into the power network of current photovoltaic plant and wind energy turbine set is still higher than traditional energy, but along with the reduction of photovoltaic apparatus and wind power equipment cost, and country is for the reinforcement of generation of electricity by new energy subsidy policy, the further reduction of the rate for incorporation into the power network of photovoltaic generation and wind power generation it is expected to.Preferentially to call new forms of energy to greatest extent for criterion in the present invention, after order subsidy, photovoltaic generation expense and wind power generation expense are lower than the online generating price of thermoelectricity and water power, and it is identical that the choosing of its meritorious cost function chooses mode with the meritorious expense of water power;

S8-4: the generating expense of reactive-load compensation equipment

With capacitor, reactor, synchronous compensator, SVC is idle expense for fixed cost expression formula (18):

$G_{\mathrm{gqj}}\left(Q_{\mathrm{gj}}\right)=\frac{C_f}{Y\&times;365\&times;24\&times;p}{Q}_{\mathrm{gi}}=f_q{Q}_{\mathrm{gj}}---\left(18\right)$

Wherein Y is the useful life of shunt capacitor, usually gets 15 years; P is average service rate, is similar to and is taken as 2/3, C _ffor the fixed cost of capacitor unit capacity, on average can be taken as 62500 yuan/MVar, calculate f with these data _q=1.97;

S9: Load flow calculation

Utilize Lagrangian method to process the equality constraint in optimization problem, thus the optimization problem with equality constraint is converted into unconfined optimization problem; Utilize the penalty function method process inequality constraints of logrithmic barrier function method, finally solve unconstrained optimization problem optimal solution by Newton method;

The following mathematical formulae of nonlinear problem is represented:

obj min.f(x)

s.t.h(x)＝0 (19)

$\underset{\&OverBar;}{g}\&le;g\left(x\right)\&le;\stackrel{\&OverBar;}{g}$

Wherein: min.f (x) is target function, be a nonlinear function; H (x)=[h ₁(x) ..., h _m(x)] ^tfor Nonlinear Equality Constrained condition, g (x)=[g ₁(x) ..., g _r(x)] ^tfor nonlinear complementary problem.Suppose total k variable in above model, m equality constraint, r inequality constraints.During with interior point method Solve problems (19), first inequality constraints is converted into equality constraint, constructs barrier function simultaneously.First introduce slack variable l > 0, u > 0, l ∈ R for this reason ^r, u ∈ R ^r, the inequality constraints of formula (19) is converted into equality constraint, and target function is transformed into barrier function, following optimization problem A can be obtained:

$\begin{array}{cc}\mathrm{obj}& \min .f\left(x\right)-\mathrm{\&mu;}(\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln k_j+\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln u_j)\end{array}$

s.t. h(x)＝0 (11)

$g\left(x\right)+u-\stackrel{\&OverBar;}{g}=0$

$g\left(x\right)-l-\underset{\&OverBar;}{g}=0$

Wherein Discontinuous Factors u > 0.Work as l _ior u _iwhen border, be tending towards infinitely great with superior function, the minimal solution therefore meeting above obstacle target function can not find on border, just can only may obtain optimal solution when meeting l > 0, u>0; Like this, just by the conversion of target function, the optimization problem limited containing inequality is become only containing the optimization problem A of equality constraint restriction, therefore can solve by direct method of Lagrange multipliers.

The Lagrangian of Optimized model A is:

$L=f\left(x\right)-y^{T}h\left(x\right)-z^T[g\left(x\right)-1-\underset{\&OverBar;}{g}]-w^T[g\left(x\right)+u-\stackrel{\&OverBar;}{g}]-\mathrm{\&mu;}(\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln l_j+\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln u_j)---\left(12\right)$

In formula: y=[y ₁..., y _m] ^t, z=[z ₁..., z _r] ^t, w=[w ₁..., w _r] ^tbe Lagrange multiplier; The necessary condition that this problem minimum exists is the partial derivative of Lagrangian to all variablees and multiplier is 0, thus constrained optimization is converted into unconstrained optimization, next can use Newton Algorithm of the prior art;

S10: record the data such as the n-th group node voltage, branch power and cost of electricity-generating;

S11: carry out next round Load flow calculation, t=t+1, turn S5.

The present invention can adopt IEEE30 node instance to carry out emulating and verifying.Below in conjunction with embodiment accompanying drawing, technical scheme of the present invention is clearly and completely described.

This example directly use sampling scale be 512 times DN method analytical system access wind energy turbine set and photovoltaic plant after optimal load flow probabilistic statistical characteristics.

Table 1 node electricity price desired value compares

The randomness of exerting oneself in order to distributed energy is on the impact of economize on electricity electricity price, 2 kinds of situations are divided into discuss below: case 1, for 8 nodes 20,61,104,123,138,171,198 and 207, it is the wind energy turbine set of 60MW that each node all accesses an installed capacity; Case 2, for 8 nodes 20,61,104,123,138,171,198 and 207, each node all accesses wind energy turbine set that an installed capacity is 30MW and an installed capacity is the photovoltaic plant of 30MW.

Use the algorithm introduced to carry out the calculating of probability optimal load flow herein, obtain the node electricity price of distributed energy access node, as shown in table 1.

Can be found out by the result of calculation compared in two kinds of situations, the node electricity price of wind energy turbine set and photovoltaic hybrid system is low relative to only there being windfarm system.

Table 2 is the network loss desired value of system in two kinds of situations, and the network loss of wind energy turbine set and photovoltaic hybrid system is less than the network loss only having windfarm system, and visible wind energy turbine set and photovoltaic hybrid system are more conducive to systematic economy and run.

Table 2 network loss desired value compares

Distributed energy access way	Network loss/(MW)
		Case 1	232.622
Case 2	230.495

Table 3 is the expected value and standard deviation of branch road 13 ~ 20 and branch road 181 ~ 138 in two kinds of situations.As can be seen from result in table, branch power average when system accesses separately wind energy turbine set and standard deviation than wind energy turbine set and photovoltaic hybrid system large, branch power fluctuation is larger, occurs that heavy duty and out-of-limit probability are also more greatly, are unfavorable for that line security is checked.

Table 3 branch power compares

In order to weigh the photovoltaic connecting system of different capabilities to the impact of probability optimal load flow, at node 20,61,104,123,138, the photovoltaic plant that 171,198 and 207 accesses, 8 installed capacitys are equal, 8 photovoltaic plant total installed capacities are increased to 500MW (each stepping increases 50MW) successively from 50MW, carry out the calculating of probability optimal load flow respectively under often kind of capacity.

Fig. 2 and the expected value and standard deviation of node electricity price when Figure 3 shows that different capabilities photovoltaic plant connecting system.Fig. 2 result shows, along with the continuous increase of connecting system photovoltaic capacity, the node electricity price of photovoltaic node presents the trend of reduction.A part of traditional fired power generating unit can be replaced to exert oneself this is because photovoltaic is exerted oneself, thus node electricity price is reduced.Fig. 3 result shows, the randomness that photovoltaic is exerted oneself and uncertainty can bring the fluctuation of node electricity price.

Claims (1)

1. a tidal current computing method for distributed power source access electrical network, is characterized in that, comprise the following steps:

S1: the primary data reading electric power system;

S2: the dimension s determining sampling number N and input stochastic variable;

S3: according to following 3 steps, generates s × N rank sampling matrixs, to form in point range the individual point (j=1 ..., s; N=1 ...) step as follows:

S3-1: N-1 integer 2 system numbers are represented, i.e. formula (1)

N-1＝a _R-1a _R-2…a ₂a ₁(1)

Wherein a _n∈ Z _b, Z _b=0,1 ..., b-1}, R are for meeting b ^rthe maximum of the r of≤N;

S3-2: to N-1=a _r-1a _r-2a ₂a ₁sort, obtain the sequence [d after sorting ₁d ₂d _nd _r] ^tfor formula (2)

${[d_1{d}_2\&CenterDot;\&CenterDot;\&CenterDot;d_n\&CenterDot;\&CenterDot;\&CenterDot;d_R]}^T=C_{N\&times;N}^i{[a_1{a}_2\&CenterDot;\&CenterDot;\&CenterDot;a_n\&CenterDot;\&CenterDot;\&CenterDot;a_{R-1}]}^T---\left(2\right)$

Wherein, for generator matrix, 0≤d _n≤ b-1; Introduce generator matrix to reset a ₁a ₂a _na _r-1in the position of each numeral; The position of numeral is after over-reset, and the Digital size of other dimension of each peacekeeping is identical, but the difference that puts in order, thus ensure that the uniformity of result;

S3-3: through the calculating of S3-2 step, 2 binary form of formula (3) can be expressed as:

$x_n^j=0.d_1{d}_2\&CenterDot;\&CenterDot;\&CenterDot;d_R---\left(3\right)$

Finally, 2 systems are represented 10 system numbers are converted into according to formula (2);

S4: by sampling number initialization: make n=1;

S5: the size judging n and sampling number N, if n > is N, the probability statistics result of direct output variable; If n≤N, turn S6;

S6: determine that wind-powered electricity generation and photovoltaic generation are exerted oneself model, determine Stochastic Load Model;

S6-1: wind speed obeys Wei Buer distribution, active power of wind power field P _wprobability density function can be expressed as formula (4):

$f\left(P_w\right)=\frac{k}{k_{1}c}\left(\frac{P_w-k_2}{k_{1}c}\right)\exp [-{\left(\frac{P_w-k_2}{k_{1}c}\right)}^k]---\left(4\right)$

In formula: k, c are respectively form parameter and the scale parameter of Wei Buer distribution, k ₂=-k ₁v _ci, P _rfor blower fan rated power, v _r, v _cibe respectively rated wind speed and incision wind speed;

Wind-powered electricity generation is treated to PQ node, makes power of fan factor in Load flow calculation invariable, then reactive power presses following formula (5) calculating:

In formula: for power-factor angle, for grid-connected blower fan, generally be positioned at fourth quadrant, for negative value;

S6-2: photovoltaic is exerted oneself stochastic model

In certain hour section, Intensity of the sunlight can think that obeying beta distributes, then photovoltaic plant power output P _pvprobability density function be expressed as formula (6):

$f\left(P_{\mathrm{pv}}\right)=\frac{\mathrm{\&Gamma;}(\mathrm{\&alpha;}+\mathrm{\&beta;})}{\mathrm{\&Gamma;}\left(\mathrm{\&alpha;}\right)+\mathrm{\&Gamma;}\left(\mathrm{\&beta;}\right)}{\left(\frac{P_{\mathrm{pv}}}{R_{\mathrm{pv}}}\right)}^{\mathrm{\&alpha;}-1}{(1-\frac{P_{\mathrm{pv}}}{R_{\mathrm{pv}}})}^{\mathrm{\&beta;}-1}---\left(6\right)$

In formula: R _pv=A η γ _maxfor emulation peak power output, A is the solar cell emulation gross area, and η emulates total photoelectric conversion efficiency, γ _maxfor the maximum intensity of illumination in a period of time, Γ is Gamma function, and α, β are the form parameter of beta distribution;

Identical with wind-powered electricity generation, in Load flow calculation using photovoltaic plant also as PQ node;

S6-3: Stochastic Load Model

Load has time variation, and the method that a lot of related documents is proposed region load is predicted obtains its probability distribution; And as medium-term and long-term load prediction results, the probability distribution rule of load accords with normal distribution substantially; Its average and variance all can be obtained by a large amount of historical statistical datas; Like this, the probability density function that is meritorious and reactive power of load is respectively formula (7) and (8):

$f\left(P\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&delta;}}_P}\exp (-\frac{{(P-{\mathrm{\&mu;}}_P)}^2}{2{{\mathrm{\&delta;}}_P}^2})---\left(7\right)$

$f\left(P\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&delta;}}_Q}\exp (-\frac{{(Q-{\mathrm{\&mu;}}_Q)}^2}{2{{\mathrm{\&delta;}}_Q}^2})---\left(8\right)$

In formula: μ _pfor the average of active power, δ _p ²for the variance of active power, μ _qfor the average of reactive power, δ _q ²for the variance of reactive power;

S7: determine power flow algorithm

S7-1: target function

The generation optimization model built is as follows:

$\min F=\underset{i\&Element;\mathrm{Ng}}{\mathrm{\&Sigma;}}[C_{\mathrm{Gpi}}\left(P_{\mathrm{Gi}}\left(t\right)\right)+C_{\mathrm{Gqi}}\left(Q_{\mathrm{Gi}}\left(t\right)\right)]+\underset{j\&Element;\mathrm{Nq}}{\mathrm{\&Sigma;}}{C}_{\mathrm{gqj}}\left(Q_{\mathrm{gj}}\left(t\right)\right)---\left(9\right)$

C in formula (9) _gpi, C _gqifor the meritorious of unit i and reactive power generation cost function, C _gqj, C _gqjfor the reactive power generation cost function of reactive power compensator j, P _gi(t), Q _git () to be i-th generating set exert oneself at period t meritorious and idlely exert oneself, Q _gjt () is exerted oneself at the idle of period t for jth platform reactive power compensator; N _g, N _qfor generator node number and reactive-load compensation equipment number; Target function makes the generating expense of each period system minimum;

S7-2: equality constraint

Equality constraint is the node trend Constraints of Equilibrium of each period:

$P_{\mathrm{Gi}}-P_{\mathrm{Di}}-V_i\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}})=0---\left(10\right)$

$Q_{\mathrm{Gi}}+Q_{\mathrm{gi}}-Q_{\mathrm{Di}}-V_i\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{V}_j(G_{\mathrm{ij}}\sin {\mathrm{\&theta;}}_{\mathrm{ij}}+B_{\mathrm{ij}}\cos {\mathrm{\&theta;}}_{\mathrm{ij}})=0---\left(11\right)$

In formula (10), formula (11): V _i, θ _ifor node voltage and phase angle, θ _ij=θ _i-θ _j; P _di, Q _difor burden with power and load or burden without work; G _ij, B _ijfor conductance and the susceptance of node admittance matrix;

S7-3: inequality constraints formula (12)

$\left\{\begin{array}{c}\underset{\&OverBar;}{P_{\mathrm{Gi}}}\&le;P_{\mathrm{Gi}}\&le;\stackrel{\&OverBar;}{P_{\mathrm{Gi}}}\\ \underset{\&OverBar;}{Q_{\mathrm{Gi}}}\&le;Q_{\mathrm{Gi}}\&le;\stackrel{\&OverBar;}{Q_{\mathrm{Gi}}}\\ \underset{\&OverBar;}{Q_{\mathrm{gi}}}\&le;Q_{\mathrm{gi}}\&le;\stackrel{\&OverBar;}{Q_{\mathrm{gi}}}\\ \underset{\&OverBar;}{V_i}\&le;V_i\&le;\stackrel{\&OverBar;}{V_i}\\ S_i\&le;\stackrel{\&OverBar;}{S_i}\end{array}\right.,\left\{\begin{array}{c}i\&Element;N_g\\ i\&Element;N_g\\ i\&Element;N_q\\ i\&Element;N\\ i\&Element;N_b\end{array}\right.---\left(12\right)$

In formula, for the meritorious bound of exerting oneself of generator i; for the idle bound of exerting oneself of generator i; q _gifor the idle bound of exerting oneself of reactive-load compensation equipment i; for node voltage amplitude bound; for circuit ij continues the transmission capacity limit (MVA); N, N _bfor set of node, set of fingers;

P _GT,i(t+1)-P _GT,i≤R _i,up(13)

P _GT,i(t)-P _GT,i(t+1)≤R _i,down(14)

In formula (13), R _{i, up}it is the upwards creep speed of i-th fired power generating unit; In formula (14), R _{i, down}it is the downward creep speed of i-th fired power generating unit;

S8: determine optimal economic model

S8-1: the generating expense of thermal power plant

Meritorious the exerting oneself with coal consumption amount of Coal-fired group is that standard carries out charging, the meritorious cost function C that exerts oneself of unit i _gpicalculate with formula (15); A in formula _i, b _i, c _iit is the coal consumption cost coefficient of i-th fired power generating unit;

$C_{\mathrm{Gpi}}=a_i{P}_{\mathrm{Gi}}^2+b_i{P}_{\mathrm{Gi}}+c_i---\left(15\right)$

Idle opportunity cost is the profit corresponding to active power generate output that this generator loses because of output reactive power; If ignore the Power generation limits of prime mover, and suppose this is idle opportunity cost C _op(Q _gi) can represent such as formula (16);

$C_{\mathrm{Gqi}}\left(Q_{\mathrm{Gi}}\right)=C_{\mathrm{op}}\left(Q_{\mathrm{Gi}}\right)=k[C_{\mathrm{Gpi}}\left(S_{\mathrm{Gi},\max }\right)-C_{\mathrm{Gpi}}\left(\sqrt{S_{\mathrm{Gi},\max }^2-Q_{\mathrm{Gi}}^2}\right)]---\left(16\right)$

Formula (15) is updated in formula (16), and carries out Taylor expansion, remain into , also arrange after ignoring high-order term and obtain formula (17)

$C_{\mathrm{Gqi}}\left(Q_{\mathrm{Gi}}\right)=d{Q}_{\mathrm{Gi}}^2+e---\left(17\right)$

C _gqi(Q _gi) be the idle cost function of exerting oneself of generating set i, S _{gi, max}for the specified apparent power of generating set i, Q _gifor generating set i idle go out force value, k is the profit margin of power plant production active power, is generally 5%-10%;

S8-2: the generating expense of hydroelectric plant

Adopt the form of water power Active Generation cost formula (15) to carry out charging, and the value of its design parameter is approximate and a in thermoelectricity _i, b _i, c _idoubly, m is the ratio of thermoelectricity operating cost and water power operating cost electricity price to difference m, a _i, b _i, c _ivalue has fine setting change, to distinguish the generating expense in similar power station; The reactive power generation expense of hydroelectric plant also adopts the idle expense of exerting oneself of similar thermal power plant, and according to the charging way of formula (16), C wherein _gpiget the Active Generation cost function of corresponding hydroelectric plant;

S8-3: photovoltaic plant and wind energy turbine set generating expense

Preferentially to call new forms of energy to greatest extent for criterion, after order subsidy, photovoltaic generation expense and wind power generation expense are lower than the online generating price of thermoelectricity and water power, and it is identical that the choosing of its meritorious cost function chooses mode with the meritorious expense of water power;

S8-4: the generating expense of reactive-load compensation equipment

With capacitor, reactor, synchronous compensator, SVC is idle expense for fixed cost expression formula (18):

$C_{\mathrm{gqj}}\left(Q_{\mathrm{gj}}\right)=\frac{C_f}{Y\&times;365\&times;24\&times;p}{Q}_{\mathrm{gj}}=f_q{Q}_{\mathrm{gj}}---\left(18\right)$

Wherein Y is the useful life of shunt capacitor, usually gets 15 years; P is average service rate, is similar to and is taken as 2/3, C _ffor the fixed cost of capacitor unit capacity, on average can be taken as 62500 yuan/MVar, calculate f with these data _q=1.97;

S9: Load flow calculation

Utilize Lagrangian method to process the equality constraint in optimization problem, thus the optimization problem with equality constraint is converted into unconfined optimization problem; Utilize the penalty function method process inequality constraints of logrithmic barrier function method, finally solve unconstrained optimization problem optimal solution by Newton method;

The following mathematical formulae of nonlinear problem is represented:

obj min.f(x)

s.t.h(x)＝0 (19)

$\underset{\&OverBar;}{g}\&le;g\left(x\right)\&le;\stackrel{\&OverBar;}{g}$

Wherein: min.f (x) is target function, be a nonlinear function; H (x)=[h ₁(x) ..., h _m(x)] ^tfor Nonlinear Equality Constrained condition, g (x)=[g ₁(x) ..., g _r(x)] ^tfor nonlinear complementary problem; Suppose total k variable in above model, m equality constraint, r inequality constraints; During with interior point method Solve problems (19), first inequality constraints is converted into equality constraint, constructs barrier function simultaneously; First introduce slack variable l > 0, u > 0, l ∈ R for this reason ^r, u ∈ R ^r, the inequality constraints of formula (19) is converted into equality constraint, and target function is transformed into barrier function, following optimization problem A can be obtained:

$\mathrm{obj},\min .f\left(x\right)-\mathrm{\&mu;}(\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln l_j+\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln u_j)$

s.t.h(x)＝0 (11)

$g\left(x\right)+u-\stackrel{\&OverBar;}{g}=0$

$g\left(x\right)-l-\underset{\&OverBar;}{g}=0$

Wherein Discontinuous Factors u > 0; Work as l _ior u _iwhen border, be tending towards infinitely great with superior function, the minimal solution therefore meeting above obstacle target function can not find on border, just can only may obtain optimal solution when meeting l > 0, u>0; Like this, just by the conversion of target function, the optimization problem limited containing inequality is become only containing the optimization problem A of equality constraint restriction, therefore can solve by direct method of Lagrange multipliers;

The Lagrangian of Optimized model A is:

$L=f\left(x\right)-y^{T}h\left(x\right)-z^T[g\left(x\right)-1-\underset{\&OverBar;}{g}]-w^T[g\left(x\right)+u-\stackrel{\&OverBar;}{g}]-\mathrm{\&mu;}(\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln l_j+\underset{j=1}{\overset{r}{\mathrm{\&Sigma;}}}\ln u_j)---\left(12\right)$

In formula: y=[y ₁..., y _m] ^t, z=[z ₁..., z _r] ^t, w=[w ₁..., w _r] ^tbe Lagrange multiplier; The necessary condition that this problem minimum exists is the partial derivative of Lagrangian to all variablees and multiplier is 0, thus constrained optimization is converted into unconstrained optimization, next can use Newton Algorithm of the prior art;

S10: record the data such as the n-th group node voltage, branch power and cost of electricity-generating;

S11: carry out next round Load flow calculation, t=t+1, turn S5.