CN104578157B - 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

Load flow calculation method for distributed power supply connected to power grid

Technical Field

The invention relates to the field of application of a power distribution network of a power system, in particular to a load flow calculation method for a distributed power supply connected to a power grid.

Background

Wind energy and solar energy are green clean energy, and the vigorous development of wind power and photovoltaic is beneficial to reducing fossil fuel consumption and reducing carbon emission level. But because the device has the characteristics of intermittency and randomness, higher requirements are put forward on the operation control of the power system. In recent years, wind power and photovoltaic power generation in China are developed rapidly, installed capacity is increased rapidly, and absorption is difficult, characteristics of wind power plants and photovoltaic power stations are considered more comprehensively in planning design and operation control of power systems, fluctuation rules of the wind power plants and the photovoltaic power stations are mastered, and the method has important significance for improving safety and economy of the systems.

The output of wind power and photovoltaic power is greatly influenced by natural weather conditions, and when large-scale wind power and photovoltaic power are connected into the system, the output fluctuation of the system is different from the output scheduling of the traditional thermal power and hydropower. The method includes the steps that firstly, new energy is scheduled under the condition that the balance of system power supply and demand is met, and a thermal power generating unit bears basic load under the condition that the fluctuation of the new energy is considered, so that the output change in different time periods is small; the peak load can be adjusted by hydroelectric power, and large fluctuation can be realized in different periods; meanwhile, active power output optimization and reactive power output optimization are considered, the total power generation cost of the system is the lowest, and higher requirements are provided for modeling of the optimal power flow. The optimal power flow of the power grid has high practical value, the economy is combined with the safety, the active power optimization and the reactive power optimization nearly perfectly for the first time, and the requirements of system planning designers and operation dispatchers after the interconnection of large systems and the scale expansion of the power grid are met.

The probability-optimal power flow calculation also becomes more complex due to uncertainty of new energy output. At present, relevant documents are available for researching the optimal power flow containing new energy. The literature, namely a random optimal power flow method considering injection power distribution, considers uncertainty of fan output and establishes an optimal power flow model with opportunity constraint. The document wind power access capability analysis based on the probability optimal power flow uses a particle swarm optimization algorithm of a random technology to solve a probability optimal power flow model, and feasibility and effectiveness of the wind power access capability are evaluated. However, the above research only considers the wind power plant generally, and the influence of the wind power plant and the photovoltaic power plant simultaneously accessed into the system on the optimal power flow is rarely researched. The output of the wind power station and the output of the photovoltaic power station are random and have different output characteristics, so that uncertain factors in the power market are increased.

At present, the optimal power flow calculation method considering the randomness influence mainly comprises a monte carlo method, an accumulative method, a point estimation method, an ant colony algorithm and the like. The Monte Carlo method can well research the influence of the randomness factors on the optimal power flow of the system, but the method needs thousands of times of simulation of different operation states of the system to obtain a reasonable result, and has long calculation time and large occupied memory. On the premise that the input random variables are mutually independent or satisfy a linear relation, the cumulative quantity method is fitted by using Gram-Charlie expansion series, Cornish-Fisher expansion series and the like, so that a probability density function of the output random variables is obtained, and the calculation efficiency is improved. Although the point estimation method has a fast calculation speed, the high-order moment error of the output random variable is large. The ant colony algorithm has a large amount of calculation and directly influences the calculation speed.

Disclosure of Invention

In order to solve the above problems, the present invention provides a load flow calculation method for a distributed power supply to be connected to a power grid, including the following steps:

s1: reading initial data of the power system;

s2: determining the sampling times N and the dimension s of an input random variable;

s3, generating an S × N-order sampling matrix according to the following 3 steps to form the first order in the point rowThe procedure for point (j 1, …, s; n 1, …) is as follows:

s3-1: the N-1 integer is expressed by 2-system number, namely formula (1)

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

Wherein a is_n∈Z_b，Z_b(0, 1, …, b-1), and R is B^rThe maximum value of r is less than or equal to N;

s3-2: for N-1 ═ a_R-1a_R-2…a₂a₁Sequencing to obtain a sequenced sequence [ d ]₁d₂…d_n…d_R]^TIs (2)

Wherein,to generate a matrix, d is 0 ≦ d_nB-1 is less than or equal to; introducing a generator matrixIs to reset a₁a₂···a_n···a_R-1The position of each digit in; position of the digitAfter the resetting, the number of each dimension is the same as that of other dimensions, but the arrangement sequence is different, thereby ensuring the uniformity of the result;

s3-3: through the calculation of step S3-2,can be expressed as a 2-ary form of equation (3):

finally, 2 is expressedConverting into 10-system number according to formula (2);

s4: initializing the sampling times: let n equal to 1;

s5: judging the N and the sampling times N, and if N is larger than N, directly outputting the probability statistical result of the variable; if N is less than or equal to N, turning to S6;

