CN104156542A - Implicit-projection-based method for simulating stability of active power distribution system - Google Patents

Abstract

A stability simulation method for active power distribution systems based on implicit projection. Aiming at the stability simulation model of active power distribution systems with rigid characteristics, based on system simulation parameters and power flow calculation results, firstly use the internal integrator to perform several small steps Step-size integral calculation, the step size is h, and any explicit integration algorithm with second-order or higher precision can be used; then, according to the calculation result of the internal integrator, an external integrator based on the implicit prediction-correction method is used to perform a step-by-step calculation with a step size of Mh Integral calculation with large step size. The method of the invention can realize the simulation calculation of system faults, is a second-order precision algorithm, has good numerical stability, and its numerical stability domain hardly changes with the change of the step size multiple of the external integrator, and the algorithm performance is better than explicit projection integration The algorithm and the traditional implicit trapezoidal method are suitable for the rapid solution of the stability simulation problems of active power distribution systems with multi-time scale characteristics, and lay the foundation for the development of efficient and reliable active power distribution system simulation programs.

Description

A Stability Simulation Method for Active Power Distribution System Based on Implicit Projection

technical field

The invention relates to a stability simulation method of an active power distribution system. In particular, it relates to an implicit projection-based active power distribution system stability simulation method suitable for the stability simulation application of the active power distribution system including distributed power supply and energy storage.

Background technique

The large-scale and extensive access of distributed generation (DG) and the change of system load characteristics after the implementation of demand-side response technology all bring new challenges to the planning and operation of power distribution systems. The distribution network without DG is "passive", and the power used by users is provided by the upper-level transmission network. When the distribution network is connected to DG to generate bidirectional power flow, the system is called "active distribution system". ". Active power distribution system is a complex power distribution system capable of combined control of various distributed energy resources (DER, such as DG, controllable load, energy storage, etc.). In the future active power distribution system, the access capacity of DG can easily exceed (at least in a certain period of time) the total load in the power distribution system. At this time, the active power distribution system is used as an external power source to transmit electrical energy. Even if the total capacity of DG does not exceed the total load, the large-scale access of DG will still lead to changes in the dynamic response characteristics of the distribution network and affect the dynamic characteristics of the entire power system, especially the dynamic characteristics of large disturbances. At the system level, the analysis and research of relevant issues often cannot be tested directly on the actual system, so effective digital simulation tools must be used as an important research method.

Traditional power system time domain simulation has developed three power system digital simulation methods for different time scales of the system dynamic process: electromagnetic transient simulation, electromechanical transient simulation and medium and long-term dynamic simulation. have distinct characteristics. The electromagnetic transient simulation of the power system focuses on the fast dynamic change process of the voltage and current generated by the interaction between the electric field and the magnetic field in the system; etc.) dynamic behavior and the ability to maintain synchronous and stable operation, that is, transient stability, the time range of concern is usually from a few seconds to tens of seconds, so it is also called stability simulation; mid- and long-term dynamic process simulation is a power system The dynamic simulation of a long process after being disturbed, that is, the usual long process dynamic stability calculation of the power system.

In addition to focusing on the transient stability of traditional power systems, power system stability simulation has also focused on the analysis of power frequency electrical quantities in the system when active power distribution systems with various distributed power sources and energy storage devices are in operation. The dynamic response characteristics under (switching operation, fault, distributed power supply and load fluctuation, etc.) can be called active distribution system stability simulation at this time. The stability simulation of active power distribution system can essentially be attributed to the calculation of the time domain response of the dynamic system, which is divided into two parts: mathematical modeling and model solving. First, according to the topological relationship between the components, the characteristic equations of the components of the active power distribution system are constructed to form a stability simulation model of the whole system, forming a set of simultaneous differential-algebraic equations, and then the steady-state working condition or the power flow solution is used as the initial Value, to solve the numerical solution under the disturbance, that is, to gradually obtain the change curve of the system state quantity and algebraic quantity with time.

Active distribution system stability simulation modeling is the process of abstracting the mathematical model from the physical prototype according to the time scale range concerned by the system simulation. The mathematical model in the stability simulation of active power distribution system consists of two parts: the differential equation describing the dynamic characteristics of the equipment and the algebraic equation describing the electrical connection between the equipment. Among them, the electrical connection relationship between equipment may change during operation, such as load switching, unit start and stop, line opening and reclosing, etc. If the relay protection device is taken into account, it should also include a large number of continuous and (or) discrete logical time-varying parameters.

Generally, the mathematical model of active power distribution system can be described by a high-dimensional nonlinear and continuous autonomous differential-algebraic equations, as shown in formula (1).

