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

Abstract

A power flow calculation method for distributed power sources connected to the power grid, comprising the following steps: S1: read the initial data of the power system; S2: determine the number of samples N and the dimension s of the input random variable; S3: generate s×N order samples Matrix; S4: Initialize the number of samples: let n=1; S5: Determine the size of n and the number of samples N, if n>N, directly output the probability statistics of variables; if n≤N, go to S6; S6: Determine the wind power and photovoltaic power generation output model to determine the load random model; S7: determine the power flow calculation model; S8: determine the optimal economic model; S9: power flow calculation S10: record the nth group of node voltage, branch power and power generation costs and other data; S11: Carry out the next round of power flow calculation, t=t+1, go to S5. The invention can better estimate the probability distribution of the output random variable, and can effectively deal with the uncertainty problem in the electric power market. The invention saves manpower and material resources for debugging, and reduces the production cost.

Description

A Power Flow Calculation Method for Distributed Power Supply Connected to Power Grid

technical field

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

Background technique

Wind energy and solar energy are green and clean energies. Vigorously developing wind power and photovoltaics will help reduce fossil fuel consumption and reduce carbon emissions. But because of its intermittent and random characteristics, higher requirements are put forward for the operation control of power system. In recent years, my country's wind power and photovoltaic power generation have developed rapidly, the installed capacity has increased rapidly, and it is difficult to accommodate them. In the planning, design and operation control of the power system, the characteristics of wind farms and photovoltaic power stations should be considered more comprehensively, and their fluctuation laws can be mastered. Safety and economy are of great significance.

The output of wind power and photovoltaics is greatly affected by natural weather conditions. When large-scale wind energy and photovoltaics are connected to the system, the fluctuation of its output will be different from the previous output scheduling of thermal power and hydropower. How to dispatch new energy first under the condition of satisfying the balance of system power supply and demand, and let thermal power units bear the basic load under the condition of considering the volatility of new energy, with little change in output at different times; There are large fluctuations; considering the optimization of active power output and reactive power output at the same time and making the total power generation cost of the system the lowest, these all put forward higher requirements for the modeling of optimal power flow. The optimal power flow of power grid has high practical value. For the first time, it combines economy and safety, active power and reactive power optimization together almost perfectly. Run the scheduler's request.

Due to the uncertainty of new energy output, the calculation of probabilistic optimal power flow becomes more complicated. At present, there are relevant literatures on the optimal power flow with new energy. The literature "Stochastic Optimal Power Flow Method Considering Injected Power Distribution" considers the uncertainty of wind turbine output and establishes a chance-constrained optimal power flow model. The document "Analysis of Wind Power Access Capability Based on Probabilistic Optimal Power Flow" uses the particle swarm optimization algorithm of stochastic technology to solve the probability optimal power flow model, and evaluates the feasibility and effectiveness of wind power access capacity. However, the above studies generally only consider wind farms, and rarely study the influence of wind farms and photovoltaic power stations connected to the system at the same time on the optimal power flow. The output of wind farms and photovoltaic power plants is random and has different output characteristics, which increases the uncertainties in the electricity market.

At present, the optimal power flow calculation methods considering the influence of randomness mainly include Monte Carlo method, cumulant method, point estimation method, ant colony algorithm, etc. The Monte Carlo method can well study the influence of random factors on the optimal power flow of the system, but this method requires thousands of simulations of different operating states of the system to obtain reasonable results, which takes a long time to calculate and takes up a lot of memory. On the premise that the input random variables are independent of each other or satisfy a linear relationship, the cumulant method uses Gram-Charlier expansion series, Cornish-Fisher expansion series, etc. to fit, so as to obtain the probability density function of the output random variable, which improves the calculation efficiency . Although the point estimation method has a faster calculation speed, the error of the high-order moment of the output random variable is relatively large. The ant colony algorithm has a large amount of calculation, which directly affects the calculation speed.

Contents of the invention

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

S1: read the initial data of the power system;

S2: Determine the number of sampling N and the dimension s of the input random variable;

S3: According to the following 3 steps, generate the s×N order sampling matrix to form the points (j=1,...,s; n=1,...) The steps are as follows:

S3-1: Express the N-1th integer in binary, that is, formula (1)

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

Where a _n ∈ Z _b , Z _b ={0,1,…,b-1}, R is the maximum value of r satisfying b ^r ≤ N;