s6: determining a wind power and photovoltaic power generation output model and determining a load random model

S6-1: wind speed follows Weibull distribution, and active power P of wind power plant_wCan be expressed as formula (4):

in the formula: k and c are respectively the shape parameter and the scale parameter of the Weibull distribution,P_rrated power of the fan, v_r,v_ciRated wind speed and cut-in wind speed respectively;

wind power is processed into PQ nodes, the wind power factor in the load flow calculation is constant, and then reactive power is calculated according to the following formula (5):

in the formula:for power factor angle, for a grid-connected fan,is generally located in the fourth quadrant of the device,is negative.

S6-2: photovoltaic output stochastic model

Within a certain time period, the solar illumination intensity can be considered to obey the beta distribution, and then the output power P of the photovoltaic power station_pvIs expressed as formula (6):

in the formula: r_pv＝Aηγ_maxTo simulate maximum output power, A is the simulated total area of the solar cell, η is the simulated total photoelectric conversion efficiency, γ_maxα are all shape parameters of beta distribution, which is the maximum illumination intensity in a period of time and is a Gamma function;

the method is the same as wind power, and a photovoltaic power station is also used as a PQ node in load flow calculation;

s6-3: load stochastic model

The load has time-varying property, and many related documents propose a method for predicting the regional load to obtain the probability distribution of the regional load; as the load forecasting result of medium and long periods, the probability distribution rule of the load basically conforms to normal distribution; the mean and variance can be obtained from a large amount of historical statistical data; thus, the probability density functions of the active and reactive power of the load are respectively equations (7) and (8):

in the formula: mu.s_pIs the average value of the active power,_p ²is the variance of active power, mu_QIs the average value of the reactive power,_Q ²is the variance of the reactive power;

s7: determining a load flow calculation model

The method establishes an optimal power flow model with the minimum total cost of active power and reactive power generation of various energy sources as a target function, adjusts the output of the generator and the output of the reactive power source as much as possible to meet load requirements and system operation constraints, and searches feasible power output arrangement and power flow distribution state of a power grid with the minimum total power generation cost under the condition of ensuring the current load requirements and meeting the upper and lower limits of the voltage of each node and the transmission power limit of a transmission line;

s7-1: objective function

The power generation optimization model constructed by the invention is as follows:

c in formula (9)_Gpi、C_GqiAs a function of the active and reactive power generation costs of the unit i, C_gqj、C_gqjAs a function of the reactive power generation cost of the reactive power compensation device j, P_Gi(t)、Q_Gi(t) is the ith generating setActive and reactive power, Q, at time t_gj(t) the reactive power output of the jth reactive power compensation device in the time period t; n is a radical of_g、N_qThe number of the generator nodes and the number of the reactive compensation equipment are calculated; the objective function is that the power generation cost of the system is minimum in each period;

s7-2: constraint of equality

The equality constraint is a node power flow balance constraint of each time interval:

in formulas (10) and (11): v_i、θ_iIs the node voltage and phase angle, θ_ij＝θ_i-θ_j；P_Di、Q_DiActive load and reactive load; g_ij、B_ijConductance and susceptance of a node admittance matrix;

s7-3: inequality constraint type (12)

In the formula,the upper and lower active output limits of the generator i are set;the upper and lower limit of reactive power output of the generator i;Q_giupper and lower reactive power output limits for the reactive power compensation equipment i;the node voltage amplitude upper and lower limits;continuously delivering a capacity limit (MVA) for line ij; n, N_bNode set and branch set;

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 the formula (13), R_i，upThe upward climbing speed of the ith thermal power generating unit is obtained; in the formula (14), R_i，downThe downward climbing speed of the ith thermal power generating unit is obtained;

s8: determining an optimal economic model

S8-1: electricity generation cost of thermal power plant

The active output of the thermal power coal-fired unit is charged by taking the coal consumption as a standard, and the unit i has an active output cost function C_GpiThe calculation was performed by equation (15). In the formula a_i、b_i、c_iThe coal consumption cost coefficient of the ith thermal power generating unit is obtained;

the reactive power price of the power generation side is divided into two parts: reactive capacity electricity prices and reactive power electricity quantity electricity prices. The reactive power electricity price mainly relates to the reactive power opportunity cost and the active loss cost of the generator, and the reactive power opportunity cost is used as the total reactive power generation cost on the generator side;

the reactive opportunity cost is the profit corresponding to the active power generating capacity lost by the generator due to the output of reactive power; if the output limit of the prime mover is ignored, and assumeThe cost of the reactive opportunity C_op(Q_Gi) Can be represented as formula (16);

substituting the formula (15) into the formula (16), performing Taylor expansion, and retainingNeglecting higher-order terms and then finishing to obtain a formula (17);