$\left\{\begin{array}{c}\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>\\ \mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; g</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>\end{array}\right.<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

In the formula, is the differential equation, is an algebraic equation, is the system state variable, representing the motor rotor speed, power electronic device control system and load dynamic parameters, etc. is an algebraic variable, representing the bus voltage amplitude and phase angle. The solution of the mathematical model is generally realized through specific numerical algorithms and corresponding simulation programs. The active power distribution system is connected to a wide variety of distributed power sources and a large number of power electronic devices, including rotating motors and various static DC distributed power sources, which have obvious multi-time scale characteristics and can be attributed to rigid problems in mathematics . Therefore, the stability simulation of active power distribution system can be attributed to solving the initial value problem of a rigid differential-algebraic equation system mathematically, which requires higher accuracy and stability of the numerical algorithm used.

Active distribution system stability simulation algorithms can be divided into two categories: alternate solution method and simultaneous solution method according to the different solution forms of differential equation and algebraic equation in formula (1). Alternate solution method first adopts a specific numerical integration algorithm, solves the differential equation according to the initial calculation results, obtains the value of the state variable at this time step, and then substitutes it into the algebraic equation to obtain the value of the algebraic variable at this time step, and finally The algebraic variable is substituted into the differential equation to solve the state variable of the next time step, and so on to realize the alternate solution of the differential-algebraic equation system; the simultaneous solution method is to combine the differential equation with the algebraic equation to form a complete algebraic equation system. Simultaneously solve for state variables and algebraic variables.

For the differential equation in formula (1), except for a few cases where analytical solutions can be obtained, in most cases only numerical solutions can be used to solve it. Among them, the differential method is widely used in the stability simulation of active power distribution systems. Can be divided into single-step method (one step method) and linear multi-step method (linear multistep method). According to different solving processes, the single-step method can be divided into explicit integration method and implicit integration method. The explicit integration method can directly calculate the state variable at the next moment according to the state variable at the current moment, while the implicit integration method needs to Only by solving the equation of the state variable at the next moment can the state variable at the next moment be obtained. Common explicit integration methods include Euler method, improved Euler method and Runge-Kutta method, while implicit integration methods mainly include backward Euler method and implicit trapezoidal method. The active power distribution system has obvious rigid characteristics. For the explicit integration method, the calculation amount in each time step is small, but because of its poor numerical stability, only a small integration step can be used to solve the rigid problem. , which will seriously limit the simulation speed in the stability calculation; and for the implicit integration method, although it needs to iteratively solve the equation system in each time step, the calculation and programming work is more complicated than the explicit integration algorithm, but its numerical stability In the process of solving rigid problems, the numerical stability can be guaranteed, and the simulation speed can be improved by adopting a larger integration step size.

The explicit projection algorithm is a numerical integral solution algorithm proposed for the initial value problem of the ordinary differential equation (ODE) with rigid characteristics shown in the following formula (2). Integral calculation, the calculation step size corresponds to the time constant of the fast dynamic process of the system; then, according to the calculation result of the small step size, use the following formula (3) or the explicit prediction-correction process based on the improved Euler method to perform a projection step For integral calculation, the projection step corresponds to the time constant of the slow dynamic process of the system. Among them, the small-step integral calculation process is called an internal integrator, and an explicit four-order Runge-Kutta method with better numerical stability and higher precision is used to improve the stability of the algorithm and numerical accuracy; the projective integration process with large steps is called an external integrator. The explicit projection integration algorithm can further improve the simulation calculation speed while realizing the improvement of the numerical stability of the traditional explicit integration algorithm. Nevertheless, for the intelligent power distribution system with obvious multi-time scale characteristics, its numerical stability is still greatly limited, and it is very difficult to further improve its calculation speed.

$\left\{\begin{array}{c}\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>\\ \mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}\end{array}\right.<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

x(t _n+k+1+M )=(M+1)x(t _n+k+1 )-Mx(t _n+k ) (3)

It can be seen that it is very important to propose a stability simulation method with high numerical accuracy, good numerical stability, high computational efficiency, and suitable for active distribution systems with rigid characteristics.

Contents of the invention

The technical problem to be solved by the present invention is to provide a stability simulation method for an active power distribution system based on implicit projection that has high numerical accuracy, good numerical stability, high computational efficiency and is suitable for rigid characteristics.

The technical solution adopted in the present invention is: a method for simulation of the stability of an active power distribution system based on implicit projection, comprising the following steps:

1) According to the topological connection relationship of the system and the dynamic equation of the components, the transient stability simulation model of the active power distribution system is established, which is in the form of $\left\{\begin{array}{c}\stackrel{\&Center\; Dot;}{x}=f(x,\mathrm{the\; y})\\ 0=g(x,\mathrm{the\; y})\end{array}\right.,$ In the formula, is the differential equation, is an algebraic equation, is the system state variable, is the system algebraic variable;

2) Perform power flow calculation on the active power distribution system to obtain system power flow data;

3) Read system parameters and simulation calculation parameters, including simulation termination time T, simulation step size h, integration step number k of the internal integrator of the implicit projection algorithm, and the integration step size of the external integrator of the implicit projection algorithm relative to the implicit projection algorithm The multiple M of the integral step size of the internal integrator, k and M are all positive integers, and the simulated fault and operation event information is set, including fault occurrence and clearing time, fault location and fault type;