S3-2: Sort N-1=a _R-1 a _R-2 ...a ₂ a ₁ , and obtain the sorted sequence [d ₁ d ₂ ...d _n ...d _R ] ^T as formula (2)

in, is the generator matrix, 0≤d _n ≤b-1; introduce the generator matrix It is to reset the position of each number in a ₁ a ₂ ···a _n ···a _R-1 ; after the position of the number is reset, the size of each dimension is the same as that of other dimensions, but the arrangement order is different, so that The uniformity of results is guaranteed;

S3-3: After the calculation of step S3-2, It can be expressed as the binary form of formula (3):

Finally, the binary representation According to the formula (2), it can be converted into a decimal number;

S4: Initialize the sampling times: make n=1;

S5: Determine the size of n and the number of samples N, if n>N, directly output the probability statistics of variables; if n≤N, go to S6;

S6: Determine the output model of wind power and photovoltaic power generation, and determine the random load model

S6-1: The wind speed obeys the Weibull distribution, and the probability density function of the active power P _w of the wind farm can be expressed as formula (4):

In the formula: k and c are the shape parameter and scale parameter of the Weibull distribution respectively, P _r is the rated power of the fan, v _r and v _ci are the rated wind speed and cut-in wind speed respectively;

The wind power processing is a PQ node, so that the power factor of the wind turbine in the power flow calculation is constant, and the reactive power is calculated according to formula (5):

In the formula: is the power factor angle, for grid-connected wind turbines, Generally located in the fourth quadrant, is a negative value.

S6-2: Photovoltaic output stochastic model

In a certain period of time, the intensity of sunlight can be considered to obey the Beta distribution, then the probability density function of the output power _Ppv of the photovoltaic power station is expressed as formula (6):

In the formula: R _pv = Aηγ _max is the maximum output power of the simulation, A is the total area of solar cell simulation, η is the total photoelectric conversion efficiency of the simulation, γ _max is the maximum light intensity in a period of time, Γ is the Gamma function, α, β Both are shape parameters of the Beta distribution;

Like wind power, photovoltaic power plants are also used as PQ nodes in power flow calculations;

S6-3: Load Stochastic Model

The load is time-varying, and many relevant literatures have proposed methods to predict the regional load to obtain its probability distribution; as the medium and long-term load forecast results, the probability distribution of the load basically conforms to the normal distribution; its 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 formulas (7) and (8) respectively:

In the formula: μ _p is the mean value of active power, δ _p ² is the variance of active power, μ _Q is the mean value of reactive power, and δ _Q ² is the variance of reactive power;

S7: Determine the power flow calculation model

The present invention establishes the optimal power flow model with the minimum total cost of active power and reactive power generation of various energy sources as the objective function, and adjusts the output of generators and reactive power sources as much as possible to meet the load requirements and system operation constraints. Under the upper and lower limits of the voltage of each node and the transmission power limit of the transmission line, find the feasible arrangement of the generator output and the power flow distribution state of the power grid with the minimum total power generation cost;

S7-1: Objective function

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

In formula (9), C _Gpi and C _Gqi are the active and reactive power generation cost functions of unit i, C _gqj and C _gqj are the reactive power generation cost functions of reactive compensation device j, P _Gi (t), Q _Gi (t ) is the active output and reactive output of the i-th generator set in the period t, Q _gj (t) is the reactive output of the j-th reactive compensation device in the period t; N _g and N _q are the number of generator nodes and the number of reactive power compensation equipment; the objective function is to minimize the power generation cost of the system in each period;

S7-2: Equality constraints

The equation constraint is the node power flow balance constraint in each time period:

In formula (10) and formula (11): V _i , θ _i are node voltage and phase angle, θ _ij = θ _i - θ _j ; P _Di , Q _Di are active load and reactive load; G _ij , B _ij is the conductance and susceptance of the nodal admittance matrix;

S7-3: Inequality constraints (12)

In the formula, The upper and lower limits of the active output of generator i; The upper and lower limits of reactive power output of generator i; Q _gi is the upper and lower limit of reactive power output of reactive power compensation equipment i; is the upper and lower limits of the node voltage amplitude; is the continuous transmission capacity limit (MVA) of the line ij; N, N _b are node sets and branch sets;

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 is the upward climbing rate of the i-th thermal power unit; in formula (14), R _i,down is the downward climbing rate of the i-th thermal power unit;

S8: Determine the optimal economic model

S8-1: Power generation costs of thermal power plants