C_Gqi(Q_Gi) As a function of the reactive power contribution cost of the generator set i, S_Gi,maxRated apparent power, Q, of the generator set i_GiThe value k is the reactive output value of the generator set i, and the profit margin of the active power produced by the power plant is generally 5% -10%;

s8-2: electricity generation cost of hydraulic power plant

At present, the running cost of hydropower in China is generally 4-9 minutes/kilowatt hour, and the running cost of thermal power in China is about 0.09-0.19 yuan/kilowatt hour, the invention adopts the form of a hydropower active power generation cost formula (15) to charge, and the value of the specific parameter is similar to that of a in thermal power_i，b_i，c_iThe difference is m times, m is the ratio of the thermal power running cost to the hydroelectric power running cost, a_i，b_i，c_iThe value is changed in a fine adjustment way so as to distinguish the power generation cost of the same power station; the reactive power generation cost of the hydraulic power plant is similar to the reactive power output cost of the thermal power plant, and the charging mode of the formula (16) is adopted, wherein C_GpiTaking an active power generation cost function of a corresponding hydraulic power plant;

s8-3: generating cost of photovoltaic power station and wind power plant

At present, the grid-surfing electricity price of a photovoltaic power station and a wind power plant is still higher than that of the traditional energy, but with the reduction of the cost of photovoltaic equipment and wind power equipment and the enhancement of the national subsidy policy for new energy power generation, the further reduction of the grid-surfing electricity price of photovoltaic power generation and wind power generation can be expected. According to the method, new energy is preferentially called to the maximum extent as a criterion, the subsidized photovoltaic power generation cost and wind power generation cost are lower than the online power generation price of thermal power and hydropower, and the active cost function is selected in the same way as the active cost of the hydropower;

s8-4: cost of power generation of reactive power compensation equipment

The reactive cost of a capacitor, a reactor, a synchronous phase modulator and SVC is taken as a fixed cost expression (18):

wherein Y is the service life of the parallel capacitor, usually 15 years; p is the average usage, taken approximately as 2/3, C_fFor a fixed cost per unit capacity of capacitor, it is preferably 62500 yuan/MVar on average, from which f is calculated_q＝1.97；

S9: load flow calculation

Processing equality constraint in the optimization problem by using a Lagrange function method, thereby converting the optimization problem with equality constraint into an unconstrained optimization problem; processing inequality constraints by using a penalty function method of a logarithmic barrier function method, and finally solving an optimal solution of an unconstrained optimization problem by using a Newton method;

the non-linear problem is expressed by the following mathematical formula:

obj min.f(x)

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

wherein: f (x) is an objective function, which is a nonlinear function; h (x) ═ h₁(x),...,h_m(x)]^TFor non-linear equation constraints, g (x) ═ g₁(x),...,g_r(x)]^TWhen solving the problem (19) by using an interior point method, firstly, the inequality constraints are converted into the equality constraints, and simultaneously, barrier functions are constructed, so that firstly, relaxed variables l & gt 0, u & gt 0 and l ∈ R are introduced^r，u∈R^rConverting the inequality constraint of equation (19) into an equality constraint, and transforming the objective function into a barrier function, the following optimization problem a can be obtained:

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

wherein the perturbation factor u is greater than 0; when l is_iOr u_iWhen the boundary is approached, the function tends to be infinite, so that a minimal solution satisfying the barrier objective function cannot be found on the boundary, and only the condition that l is greater than 0 and u is satisfied>An optimal solution is possible to obtain when the value is 0; therefore, the optimization problem with inequality constraint is changed into the optimization problem A with equality constraint only through the transformation of the objective function, and therefore the optimization problem A can be directly solved by the Lagrange multiplier method.

The lagrangian function of the optimization model a is:

in the formula: y ═ y₁,...,y_m]^T，z＝[z₁,...,z_r]^T，w＝[w₁,...,w_r]^TAre all lagrange multipliers; the minimum value of the problem has the necessary condition that partial derivatives of the Lagrangian function to all variables and multipliers are 0, so that constrained optimization is converted into unconstrained optimization, and then a Newton method in the prior art can be used for solving;

s10: recording data such as nth group node voltage, branch power, power generation cost and the like;

s11: performing next round of load flow calculation, wherein t is t +1, and turning to S5;

compared with the prior art, the invention has the following advantages and beneficial effects:

1. the method has the advantages of high calculation speed and high accuracy, the obtained probability statistical information can comprehensively reflect the running condition of the electric power market, the uncertainty problem in the electric power market can be effectively processed, and the method has good engineering practical value;

2. through simulation and verification of the embodiment, the safe and economic operation of the distributed power supply connected to the power grid can be improved, the grid loss is reduced, the node voltage level is improved, and the economy and the practicability are effectively improved;