4) According to the calculation results of the system power flow, the simulation initialization calculation of the dynamic components of the whole system is carried out;

5) Set simulation time t=0;

6) The number of integration steps s=1 of the internal integrator of the current implicit projection algorithm is set, and s is a positive integer;

7) Set the simulation time t=t+h, h is the integral step size of the internal integrator of the implicit projection algorithm, and use the internal integrator of the implicit projection algorithm to calculate a step size for the model of the active power distribution system to obtain the state of the system at this moment Variable x _n+s and algebraic variable y _n+s , and set s=s+1;

8) according to step 3) simulated failure and operation event information that is set, judge whether system breaks down or operate at this moment, if the occurrence time T _event =t of failure and operation event, then return to step 6), otherwise proceed to the next step;

9) Judging whether the simulation time t reaches the simulation termination time T, if t=T, then the simulation ends, otherwise proceed to the next step;

10) Determine whether the number of integration steps s of the internal integrator of the implicit projection algorithm is greater than the number of integration steps k+1 of the internal integrator of the implicit projection algorithm input by the user in step 3), if s≤k+1, return to step 7 ), otherwise proceed to the next step;

11) According to the simulated failure and operation event information set in step 3), it is judged whether a failure or operation occurs within t～t+Mh, if t<T _event <t+Mh, then enter step 13), otherwise proceed to the next step;

12) Set the integral step size of the external integrator of the implicit projection algorithm H=Mh, set the simulation time t=t+Mh, and use the external integrator of the implicit projection algorithm to obtain the state variables x _n+k+1+H and Algebraic variable y _n+k+1+H , then directly enter step 14);

13) Set the integral step size of the external integrator of the implicit projection algorithm H=T _event -t, set the simulation time t=T _event , and use the external integrator of the implicit projection algorithm to obtain the state variable x _n+k of the system before the fault or operation occurs _+1+H and the algebraic variable y _n+k+1+H ;

14) Determine whether the simulation time t reaches the simulation termination time T, if t=T, then the simulation ends, otherwise return to step 6), and repeat steps 6) to 14) until the simulation ends.

The internal integrator of the implicit projection algorithm described in step 3) uses an explicit alternate solution method to solve the active power distribution system model, and selects a method with second-order or higher precision for solving the differential equation in the active power distribution system model. Arbitrary explicit numerical integration method.

Step 3) the described implicit projection algorithm external integrator, within an implicit projection algorithm external integrator integration step size H, the specific steps for solving the differential-algebraic equation describing the active power distribution system are as follows:

(1) Let the current simulation time be t _n+k+1 , where n is the total number of simulation steps at the current moment, the system state variable at the current moment is x _n+k+1 , and the algebraic variable is y _n+k+1 , After the step size H, the simulation time of the system is t _n+k+1+H , and the state variables and algebraic variables of the system are x _n+k+1+H and y _n+k+1+H respectively, and the description has The differential-algebraic equation of the source distribution system model is implicitly differentiated, and the following formula is obtained:

$\left\{\begin{array}{c}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span>\\ \mathrm{<span\; class="notranslate">\; g</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right.<span\; class="notranslate">\; ;</span>$

(2) Use the forward Euler method to obtain the predicted value of the initial estimated value of x _n+k+1+H as shown below

${\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>$

Then substitute into the equation Get the predicted value of the initial estimate of y _n+k+1+H

(3) Use the following formula to correct the predicted value to obtain the initial estimated value of x _n+k+1+H

${\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span>$

Then substitute into the algebraic equation Get the initial estimate of y _n+k+1+H in

(4) Obtain the correction value of x _n+k+1+H by the following formula

${\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span>$

Then Substitute the value of into the algebraic equation In the solution, the corrected value of y _n+k+1+H is obtained

(5) respectively and Substitute into the following formula to judge whether the convergence condition is satisfied,

$<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; \&xi;</span>}$

In the formula, ξ is the allowable error value set by the user. If the convergence condition is satisfied, the calculation step of the external integrator of the implicit projection algorithm ends; otherwise, the and replace and Return to step (4), repeat steps (4) and (5) until the convergence condition is satisfied.

An implicit projection-based active power distribution system stability simulation method of the present invention fully considers the rigid characteristics of the active power distribution system, and adopts an alternate solution method to describe the differential-algebraic equation of the active power distribution system simulation model Alternately solve the group, and use the implicit projection algorithm to solve the differential equations. The method of the invention is a second-order precision algorithm, and the performance of the algorithm is better than that of the traditional implicit trapezoidal method. At the same time, this method has good numerical stability, and its numerical stability domain hardly changes with the change of the step size multiple of the external integrator. It has obvious advantages over the explicit projection integration algorithm in terms of numerical stability and calculation speed, and is suitable for The rapid solution of the stability simulation problem of active distribution system with obvious multi-time scale characteristics has laid a good foundation for the development of efficient and reliable simulation program of active distribution system.

Description of drawings

Fig. 1 is the overall flowchart of the inventive method;

Fig. 2 is the numerical stability domain of method of the present invention and explicit projection algorithm and traditional explicit 4th order Runge-Kutta method;

Fig. 3 is a partially enlarged schematic diagram of A in Fig. 2;

Figure 4 is a structural diagram of a low-voltage active power distribution system example;

In the figure 1: the first storage battery; 2: the first photovoltaic cell; 3: the second photovoltaic cell; 4: the second storage battery; M1: medium voltage bus; S1: switch; L1~L19: low voltage bus; Load1~Load7: load ;

Figure 5 is the L17 bus voltage simulation results and a partial enlarged view;

Fig. 6 is the second photovoltaic cell active power output simulation result and a partial enlarged view;

Fig. 7 is a structure diagram of IEEE123 node active power distribution system calculation example;

Fig. 8 is the simulation result of the active power output of the photovoltaic cell at node 56 and a partial enlarged view;

Fig. 9 is the simulation result of grid-connected voltage output of photovoltaic cells at node 56 and a partial enlarged view;

Figure 10 is a comparison of the numerical accuracy of the implicit projection algorithm and the implicit trapezoidal method with a step size of 0.5ms (logarithmic coordinate system);

Figure 11 is a comparison of the numerical accuracy of the implicit projection algorithm and the implicit trapezoidal method with a step size of 0.02s (logarithmic coordinate system);

Figure 12 is a comparison of the numerical accuracy of the implicit projection algorithm and the implicit trapezoidal method with a step size of 0.0025s (logarithmic coordinate system).

Detailed ways

A method for simulating the stability of an active power distribution system based on implicit projection according to the present invention will be described in detail below in conjunction with the embodiments and drawings.

The invention relates to an implicit projection-based active power distribution system stability simulation method, belonging to an explicit and implicit mixed integral method. The active power distribution system is connected to a wide variety of distributed power sources with large differences in dynamic response characteristics and a large number of power electronic devices, which include both rotating AC motors and static DC power sources such as photovoltaic cells and batteries. The active power distribution system has obvious multi-time scale characteristics, so its digital simulation problem shows strong rigidity, and it is necessary to use a numerical integration algorithm with good numerical accuracy and numerical stability to realize its simulation calculation. An implicit projection-based active power distribution system stability simulation method proposed by the present invention fully considers the rigid characteristics of the active power distribution system, and uses an alternate solution method to describe the differential-algebraic equations of the active power distribution system model The group is solved alternately, and the differential equations are solved by the implicit projection algorithm. The method of the invention is a second-order precision algorithm, and the performance of the algorithm is better than that of the traditional implicit trapezoidal method. At the same time, this method has good numerical stability, and its numerical stability domain hardly changes with the change of the step size multiple of the external integrator. It has obvious advantages over the explicit projection integration algorithm in terms of numerical stability and calculation speed, and is suitable for The rapid solution of the stability simulation problem of active distribution system with obvious multi-time scale characteristics has laid a good foundation for the development of efficient and reliable simulation program of active distribution system.

The present invention uses an alternate solution algorithm to realize the calculation of the mathematical model of the active distribution system based on the differential-algebraic equation description, and uses the implicit projection algorithm to solve the differential equation. Integral calculation, the simulation step size corresponds to the fast dynamic process of the system; then, according to the result of the small step integral calculation, the implicit prediction-correction method is used to carry out a projection integration step with a large step size, and the step size corresponds to the slow dynamic process of the system correspond. Among them, the small-step integral calculation process is called an internal integrator, and any explicit numerical integration algorithm with second-order or higher precision can be used; the large-step integral calculation process is called an external integrator.

As shown in Figure 1, a kind of implicit projection based active power distribution system stability simulation method of the present invention is characterized in that, comprises the following steps:

1) According to the topological connection relationship of the system and the dynamic equation of the components, the transient stability simulation model of the active power distribution system is established, which is in the form of $\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=f(x,\mathrm{the\; y})\\ 0=g(x,\mathrm{the\; y})\end{array}\right.,$ In the formula, is the differential equation, is an algebraic equation, is the system state variable, is the system algebraic variable;

2) Perform power flow calculation on the active power distribution system, and the power flow calculation is calculated by using Gauss method or Newton-Raphson method. The system power flow data is obtained through power flow calculation, including node voltage, current, load active power and reactive power, power supply active power output, reactive power output and injection current, etc.;

3) Read system parameters and simulation calculation parameters, including simulation termination time T, simulation step size h, integration step number k of the internal integrator of the implicit projection algorithm, and the integration step size of the external integrator of the implicit projection algorithm relative to the implicit projection algorithm The multiple M of the integral step size of the internal integrator, k and M are all positive integers, and the simulated fault and operation event information is set, including fault occurrence and clearing time, fault location and fault type;

The internal integrator of the implicit projection algorithm is to use the explicit alternate solution method to solve the active power distribution system model, and to solve the differential equation in the active power distribution system model, select any explicit Numerical integration method.

The specific steps for solving the differential-algebraic equation describing the active power distribution system within an implicit projection algorithm external integrator integration step size H are as follows:

(1) Let the current simulation time be t _n+k+1 , where n is the total number of simulation steps at the current moment, the system state variable at the current moment is x _n+k+1 , and the algebraic variable is y _n+k+1 , After the step size H, the simulation time of the system is t _n+k+1+H , and the state variables and algebraic variables of the system are x _n+k+1+H and y _n+k+1+H respectively, and the description has The differential-algebraic equation of the source distribution system model is implicitly differentiated, and the following formula is obtained:

$\left\{\begin{array}{c}{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span>\\ \mathrm{<span\; class="notranslate">\; g</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right.<span\; class="notranslate">\; ;</span>$

(2) Use the forward Euler method to obtain the predicted value of the initial estimated value of x _n+k+1+H as shown below

${\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>$

Then substitute into the equation Get the predicted value of the initial estimate of y _n+k+1+H

(3) Use the following formula to correct the predicted value to obtain the initial estimated value of x _n+k+1+H

${\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span>$

Then substitute into the algebraic equation Get the initial estimate of y _n+k+1+H in

(4) Obtain the correction value of x _n+k+1+H by the following formula

${\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span>$

Then Substitute the value of into the algebraic equation In the solution, the corrected value of y _n+k+1+H is obtained

(5) respectively and Substitute into the following formula to judge whether the convergence condition is satisfied

$<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; \&xi;</span>}$

In the formula, ξ is the allowable error value set by the user. If the convergence condition is satisfied, the calculation step of the external integrator of the implicit projection algorithm ends; otherwise, the and replace and Return to step (4), repeat steps (4) and (5) until the convergence condition is satisfied.

4) According to the calculation results of the system power flow, the simulation initialization calculation of the dynamic components of the whole system is carried out;

5) Set simulation time t=0;

6) The number of integration steps s=1 of the internal integrator of the current implicit projection algorithm is set, and s is a positive integer;

7) Set the simulation time t=t+h, h is the integral step size of the internal integrator of the implicit projection algorithm, and use the internal integrator of the implicit projection algorithm to calculate a step size for the model of the active power distribution system to obtain the state of the system at this moment Variable x _n+s and algebraic variable y _n+s , and set s=s+1;

8) according to step 3) simulated failure and operation event information that is set, judge whether system breaks down or operate at this moment, if the occurrence time T _event =t of failure and operation event, then return to step 6), otherwise proceed to the next step;