The active output of thermal power coal-fired units is billed based on coal consumption, and the active output cost function C _Gpi of unit i is calculated according to formula (15). In the formula, a _i , b _i , and c _i are the coal consumption coefficients of the i-th thermal power unit;

The reactive power price on the power generation side is divided into two parts: reactive capacity power price and reactive power power price. The reactive electricity price mainly involves the reactive opportunity cost and active power loss expense of the generator, and the present invention uses the reactive opportunity cost as the total reactive power generation expense on the generator side;

The reactive opportunity cost is the profit corresponding to the active power generation capacity lost by the generator due to the output of reactive power; if the output limit of the prime mover is ignored, and assuming The reactive opportunity cost C _op (Q _Gi ) can be expressed as formula (16);

Substitute Equation (15) into Equation (16), and carry out Taylor expansion, retaining to , after ignoring the high-order terms and sorting out the formula (17);

C _Gqi (Q _Gi ) is the reactive output cost function of generator i, S _Gi,max is the rated apparent power of generator i, Q _Gi is the reactive output value of generator i, and k is the active power produced by the power plant The profit margin is generally 5%-10%;

S8-2: Electricity generation costs of hydropower plants

At present, the operating cost of hydropower in China is generally 4 to 9 cents/kWh, while the operating cost of thermal power in China is about 0.09-0.19 yuan/kWh. The values of the parameters are approximately m times different from a _i , b _i , and c _i in thermal power _{. m} is the ratio of the operating cost of thermal power to the electricity price of hydropower _. The power generation cost of the power station; the reactive power generation cost of the hydropower plant also adopts the reactive power output cost similar to that of the thermal power plant, and according to the billing method of formula (16), wherein C _Gpi takes the active power generation cost function of the corresponding hydropower plant;

S8-3: Power generation costs of photovoltaic power plants and wind farms

At present, the on-grid electricity price of photovoltaic power plants and wind farms is still higher than that of traditional energy sources, but with the reduction of the cost of photovoltaic equipment and wind power equipment, and the strengthening of the national subsidy policy for new energy power generation, the further reduction of on-grid electricity prices for photovoltaic power generation and wind power generation is can be expected. In the present invention, the maximum priority is given to the use of new energy as a criterion, so that the subsidized photovoltaic power generation cost and wind power generation cost are lower than the grid-connected power generation prices of thermal power and hydropower, and the selection of the active cost function is the same as that of hydropower;

S8-4: Power generation cost of reactive power compensation equipment

Taking capacitors, reactors, synchronous condensers, and SVC reactive costs as fixed cost expressions (18):

Among them, Y is the service life of parallel capacitors, which is usually 15 years; p is the average utilization rate, which is approximately 2/3; C _f is the fixed cost of capacitor unit capacity, which can be 62500 yuan/MVar on average, and is calculated from this data out f _q = 1.97;

S9: Power Flow Calculation

Use the Lagrange function method to deal with the equality constraints in the optimization problem, so as to transform the optimization problem with the equality constraint into an unconstrained optimization problem; use the penalty function method of the logarithmic barrier function method to deal with the inequality constraints, and finally use Newton's method to solve the optimal solution of unconstrained optimization problems;

Express the nonlinear problem with the following mathematical formula:

obj min.f(x)

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

Among them: min.f(x) is the objective function, which is a nonlinear function; h(x)=[h ₁ (x),...,h _m (x)] ^T is the nonlinear equality constraint, g (x)=[g ₁ (x),...,g _r (x)] ^T is a nonlinear inequality constraint. Assume that there are k variables in the above model, m equality constraints, and r inequality constraints. When using the interior point method to solve problem (19), the inequality constraints are converted into equality constraints first, and the barrier function is constructed at the same time; for this purpose, the slack variables l>0, u>0, l∈R ^r , u∈R ^r , Transforming the inequality constraints of formula (19) into equality constraints, and transforming the objective function into an obstacle function, the following optimization problem A can be obtained:

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

Among them, the disturbance factor u>0; when l _i or u _i is close to the boundary, the above function tends to infinity, so the minimal solution satisfying the above obstacle objective function cannot be found on the boundary, only when l>0, u> In this way, the optimization problem containing inequality constraints is transformed into an optimization problem A containing only equality constraints through the transformation of the objective function, so the Lagrange multiplier method can be used directly solve.

The Lagrangian function of optimizing model A is:

In the formula: y=[y ₁ ,...,y _m ] ^T , z=[z ₁ ,...,z _r ] ^T , w=[w ₁ ,...,w _r ] ^T are pull Grangian multipliers; the necessary condition for the existence of the minimum value of this problem is that the partial derivatives of the Lagrangian function to all variables and multipliers are 0, so that the constrained optimization can be transformed into unconstrained optimization, and then the existing Newton's method solution in technology;

S10: Record data such as node voltage, branch power and power generation cost of the nth group;

S11: Carry out the next round of power flow calculation, t=t+1, go to S5;

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

1. The present invention has fast calculation speed and high accuracy, and the obtained probability and statistical information can fully reflect the operation status of the electric power market, can effectively deal with uncertainties in the electric power market, and has good engineering practical value;

2. Through the simulation and verification of the embodiment, the present invention can improve the safe and economical operation of the distributed power supply connected to the power grid, reduce the network loss, improve the node voltage level, and effectively improve the economy and practicability;

3. The node power price, network loss and branch power fluctuations of the wind farm and photovoltaic power plant hybrid system adopted in the present invention are smaller than those of the independent wind farm system. The greater the photovoltaic capacity connected to the system, the lower the node power price, which can be more comprehensive , Effectively and quickly give full play to the functions of power flow calculation and control, the present invention saves manpower and material resources for commissioning, reduces production costs, and has certain economic benefits.

Description of drawings

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

Fig. 2 is the expected value of the node electricity price when the photovoltaic of the embodiment is connected to the grid;

Fig. 3 is the standard deviation of the node electricity price when the photovoltaic of the embodiment is connected to the grid.

detailed description

A power flow calculation method for connecting distributed power sources to a power grid, comprising the following steps:

S1: read the initial data of the power system;

S2: Determine the number of sampling N and the dimension s of the input random variable;

S3: According to the following 3 steps, generate the s×N order sampling matrix to form the points (j=1,...,s; n=1,...) the steps are as follows:

S3-1: Express the N-1th integer in binary, that is, formula (1)

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

Where a _n ∈ Z _b , Z _b = {0,1,...,b-1}, R is the maximum value of r satisfying b ^r ≤ N;

S3-2: Sort N-1=a _R-1 a _R-2 ···a ₂ a ₁ , and obtain the sorted sequence [d ₁ d ₂ ···d _n ···d _R ] ^T as Formula (2)

in, is the generator matrix, 0≤d _n ≤b-1; introduce the generator matrix It is to reset the position of each number in a ₁ a ₂ ···a _n ···a _R-1 ; after the position of the number is reset, the size of each dimension is the same as that of other dimensions, but the arrangement order is different, so that The uniformity of results is guaranteed;

S3-3: After the calculation of step S3-2, It can be expressed as the binary form of formula (3):

Finally, the binary representation According to the formula (2), it can be converted into a decimal number;

S4: Initialize the sampling times: make n=1;

S5: Determine the size of n and the number of samples N, if n>N, directly output the probability statistics of variables; if n≤N, go to S6;

S6: Determine the output model of wind power and photovoltaic power generation, and determine the random load model

S6-1: The wind speed obeys the Weibull distribution, and the probability density function of the active power P _w of the wind farm can be expressed as formula (4):

In the formula: k and c are the shape parameter and scale parameter of the Weibull distribution respectively, P _r is the rated power of the fan, v _r and v _ci are the rated wind speed and cut-in wind speed respectively;

The wind power processing is a PQ node, so that the power factor of the wind turbine in the power flow calculation is constant, and the reactive power is calculated according to formula (5):

In the formula: is the power factor angle, for grid-connected wind turbines, Generally located in the fourth quadrant, is a negative value;

S6-2: Photovoltaic output stochastic model

In a certain period of time, the intensity of sunlight can be considered to obey the Beta distribution, then the probability density function of the output power _Ppv of the photovoltaic power station is expressed as formula (6):

In the formula: R _pv = Aηγ _max is the maximum output power of the simulation, A is the total area of solar cell simulation, η is the total photoelectric conversion efficiency of the simulation, γ _max is the maximum light intensity in a period of time, Γ is the Gamma function, α, β Both are shape parameters of the Beta distribution;

Like wind power, photovoltaic power plants are also used as PQ nodes in power flow calculations;

S6-3: Load Stochastic Model