3. the node power price, the grid loss and the branch power fluctuation conditions of the wind power plant and photovoltaic power station hybrid system adopted by the invention are smaller than those of a single wind power plant system, the higher the photovoltaic capacity of the access system is, the lower the node power price is, and the functions of load flow calculation and control can be fully, effectively and quickly exerted.

Drawings

FIG. 1 is a flow chart of the steps of the present invention;

FIG. 2 is an expected value of node electricity prices when the photovoltaic power of the embodiment is accessed to a power grid;

FIG. 3 is a standard deviation of node electricity prices when the embodiment photovoltaic is connected into a power grid.

Detailed Description

A load flow calculation method for a distributed power supply to be connected into a power grid comprises the following steps:

s1: reading initial data of the power system;

s2: determining the sampling times N and the dimension s of an input random variable;

s3, generating an S × N-order sampling matrix according to the following 3 steps to form the first order in the point rowThe point (j 1, s, n 1) is as follows:

s3-1: the N-1 integer is expressed by 2-system number, namely formula (1)

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

Wherein a is_n∈Z_b，Z_b0,1, b-1, R is B^rThe maximum value of r is less than or equal to N;

s3-2: for N-1 ═ a_R-1a_R-2···a₂a₁Sequencing to obtain a sequenced sequence [ d ]₁d₂···d_n···d_R]^TIs (2)

Wherein,to generate a matrix, d is 0 ≦ d_nB-1 is less than or equal to; introducing a generator matrixIs to reset a₁a₂···a_n···a_R-1The position of each digit in; after the positions of the numbers are reset, the numbers of each dimension and other dimensions have the same size but different arrangement sequences, so that the uniformity of the result is ensured;

s3-3: through the calculation of step S3-2,can be expressed as a 2-ary form of equation (3):

finally, 2 is expressedConverting into 10-system number according to formula (2);

s4: initializing the sampling times: let n equal to 1;

s5: judging the N and the sampling times N, and if N is larger than N, directly outputting the probability statistical result of the variable; if N is less than or equal to N, turning to S6;

s6: determining a wind power and photovoltaic power generation output model and determining a load random model

S6-1: wind speed follows Weibull distribution, and active power P of wind power plant_wCan be expressed as formula (4):

in the formula: k and c are respectively the shape parameter and the scale parameter of the Weibull distribution,P_rrated power of the fan, v_r,v_ciRated wind speed and cut-in wind speed respectively;

wind power is processed into PQ nodes, the wind power factor in the load flow calculation is constant, and then reactive power is calculated according to the following formula (5):

in the formula:for power factor angle, for a grid-connected fan,is generally located in the fourth quadrant of the device,is a negative value;

s6-2: photovoltaic output stochastic model

Within a certain time period, the solar illumination intensity can be considered to obey the beta distribution, and then the output power P of the photovoltaic power station_pvIs expressed as formula (6):

in the formula: r_pv＝Aηγ_maxIn order to simulate the maximum output power, A is the total simulated area of the solar cell,η is a simulation of the total photoelectric conversion efficiency, γ_maxα are all shape parameters of beta distribution, which is the maximum illumination intensity in a period of time and is a Gamma function;

the method is the same as wind power, and a photovoltaic power station is also used as a PQ node in load flow calculation;

s6-3: load stochastic model

The load has time-varying property, and many related documents propose a method for predicting the regional load to obtain the probability distribution of the regional load; and as the load forecasting result of medium and long periods, the probability distribution rule of the load basically conforms to normal distribution. The mean and variance can be obtained from a large amount of historical statistical data; thus, the probability density functions of the active and reactive power of the load are respectively equations (7) and (8):

in the formula: mu.s_pIs the average value of the active power,_p ²is the variance of active power, mu_QIs the average value of the reactive power,_Q ²is the variance of the reactive power;

s7: determining a load flow calculation model

The method establishes an optimal power flow model with the minimum total cost of active power and reactive power generation of various energy sources as a target function, adjusts the output of the generator and the output of the reactive power source as much as possible to meet load requirements and system operation constraints, and searches feasible power output arrangement and power flow distribution state of a power grid with the minimum total power generation cost under the condition of ensuring the current load requirements and meeting the upper and lower limits of the voltage of each node and the transmission power limit of a transmission line;

s7-1: objective function

The power generation optimization model constructed by the invention is as follows:

c in formula (9)_Gpi、C_GqiAs a function of the active and reactive power generation costs of the unit i, C_gqj、C_gqjAs a function of the reactive power generation cost of the reactive power compensation device j, P_Gi(t)、Q_Gi(t) is the active and reactive power of the ith generator set in time period t, Q_gj(t) the reactive power output of the jth reactive power compensation device in the time period t; n is a radical of_g、N_qThe number of the generator nodes and the number of the reactive compensation equipment are calculated; the objective function is that the power generation cost of the system is minimum in each period;

s7-2: constraint of equality

The equality constraint is a node power flow balance constraint of each time interval:

in formulas (10) and (11): v_i、θ_iIs the node voltage and phase angle, θ_ij＝θ_i-θ_j；P_Di、Q_DiActive load and reactive load; g_ij、B_ijConductance and susceptance of a node admittance matrix;

s7-3: inequality constraint type (12)

In the formula,the upper and lower active output limits of the generator i are set;the upper and lower limit of reactive power output of the generator i;Q_giupper and lower reactive power output limits for the reactive power compensation equipment i;the node voltage amplitude upper and lower limits;continuously delivering a capacity limit (MVA) for line ij; n, N_bNode set and branch set;

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 the formula (13), R_i，upThe upward climbing speed of the ith thermal power generating unit is obtained; in the formula (14), R_i，downThe downward climbing speed of the ith thermal power generating unit is obtained;

s8: determining an optimal economic model

S8-1: electricity generation cost of thermal power plant

The active output of the thermal power coal-fired unit is charged by taking the coal consumption as a standard, and the unit i has an active output cost function C_GpiThe calculation was performed by equation (15). In the formula a_i、b_i、c_iThe coal consumption cost coefficient of the ith thermal power generating unit is obtained;

the reactive power price of the power generation side is divided into two parts: reactive capacity electricity prices and reactive power electricity quantity electricity prices. The reactive power electricity price mainly relates to the reactive power opportunity cost and the active loss cost of the generator, and the reactive power opportunity cost is used as the total reactive power generation cost on the generator side;

the reactive opportunity cost is the profit corresponding to the active power generating capacity lost by the generator due to the output of reactive power; if the output limit of the prime mover is ignored, and assumeThe cost of the reactive opportunity C_op(Q_Gi) Can be represented as formula (16);

substituting the formula (15) into the formula (16), performing Taylor expansion, and retainingNeglecting higher-order terms and then finishing to obtain a formula (17);

C_Gqi(Q_Gi) As a function of the reactive power contribution cost of the generator set i, S_Gi,maxRated apparent power, Q, of the generator set i_GiThe value k is the reactive output value of the generator set i, and the profit margin of the active power produced by the power plant is generally 5% -10%;

s8-2: electricity generation cost of hydraulic power plant

At present, the running cost of hydropower in China is generally 4-9 minutes/kilowatt hour, and the running cost of thermal power in China is about 0.09-0.19 yuan/kilowatt hour, the method adoptsThe water and electricity active power generation cost formula (15) is used for charging, and the value of the specific parameter is similar to that of a in thermal power_i，b_i，c_iThe difference is m times, m is the ratio of the thermal power running cost to the hydroelectric power running cost, a_i，b_i，c_iThe value is changed in a fine adjustment mode so as to distinguish the power generation cost of the same power station. The reactive power generation cost of the hydraulic power plant is similar to the reactive power output cost of the thermal power plant, and the charging mode of the formula (16) is adopted, wherein C_GpiTaking an active power generation cost function of a corresponding hydraulic power plant;

s8-3: generating cost of photovoltaic power station and wind power plant

At present, the grid-surfing electricity price of a photovoltaic power station and a wind power plant is still higher than that of the traditional energy, but with the reduction of the cost of photovoltaic equipment and wind power equipment and the enhancement of the national subsidy policy for new energy power generation, the further reduction of the grid-surfing electricity price of photovoltaic power generation and wind power generation can be expected. According to the method, new energy is preferentially called to the maximum extent as a criterion, the subsidized photovoltaic power generation cost and wind power generation cost are lower than the online power generation price of thermal power and hydropower, and the active cost function is selected in the same way as the active cost of the hydropower;

s8-4: cost of power generation of reactive power compensation equipment

The reactive cost of a capacitor, a reactor, a synchronous phase modulator and SVC is taken as a fixed cost expression (18):

wherein Y is the service life of the parallel capacitor, usually 15 years; p is the average usage, taken approximately as 2/3, C_fFor a fixed cost per unit capacity of capacitor, it is preferably 62500 yuan/MVar on average, from which f is calculated_q＝1.97；

S9: load flow calculation

Processing equality constraint in the optimization problem by using a Lagrange function method, thereby converting the optimization problem with equality constraint into an unconstrained optimization problem; processing inequality constraints by using a penalty function method of a logarithmic barrier function method, and finally solving an optimal solution of an unconstrained optimization problem by using a Newton method;

the non-linear problem is expressed by the following mathematical formula:

obj min.f(x)

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

wherein: f (x) is an objective function, which is a nonlinear function; h (x) ═ h₁(x),...,h_m(x)]^TFor non-linear equation constraints, g (x) ═ g₁(x),...,g_r(x)]^TWhen solving the problem (19) by the interior point method, the inequality constraints are firstly converted into the equality constraints, and simultaneously barrier functions are constructed, so that firstly, relaxed variables l & gt 0, u & gt 0 and l ∈ R are introduced^r，u∈R^rConverting the inequality constraint of equation (19) into an equality constraint, and transforming the objective function into a barrier function, the following optimization problem a can be obtained:

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