9) Judging whether the simulation time t reaches the simulation termination time T, if t=T, then the simulation ends, otherwise proceed to the next step;

10) Determine whether the number of integration steps s of the internal integrator of the implicit projection algorithm is greater than the number of integration steps k+1 of the internal integrator of the implicit projection algorithm input by the user in step 3), if s≤k+1, return to step 7 ), otherwise proceed to the next step;

11) According to the simulated failure and operation event information set in step 3), it is judged whether a failure or operation occurs within t～t+Mh, if t<T _event <t+Mh, then enter step 13), otherwise proceed to the next step;

12) Set the integral step size of the external integrator of the implicit projection algorithm H=Mh, set the simulation time t=t+Mh, and use the external integrator of the implicit projection algorithm to obtain the state variables x _n+k+1+H and Algebraic variable y _n+k+1+H , then directly enter step 14);

13) Set the integral step size of the external integrator of the implicit projection algorithm H=T _event -t, set the simulation time t=T _event , and use the external integrator of the implicit projection algorithm to obtain the state variable x _n+k of the system before the fault or operation occurs _+1+H and the algebraic variable y _n+k+1+H ;

14) Determine whether the simulation time t reaches the simulation termination time T, if t=T, then the simulation ends, otherwise return to step 6), and repeat steps 6) to 14) until the simulation ends.

Specific examples are given below:

For the linear integration algorithm, the linear constant coefficient differential equation The solution of and the scalar equation

$\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&lambda;x</span>}$

The solutions are equivalent, where λ is the characteristic root of the matrix A, so the above formula is also called the scalar test equation. For a given numerical integration algorithm, its numerical stability domain refers to the s domain, and the algorithm satisfies the numerical stability condition when solving the scalar test equation

|σ(hλ)|≤1

The set of hλ.

This example takes the explicit fourth-order Runge-Kutta method (hereinafter referred to as RK4 algorithm) as an example for the internal integrator of the implicit projection algorithm. According to the numerical stability conditions, the implicit projection algorithm, the traditional RK4 algorithm and the explicit The numerical stability domain of the formula projection algorithm in the hλ plane, as shown in Figure 2 and Figure 3, in which the internal integrator of the explicit projection algorithm adopts the RK4 algorithm, and the external integrator adopts the explicit prediction-correction method based on the improved Euler method . It can be seen from the figure that as the multiple M of the step size of the external integrator increases, the numerical stability domain of the explicit projection algorithm decreases significantly, and gradually splits from a large area into several smaller areas. However, the numerical stability domain of the implicit projection algorithm is basically the same as that of the traditional RK4 algorithm, and hardly changes with the change of M. The numerical stability domains of the RK4 algorithm coincide. Therefore, the numerical stability of the implicit projection integral algorithm proposed by the present invention is better than that of the explicit projection integral algorithm.

