CN104375876A - 0+ error immune electromagnetic transient simulation algorithm under input quantity mutation condition - Google Patents

Abstract

The invention relates to a 0+ error immune electromagnetic transient simulation algorithm under the condition of sudden change of input quantity, belonging to the technical field of electromagnetic transient analysis of a power system. The method utilizes the impulse response invariant principle, the numerical sampling of the system continuous signal corresponds to the Z transformation of the discrete signal, and the basic form of the element based on the operation conductance and the historical current item is given. Based on the form, a system network equation is obtained according to a node analysis method, and a system electromagnetic transient simulation result is obtained according to the same steps. The method has the advantages that the truncation error is smaller than that of the currently adopted method, and the 0+ error is caused under the condition of mutation of the immune input quantity. The method solves the problems of nonlinear load modeling of the power system and jump numerical oscillation of the state quantity of the power electronic switch model. The 0+ error immune electromagnetic transient simulation algorithm under the condition of sudden change of the input quantity provides a new means for electromagnetic transient branch level modeling of the power system.

Description

0+ error immunity electromagnetic transient simulation algorithm under a kind of input quantity catastrophe

Technical field

The present invention relates to the 0+ error immunity electromagnetic transient simulation algorithm under a kind of input quantity catastrophe, belong to Electromagnetic Transient Analysis of Power System technical field.

Background technology

Along with the continuous expansion, interconnected of electric system scale, power electronic equipment increasing, electromagnetic transient analysis technology, becomes the gordian technique of power network safety operation.But when carrying out electromagnetic transient in power system emulation, if network has switch motion, network structure change, fault etc. will cause the transition of state variable, if at this moment still adopt the trapezoidal algorithm of traditional implicit expression will produce the numerical oscillation of non-prototype, in addition when emulating control system, a step time lag of control system and main system also can cause numerical value instability problem.

Tradition electromagnetic transient analysis method carries out modeling based on implicit trapezoid method to branch road.Result of study shows that the method there will be non-prototype numerical value oscillatory occurences under quantity of state saltus step situation.Even if adopt method of interpolation and damping resistance measure, non-prototype can only be suppressed to vibrate, and cannot fundamentally eliminate, the affix error that the use of damping resistance can be artificial in addition.Especially at failure situations and network structure changes, now fault current is larger, improves one's methods and can bring obvious error.Therefore, the electro-magnetic transient numerical algorithm of quantity of state jump error immunity has important theory and practical meaning in engineering at relay protection setting calculation, system failure Current calculation.

Summary of the invention

The object of this invention is to provide the 0+ error immunity electromagnetic transient simulation algorithm under a kind of input quantity catastrophe, during to realize quantity of state saltus step, the immunity of quantity of state jump error, avoid occurring non-prototype oscillatory occurences, quantity of state saltus step numerical value oscillation problem during solution electric system nonlinear-load modeling.

0+ error immunity electromagnetic transient simulation algorithm under the input quantity catastrophe that the present invention proposes, comprises the following steps:

(1) set up the branch road level Equivalent Model of an electric system, detailed process is as follows:

(1-1) the voltage-current characteristic equation setting electric system branch road is:

$\frac{\mathrm{di}\left(t\right)}{\mathrm{dt}}=f(i\left(t\right),u\left(t\right))$

Wherein, i (t), u (t) are respectively voltage time-domain signal and the electric current time-domain signal of port corresponding with each branch road, and f (i (t), u (t)) is a smooth function in real number field;

(1-2) according to above-mentioned voltage-current characteristic equation, the frequency domain equation of each branch road under electric current and voltage Laplace transform is obtained:

I(s)＝H(s)U(s)

Wherein, I (s), U (s) are respectively voltage frequency-region signal and the electric current frequency-region signal of port corresponding with each branch road, and H (s) is the transport function between the electric current and voltage corresponding with each branch road;

(1-3) transfer function H (s) of above-mentioned steps (1-2) is rewritten into zero pole point expression formula:

$H_a\left(s\right)=A\frac{\underset{s=1}{\overset{N}{\mathrm{\Π}}}(s-z_s)}{\underset{t=1}{\overset{M}{\mathrm{\Π}}}(s-p_t)},$

Wherein: A is the first proportional gain factor, z _sfor the zero point of transport function, p _tfor the limit of transport function, N is the zero number in transport function, s=1,2,3 ... N, M are the zero number in transport function, t=1,2,3 ... M;

(1-4) according to the expression formula of above-mentioned steps (1-3), according to Impulse invariance procedure principle, the transport function of each branch road under discrete domain is obtained as follows:

$H_d\left(z\right)=B\frac{\underset{s=1}{\overset{N}{\mathrm{\Π}}}(z-e^{z_{s}h})}{\underset{t=1}{\overset{M}{\mathrm{\Π}}}(z-e^{p_{t}h})}$

Wherein: B is the second proportional gain factor, z _sfor the zero point of transport function, p _tfor the limit of transport function, h is electromagnetic transient simulation material calculation, and N is the zero number in transport function, s=1,2,3 ... N, M are the zero number in transport function, t=1,2,3 ... M;

(1-5) transport function under the discrete domain of above-mentioned steps (1-4) is converted into the difference equation under time domain:

$i\left(n\right)=\mathrm{Cv}\left(n\right)+\underset{k=1}{\overset{s}{\mathrm{\Σ}}}{D}_{k}v(n-k)+\underset{k=1}{\overset{t}{\mathrm{\Σ}}}{E}_{k}i(n-k),$

Obtain the coefficient E of historical current item in the coefficient D of history voltage item in the coefficient C of current voltage item in difference equation, difference equation and difference equation,

Wherein, n is the n-th electromagnetic transient simulation step-length, and k is the Difference Terms number of difference equation; V (n) is the voltage item of the n-th electromagnetic transient simulation; I (n) is the electric current item of the n-th electromagnetic transient simulation, and v (n-k) is the history voltage item of the n-th-k electromagnetic transient simulation; I (n-k) is the historical current item of the n-th-k electromagnetic transient simulation;

(1-6) difference equation of above-mentioned steps (1-5) is rewritten into computing conductance in parallel and historical current item form, the Equivalent Model obtaining each branch road of electric system is: i (n)=gv (n)+i _hist(n-1),

Wherein, g is computing specific conductance in parallel, i _hist(n-1) be the historical current item of (n-1)th electromagnetic transient simulation,

$\left\{\begin{array}{c}g=C\\ i_{\mathrm{hist}}(n-1)=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}{D}_{k}v(n-k)+\underset{k=1}{\overset{t}{\mathrm{\Σ}}}{E}_{k}i(n-k)\end{array}\right.;$

(2) load in electric system is judged, if without load in electric system, then carries out step (3), if having nonlinear-load in electric system, then carry out following steps:

(2-1) constant impedance setting electric system is Z _abc, Z _abcbe 3 × 3 matrixes, the steady current of electric system is I _abc, I _abcbe 3 × 1 column vectors, the firm power of electric system is , wherein P is the active power of firm power, and Q is the reactive power of firm power;

The equivalent impedance calculating firm power load is as follows:

Wherein for the complex power of firm power load under electric power system tide stable state, for conjugation, for firm power load bus voltage under electric power system tide stable state,

Above-mentioned equivalent impedance is converted to three phase of impedance: i ₃for 3*3 unit matrix, for the tensor product of matrix;

The offset current calculating firm power load is as follows:

(2-2) according to voltage time-domain signal and the electric current time-domain signal of electric system, the power difference of electric system firm power load is obtained $\mathrm{\Δ}\stackrel{\·}{S}\left(t\right):\mathrm{\Δ}\stackrel{\·}{S}\left(t\right)=\stackrel{\·}{S}-(p\left(t\right)+\mathrm{jq}\left(t\right))$

p(t)＝v _a(t)i _a(t)+v _b(t)i _b(t)+v _c(t)i _c(t)

Wherein, $q\left(t\right)=\frac{1}{\sqrt{3}}[(v_b\left(t\right)-v_c\left(t\right))i_a\left(t\right)+(v_c\left(t\right)-v_a\left(t\right))i_b\left(t\right)+(v_a\left(t\right)-v_b\left(t\right))i_c\left(t\right)]$

V _a(t), v _b(t), v _ct () is respectively the three-phase voltage time-domain signal of electric system firm power load, i _a(t), i _b(t), i _ct () is respectively the three-phase current time-domain signal of electric system firm power load, for the complex power of firm power load under electric power system tide stable state, j ²=-1;

(2-3) according to the power difference of above-mentioned electric system firm power load obtain the amplitude I of electric system firm power load compensation electric current _m(t):

$I_m\left(t\right)=\frac{2\left|\right|\mathrm{\Δ}\stackrel{\·}{S}\left(t\right)\left|\right|}{{3V}_m\left(t\right)}$

Wherein V _mt () is the voltage magnitude of firm power load bus:

$V_m\left(t\right)=\sqrt{\frac{2(v_a^2\left(t\right)+v_b^2\left(t\right)+v_c^2\left(t\right))}{3}}$

V _a(t), v _b(t), v _ct () is respectively the three-phase voltage time-domain signal of electric system firm power load;

(2-4) according to the result of calculation of above-mentioned steps (2-2) and step (2-3), the three-phase offset current of firm power load is obtained:

ω is electric system angular frequency, be respectively the initial phase of firm power load three-phase current time-domain signal, ω, with obtained by the Load flow calculation of electric system;

(2-5) according to the result of calculation of above-mentioned steps (2-4), obtain computing conductance in parallel and the historical current item of nonlinear-load in electric system, wherein computing impedance in parallel is by Z _abcand Z' _abcparallel connection obtains, and computing conductance in parallel is the inverse of computing impedance in parallel.The electric current of historical current item is by I _abcwith Δ I _abcparallel connection obtains, wherein Δ I _abcfor constant impedance offset current item, Δ I _abc=[Δ i _aΔ i _bΔ i _c] ^t, T is matrix transpose, Z _abcfor the constant impedance of electric system, I _abcfor the steady current of electric system;

(3) set up the nodal voltage equation of an electric system, comprise the following steps:

(3-1) according to the computing specific conductance g in parallel of above-mentioned steps (1-6), the branch road conductance matrix G of electric system is obtained,

$G=\tilde{A}\mathrm{diag}\left(g_i\right){\tilde{A}}^T$

Wherein for the topological correlation matrix of electric system, the element in topological correlation matrix is ± 1 or 0, g _ifor the computing specific conductance in parallel of bar branch road every in electric system, diag (g _i) diagonal matrix that forms for the computing specific conductance in parallel of every bar branch road;