wherein the perturbation factor u > 0. When l is_iOr u_iWhen the boundary is approached, the function tends to be infinite, so that a minimal solution satisfying the barrier objective function cannot be found on the boundary, and only the condition that l is greater than 0 and u is satisfied>An optimal solution is possible to obtain when the value is 0; therefore, the optimization problem with inequality constraint is changed into the optimization problem A with equality constraint only through the transformation of the objective function, and therefore the optimization problem A can be directly solved by the Lagrange multiplier method.

The lagrangian function of the optimization model a is:

in the formula: y ═ y₁,...,y_m]^T，z＝[z₁,...,z_r]^T，w＝[w₁,...,w_r]^TAre all lagrange multipliers; the minimum value of the problem has the necessary condition that partial derivatives of the Lagrangian function to all variables and multipliers are 0, so that constrained optimization is converted into unconstrained optimization, and then a Newton method in the prior art can be used for solving;

s10: recording data such as nth group node voltage, branch power, power generation cost and the like;

s11: and (5) performing next round of load flow calculation, wherein t is t +1, and turning to S5.

The present invention may be simulated and verified using IEEE30 node instances. The technical scheme of the invention is clearly and completely described below by combining the embodiment drawings.

The optimal power flow probability statistical characteristic of the system after the system is connected into the wind power plant and the photovoltaic power station is directly analyzed by using a DN method with the sampling scale of 512 times.

TABLE 1 node Electricity price expected value comparison

In case 1, for 8 nodes 20, 61, 104, 123, 138, 171, 198 and 207, each node is accessed to a wind power plant with installed capacity of 60 MW; case 2, for 8 nodes 20, 61, 104, 123, 138, 171, 198 and 207, each node has access to one wind farm with 30MW installed capacity and one photovoltaic power plant with 30MW installed capacity.

And performing probability optimal power flow calculation by using the algorithm introduced herein to obtain the node electricity price of the distributed energy access node, as shown in table 1.

By comparing the calculation results in the two cases, the node electricity prices of the wind power plant and the photovoltaic hybrid system are lower than those of the wind power plant system only.

Table 2 shows expected grid loss values of the system under the two conditions, where the grid loss of the wind farm and the photovoltaic hybrid system is smaller than that of the wind farm system alone, and it can be seen that the wind farm and the photovoltaic hybrid system are more favorable for economic operation of the system.

Table 2 expected loss comparison

Distributed energy access mode	Loss of network/(MW)
		Case 1	232.622
Case 2	230.495

Table 3 shows the expected values and standard deviations for the arms 13-20 and arms 181-138 for the two cases. The results in the table show that the mean value and the standard deviation of the branch power when the system is independently accessed into the wind power plant are larger than those of the wind power plant and the photovoltaic hybrid system, the branch power fluctuation is larger, the probability of occurrence of heavy load and out-of-limit is also larger, and the line safety check is not facilitated.

TABLE 3 Branch Power comparison

In order to measure the influence of photovoltaic access systems with different capacities on the probability-optimal power flow, 8 photovoltaic power stations with the same installed capacity are accessed to the nodes 20, 61, 104, 123, 138, 171, 198 and 207, the total installed capacity of the 8 photovoltaic power stations is sequentially increased to 500MW (50 MW is increased in each step), and the probability-optimal power flow calculation is respectively carried out at each capacity.

Fig. 2 and 3 show the expected value and standard deviation of the node electricity prices when the photovoltaic power stations with different capacities are accessed into the system. The results of fig. 2 show that as the photovoltaic capacity of the access system is increased, the node electricity prices of the photovoltaic nodes show a trend of decreasing. This is because the photovoltaic output can replace a part of the output of the traditional thermal power generating unit, so that the node electricity price is reduced. The results in fig. 3 show that the randomness and uncertainty of the photovoltaic output can bring about the fluctuation of the node electricity price.

Claims (1)

1. A load flow calculation method for a distributed power supply to be connected into a power grid is characterized by comprising the following steps:

s1: reading initial data of the power system;

s2: determining the sampling times N and the dimension s of an input random variable;

s3, generating an S × N-order sampling matrix according to the following 3 steps to form the first order in the point rowAt a point (j-1, …, s; n-1, …)The method comprises the following steps:

s3-1: the N-1 integer is expressed by 2-system number, namely formula (1)

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

Wherein a is_n∈Z_b，Z_b(0, 1, …, b-1), and R is B^rThe maximum value of r is less than or equal to N;

s3-2: for N-1 ═ a_R-1a_R-2…a₂a₁Sequencing to obtain a sequenced sequence [ d ]₁d₂…d_n…d_R]^TIs (2)

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