Digital simulation and grid calculation program (DIgSILENT PowerFactory) is a commercial power system simulation software developed by DIgSLENT GmbH in Germany. This example realizes a kind of implicit projection-based active power distribution system stability simulation method proposed by the present invention in the C++ programming environment, and compares the simulation results and calculation performance of the implicit projection integral method with commercial software DIgSILENTPowerFactory To verify the correctness and effectiveness of the method, the hardware platform for the simulation test is a PC with Intel(R) Core(TM) i5-3470 CPU 3.20GHz and 4GB RAM; the software environment is a 32-bit Windows 7 operating system. This example tests the method of the present invention by selecting different algorithm parameters. The algorithm parameters can be arbitrarily selected under the condition of satisfying the numerical accuracy according to the actual situation during specific implementation, and the implementation of the present invention is not limited to this.

First, this example uses a low-voltage active power distribution system example with distributed power sources to test and verify the method of the present invention, as shown in Figure 4. The voltage level of the low-voltage active power distribution system example is 400V. The main feeder is connected to the medium-voltage bus M1 through a 0.4/10kV transformer. 50m, using three-phase symmetrical lines and loads. In addition, various types of distributed power sources are connected in the calculation example, including: photovoltaic power generation system and battery energy storage system with maximum power tracking control. The control methods, access capacity and active power output of each distributed power source are shown in Table 1 shown.

Table 1 Distributed power supply control mode, access capacity and output power

An implicit projection-based active power distribution system stability simulation method of the present invention is used to perform stability simulation calculations on test examples, and the simulation time is set to 9s, and the simulation step size is 0.5ms. At 2.0s, the switch S1 of the low-voltage active power distribution system is turned off, and the system switches from the grid-connected operation mode to the island operation mode; at 4.7s, the S1 switch is closed, and the system switches from the island operation mode to the grid-connected operation mode.

A kind of active power distribution system stability simulation method based on implicit projection of the present invention, wherein the algorithm parameters take k=3, M=4, compare with the fixed-step simulation results of DIgSILENT PowerFactory and traditional implicit trapezoidal method , where the internal integration step of the implicit projection algorithm is 0.5ms, the external integration step H=Mh=0.002s, the simulation step of DIgSILENT is the same as the internal integration step of the implicit projection algorithm, and the simulation step of the traditional implicit trapezoidal method The length is the same as the outer integration step size. The simulation results of the L17 bus voltage and the active power output of the No. 2 battery and the partial enlarged diagrams are shown in attached drawings 5 and 6. It can be seen that the simulation results of the implicit projection algorithm are basically consistent with DIgSILENT. In addition, by comparing with the traditional implicit trapezoidal method, it can be seen that the numerical accuracy of the method of the present invention is better than that of the traditional implicit trapezoidal method when the step size is taken as the external integral step size due to the influence of the integral calculation of the small step size of the internal integrator.

In order to verify the adaptability of the implicit projection integral algorithm of the present invention to large-scale active power distribution systems with rigid characteristics, this example is based on the IEEE123 node distribution network standard example (as shown in Figure 7), from the numerical value The algorithm is comprehensively tested in terms of accuracy and computing performance. The IEEE123 node calculation example describes a radial power distribution network with a complex structure, a total of 123 nodes, and a voltage level of 4.16kV. Various types of loads are considered internally, and the node 150 is connected to the external network. In this example, a total of 50 photovoltaic power generation systems with a capacity of 30kWp and an active power output of 20.4kW are connected to the nodes in the dotted line box in Figure 7.

Using the implicit projection-based active power distribution system stability simulation method of the present invention to perform stability simulation calculations on the IEEE123 node example, the simulation time is set to 9s, the simulation step is 0.5ms, and the environment of the example is initially The light intensity is set to 1000W/m ² , and the light intensity changes to 1025W/m ² , 1010W/m ² and 1000W/m ² at 1.5s, 3.5s and 6s, respectively.

In this example, the parameters of the implicit projection algorithm are k=3, M=4, the internal integration step is 0.5ms, and the external integration step is H=Mh=0.002s. Comparing the fixed-step simulation results of the implicit projection algorithm with DIgSILENT whose step size is the internal integration step and the traditional implicit trapezoidal method whose step size is the external integration step, the photovoltaic active power output at node 56 and the grid-connected bus The simulation results of the voltage and the partial enlarged diagrams are shown in Figure 8 and Figure 9 . It can be seen from the simulation results that the simulation results of the implicit projection algorithm are still basically consistent with DIgSILENT. At the same time, it can be seen from the comparison with the traditional implicit trapezoidal method that the numerical accuracy of the implicit projection algorithm of the present invention is better than that of the traditional implicit trapezoidal method in which the step size is taken as the external integral step size.