(3-2) the historical current item of each branch road in the electric system obtained according to above-mentioned steps (1-6), obtains branch node historical current vector I _hist(n):

$I_{\mathrm{hist}}(n-1)=\tilde{A}\mathrm{diag}\left\{i_{\mathrm{hist}}(n-1)\right\};$

Wherein, diag{i _hist(n-1) } be diagonal matrix that the historical current item of every bar branch road is formed;

(3-3) according to the result of calculation of above-mentioned steps (3-1) and step (3-2), the nodal voltage equation of electric system is obtained:

GV(n)＝I(n)+I _hist(n-1)

V (n) is the voltage vector of electric system node in the n-th electromagnetic transient simulation step-length, solves this nodal voltage equation and obtains V (n), and I (n) is the current vector of electric system node in the n-th electromagnetic transient simulation;

(4) implicit trapezoid method is utilized, the inertia branch model of the Synchronous Machine Models of electric system, electric system, generator excitation, speed governing are incorporated to together with prime mover model in the Electric Power System Node Voltage equation of above-mentioned steps (3), obtain an expanding node voltage equation;

(5) solve the expanding node voltage equation of above-mentioned steps (4), obtain a, b, c three-phase node voltage of electric system, and according to three-phase node voltage, calculate the Current calculation value in electric system generator model;

(6) the current error threshold value of a generator model is set, the Current calculation value of the electric current predictand in electric system generator model and above-mentioned steps (5) is compared, if error is less than or equal to error threshold, then obtain the electromagnetic transient simulation result of current state simulation step length, carry out step (7), if error is greater than error threshold, then repeat step (4)-step (6), when returning step (4), keep the terminal voltage of the node be connected with generator constant;

(7) repeat step (2)-step (6), obtain the electromagnetic transient simulation result of the full step-length of electric system.

0+ error immunity electromagnetic transient simulation algorithm under the input quantity catastrophe that the present invention proposes, its advantage is:

The inventive method under Electromagnetic Transient Analysis of Power System, based on impulse response invariance principle, can the non-prototype oscillatory occurences brought of immune state amount saltus step amount, have essential distinction with classic method.The electro-magnetic transient numerical algorithm of quantity of state jump error immunity is that the modeling of electromagnetic transient in power system branch road level provides new means.This method is A stable algorithm for typical branch and electronic power switch model.This method can the non-prototype oscillatory occurences that causes of immune state amount saltus step, accurately obtains the Current calculation value under electric power system fault, contributes to relay protection setting calculation, system failure Current calculation.Be particularly useful for the occasion of many electronic power switches, avoid switch motion to cause quantity of state jump error.In addition, this method is in stable state calculating, higher than classic method precision, and truncation error is less, contributes to the precision improving electromagnetic transient in power system simulation calculation further.Further, present method solves the modeling of electric system nonlinear-load and electronic power switch model state amount saltus step numerical value oscillation problem.

Accompanying drawing explanation

Fig. 1 is the FB(flow block) of the inventive method.

Fig. 2 is the schematic diagram of the resistance capacitance series loop system that electromagnetic transient in power system emulation mode of the present invention relates to.

Fig. 3 is the effect schematic diagram of the inventive method compared with the capacitance voltage curves (implicit trapezoid method, this method, backward Euler method, theoretical curve) of prior art.

Embodiment

0+ error immunity electromagnetic transient simulation algorithm under the input quantity catastrophe that the present invention proposes, its FB(flow block) as shown in Figure 1, comprises the following steps:

(1) set up the branch road level Equivalent Model of an electric system, detailed process is as follows:

In electromagnetic transient in power system emulation, the typical element in branch road comprises: generator, transmission line, resistance-inductance (RL) series arm, resistance capacitance (RC) series arm, regulator, load, transformer, bus, AC line, AC line, reactive-load compensator and shunt capacitance reactor etc.Quantity of state in electromagnetic transient in power system emulation generally refers to the capacitance voltage in resistance capacitance (RC) series arm and the inductive current in resistance-inductance (RL) series arm.

(1-1) the voltage-current characteristic equation setting electric system branch road is:

$\frac{\mathrm{di}\left(t\right)}{\mathrm{dt}}=f(i\left(t\right),u\left(t\right))$

Wherein, i (t), u (t) are respectively voltage time-domain signal and the electric current time-domain signal of port corresponding with each branch road, and f (i (t), u (t)) is a smooth function in real number field;

(1-2) according to above-mentioned voltage-current characteristic equation, the frequency domain equation of each branch road under electric current and voltage Laplace transform is obtained:

I(s)＝H(s)U(s)

Wherein, I (s), U (s) are respectively voltage frequency-region signal and the electric current frequency-region signal of port corresponding with each branch road, and H (s) is the transport function between the electric current and voltage corresponding with each branch road;

(1-3) transfer function H (s) of above-mentioned steps (1-2) is rewritten into zero pole point expression formula:

$H_a\left(s\right)=A\frac{\underset{s=1}{\overset{N}{\mathrm{\Π}}}(s-z_s)}{\underset{t=1}{\overset{M}{\mathrm{\Π}}}(s-p_t)},$

Wherein: A is the first proportional gain factor, z _sfor the zero point of transport function, p _tfor the limit of transport function, N is the zero number in transport function, s=1,2,3 ... N, M are the zero number in transport function, t=1,2,3 ... M;