Wherein,to generate a matrix, d is 0 ≦ d_nB-1 is less than or equal to; introducing a generator matrixIs to reset a₁a₂…a_n…a_R-1The position of each digit in; after the positions of the numbers are reset, the numbers of each dimension and other dimensions have the same size but different arrangement sequences, so that the uniformity of the result is ensured;

s3-3: through the calculation of step S3-2,can be expressed as a 2-ary form of equation (3):

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

finally, 2 is expressedConverting into 10-system number according to formula (2);

s4: initializing the sampling times: let n equal to 1;

s5: judging the N and the sampling times N, and if N is larger than N, directly outputting the probability statistical result of the variable; if N is less than or equal to N, turning to S6;

s6: determining a wind power and photovoltaic power generation output model and determining a load random model;

s6-1: wind speed follows Weibull distribution, and active power P of wind power plant_wCan be expressed as(4)：

$f\left(P_w\right)=\frac{k}{k_{1}c}\left(\frac{P_w-k_2}{k_{1}c}\right)\exp \&lsqb;-{\left(\frac{P_w-k_2}{k_{1}c}\right)}^k\&rsqb;---\left(4\right)$

In the formula: k and c are respectively the shape parameter and the scale parameter of the Weibull distribution,k₂＝-k₁v_ci，P_rrated power of the fan, v_r,v_ciRated wind speed and cut-in wind speed respectively;

wind power is processed into PQ nodes, the wind power factor in the load flow calculation is constant, and then reactive power is calculated according to the following formula (5):

in the formula:for power factor angle, for a grid-connected fan,is positioned in the fourth quadrant of the device,is a negative value;

s6-2: photovoltaic output stochastic model

Within a certain time period, the solar illumination intensity can be considered to obey the beta distribution, and then the output power P of the photovoltaic power station_pvIs expressed as formula (6):

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

in the formula: r_pv＝Aηγ_maxTo simulate maximum output power, A is the simulated total area of the solar cell, η is the simulated total photoelectric conversion efficiency, γ_maxα are all shape parameters of beta distribution, which is the maximum illumination intensity in a period of time and is a Gamma function;

the method is the same as wind power, and a photovoltaic power station is also used as a PQ node in load flow calculation;

s6-3: load stochastic model