In order to compare the numerical accuracy of the implicit projection algorithm and the traditional implicit trapezoidal method under different step lengths, in this example, the parameters of the projection algorithm are set to k=6, M=40, the internal integration step is 0.5ms, and the external integration step is H=Mh = 0.02s. Assuming that the simulation step size of the traditional implicit trapezoidal method is h _TR , the step sizes are respectively h _TR =0.5ms, h _TR =0.02s and h _TR =0.0025s, and the explicit fourth-order Runge- Based on the Kutta method, the absolute error of the simulation results of the implicit projection algorithm and the implicit trapezoidal method under different step lengths relative to the simulation results of the RK4 algorithm was compared in the logarithmic coordinate system, as shown in Figures 10 to 12. It can be seen from the figure that when h _TR is the same as the internal integration step of the implicit projection algorithm, the absolute error of the simulation result of the trapezoidal method is smaller than that of the implicit projection algorithm, and the numerical accuracy is higher; when h _TR is the same as the external integration step, the projection The numerical precision of the algorithm is higher than that of the trapezoidal method; when h _TR is 0.0025s, the absolute error of the projection algorithm is slightly lower than that of the trapezoidal method.

It can be seen from accompanying drawings 2 and 3 that a kind of implicit projection-based active power distribution system stability simulation method of the present invention has greater advantages in numerical stability than the explicit projection-integration algorithm, so it can be obtained by A larger M value is used to further improve the simulation calculation speed. This example takes the IEEE123 node active power distribution system as a test example, and uses the explicit fourth-order Runge-Kutta method with a step size of 0.5ms as the benchmark, selects different simulation step sizes and algorithm parameters, and compares the implicit projections The calculation efficiency of algorithm, explicit projection algorithm, DIgSILENT and traditional implicit trapezoidal legal step simulation are shown in Table 2.

Table 2 Algorithm performance comparison

It can be seen from Table 2 that as the value of k decreases or the value of M increases, the calculation time of the implicit projection algorithm decreases gradually. When the parameters of the projection algorithm are the same, the computational efficiency of the implicit projection algorithm is slightly lower than that of the explicit projection algorithm because the calculation process of the external integrator of the implicit projection algorithm is slightly more complicated than that of the explicit projection algorithm and consumes more computing resources. However, when k=3 and M=8, the explicit projection algorithm has not converged numerically, while the implicit projection algorithm can still maintain numerical stability even when M is 60. At this time, the simulation time of the implicit projection algorithm is much shorter than that of the traditional RK4 Algorithm and DIgSILENT, the speedup ratio can reach more than 7 times. In addition, the implicit projection algorithm with algorithm parameters k=6, M=40 is compared with the implicit trapezoidal method when h _TR is equal to the inner integration step size, outer integration step size and 0.0025s. It can be seen that when h _TR takes the internal integral step of the implicit projection algorithm, the computational efficiency of the implicit projection algorithm is much higher than that of the trapezoidal method; when h _TR takes the external integral step, the computational efficiency of the implicit projection algorithm is slightly lower Compared with the trapezoidal method; when h _TR is 0.0025s, the calculation efficiency of the projection algorithm is higher than that of the trapezoidal method. From the analysis of accompanying drawings 10 to 12 and Table 2, it can be seen that for the traditional implicit trapezoidal method under different simulation step sizes, the implicit projection algorithm has certain advantages in terms of numerical accuracy or calculation efficiency compared with the implicit trapezoidal method. , and when the implicit trapezoidal method takes a certain intermediate step size, the numerical accuracy and computational efficiency of the implicit projection algorithm are superior to the trapezoidal method at the same time. Therefore, the implicit projection-based active power distribution system stability simulation method proposed by the present invention has better algorithm performance than the traditional implicit trapezoidal method.

In summary, the implicit projection-based active power distribution system stability simulation method of the present invention has good numerical accuracy and numerical stability, and can greatly improve the calculation efficiency, especially suitable for rigid characteristics. The large-scale active power distribution system stability simulation calculation has laid a good foundation for the development of efficient and reliable active power distribution system simulation programs.

Claims (3)

1. the active distribution system Simulation of stability method based on implicit expression projection, is characterized in that, comprises the steps:

1) according to the dynamic equation of the topological connection relation of system and element, set up active distribution system transient stability realistic model, shape as $\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=f(x,y)\\ 0=g(x,y)\end{array}\right.,$ In formula, for the differential equation, for algebraic equation, for system state variables, for system algebraically variable;

2) active distribution system is carried out to trend calculating, obtain system load flow data;

3) reading system parameter and simulation calculation parameter, comprise emulation termination time T, simulation step length h, the outside integrator integration step of the integration step number k of implicit expression projection algorithm internal integral device and implicit expression projection algorithm is with respect to the multiple M of implicit expression projection algorithm internal integral device integration step, k and M are positive integer, emulation fault and operation event information are set, comprise fault generation and checkout time, abort situation and fault type;

4), according to system load flow result of calculation, system-wide dynamic element is carried out to simulation initialisation calculating;

5) simulation time t=0 is set;