The load is time-varying, and many relevant literatures have proposed the method of forecasting the regional load to obtain its probability distribution; as the medium and long-term load forecast results, the probability distribution of the load basically conforms to the normal distribution. Both its 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 formulas (7) and (8) respectively:

In the formula: μ _p is the mean value of active power, δ _p ² is the variance of active power, μ _Q is the mean value of reactive power, and δ _Q ² is the variance of reactive power;

S7: Determine the power flow calculation model

The present invention establishes the optimal power flow model with the minimum total cost of active power and reactive power generation of various energy sources as the objective function, and adjusts the output of generators and reactive power sources as much as possible to meet the load requirements and system operation constraints. Under the upper and lower limits of the voltage of each node and the transmission power limit of the transmission line, find the feasible arrangement of the generator output and the power flow distribution state of the power grid with the minimum total power generation cost;

S7-1: Objective function

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

In formula (9), C _Gpi and C _Gqi are the active and reactive power generation cost functions of unit i, C _gqj and C _gqj are the reactive power generation cost functions of reactive compensation device j, P _Gi (t), Q _Gi (t ) is the active output and reactive output of the i-th generator set in the period t, Q _gj (t) is the reactive output of the j-th reactive compensation device in the period t; N _g and N _q are the number of generator nodes and the number of reactive power compensation equipment; the objective function is to minimize the power generation cost of the system in each period;

S7-2: Equality constraints

The equation constraint is the node power flow balance constraint in each time period:

In formula (10) and formula (11): V _i , θ _i are node voltage and phase angle, θ _ij = θ _i - θ _j ; P _Di , Q _Di are active load and reactive load; G _ij , B _ij is the conductance and susceptance of the nodal admittance matrix;

S7-3: Inequality constraints (12)

In the formula, The upper and lower limits of the active output of generator i; The upper and lower limits of reactive power output of generator i; Q _gi is the upper and lower limit of reactive power output of reactive power compensation equipment i; is the upper and lower limits of the node voltage amplitude; is the continuous transmission capacity limit (MVA) of the line ij; N, N _b are node sets and branch sets;

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 is the upward climbing rate of the i-th thermal power unit; in formula (14), R _i,down is the downward climbing rate of the i-th thermal power unit;

S8: Determine the optimal economic model

S8-1: Power generation costs of thermal power plants

The active output of thermal power coal-fired units is billed based on coal consumption, and the active output cost function C _Gpi of unit i is calculated according to formula (15). In the formula, a _i , b _i , and c _i are the coal consumption coefficients of the i-th thermal power unit;

The reactive power price on the power generation side is divided into two parts: reactive capacity power price and reactive power power price. The reactive electricity price mainly involves the reactive opportunity cost and active power loss expense of the generator, and the present invention uses the reactive opportunity cost as the total reactive power generation expense on the generator side;

The reactive opportunity cost is the profit corresponding to the active power generation capacity lost by the generator due to the output of reactive power; if the output limit of the prime mover is ignored, and assuming The reactive opportunity cost C _op (Q _Gi ) can be expressed as formula (16);

Substitute Equation (15) into Equation (16), and carry out Taylor expansion, retaining to , after ignoring the high-order terms and sorting out the formula (17);

C _Gqi (Q _Gi ) is the reactive output cost function of generator i, S _Gi,max is the rated apparent power of generator i, Q _Gi is the reactive output value of generator i, and k is the active power produced by the power plant The profit margin is generally 5%-10%;

S8-2: Electricity generation costs of hydropower plants

At present, the operating cost of hydropower in China is generally 4 to 9 cents/kWh, while the operating cost of thermal power in China is about 0.09-0.19 yuan/kWh. The values of the parameters are approximately m times different from a _i , b _i , and c _i in thermal power _{. m} is the ratio of the operating cost of thermal power to the electricity price of hydropower _. power generation costs of the power station. The reactive power generation cost of the hydropower plant is also similar to the reactive power output cost of the thermal power plant, and is billed according to the formula (16), where C _Gpi takes the active power generation cost function of the corresponding hydropower plant;

S8-3: Power generation costs of photovoltaic power plants and wind farms

At present, the on-grid electricity price of photovoltaic power plants and wind farms is still higher than that of traditional energy sources, but with the reduction of the cost of photovoltaic equipment and wind power equipment, and the strengthening of the national subsidy policy for new energy power generation, the further reduction of on-grid electricity prices for photovoltaic power generation and wind power generation is can be expected. In the present invention, the maximum priority is given to the use of new energy as a criterion, so that the subsidized photovoltaic power generation cost and wind power generation cost are lower than the grid-connected power generation prices of thermal power and hydropower, and the selection of the active cost function is the same as that of hydropower;