(1-4) according to the expression formula of above-mentioned steps (1-3), according to Impulse invariance procedure principle, the transport function of each branch road under discrete domain is obtained as follows:

$H_d\left(z\right)=B\frac{\underset{s=1}{\overset{N}{\mathrm{\Π}}}(z-e^{z_{s}h})}{\underset{t=1}{\overset{M}{\mathrm{\Π}}}(z-e^{p_{t}h})}$

Wherein: B is the second proportional gain factor, z _sfor the zero point of transport function, p _tfor the limit of transport function, h is electromagnetic transient simulation material calculation, and N is the zero number in transport function, s=1,2,3 ... N, M are the zero number in transport function, t=1,2,3 ... M;

(1-5) transport function under the discrete domain of above-mentioned steps (1-4) is converted into the difference equation under time domain:

$i\left(n\right)=\mathrm{Cv}\left(n\right)+\underset{k=1}{\overset{s}{\mathrm{\Σ}}}{D}_{k}v(n-k)+\underset{k=1}{\overset{t}{\mathrm{\Σ}}}{E}_{k}i(n-k),$

Obtain the coefficient E of historical current item in the coefficient D of history voltage item in the coefficient C of current voltage item in difference equation, difference equation and difference equation,

Wherein, n is the n-th electromagnetic transient simulation step-length, k is the Difference Terms number of difference equation, v (n) is the voltage item of the n-th electromagnetic transient simulation, i (n) is the electric current item of the n-th electromagnetic transient simulation, v (n-k) is the history voltage item of the n-th-k electromagnetic transient simulation, and i (n-k) is the historical current item of the n-th-k electromagnetic transient simulation;

(1-6) difference equation of above-mentioned steps (1-5) is rewritten into computing conductance in parallel and historical current item form, the Equivalent Model obtaining each branch road of electric system is: i (n)=gv (n)+i _hist(n-1),

Wherein, g is computing specific conductance in parallel, i _hist(n-1) be the historical current item of (n-1)th electromagnetic transient simulation,

$\left\{\begin{array}{c}g=C\\ i_{\mathrm{hist}}(n-1)=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}{D}_{k}v(n-k)+\underset{k=1}{\overset{t}{\mathrm{\Σ}}}{E}_{k}i(n-k)\end{array}\right.;$

(2) load in electric system is judged, if without load in electric system, then carry out step (3), if there is nonlinear-load (comprising typical constant impedance load, steady current load and firm power load) in electric system, then carry out following steps:

(2-1) constant impedance setting electric system is Z _abc, Z _abcbe 3 × 3 matrixes, the steady current of electric system is I _abc, I _abcbe 3 × 1 column vectors, the firm power of electric system is , wherein P is the active power of firm power, and Q is the reactive power of firm power;

Determine power load and carry out special processing by following thinking.Being equivalent to constant impedance model with trend steady result, determining power load under disturbance situation, in order to keep power invariability, adopting time-dependent current model to compensate.

The equivalent impedance calculating firm power load is as follows:

Wherein for the complex power of firm power load under electric power system tide stable state, for conjugation, for firm power load bus voltage under electric power system tide stable state,

Above-mentioned equivalent impedance is converted to three phase of impedance: i ₃for 3*3 unit matrix, for the tensor (knoneck) of matrix amasss;

The offset current calculating firm power load is as follows:

(2-2) according to voltage time-domain signal and the electric current time-domain signal of electric system, the power difference of electric system firm power load is obtained $\mathrm{\Δ}\stackrel{\·}{S}\left(t\right):\mathrm{\Δ}\stackrel{\·}{S}\left(t\right)=\stackrel{\·}{S}-(p\left(t\right)+\mathrm{jq}\left(t\right))$

p(t)＝v _a(t)i _a(t)+v _b(t)i _b(t)+v _c(t)i _c(t)

Wherein, $q\left(t\right)=\frac{1}{\sqrt{3}}[(v_b\left(t\right)-v_c\left(t\right))i_a\left(t\right)+(v_c\left(t\right)-v_a\left(t\right))i_b\left(t\right)+(v_a\left(t\right)-v_b\left(t\right))i_c\left(t\right)]$

V _a(t), v _b(t), v _ct () is respectively the three-phase voltage time-domain signal of electric system firm power load, i _a(t), i _b(t), i _ct () is respectively the three-phase current time-domain signal of electric system firm power load, for the complex power of firm power load under electric power system tide stable state, j ²=-1;

(2-3) according to the power difference of above-mentioned electric system firm power load obtain the amplitude I of electric system firm power load compensation electric current _m(t):

$I_m\left(t\right)=\frac{2\left|\right|\mathrm{\Δ}\stackrel{\·}{S}\left(t\right)\left|\right|}{{3V}_m\left(t\right)}$

Wherein V _mt () is the voltage magnitude of firm power load bus:

$V_m\left(t\right)=\sqrt{\frac{2(v_a^2\left(t\right)+v_b^2\left(t\right)+v_c^2\left(t\right))}{3}}$

V _a(t), v _b(t), v _ct () is respectively the three-phase voltage time-domain signal of electric system firm power load;

(2-4) according to the result of calculation of above-mentioned steps (2-2) and step (2-3), the three-phase offset current of firm power load is obtained:

ω is electric system angular frequency, be respectively the initial phase of firm power load three-phase current time-domain signal, ω, with can be obtained by the Load flow calculation of electric system;

(2-5) according to the result of calculation of above-mentioned steps (2-4), obtain computing conductance in parallel and the historical current item of nonlinear-load in electric system, wherein computing impedance in parallel is by Z _abcand Z' _abcparallel connection obtains, and computing conductance in parallel is the inverse of computing impedance in parallel.The electric current of historical current item is by I _abcwith Δ I _abcparallel connection obtains, wherein Δ I _abcfor constant impedance offset current item, Δ I _abc=[Δ i _aΔ i _bΔ i _c] ^t, T is matrix transpose, Z _abcfor the constant impedance of electric system, I _abcfor the steady current of electric system;

(3) set up the nodal voltage equation of an electric system, comprise the following steps:

(3-1) according to the computing specific conductance g in parallel of above-mentioned steps (1-6), the branch road conductance matrix G of electric system is obtained,

$G=\tilde{A}\mathrm{diag}\left(g_i\right){\tilde{A}}^T$

Wherein for the topological correlation matrix of electric system, the element in topological correlation matrix is ± 1 or 0, g _ifor the computing specific conductance in parallel of bar branch road every in electric system, diag (g _i) diagonal matrix that forms for the computing specific conductance in parallel of every bar branch road;

(3-2) the historical current item of each branch road in the electric system obtained according to above-mentioned steps (1-6), obtains branch node historical current vector I _hist(n):

$I_{\mathrm{hist}}(n-1)=\tilde{A}\mathrm{diag}\left\{i_{\mathrm{hist}}(n-1)\right\};$

Wherein, diag{i _hist(n-1) } be diagonal matrix that the historical current item of every bar branch road is formed;

(3-3) according to the result of calculation of above-mentioned steps (3-1) and step (3-2), the nodal voltage equation of electric system is obtained:

GV(n)＝I(n)+I _hist(n-1)

V (n) is the voltage vector of electric system node in the n-th electromagnetic transient simulation step-length, solve this nodal voltage equation and obtain V (n), I (n) is the current vector of electric system node in the n-th electromagnetic transient simulation, this current vector is a known quantity, can obtain from the operational factor of electric system;

(4) implicit trapezoid method is utilized, the inertia branch model (comprising generator, motor etc. in inertia branch model) of the Synchronous Machine Models of electric system, electric system and generator excitation, speed governing are incorporated in the Electric Power System Node Voltage equation of above-mentioned steps (3) together with prime mover model, obtain an expanding node voltage equation;

(5) solve the expanding node voltage equation of above-mentioned steps (4), obtain a, b, c three-phase node voltage of electric system, and according to three-phase node voltage, calculate the Current calculation value in electric system generator model;

(6) the current error threshold value of a generator model is set, the Current calculation value of the electric current predictand in electric system generator model and above-mentioned steps (5) is compared, if error is less than or equal to error threshold, then obtain the electromagnetic transient simulation result of current state simulation step length, carry out step (7), if error is greater than error threshold, then repeat step (4)-step (6), when returning step (4), keep the terminal voltage of the node be connected with generator constant;

(7) repeat step (2)-step (6), obtain the electromagnetic transient simulation result of the full step-length of electric system.

0+ error immunity electromagnetic transient simulation algorithm under the input quantity catastrophe that the present invention proposes, the typical branch for electric system is A stable algorithm.For RC series arm, Equivalent Model equation is as follows:

$\left\{\begin{array}{c}G=\frac{C(1-e^{-h/\mathrm{RC}})}{h}\\ i_{\mathrm{hist}}(n-1)=e^{-h/\mathrm{RC}}i(n-1)-\frac{C(1-e^{-h/\mathrm{RC}})}{h}u(n-1)\end{array}\right.$

Wherein: h is electromagnetic transient simulation material calculation, u (n-1) is the voltage item of (n-1)th electromagnetic transient simulation, and i (n-1) is the electric current item of (n-1)th electromagnetic transient simulation.

Can obtain corresponding equation Stable Polynomials is:

$p\left(r\right)=r-\frac{C(1-e^{-h/\mathrm{RC}})}{h}\mathrm{r\λ}+e^{-h/\mathrm{RC}}-\frac{C(1-e^{-h/\mathrm{RC}})}{h}\mathrm{\λ}$

Characteristic root is:

$r=\frac{e^{-h/\mathrm{RC}}-\frac{C(1-e^{-h/\mathrm{RC}})}{h}\mathrm{\λ}}{1-\frac{C(1-e^{-h/\mathrm{RC}})}{h}\mathrm{\λ}}$

According to stability distinguishing condition | r| < 1 (as λ < 0) this algorithm known is A stable algorithm.

In the electromagnetic transient in power system emulation mode of the quantity of state jump error immunity that the present invention proposes, steady-state error is much less than traditional trapezoidal integration, and error coefficient is multiplied by Ah, and wherein h is electromagnetic transient simulation material calculation.

For RL series arm, if the Fourier expansion of voltage and current signal is as follows:

Wherein i (t), u (t) are respectively voltage time-domain signal and the electric current time-domain signal of corresponding ports, I _k, U _kbe respectively the voltage time-domain signal of corresponding ports and the Fourier series coefficient of electric current time-domain signal; be respectively the initial phase of the voltage time-domain signal of corresponding ports and the Fourier series coefficient of electric current time-domain signal.

The ratio of defining ideal calculated value and numerical computation is error coefficient e (ω).Can prove that the error coefficient of implicit trapezoid method and this method is as follows respectively:

$e_{\mathrm{ti}}\left(\mathrm{\&omega;}\right)=\frac{\sin \frac{\mathrm{\&omega;h}}{2}}{\frac{\mathrm{\&omega;h}}{2}}-\cos \frac{\mathrm{\&omega;h}}{2}$

$e_{\mathrm{ch}}\left(\mathrm{\&omega;}\right)=\frac{h}{L}(\frac{\sin \frac{\mathrm{\&omega;h}}{2}}{\frac{\mathrm{\&omega;h}}{2}}-\cos \frac{\mathrm{\&omega;h}}{2})+O\left(h^2\right)$

Wherein ω is electric system angular frequency, and h is electromagnetic transient simulation material calculation, and L is this shunt inductance parameter size.O (h ²) represent h ²higher-order shear deformation.

0+ error immunity electromagnetic transient simulation algorithm under the input quantity catastrophe that the present invention proposes, can the non-prototype oscillatory occurences brought of immune state amount saltus step amount.Can prove that implicit trapezoid method can produce the non-prototype vibration of decay, even if adopt the measures such as interpolation, this oscillatory occurences can only suppress to eliminate, and this is the advantage of the inventive method.

For resistance capacitance series loop system, suppose under the n-th simulation step length, quantity of state has a sudden change u' _n=u _n+ δ.Wherein u _nbe the capacitance voltage of the n-th electromagnetic transient simulation, δ is the capacitance voltage saltus step amount of the n-th electromagnetic transient simulation, definition status amount jump error be respectively capacitance voltage and the capacitance current quantity of state error of the n-th electromagnetic transient simulation.Can prove that the error recursion formula of this method is as follows:

${\mathrm{\&sigma;}}_{n+2}^u={\mathrm{\&sigma;}}_{n+1}^u+\frac{h}{C(1-e^{-\frac{h}{\mathrm{RC}}})}[{\mathrm{\&sigma;}}_{n+2}^i-e^{-\frac{h}{\mathrm{RC}}}{\mathrm{\&sigma;}}_{n+1}^i]=0$

${\mathrm{\&sigma;}}_{n+2}^i={i^{\&prime;}}_{n+2}-i_{n+2}=e^{-\frac{h}{\mathrm{RC}}}{\mathrm{\&sigma;}}_{n+1}^i$

As shown from the above formula, state variable walked in lower a period of time, and quantity of state error is 0, and non-state variable decays rapidly.

In an embodiment of the inventive method, electromagnetic transient analysis carries out to the resistance-capacitance network node of electric system as follows:

Resistance capacitance series loop system simulation capacitance state amount saltus step situation as shown in Figure 2, adopts the inventive method to simulate electromagnetic transient, analyzes the method and the difference of classic method in quantity of state error and truncation error.

Step 1: the branch road level Equivalent Model setting up resistance capacitance series loop system.

The concrete time domain external characteristics equation first writing out the electric current and voltage of branch road, does branch model equivalence afterwards.

The voltage-current characteristic equation writing out branch road is:

$i=C\frac{\mathrm{du}}{\mathrm{dt}}-\mathrm{CR}\frac{\mathrm{di}}{\mathrm{dt}}$

According to above-mentioned voltage-current characteristic equation, obtain the frequency domain equation of branch road under electric current and voltage Laplace transform:

Transform to s domain equation corresponding as follows:

$I\left(s\right)=\frac{\mathrm{Cs}}{1+\mathrm{CRs}}U\left(s\right)$

Can finally push away RC branch road time-domain expression is as follows according to zero pole point correspondence rule:

i(n)＝Gv(n)+i _hist(n-1)

Wherein:

$\left\{\begin{array}{c}G=\frac{C(1-e^{-h/\mathrm{RC}})}{h}\\ i_{\mathrm{hist}}(n-1)=e^{-h/\mathrm{RC}}i(n-1)-\frac{C(1-e^{-h/\mathrm{RC}})}{h}u(n-1)\end{array}\right.$

Wherein: h is electromagnetic transient simulation material calculation, u (n-1) is the voltage item of (n-1)th electromagnetic transient simulation, and i (n-1) is the electric current item of (n-1)th electromagnetic transient simulation.

Step 2: because without load is skipped.

Step 3: system node establishing equation: adopt implicit trapezoid method and this method (abbreviation index method) and theoretical curve comparative result respectively.Result is as Fig. 3.Draw various algorithm in figure respectively and calculate capacitance voltage curves result.

In figure, implicit trapezoid method is the practical approach curve simultaneously using damping resistance and method of interpolation.As can be seen from the figure: improve one's methods and can suppress but fundamentally can not eliminate numerical oscillation.This method has great advantage at the party's mask.

Error analysis 1: truncation error analysis.

If the Fourier expansion of voltage and current signal is as follows:

Wherein i (t), u (t) are respectively voltage time-domain signal and the electric current time-domain signal of corresponding ports, I _k, U _kbe respectively the voltage time-domain signal of corresponding ports and the Fourier series coefficient of electric current time-domain signal; be respectively the initial phase of the voltage time-domain signal of corresponding ports and the Fourier series coefficient of electric current time-domain signal.

The ratio of defining ideal calculated value and numerical computation is error coefficient, is designated as: e (ω).Can prove that the error coefficient of implicit trapezoid method and this method is as follows respectively:

$e_{\mathrm{ti}}\left(\mathrm{\&omega;}\right)=\frac{\sin \frac{\mathrm{\&omega;h}}{2}}{\frac{\mathrm{\&omega;h}}{2}}-\cos \frac{\mathrm{\&omega;h}}{2}$

$e_{\mathrm{ch}}\left(\mathrm{\&omega;}\right)=\frac{h}{L}(\frac{\sin \frac{\mathrm{\&omega;h}}{2}}{\frac{\mathrm{\&omega;h}}{2}}-\cos \frac{\mathrm{\&omega;h}}{2})+O\left(h^2\right)$

Wherein ω is electric system angular frequency, and h is electromagnetic transient simulation material calculation, and L is this shunt inductance parameter size.O (h ²) represent h ²higher-order shear deformation.

Error analysis 2: quantity of state error analysis.

For resistance capacitance series loop system, suppose under the n-th simulation step length, quantity of state has a sudden change u' _n=u _n+ δ.Wherein u _nbe the capacitance voltage of the n-th electromagnetic transient simulation, δ is the capacitance voltage saltus step amount of the n-th electromagnetic transient simulation, definition status amount jump error be respectively capacitance voltage and the capacitance current quantity of state error of the n-th electromagnetic transient simulation.Can prove that the error recursion formula of this method is as follows:

${\mathrm{\&sigma;}}_{n+2}^u={\mathrm{\&sigma;}}_{n+1}^u+\frac{h}{C(1-e^{-\frac{h}{\mathrm{RC}}})}[{\mathrm{\&sigma;}}_{n+2}^i-e^{-\frac{h}{\mathrm{RC}}}{\mathrm{\&sigma;}}_{n+1}^i]=0$

${\mathrm{\&sigma;}}_{n+2}^i={i^{\&prime;}}_{n+2}-i_{n+2}=e^{-\frac{h}{\mathrm{RC}}}{\mathrm{\&sigma;}}_{n+1}^i$

As shown from the above formula, state variable walked in lower a period of time, and quantity of state error is 0, and non-state variable decays rapidly.

Concrete matlab numerical result and theoretical value contrast as follows:

Can find out do not have quantity of state error, error is the basic truncation error of algorithm.

Can prove that implicit trapezoid method quantity of state error recurrence formula is as follows:

${\mathrm{\&sigma;}}_{k+2}^u={\mathrm{\&sigma;}}_{k+1}^u+\frac{\mathrm{\&Delta;t}}{2C}[{\mathrm{\&sigma;}}_{k+2}^i+{\mathrm{\&sigma;}}_{k+1}^i]$

${\mathrm{\&sigma;}}_{k+2}^i=\frac{(R-\frac{h}{2C})}{(R+\frac{h}{2C})}{\mathrm{\&sigma;}}_{k+1}^i$

Actual matlab calculate data and theoretical calculation data as follows: (starting during N=2000 to break down)

Voltage error theoretical calculation formula compares with MATLAB numerical result:

Current error theoretical calculation formula compares with MATLAB numerical result:

Can find out that the calculated results and maltab numerical result are identical.The method has the advantage of immune state amount saltus step than implicit trapezoid method.

Claims (1)

1. the 0+ error immunity electromagnetic transient simulation algorithm under input quantity catastrophe, it is characterized in that, the method comprises the following steps:

(1) set up the branch road level Equivalent Model of an electric system, detailed process is as follows:

(1-1) the voltage-current characteristic equation setting electric system branch road is:

Wherein, i (t), u (t) are respectively voltage time-domain signal and the electric current time-domain signal of port corresponding with each branch road, and f (i (t), u (t)) is a smooth function in real number field;

(1-2) according to above-mentioned voltage-current characteristic equation, the frequency domain equation of each branch road under electric current and voltage Laplace transform is obtained:

I(s)＝H(s)U(s)

Wherein, I (s), U (s) are respectively voltage frequency-region signal and the electric current frequency-region signal of port corresponding with each branch road, and H (s) is the transport function between the electric current and voltage corresponding with each branch road;

(1-3) transfer function H (s) of above-mentioned steps (1-2) is rewritten into zero pole point expression formula:

Wherein: A is the first proportional gain factor, z _sfor the zero point of transport function, p _tfor the limit of transport function, N is the zero number in transport function, s=1,2,3 ... N, M are the zero number in transport function, t=1,2,3 ... M;

(1-4) according to the expression formula of above-mentioned steps (1-3), according to Impulse invariance procedure principle, the transport function of each branch road under discrete domain is obtained as follows:

Wherein: B is the second proportional gain factor, z _sfor the zero point of transport function, p _tfor the limit of transport function, h is electromagnetic transient simulation material calculation, and N is the zero number in transport function, s=1,2,3 ... N, M are the zero number in transport function, t=1,2,3 ... M;

(1-5) transport function under the discrete domain of above-mentioned steps (1-4) is converted into the difference equation under time domain:

Obtain the coefficient E of historical current item in the coefficient D of history voltage item in the coefficient C of current voltage item in difference equation, difference equation and difference equation;

Wherein, n is the n-th electromagnetic transient simulation step-length, and k is the Difference Terms number of difference equation; V (n) is the voltage item of the n-th electromagnetic transient simulation; I (n) is the electric current item of the n-th electromagnetic transient simulation, and v (n-k) is the history voltage item of the n-th-k electromagnetic transient simulation; I (n-k) is the historical current item of the n-th-k electromagnetic transient simulation;

(1-6) difference equation of above-mentioned steps (1-5) is rewritten into computing conductance in parallel and historical current item form, the Equivalent Model obtaining each branch road of electric system is: i (n)=gv (n)+i _hist(n-1),

Wherein, g is computing specific conductance in parallel, i _hist(n-1) be the historical current item of (n-1)th electromagnetic transient simulation,

(2) load in electric system is judged, if without load in electric system, then carries out step (3), if having nonlinear-load in electric system, then carry out following steps:

(2-1) constant impedance setting electric system is Z _abc, Z _abcbe 3 × 3 matrixes, the steady current of electric system is I _abc, I _abcbe 3 × 1 column vectors, the firm power of electric system is wherein P is the active power of firm power, and Q is the reactive power of firm power;

The equivalent impedance calculating firm power load is as follows:

Wherein for the complex power of firm power load under electric power system tide stable state, for conjugation, for firm power load bus voltage under electric power system tide stable state,

Above-mentioned equivalent impedance is converted to three phase of impedance: i ₃for 3*3 unit matrix, for the tensor product of matrix;

The offset current calculating firm power load is as follows:

(2-2) according to voltage time-domain signal and the electric current time-domain signal of electric system, the power difference of electric system firm power load is obtained

Wherein,

V _a(t), v _b(t), v _ct () is respectively the three-phase voltage time-domain signal of electric system firm power load, i _a(t), i _b(t), i _ct () is respectively the three-phase current time-domain signal of electric system firm power load, for the complex power of firm power load under electric power system tide stable state, j ²=-1;

(2-3) according to the power difference of above-mentioned electric system firm power load obtain the amplitude I of electric system firm power load compensation electric current _m(t):

Wherein V _mt () is the voltage magnitude of firm power load bus:

V _a(t), v _b(t), v _ct () is respectively the three-phase voltage time-domain signal of electric system firm power load;

(2-4) according to the result of calculation of above-mentioned steps (2-2) and step (2-3), the three-phase offset current of firm power load is obtained:

ω is electric system angular frequency, be respectively the initial phase of firm power load three-phase current time-domain signal, ω, with obtained by the Load flow calculation of electric system;

(2-5) according to the result of calculation of above-mentioned steps (2-4), obtain computing conductance in parallel and the historical current item of nonlinear-load in electric system, wherein computing impedance in parallel is by Z _abcand Z' _abcparallel connection obtains, and computing conductance in parallel is the inverse of computing impedance in parallel.The electric current of historical current item is by I _abcwith Δ I _abcparallel connection obtains, wherein Δ I _abcfor constant impedance offset current item, Δ I _abc=[Δ i _aΔ i _bΔ i _c] ^t, T is matrix transpose, Z _abcfor the constant impedance of electric system, I _abcfor the steady current of electric system;

(3) set up the nodal voltage equation of an electric system, comprise the following steps:

(3-1) according to the computing specific conductance g in parallel of above-mentioned steps (1-6), the branch road conductance matrix G of electric system is obtained,

Wherein for the topological correlation matrix of electric system, the element in topological correlation matrix is ± 1 or 0, g _ifor the computing specific conductance in parallel of bar branch road every in electric system, diag (g _i) diagonal matrix that forms for the computing specific conductance in parallel of every bar branch road;

(3-2) the historical current item of each branch road in the electric system obtained according to above-mentioned steps (1-6), obtains branch node historical current vector I _hist(n):

Wherein, diag{i _hist(n-1) } be diagonal matrix that the historical current item of every bar branch road is formed;

(3-3) according to the result of calculation of above-mentioned steps (3-1) and step (3-2), the nodal voltage equation of electric system is obtained:

GV(n)＝I(n)+I _hist(n-1)

V (n) is the voltage vector of electric system node in the n-th electromagnetic transient simulation step-length, solves this nodal voltage equation and obtains V (n), and I (n) is the current vector of electric system node in the n-th electromagnetic transient simulation;

(4) implicit trapezoid method is utilized, the inertia branch model of the Synchronous Machine Models of electric system, electric system, generator excitation, speed governing are incorporated to together with prime mover model in the Electric Power System Node Voltage equation of above-mentioned steps (3), obtain an expanding node voltage equation;

(5) solve the expanding node voltage equation of above-mentioned steps (4), obtain a, b, c three-phase node voltage of electric system, and according to three-phase node voltage, calculate the Current calculation value in electric system generator model;

(6) the current error threshold value of a generator model is set, the Current calculation value of the electric current predictand in electric system generator model and above-mentioned steps (5) is compared, if error is less than or equal to error threshold, then obtain the electromagnetic transient simulation result of current state simulation step length, carry out step (7), if error is greater than error threshold, then repeat step (4)-step (6), when returning step (4), keep the terminal voltage of the node be connected with generator constant;

(7) repeat step (2)-step (6), obtain the electromagnetic transient simulation result of the full step-length of electric system.