As the load forecasting result of medium and long periods, the probability distribution rule of the load basically conforms to normal distribution; the mean and variance can be obtained from a large amount of historical statistical data; thus, the probability density functions of the active and reactive power of the load are respectively equations (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 the formula: mu.s_pIs the average value of the active power,_p ²is the variance of active power, mu_QIs the average value of the reactive power,_Q ²is the variance of the reactive power;

s7: determining a load flow calculation model

S7-1: objective function

The power generation optimization model is constructed as follows:

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

c in formula (9)_Gpi、C_GqiAs a function of the active and reactive power generation costs of the unit i, C_gqjAs a function of the reactive power generation cost of the reactive power compensation device j, P_Gi(t)、Q_Gi(t) is the active and reactive power of the ith generator set in time period t, Q_gj(t) the reactive power output of the jth reactive power compensation device in the time period t; n is a radical of_g、N_qThe number of the generator nodes and the number of the reactive compensation equipment are calculated; the objective function enables the power generation cost of the system in each time period to be minimum;

s7-2: constraint of equality

The equality constraint is a node power flow balance constraint of each time interval:

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

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

in formulas (10) and (11): v_i、θ_iIs the node voltage and phase angle, θ_ij＝θ_i-θ_j；P_Di、Q_DiActive load and reactive load; g_ij、B_ijConductance and susceptance of a node admittance matrix;

s7-3: inequality constraint type (12)

$\left\{\begin{array}{c}\underset{\&OverBar;}{P_{Gi}}\&le;P_{Gi}\&le;\stackrel{\&OverBar;}{P_{Gi}}\\ \underset{\&OverBar;}{Q_{Gi}}\&le;Q_{Gi}\&le;\stackrel{\&OverBar;}{Q_{Gi}}\\ \underset{\&OverBar;}{Q_{gi}}\&le;Q_{gi}\&le;\stackrel{\&OverBar;}{Q_{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 the formula, P _Gi the upper and lower active output limits of the generator i are set; Q _Gi the upper and lower limit of reactive power output of the generator i;Q_giupper and lower reactive power output limits for the reactive power compensation equipment i; V _i the node voltage amplitude upper and lower limits;continuously delivering a capacity limit (MVA) for line i; n, N_bNode set and branch set;

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 the formula (13), R_i，upThe upward climbing speed of the ith thermal power generating unit is obtained; in the formula (14), R_i，downThe downward climbing speed of the ith thermal power generating unit is obtained;

s8: determining an optimal economic model

S8-1: electricity generation cost of thermal power plant

The active output of the thermal power coal-fired unit is charged by taking the coal consumption as a standard, and the unit i has an active output cost function C_GpiCalculating by the formula (15); in the formula a_i、b_i、c_iThe coal consumption cost coefficient of the ith thermal power generating unit is obtained;

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

the reactive opportunity cost is the profit corresponding to the active power generating capacity lost by the generator due to the output of reactive power; if the output limit of the prime mover is ignored, and assumeThe cost of the reactive opportunity C_op(Q_Gi) Can be represented as formula (16);

$C_{Gqi}\left(Q_{Gi}\right)=C_{op}\left(Q_{Gi}\right)=k\&lsqb;C_{Gpi}\left(S_{Gi,\max }\right)-C_{Gpi}\left(\sqrt{S_{Gi,\max }^2-Q_{Gi}^2}\right)\&rsqb;---\left(16\right)$

substituting the formula (15) into the formula (16), performing Taylor expansion, and retainingNeglecting higher order terms and then arranging to obtain formula (17)

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

C_Gqi(Q_Gi) As a function of the reactive power contribution cost of the generator set i, S_Gi,maxRated apparent power, Q, of the generator set i_GiThe value k is the reactive output value of the generator set i, and the profit margin of the active power produced by the power plant is generally 5% -10%;

s8-2: electricity generation cost of hydraulic power plant

The method adopts a hydropower active power generation cost formula (15) for charging, and the value of the specific parameter is similar to that of a in thermal power_i，b_i，c_iThe difference is m times, m is the ratio of the thermal power running cost to the hydroelectric power running cost, a_i，b_i，c_iThe value is changed in a fine adjustment way so as to distinguish the power generation cost of the same power station; the reactive power generation cost of the hydraulic power plant is similar to the reactive power output cost of the thermal power plant, and the charging mode of the formula (16) is adopted, wherein C_GpiTaking an active power generation cost function of a corresponding hydraulic power plant;

s8-3: generating cost of photovoltaic power station and wind power plant

Taking priority to calling new energy to the maximum extent as a criterion, ensuring that the subsidized photovoltaic power generation cost and wind power generation cost are lower than the online power generation price of thermal power and hydropower, and selecting an active cost function in the same way as the active cost of the hydropower;

s8-4: cost of power generation of reactive power compensation equipment

The reactive cost of a capacitor, a reactor, a synchronous phase modulator and SVC is taken as a fixed cost expression (18):

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

wherein Y is the service life of the parallel capacitor, and is 15 years; p is the average usage, taken approximately as 2/3, C_fFor a fixed cost per unit capacity of capacitor, it is preferably 62500 yuan/MVar on average, from which f is calculated_q＝1.97；

S9: load flow calculation

Processing equality constraint in the optimization problem by using a Lagrange function method, thereby converting the optimization problem with equality constraint into an unconstrained optimization problem; processing inequality constraints by using a penalty function method of a logarithmic barrier function method, and finally solving an optimal solution of an unconstrained optimization problem by using a Newton method;

the non-linear problem is expressed by the following mathematical formula:

obj min.f(x)

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

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

wherein: f (x) is an objective function, which is a nonlinear function; h (x) ═ h₁(x),...,h_m(x)]^TFor non-linear equation constraints, g (x) ═ g₁(x),...,g_r(x)]^TFor nonlinear inequality constraint, assuming that the above model has k variables, m equality constraints and R inequality constraints, when solving problem (19) by interior point method, firstly converting inequality constraints into equality constraints and simultaneously constructing barrier function, firstly introducing relaxed variables l > 0, u > 0 and l ∈ R^r，u∈R^rConverting the inequality constraint of equation (19) into an equality constraint, and transforming the objective function into a barrier function, the following optimization problem a can be obtained:

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

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

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

g(x)-l-g＝0

wherein the perturbation factor u is greater than 0; when l is_iOr u_iWhen the boundary is approached, the function tends to be infinite, so that a minimal solution satisfying the barrier objective function cannot be found on the boundary, and only the condition that l is greater than 0 and u is satisfied>An optimal solution is possible to obtain when the value is 0; therefore, the optimization problem with inequality limitation is changed into the optimization problem A with equality constraint limitation only through the transformation of the objective function, and therefore the Lagrange multiplier method can be directly used for solving;

the lagrangian function of the optimization model a is:

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

in the formula: y ═ y₁,...,y_m]^T，z＝[z₁,...,z_r]^T，w＝[w₁,...,w_r]^TAre all lagrange multipliers; the minimum value of the problem has the necessary condition that the partial derivatives of the Lagrangian function to all variables and multipliers are 0, so that the constrained optimization is converted into the non-constrained optimizationConstraint optimization, which can then be solved using the newton method in the prior art;

s10: recording data such as nth group node voltage, branch power, power generation cost and the like;

s11: and (5) performing next round of load flow calculation, wherein t is t +1, and turning to S5.