S8-4: Power generation cost of reactive power compensation equipment

Taking capacitors, reactors, synchronous condensers, and SVC reactive costs as fixed cost expressions (18):

Among them, Y is the service life of parallel capacitors, which is usually 15 years; p is the average utilization rate, which is approximately 2/3; C _f is the fixed cost of capacitor unit capacity, which can be 62500 yuan/MVar on average, and is calculated from this data out f _q = 1.97;

S9: Power Flow Calculation

Use the Lagrange function method to deal with the equality constraints in the optimization problem, so as to transform the optimization problem with the equality constraint into an unconstrained optimization problem; use the penalty function method of the logarithmic barrier function method to deal with the inequality constraints, and finally use Newton's method to solve the optimal solution of unconstrained optimization problems;

Express the nonlinear problem with the following mathematical formula:

obj min.f(x)

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

Among them: min.f(x) is the objective function, which is a nonlinear function; h(x)=[h ₁ (x),...,h _m (x)] ^T is the nonlinear equality constraint, g (x)=[g ₁ (x),...,g _r (x)] ^T is a nonlinear inequality constraint. Assume that there are k variables in the above model, m equality constraints, and r inequality constraints. When using the interior point method to solve the problem (19), the inequality constraints are converted into equality constraints first, and the barrier function is constructed at the same time. To this end, first introduce the slack variables l>0, u>0, l∈R ^r , u∈R ^r , transform the inequality constraints in equation (19) into equality constraints, and transform the objective function into an obstacle function, the following can be obtained Optimization problem A:

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

Wherein the disturbance factor u>0. When l _i or u _i is close to the boundary, the above function tends to infinity, so the minimal solution satisfying the above obstacle objective function cannot be found on the boundary, and the optimal solution can only be obtained when l>0, u>0 In this way, through the transformation of the objective function, the optimization problem containing inequality constraints is transformed into an optimization problem A containing only equality constraints, so it can be solved directly by the Lagrange multiplier method.

The Lagrangian function of optimizing model A is:

In the formula: y=[y ₁ ,...,y _m ] ^T , z=[z ₁ ,...,z _r ] ^T , w=[w ₁ ,...,w _r ] ^T are pull Grangian multipliers; the necessary condition for the existence of the minimum value of this problem is that the partial derivatives of the Lagrangian function to all variables and multipliers are 0, so that the constrained optimization can be transformed into unconstrained optimization, and then the existing Newton's method solution in technology;

S10: Record data such as node voltage, branch power and power generation cost of the nth group;

S11: Carry out the next round of power flow calculation, t=t+1, go to S5.

The present invention can adopt the IEEE30 node instance for simulation and verification. The technical solutions of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the embodiments.

This calculation example directly uses the DN method with a sampling scale of 512 times to analyze the probability and statistical characteristics of the optimal power flow after the system is connected to the wind farm and photovoltaic power station.

Table 1 Comparison of node electricity price expectations

In order to influence the randomness of distributed energy output on power-saving electricity price, the following is divided into two cases for discussion: Case 1, for 8 nodes 20, 61, 104, 123, 138, 171, 198 and 207, each node Both are connected to a wind farm with an installed capacity of 60MW; Case 2, for eight nodes 20, 61, 104, 123, 138, 171, 198 and 207, each node is connected to a wind farm with an installed capacity of 30MW and A photovoltaic power station with an installed capacity of 30MW.

Using the algorithm introduced in this paper to calculate the probabilistic optimal power flow, the node electricity price of the distributed energy access node is obtained, as shown in Table 1.

By comparing the calculation results in the two cases, it can be seen that the node electricity price of the wind farm and the photovoltaic hybrid system is lower than that of the wind farm system only.

Table 2 shows the expected value of the grid loss of the system in the two cases. The grid loss of the wind farm and the photovoltaic hybrid system is smaller than that of the wind farm system alone. It can be seen that the wind farm and the photovoltaic hybrid system are more conducive to the economic operation of the system.

Table 2 Comparison of network loss expectation

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

Table 3 shows the expected value and standard deviation of branches 13-20 and branches 181-138 in two cases. It can be seen from the results in the table that the mean value and standard deviation of the branch power when the system is connected to the wind farm alone are larger than those of the wind farm and the photovoltaic hybrid system, the branch power fluctuates more, and the probability of overloading and exceeding the limit is also higher. Large, which is not conducive to line safety check.