6) the integration step number s=1 of current implicit expression projection algorithm internal integral device is set, s is positive integer;

7) simulation time t=t+h is set, h is implicit expression projection algorithm internal integral device integration step, adopts implicit expression projection algorithm internal integral device to calculate a step-length to active distribution system model, obtains the state variable x in this moment of system _n+swith algebraically variable y _n+s, and s=s+1 is set;

8) according to step 3) the emulation fault and the operation event information that arrange judge whether system now breaks down or operate, if the time of origin T of fault and Action Events _event=t, returns to step 6), otherwise carry out next step;

9) judge whether simulation time t reaches emulation termination time T, if t=T, emulation finishes, otherwise carries out next step;

10) whether the integration step number s that judges implicit expression projection algorithm internal integral device is greater than step 3) in the integration step number k+1 of implicit expression projection algorithm internal integral device of user's input, if s≤k+1 returns to step 7), otherwise carry out next step;

11) according to step 3) the emulation fault that arranges and operation event information judge in t～t+Mh time, whether break down or operate, if t<T _event<t+Mh, enters step 13), otherwise carry out next step;

12) the outside integrator integration step of implicit expression projection algorithm H=Mh is set, simulation time t=t+Mh is set, utilize the outside integrator of implicit expression projection algorithm to obtain the state variable x in this moment of system _n+k+1+Hwith algebraically variable y _n+k+1+H, then directly enter step 14);

13) the outside integrator integration step of implicit expression projection algorithm H=T is set _event-t, arranges simulation time t=T _event, utilize the outside integrator of implicit expression projection algorithm to obtain the state variable x of fault or the front system of operation generation _n+k+1+Hwith algebraically variable y _n+k+1+H;

14) judge whether simulation time t reaches emulation termination time T, if t=T, emulation finishes, otherwise returns to step 6), according to step 6) to 14) repeatedly carry out until emulation finishes.

2. a kind of active distribution system Simulation of stability method based on implicit expression projection according to claim 1, it is characterized in that, step 3) described implicit expression projection algorithm internal integral device, be to adopt explicit alternately method for solving to solve active distribution system model, the differential equation in active distribution system model chosen to any explicit numerical integration algorithm with the above precision of second order.

3. a kind of active distribution system Simulation of stability method based on implicit expression projection according to claim 1, it is characterized in that, step 3) the outside integrator of described implicit expression projection algorithm, as follows to describing the concrete solution procedure of differential-algebraic equation of active distribution system in the outside integrator integration step of implicit expression projection algorithm H:

(1) establishing current simulation time is t _n+k+1, wherein, the emulation total step number that n is current time, current time system state variables is x _n+k+1, algebraically variable is y _n+k+1, through step-length H, the system emulation time is t _n+k+1+H, now the state variable of system and algebraically variable are respectively x _n+k+1+Hand y _n+k+1+H, to describing differential-algebraic equation implicit expression differencing of active distribution system model, obtain following formula:

$\left\{\begin{array}{c}{x}_{n+k+1+H}=x_{n+k+1}+\frac{1}{2}H[f(x_{n+k+1},y_{n+k+1})+f(x_{n+k+1+H},y_{n+k+1+H})]\\ g(x_{n+k+1+H},y_{n+k+1+H})=0\end{array}\right.;$

(2) utilize forward direction Euler method to obtain x _n+k+1+Hthe predicted value of initial estimate be shown below

$x_{n+k+1+H}^{*}=x_{n+k+1}+\mathrm{Hf}(x_{n+k+1},y_{n+k+1})$

Substitution equation then obtain y _n+k+1+Hthe predicted value of initial estimate value

(3) utilize following formula to proofread and correct predicted value, obtain x _n+k+1+Hinitial estimate

$x_{n+k+1+H}^{\left(0\right)}=x_{n+k+1}+\frac{1}{2}H[f(x_{n+k+1},y_{n+k+1})+f(x_{n+k+1+H}^{*},y_{n+k+1+H}^{*})]$

Then be updated to algebraic equation in obtain y _n+k+1+Hinitial estimate

(4) by following formula, obtain x _n+k+1+Hmodified value

$x_{n+k+1+H}^{\left(1\right)}=x_{n+k+1}+\frac{1}{2}H[f(x_{n+k+1},y_{n+k+1})+f(x_{n+k+1+H}^{\left(0\right)},y_{n+k+1+H}^{\left(0\right)})]$

Then will value substitution algebraic equation in, solve and obtain y _n+k+1+Hmodified value

(5) respectively will with substitution following formula judges whether to meet the condition of convergence,

$\left|\right|x_{n+k+1+H}^{\left(1\right)}-x_{n+k+1+H}^{\left(0\right)}\left|\right|<\mathrm{\&xi;}$

In formula, ξ is the error permissible value being set by the user, if meet the condition of convergence, the outside integrator calculation procedure of implicit expression projection algorithm finishes; Otherwise, respectively will with replace with return to step (4), repeating step (4) and (5) are until meet the condition of convergence.