Table 3 Branch power comparison

In order to measure the impact of photovoltaic access systems with different capacities on the probabilistic optimal power flow, eight photovoltaic power stations with equal installed capacity are connected to nodes 20, 61, 104, 123, 138, 171, 198 and 207. The installed capacity increases sequentially from 50MW to 500MW (50MW per step), and the probabilistic optimal power flow calculation is performed for each capacity.

Figure 2 and Figure 3 show the expected value and standard deviation of node power prices when photovoltaic power plants with different capacities are connected to the system. The results in Figure 2 show that with the continuous increase of photovoltaic capacity of the access system, the node electricity price of photovoltaic nodes shows a downward trend. This is because photovoltaic output can replace part of the output of traditional thermal power units, thereby reducing the node electricity price. The results in Figure 3 show that the randomness and uncertainty of photovoltaic output will bring fluctuations in node electricity prices.

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)

${\[d_1{d}_2...d_n...d_R\]}^T=C_{N\×N}^i{\[a_1{a}_2...a_n...a_{R-1}\]}^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 \[-{\left(\frac{P_w-k_2}{k_{1}c}\right)}^k\]---\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{\Γ}(\mathrm{\α}+\mathrm{\β})}{\mathrm{\Γ}\left(\mathrm{\α}\right)+\mathrm{\Γ}\left(\mathrm{\β}\right)}{\left(\frac{P_{pv}}{R_{pv}}\right)}^{\mathrm{\α}-1}{(1-\frac{P_{pv}}{R_{pv}})}^{\mathrm{\β}-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{\π}}{\mathrm{\δ}}_P}\exp (-\frac{{(P-{\mathrm{\μ}}_P)}^2}{2{{\mathrm{\δ}}_P}^2})---\left(7\right)$

$f\left(P\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\δ}}_Q}\exp (-\frac{{(Q-{\mathrm{\μ}}_Q)}^2}{2{{\mathrm{\δ}}_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\∈Ng}{\Σ}\[C_{Gpi}\left(P_{Gi}\right(t\left)\right)+C_{Gqi}\left(Q_i\right(t\left)\right)\]+\underset{j\∈Nq}{\Σ}{C}_{gqj}\·\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}{\Σ}}{V}_j(G_{ij}{\mathrm{cos\θ}}_{ij}+B_{ij}{\mathrm{sin\θ}}_{ij})=0---\left(10\right)$

$Q_{Gi}+Q_{gi}-Q_{Di}-V_i\underset{j=1}{\overset{N}{\Σ}}{V}_j(G_{ij}{\mathrm{sin\θ}}_{ij}+B_{ij}{\mathrm{cos\θ}}_{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{\‾}{P_{Gi}}\≤P_{Gi}\≤\stackrel{\‾}{P_{Gi}}\\ \underset{\‾}{Q_{Gi}}\≤Q_{Gi}\≤\stackrel{\‾}{Q_{Gi}}\\ \underset{\‾}{Q_{gi}}\≤Q_{gi}\≤\stackrel{\‾}{Q_{gi}}\\ \underset{\‾}{V_i}\≤V_i\≤\stackrel{\‾}{V_i}\\ S_i\≤\stackrel{\‾}{S_i}\end{array}\right.,\left\{\begin{array}{c}i\∈N_g\\ i\∈N_g\\ i\∈N_q\\ i\∈N\\ i\∈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\[C_{Gpi}\left(S_{Gi,\max }\right)-C_{Gpi}\left(\sqrt{S_{Gi,\max }^2-Q_{Gi}^2}\right)\]---\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\×365\×24\×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{\‾}{g}\≤g\left(x\right)\≤\stackrel{\‾}{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{\μ}(\underset{j=1}{\overset{r}{\Σ}}\ln l_j+\underset{j=1}{\overset{r}{\Σ}}\ln u_j)\end{array}$

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

$g\left(x\right)+u-\stackrel{\‾}{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\[g\left(x\right)-l-\underset{\‾}{g}\]-w^T\[g\left(x\right)+u-\stackrel{\‾}{g}\]-\mathrm{\μ}(\underset{j=1}{\overset{r}{\Σ}}\ln l_j+\underset{j=1}{\overset{r}{\Σ}}\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.