CN103401238A - Method for modeling power load based on measurement-based method - Google Patents

Abstract

The invention discloses a method for modeling a power load based on a measurement-based method. The method sequentially comprises the following steps of (1) preprocessing recording data to obtain data required by load modeling; (2) performing damped least-squares iterative equation solving on the preprocessed data required by the load modeling to obtain optimal load model parameters; (3) establishing an optimal load model by adopting the load model parameters; and (4) checking the obtained load model parameters, and verifying a load model with practical values comprising higher extrapolation capability, higher interpolation capability and higher capability in the comprehensive description of different load compositions to be the optimal load model. According to the method, the modeling cost is reduced; and the load model is more consistent with the actual condition of a field, high simulation accuracy and credibility are ensured, contrary conclusions under a clinical condition are completely avoided, the load model can be coordinated with the conventional power generation and distribution model, and the analysis of a power system is accurate and credible, so that the simulation analysis accuracy and credibility of the power system are remarkably improved.

Description

A kind of power load modelling approach of distinguishing method based on total body examination

Technical field

The present invention relates to Modeling for Electric Loads, particularly relate to a kind of power load modelling approach of distinguishing method based on total body examination.

Background technology

Power system load model is the power of reflection practical power systems load ports and electric current math equation and the corresponding parameter with its port voltage and frequency variation characteristics.Load modeling is not only set up model to various concrete power consumption equipment elements, the more important thing is the relation that the power of the overall load absorption on research load bus changes with load busbar voltage and frequency variation, and describe the form of the math equation of determining this relation and parameter wherein.The Digital Simulation of modern power systems design, planning, decision-making and operation, all need the load model that accuracy is higher.Because the result that the variation of load model is calculated system transient modelling, voltage stabilization and trend has impact in various degree, the variation of matter even occurs under critical condition.And load modeling not only will obtain with one group of data fitting the structure and parameter of model, also must be with the test of echoing of other data, and comprehensive and to the robustness of noise to verify it.Load model by the measurement data identification under a certain change in voltage remains correct under larger or less change in voltage, be extrapolation and the interpolation of load model, have the load model of the ability of extrapolability, interpolation ability and comprehensive description different load composition preferably, just with practical value.

Existing load model is from basic conception mostly, and adopting Utopian model is that permanent power, constant-impedance, constant current or three's combination is as integrated load model.This coarse load model is seriously inharmonious with accurate generating and distribution model, has obviously reduced accuracy and the confidence level of emulation, under critical condition, even can draw diametrically opposite conclusion.So far there is not yet based on recorder data and belong to the load modeling method that method is distinguished in total body examination.

Summary of the invention

Technical problem to be solved by this invention is the defect that makes up above-mentioned prior art, and a kind of power load modelling approach of distinguishing method based on total body examination is provided.

Technical problem of the present invention is solved by the following technical programs.

This power load modelling approach of distinguishing method based on total body examination has following steps successively:

1) by recorder data is carried out to preliminary treatment, obtain the load modeling desired data;

Described recorder data is three-phase voltage instantaneous value and the three-phase current instantaneous value of load modeling Nodes, i.e. three-phase voltage u _a, u _b, u _cAnd three-phase current i _a, i _b, i _c, comprise fault recorder data, namely comprise the data of at least 2 cycles of presteady state process, the whole obvious disturbance transient process after fault, and the tend to be steady data of process or rear at least 5 cycles of steady-state process of disturbance;

2) the load modeling desired data of selecting levenberg-marquart algorithm as Identification of parameter, preliminary treatment to be obtained carries out the damped least squares iterative equation and solves, and obtains optimum load model parameters;

3) adopt load model parameters to set up optimum load model, load model that can identified parameters comprises the integrated load model of static load model, dynamic load model, integrated load model, consideration distribution, and the static load model of considering distribution;

4) load model parameters that obtains is carried out to the parameter check, checking has practical value and comprises the load model of the ability of extrapolability, interpolation ability and comprehensive description different load composition preferably, is optimum load model.

Technical problem of the present invention is solved by following further technical scheme.

The sample frequency of fault recorder data described step 1) is every cycle at least 20 points, namely is at least 1KHz, and sampling precision is that amplitude error is at most 1%, described three-phase voltage u _a, u _b, u _cAnd three-phase current i _a, i _b, i _cSynchronism is good, and phase error approaches zero, if there is phase error in the instantaneous value collection, can bring very large impact to subsequent calculations.

Recorder data carried out to preliminary treatment comprise described step 1):

Conversion perunit value: the ratio of the actual value of conversion three-phase voltage instantaneous value and three-phase current instantaneous value and the fixed numbers of selected commensurate;

Extract positive sequence component: utilize symmetrical component method to extract the positive sequence component of voltage, the magnitude of current, during stable state, there are following relation in the 0-1-2 component of electric weight and a-b-c component:

In formula:

Operator $\mathrm{\α}=e^{{j120}^0};$

The phasor value of a phase electric weight, x is voltage u or current i;

The phasor value of b phase electric weight, x is voltage u or current i;

The phasor value of c phase electric weight, x is voltage u or current i;

Subscript 0 expression zero-sequence component;

Subscript 1 expression positive sequence component;

Subscript 2 expression negative sequence components;

High-frequency noise in the smothing filtering data: adopt the balance filtering algorithm to step 2) positive sequence voltage and the forward-order current that obtain carry out filtering, the high-frequency noise that contains in the busbar voltage amplitude V that removal calculates, active-power P, reactive power Q, improve signal to noise ratio; Described balance filtering algorithm comprises 5 secondary balancing filtering algorithms, and the computing formula of 5 secondary balancing filtering algorithms is as follows:

$\left\{\begin{array}{c}y\left(n\right)=[-3x(n-2)+12x(n-1)+17x(n+1)-3x(n+2)]/35\\ y\left(0\right)=[69x\left(0\right)+4x\left(1\right)+4x\left(3\right)-x\left(4\right)]/60\\ y\left(1\right)=[2x\left(0\right)+27x\left(1\right)+12x\left(2\right)-8x\left(3\right)+2x\left(4\right)]/35\\ y(N-2)=[2x(N-5)-8x(N-4)+12x(N-3)+27x(N-2)]+2x(N-1)]/35\\ y(N-1)=[-x(N-5)+4x(N-4)-6x(N-3)+4x(N-2)+69x(N-1)]/70\end{array}\right.$

2≤n≤N-3

In formula:

Y (n): the signal after level and smooth;

X (n): the signal before level and smooth;

N: data length;

Quadrature Park conversion: adopt quadrature park transformation for mula to carry out conversion to positive sequence voltage and forward-order current after level and smooth, obtain d shaft voltage u _d, q shaft voltage u _q, d shaft current i _dWith q shaft current i _qFundamental formular is as follows:

$\left[\begin{array}{c}{f}_d\\ f_q\\ f_0\end{array}\right]=\sqrt{\frac{2}{3}}\left[\begin{array}{ccc}{\cos \mathrm{\θ}}_a& {\cos \mathrm{\θ}}_b& {\cos \mathrm{\θ}}_c\\ -{\sin \mathrm{\θ}}_a& {-\sin \mathrm{\θ}}_b& {-\sin }_c\\ 1/\sqrt{2}& 1/\sqrt{2}& 1/\sqrt{2}\end{array}\right]\left[\begin{array}{c}{f}_a\\ f_b\\ f_c\end{array}\right];$

In formula:

θ _a＝ωt；

f _a: phase voltage u _aWith phase current i _aThe instantaneous value sequence;

f _b: phase voltage u _bWith phase current i _bThe instantaneous value sequence;

f _c: phase voltage u _cWith phase current i _cThe instantaneous value sequence;

Fundamental formular by quadrature Park conversion calculates:

$\left\{\begin{array}{c}{u}_d=\sqrt{\frac{2}{3}}(u_a{\cos \mathrm{\θ}}_a+u_b{\cos \mathrm{\θ}}_b+u_c{\cos \mathrm{\θ}}_c)\\ u_q=\sqrt{\frac{2}{3}}(-u_a{\sin \mathrm{\θ}}_a-u_b{\sin \mathrm{\θ}}_b-u_c{\sin \mathrm{\θ}}_c)\\ i_d=\sqrt{\frac{2}{3}}(i_a{\cos \mathrm{\θ}}_a+i_b{\cos \mathrm{\θ}}_b+i_c{\cos \mathrm{\θ}}_c)\\ i_q=\sqrt{\frac{2}{3}}(-i_a{\sin \mathrm{\θ}}_a-i_b{\sin \mathrm{\θ}}_b-i_c{\sin \mathrm{\θ}}_c)\end{array}\right.$

In formula:

u _d: d axle instantaneous voltage data sequence;

u _q: q axle instantaneous voltage data sequence;

i _d: d axle transient current data sequence;

i _q: q axle transient current data sequence;

The calculated load busbar voltage: fundamental formular is as follows:

$V=\sqrt{{u_d}^2+{u_q}^2};$

In formula:

u _d: d axle instantaneous voltage data sequence;

u _q: q axle instantaneous voltage data sequence;

V: load busbar voltage amplitude data sequence;

Calculate active power: fundamental formular is as follows:

P＝u _di _d+u _qi _q；

In formula:

u _d: d axle instantaneous voltage data sequence;

i _d: d axle transient current data sequence;

u _q: q axle instantaneous voltage data sequence;

i _q: q axle transient current data sequence;

P: actual measurement active power data sequence;

Calculate reactive power: fundamental formular is as follows:

Q＝u _qi _d-u _di _q；

In formula:

u _d: d axle instantaneous voltage data sequence;

i _d: d axle transient current data sequence;

u _q: q axle instantaneous voltage data sequence;

i _q: q axle transient current data sequence;

Q: actual measurement reactive power data sequence.

Described step 2) obtain optimum load model parameters, comprising:

The measured value x of identification objects electric load model _i, y _iThere are two groups: (V, P), (V, Q), wherein V is actual measurement load busbar voltage, and P is actual measurement load active power, and Q is actual measurement reactive load power, y (x)=y (x; a ₀, a ₁..., a _M-1) be the functional relation of surveying independent variable and dependent variable of electric load model, wherein x _i, y _iFor measured value, a _jFor parameter to be identified; Utilize the N group data (x that has recorded _i, y _i), i=0, L, N-1, adopt levenberg-marquart algorithm to estimate unknown parameter a=(a ₀, a ₁, L, a _M-1) ^T, make target function J (a) (residual sum of squares (RSS)) minimum, namely

Minimum;

If y=is (y ₀, y ₁, L, y _N-1) ^TAnd linear y=Aa+b between a, directly solve a=(A with linear least square ^TA) ^-1A ^T(y-b);

If y is (x; a ₀, L, a _M-1) be the non-linear form of a, can only first by one group of initial solution of a, set out and carry out iterative computation, namely with nonlinear least square method, solve a;

If A (a) is function y (x)=y (x; a ₀, L, a _M-1) the Jacobian matrix, namely

$A\left(a\right)=\left[\begin{array}{cccc}\frac{{\∂y}_0}{{\∂a}_0}& \frac{{\∂y}_0}{{\∂a}_1}& L& \frac{{\∂y}_0}{{\∂a}_{M-1}}\\ \frac{{\∂y}_1}{{\∂a}_0}& \frac{{\∂y}_1}{{\∂a}_1}& L& \frac{{\∂y}_1}{{\∂a}_{M-1}}\\ M& M& & M\\ \frac{{\∂y}_{N-1}}{{\∂\mathrm{\α}}_0}& \frac{{\∂y}_{N-1}}{{\∂a}_1}& L& \frac{{\∂y}_{N-1}}{{\∂a}_{M-1}}\end{array}\right];$

Described levenberg-marquart algorithm namely solves following system of linear equations:

$(A^{T}A+\mathrm{\λI})\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ M\\ y_{N-1}-y(x_{N-1};a)\end{array}\right);$

In formula:

$\mathrm{\λI}\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ y_{N-1}-y(x_{N-1};a)\end{array}\right)$ It is the system of linear equations that solves of steepest descent method;

$\left(A^{T}A\right)\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ M\\ y_{N-1}-y(x_{N-1};a)\end{array}\right)$ It is the system of linear equations that solves of Gauss-Newton method;

Wherein, damping factor λ=10 ^-4～10 ^-2

When λ>0 is too large, may occur that program does steepest at the little place, the lowest point of descending grade and descend, { the a} convergence rate descends to make sequence;

When λ too hour, convergence domain is too small, initial approximation a ₀Restricted;

Levenberg-marquart algorithm uses steepest descent method away from minimum point the time, selects larger λ value;

Levenberg-marquart algorithm is switched to Gauss-Newton method gradually when near minimum value, reduces λ value.

Described step 2) damped least squares iterative equation solves, and comprising:

Input solve for parameter initial value scope;

According to the motor steady-state equation, ask the initial value of induction motor state equation;

Utilize 4 rank Runge Kuttas of fixed step size to ask the electric motor state equation, i.e. the array of the differential equation;

Utilize initial parameter value to calculate active power and the reactive power of each data point, form the model output of vector form;

Solve following damped least squares iterative equation:

$(A^{T}A+\mathrm{\λI})\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ M\\ y_{N-1}-y(x_{N-1};a)\end{array}\right),$

And judge whether to meet the condition of convergence:

If do not met convergence, the Returning utilization initial parameter value continues to calculate active power and the reactive power of each data point;

As met convergence, the result of output global convergence, and whether judgement convergence result meets the constrained condition again, namely the load model parameters of motor can not be less than 0;

If do not met constraints, return to input solve for parameter initial value scope, automatically change initial parameter value identification again;

As meet constraints, output model parameter identification value.

Described step 3) static load model means namely the load model of functional relation of busbar voltage V of actual measurement active-power P and actual measurement reactive power Q and synchronization lower node voltage magnitude.

$\left\{\begin{array}{c}P=P_0(a_p{\left(\frac{V}{V_0}\right)}^2+b_p\left(\frac{V}{V_0}\right)+c_p)\\ Q=Q_0(a_q{\left(\frac{V}{V_0}\right)}^2+b_q\left(\frac{V}{V_0}\right)+c_q)\end{array};\right.$

a _p+b _p+c _p＝1；

a _q+b _q+c _q＝1；

In formula:

Coefficient a _p: constant-impedance (Z) component occupies the percentage of merit power;

Coefficient b _p: constant current (I) component occupies the percentage of merit power;

Coefficient c _p: permanent power (P) component occupies the percentage of merit power;

Coefficient a _q: constant-impedance (Z) component accounts for the percentage of reactive power;

Coefficient b _q: constant current (I) component accounts for the percentage of reactive power;

Coefficient c _q: permanent power (P) component accounts for the percentage of reactive power.

Described dynamic load model refers to the three order induction motor model,

Its system state equation is as follows:

$\left\{\begin{array}{c}\frac{d{E}_q^{\′}}{\mathrm{dt}}=-w_{B}s{E}_d^{\′}-\frac{1}{T_0^{\′}}{E}_q^{\′}+\frac{X-X^{\′}}{T_0^{\′}}{i}_d\\ \frac{d{E}_d^{\′}}{\mathrm{dt}}=w_{B}s{E}_q^{\′}-\frac{1}{T_0^{\′}}{E}_d^{\′}-\frac{X-X^{\′}}{T_0^{\′}}{i}_q\\ \frac{{\mathrm{dw}}_m}{\mathrm{dt}}=\frac{1}{2H}(T_e-T_m)\end{array}\right.;$

In formula:

E ' _d: the transient potential of d axle;

E ' _q: the transient potential of q axle;

X: synchronous reactance;

X ': transient state reactance;

T ' ₀: the transient state open circuit time constant;

w _m: rotor speed;

w _s: synchronous speed;

i _q: the electric current in transient circuit is at the component of q axle;

i _d: the electric current in transient circuit is at the component of d axle;

T _s: be loaded into epitrochanterian electromagnetic torque;

T _m: the machine torque that rotor bears;

H: inertia time constant;

As follows through the power stage equation that derivation draws:

$\left\{\begin{array}{c}P=\frac{X_{t1}}{G}({u^2}_d-u_d{E}_d^{\′}+{u^2}_q-u_q{E}_q^{\′})+\frac{X_{t2}}{G}(u_q{E}_d^{\′}-u_d{E}_q^{\′})\\ Q=\frac{X_{t1}}{G}(u_d{E}_q^{\′}-u_q{E}_d^{\′})+\frac{X_{t2}}{G}({u^2}_d-u_d{E}_d^{\′}+{u^2}_q-u_q{E}_q^{\′})\end{array}\right.;$

In formula:

$X_{t1}=\frac{R_s(X-X^{\′})}{{R_s}^2+X^{\′2}};$

$X_{t2}=\frac{X^{\′}(X-X^{\′})}{{R_s}^2+X^{\′2}};$

G＝X-X′；

X＝X _s+X _m；

$X^{\′}=X_s+\frac{X_m{X}_r}{X_m+X_r}.$

Described step 3) integrated load model is the comprehensive load model that comprises static load model and dynamic load model.

The integrated load model of consideration distribution described step 3) is to have considered equivalent distribution network, reactive power compensation, and virtual bus is set between transformer and distribution network

Integrated load model, described virtual bus

With the actual load bus

Between be the equivalent impedance of transmission and distribution networks.

The static load model of consideration distribution described step 3) is to have considered equivalent distribution network, reactive power compensation, and virtual bus is set between transformer and distribution network

Static load model, described virtual bus

With the actual load bus

Between be the equivalent impedance of transmission and distribution networks.

Described step 4) checking extrapolability, to set up load model with the less data of change in voltage, with the test of echoing of the larger data of change in voltage, obtain the meritorious and reactive power curve of match again, by the extrapolability of the size checking load model of error of fitting;

Described step 4) checking interpolation ability, to set up load model with the larger data of change in voltage, with the test of echoing of the less data of change in voltage, obtain the meritorious and reactive power curve of match again, by the interpolation ability of the size checking load model of error of fitting;

Described step 4) ability that checking comprehensive description different load forms is can describe the ability of the part throttle characteristics of different load composition, different change in voltage amplitudes with the module verification that the data that measure under a certain change in voltage are set up.

The present invention's beneficial effect compared with prior art is:

Load modeling method of the present invention takes full advantage of the transient state recorder data resource of oscillograph record, without harvester is installed again, saved the load modeling cost, can identification obtain in real time load model parameters accurately, and can carry out parameter identification to the integrated load model of static load model, dynamic load model, integrated load model, consideration distribution, the static load model of consideration distribution, for electrical network is set up reasonable and realistic Optimal Load model.Matched curve and the measured curve of the active power of the load model of this optimum and reactive power curve are more approaching, more meet on-the-spot actual, the accuracy of emulation and confidence level are all very high, and can never draw diametrically opposite conclusion under critical condition, with the accurate generating of existing electric power system and distribution model, can mutually coordinate, the analysis that makes electric power system is accurate and believable, thereby significantly improves accuracy and the confidence level that electric system simulation is analyzed.

Data preliminary treatment in load modeling method of the present invention, utilize symmetry and the asymmetric noisy data of recorder data to carry out load modeling, by quadrature Park conversion, it is the space vector conversion, make symmetry and asymmetric noisy data in recorder data all can be used for load modeling, effectively enlarged the available data sources of load modeling; The advantages such as that the levenberg-marquart algorithm that adopts has is less demanding to initial value, fast convergence rate,

Load modeling method of the present invention can be widely used in the load modeling of transformer station, large size industrial enterprise power distribution station, and its load model parameters that accurately truly reflects part throttle characteristics can provide basis for power system operation mode selection, simulation calculation and systems organization.

Embodiment

The present invention will be described below in conjunction with embodiment.

A kind of based on recorder data and belong to the power load modelling approach that method is distinguished in total body examination, following steps are arranged successively:

1) by recorder data is carried out to preliminary treatment, obtain the load modeling desired data;

Recorder data is 2010 05 month 17: 08: 25 on the 14th SSH339 that gather, the three-phase voltage instantaneous value 11DL voltage at load modeling node 11DL place and three-phase current instantaneous value 11DL electric current, i.e. three-phase voltage u _a, u _b, u _cAnd three-phase current i _a, i _b, i _cComprise fault recorder data, namely comprise the data of at least 2 cycles of presteady state process, the whole obvious disturbance transient process after fault, and the tend to be steady data of process or rear at least 5 cycles of steady-state process of disturbance, sample frequency is 24 points of every cycle, be 1.2KHz, sampling precision is that amplitude error is at most 1%, three-phase voltage u _a, u _b, u _cAnd three-phase current i _a, i _b, i _cSynchronism is good, and phase error approaches zero, if there is phase error in the instantaneous value collection, can bring very large impact to subsequent calculations;

11DL voltage and 11DL current data are carried out to preliminary treatment to be comprised:

The conversion perunit value: power taking press fiducial value namely the fixed numbers of selected commensurate be 57.74V, the current reference value is 1A, the ratio of the actual value of conversion voltage and current and the fixed numbers of selected commensurate;

Extract positive sequence component: utilize symmetrical component method to extract the positive sequence component of voltage, the magnitude of current, during stable state, there are following relation in the 0-1-2 component of electric weight and a-b-c component:

In formula:

Operator $\mathrm{\α}=e^{{j120}^0};$

The phasor value of a phase electric weight, x is voltage u or current i;

The phasor value of b phase electric weight, x is voltage u or current i;

The phasor value of c phase electric weight, x is voltage u or current i;

Subscript 0 expression zero-sequence component;

Subscript 1 expression positive sequence component;

Subscript 2 expression negative sequence components;

High-frequency noise in the smothing filtering data: adopt 5 secondary smothing filtering algorithms to step 2) positive sequence voltage and the forward-order current that obtain carry out filtering, the high-frequency noise that contains in the busbar voltage amplitude V that removal calculates, active-power P, reactive power Q, improve signal to noise ratio; The computing formula of 5 secondary balancing filtering algorithms is as follows:

$\left\{\begin{array}{c}y\left(n\right)=[-3x(n-2)+12x(n-1)+17x(n+1)-3x(n+2)]/35\\ y\left(0\right)=[69x\left(0\right)+4x\left(1\right)+4x\left(3\right)-x\left(4\right)]/60\\ y\left(1\right)=[2x\left(0\right)+27x\left(1\right)+12x\left(2\right)-8x\left(3\right)+2x\left(4\right)]/35\\ y(N-2)=[2x(N-5)-8x(N-4)+12x(N-3)+27x(N-2)]+2x(N-1)]/35\\ y(N-1)=[-x(N-5)+4x(N-4)-6x(N-3)+4x(N-2)+69x(N-1)]/70\end{array}\right.$

2≤n≤N-3

In formula:

Y (n): the signal after level and smooth;

X (n): the signal before level and smooth;

N: data length;

Quadrature park conversion: adopt quadrature park transformation for mula to carry out conversion to positive sequence voltage and forward-order current after level and smooth, obtain d shaft voltage u _d, q shaft voltage u _q, d shaft current i _dWith q shaft current i _q

Fundamental formular is as follows:

$\left[\begin{array}{c}{f}_d\\ f_q\\ f_0\end{array}\right]=\sqrt{\frac{2}{3}}\left[\begin{array}{ccc}{\cos \mathrm{\θ}}_a& {\cos \mathrm{\θ}}_b& {\cos \mathrm{\θ}}_c\\ -{\sin \mathrm{\θ}}_a& {-\sin \mathrm{\θ}}_b& {-\sin }_c\\ 1/\sqrt{2}& 1/\sqrt{2}& 1/\sqrt{2}\end{array}\right]\left[\begin{array}{c}{f}_a\\ f_b\\ f_c\end{array}\right];$

In formula:

θ _a＝ωt；

f _a: phase voltage u _aWith phase current i _aThe instantaneous value sequence;

f _b: phase voltage u _bWith phase current i _bThe instantaneous value sequence;

f _c: phase voltage u _cWith phase current i _cThe instantaneous value sequence;

Fundamental formular by quadrature Park conversion calculates:

$\left\{\begin{array}{c}{u}_d=\sqrt{\frac{2}{3}}(u_a{\cos \mathrm{\θ}}_a+u_b{\cos \mathrm{\θ}}_b+u_c{\cos \mathrm{\θ}}_c)\\ u_q=\sqrt{\frac{2}{3}}(-u_a{\sin \mathrm{\θ}}_a-u_b{\sin \mathrm{\θ}}_b-u_c{\sin \mathrm{\θ}}_c)\\ i_d=\sqrt{\frac{2}{3}}(i_a{\cos \mathrm{\θ}}_a+i_b{\cos \mathrm{\θ}}_b+i_c{\cos \mathrm{\θ}}_c)\\ i_q=\sqrt{\frac{2}{3}}(-i_a{\sin \mathrm{\θ}}_a-i_b{\sin \mathrm{\θ}}_b-i_c{\sin \mathrm{\θ}}_c)\end{array}\right.$

In formula:

u _d: d axle instantaneous voltage data sequence;

u _q: q axle instantaneous voltage data sequence;

i _d: d axle transient current data sequence;

i _q: q axle transient current data sequence;

The calculated load busbar voltage: fundamental formular is as follows:

$V=\sqrt{{u_d}^2+{u_q}^2};$

In formula:

u _d: d axle instantaneous voltage data sequence;

u _q: q axle instantaneous voltage data sequence;

V: load busbar voltage amplitude data sequence;

Calculate active power: fundamental formular is as follows:

P＝u _di _d+u _qi _q；

In formula:

u _d: d axle instantaneous voltage data sequence;

i _d: d axle transient current data sequence;

u _q: q axle instantaneous voltage data sequence;

i _q: q axle transient current data sequence;

P: actual measurement active power data sequence;

Calculate reactive power: fundamental formular is as follows:

Q＝u _qi _d-u _di _q；

In formula:

u _d: d axle instantaneous voltage data sequence;

i _d: d axle transient current data sequence;

u _q: q axle instantaneous voltage data sequence;

i _q: q axle transient current data sequence;

Q: actual measurement reactive power data sequence.

2) the load modeling desired data of selecting levenberg-marquart algorithm as Identification of parameter, preliminary treatment to be obtained carries out the damped least squares iterative equation and solves, and obtains optimum load model parameters, the measured value x of identification objects electric load model _i, y _iThere are two groups: (V, P), (V, Q), wherein V is actual measurement load busbar voltage, and P is actual measurement load active power, and Q is actual measurement reactive load power, y (x)=y (x; a ₀, a ₁..., a _M-1) be the functional relation of surveying independent variable and dependent variable of electric load model, wherein x _i, y _iFor measured value, a _jFor parameter to be identified; Utilize the N group data (x that has recorded _i, y _i), i=0, L, N-1, adopt levenberg-marquart algorithm to estimate unknown parameter a=(a ₀, a ₁, L, a _M-1) ^T, make target function J (a) (residual sum of squares (RSS)) minimum, namely Minimum;

If y=is (y ₀, y ₁, L, y _N-1) ^TAnd linear y=Aa+b between a, directly solve a=(A with linear least square ^TA) ^-1A ^T(y-b);

If y is (x; a ₀, L, a _M-1) be the non-linear form of a, can only first by one group of initial solution of a, set out and carry out iterative computation, namely with nonlinear least square method, solve a;

If A (a) is function y (x)=y (x; a ₀, L, a _M-1) the Jacobian matrix, namely

$A\left(a\right)=\left[\begin{array}{cccc}\frac{{\∂y}_0}{{\∂a}_0}& \frac{{\∂y}_0}{{\∂a}_1}& L& \frac{{\∂y}_0}{{\∂a}_{M-1}}\\ \frac{{\∂y}_1}{{\∂a}_0}& \frac{{\∂y}_1}{{\∂a}_1}& L& \frac{{\∂y}_1}{{\∂a}_{M-1}}\\ M& M& & M\\ \frac{{\∂y}_{N-1}}{{\∂a}_0}& \frac{{\∂y}_{N-1}}{{\∂a}_1}& L& \frac{{\∂y}_{N-1}}{{\∂a}_{M-1}}\end{array}\right];$

Described levenberg-marquart algorithm namely solves following system of linear equations:

$(A^{T}A+\mathrm{\λI})\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ M\\ y_{N-1}-y(x_{N-1};a)\end{array}\right);$

In formula:

$\mathrm{\λI}\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ y_{N-1}-y(x_{N-1};a)\end{array}\right)$ It is the system of linear equations that solves of steepest descent method;

$\left(A^{T}A\right)\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ M\\ y_{N-1}-y(x_{N-1};a)\end{array}\right)$ It is the system of linear equations that solves of Gauss-Newton method;

Wherein, damping factor λ=10 ^-4～10 ^-2

When λ>0 is too large, may occur that program does steepest at the little place, the lowest point of descending grade and descend, { the a} convergence rate descends to make sequence;

When λ too hour, convergence domain is too small, initial approximation a ₀Restricted;

Levenberg-marquart algorithm uses steepest descent method away from minimum point the time, selects larger λ value;

Levenberg-marquart algorithm is switched to Gauss-Newton method gradually when near minimum value, reduces λ value.

Described damped least squares iterative comprises:

Input solve for parameter initial value scope;

According to the motor steady-state equation, ask the initial value of induction motor state equation;

Utilize 4 rank Runge Kuttas of fixed step size to ask the electric motor state equation, i.e. the array of the differential equation;

Utilize initial parameter value to calculate active power and the reactive power of each data point, form the model output of vector form;

Solve the damped least squares iterative equation, and judge whether to meet the condition of convergence:

If do not met convergence, the Returning utilization initial parameter value continues to calculate active power and the reactive power of each data point;

As met convergence, the result of output global convergence, and whether judgement convergence result meets the constrained condition again, namely the load model parameters of motor can not be less than 0;

If do not met constraints, return to input solve for parameter initial value scope, automatically change initial parameter value identification again;

As meet constraints, output model parameter identification value;

3) adopting load model parameters to set up optimum load model is the integrated load model of considering distribution, has namely considered equivalent distribution network, reactive power compensation, and between transformer and distribution network, virtual bus has been set

integrated load model, virtual bus

with the actual load bus

between be the equivalent impedance of transmission and distribution networks, this embodiment is set up to integrated load model and the permanent power of the consideration distribution that obtains based on recorder data on-line identification, constant-impedance, constant current or three's combination compares as the power matched curve of the idealized load model of integrated load model, matched curve and the measured curve of the active power of this embodiment and reactive power curve are more approaching, more meet on-the-spot actual, the accuracy of emulation and confidence level are all very high, and can never draw diametrically opposite conclusion under critical condition, with the accurate generating of existing electric power system and distribution model, can mutually coordinate, the analysis that makes electric power system is accurate and believable, thereby significantly improve accuracy and confidence level that electric system simulation is analyzed,

4) load model parameters that obtains is carried out to the parameter check, checking has practical value and comprises the load model of the ability of extrapolability, interpolation ability and comprehensive description different load composition preferably, is optimum load model.

Checking extrapolability: set up load model with the data that change in voltage is less, then with the test of echoing of the larger data of change in voltage, obtain the meritorious and reactive power curve of match, verified the extrapolability of load model by the size of error of fitting;

Checking interpolation ability: set up load model with the data that change in voltage is larger, then with the test of echoing of the less data of change in voltage, obtain the meritorious and reactive power curve of match, verified the interpolation ability of load model by the size of error of fitting;

The ability that checking comprehensive description different load forms: the ability that can describe the part throttle characteristics of different load composition, different change in voltage amplitudes with the module verification that the data that measure under a certain change in voltage are set up.

Above content is in conjunction with concrete preferred implementation further description made for the present invention, can not assert that specific embodiment of the invention is confined to these explanations.For the general technical staff of the technical field of the invention; make without departing from the inventive concept of the premise some alternative or obvious modification that are equal to; and performance or purposes identical, all should be considered as belonging to the scope of patent protection that the present invention is determined by claims of submitting to.

Claims (10)

1. power load modelling approach of distinguishing method based on total body examination is characterized in that:

Following steps are arranged successively:

1) by recorder data is carried out to preliminary treatment, obtain the load modeling desired data;

Described recorder data is three-phase voltage instantaneous value and the three-phase current instantaneous value of load modeling Nodes, i.e. three-phase voltage u _a, u _b, u _cAnd three-phase current i _a, i _b, i _c, comprise fault recorder data, namely comprise the data of at least 2 cycles of presteady state process, the whole obvious disturbance transient process after fault, and the tend to be steady data of process or rear at least 5 cycles of steady-state process of disturbance;

2) the load modeling desired data of selecting levenberg-marquart algorithm as Identification of parameter, preliminary treatment to be obtained carries out the damped least squares iterative equation and solves, and obtains optimum load model parameters;

3) adopt load model parameters to set up optimum load model, load model that can identified parameters comprises the integrated load model of static load model, dynamic load model, integrated load model, consideration distribution, and the static load model of considering distribution;

4) load model parameters that obtains is carried out to the parameter check, checking has practical value and comprises the load model of the ability of extrapolability, interpolation ability and comprehensive description different load composition preferably, is optimum load model.

2. the power load modelling approach of distinguishing method based on total body examination as claimed in claim 1 is characterized in that:

The sample frequency of fault recorder data described step 1) is every cycle at least 20 points, namely is at least 1KHz, and sampling precision is that amplitude error is at most 1%, described three-phase voltage u _a, u _b, u _cAnd three-phase current i _a, i _b, i _cSynchronism is good, and phase error approaches zero, if there is phase error in the instantaneous value collection, can bring very large impact to subsequent calculations.

3. the power load modelling approach of distinguishing method based on total body examination as claimed in claim 2 is characterized in that:

Recorder data carried out to preliminary treatment comprise described step 1):

Conversion perunit value: the ratio of the actual value of conversion three-phase voltage instantaneous value and three-phase current instantaneous value and the fixed numbers of selected commensurate;

Extract positive sequence component: utilize symmetrical component method to extract the positive sequence component of voltage, the magnitude of current, during stable state, there are following relation in the 0-1-2 component of electric weight and a-b-c component:

In formula:

Operator $\mathrm{\α}=e^{{j120}^0};$

The phasor value of a phase electric weight, x is voltage u or current i;

The phasor value of b phase electric weight, x is voltage u or current i;

The phasor value of c phase electric weight, x is voltage u or current i;

Subscript 0 expression zero-sequence component;

Subscript 1 expression positive sequence component;

Subscript 2 expression negative sequence components;

High-frequency noise in the smothing filtering data: adopt the balance filtering algorithm to step 2) positive sequence voltage and the forward-order current that obtain carry out filtering, the high-frequency noise that contains in the busbar voltage amplitude V that removal calculates, active-power P, reactive power Q, improve signal to noise ratio; Described balance filtering algorithm comprises 5 secondary balancing filtering algorithms, and the computing formula of 5 secondary balancing filtering algorithms is as follows:

$\left\{\begin{array}{c}y\left(n\right)=[-3x(n-2)+12x(n-1)+17x(n+1)-3x(n+2)]/35\\ y\left(0\right)=[69x\left(0\right)+4x\left(1\right)+4x\left(3\right)-x\left(4\right)]/60\\ y\left(1\right)=[2x\left(0\right)+27x\left(1\right)+12x\left(2\right)-8x\left(3\right)+2x\left(4\right)]/35\\ y(N-2)=[2x(N-5)-8x(N-4)+12x(N-3)+27x(N-2)]+2x(N-1)]/35\\ y(N-1)=[-x(N-5)+4x(N-4)-6x(N-3)+4x(N-2)+69x(N-1)]/70\end{array}\right.$

2≤n≤N-3

In formula:

Y (n): the signal after level and smooth;

X (n): the signal before level and smooth;

N: data length;

Quadrature Park conversion: adopt quadrature park transformation for mula to carry out conversion to positive sequence voltage and forward-order current after level and smooth, obtain d shaft voltage u _d, q shaft voltage u _q, d shaft current i _dWith q shaft current i _qFundamental formular is as follows:

$\left[\begin{array}{c}{f}_d\\ f_q\\ f_0\end{array}\right]=\sqrt{\frac{2}{3}}\left[\begin{array}{ccc}{\cos \mathrm{\θ}}_a& {\cos \mathrm{\θ}}_b& {\cos \mathrm{\θ}}_c\\ -{\sin \mathrm{\θ}}_a& {-\sin \mathrm{\θ}}_b& {-\sin }_c\\ 1/\sqrt{2}& 1/\sqrt{2}& 1/\sqrt{2}\end{array}\right]\left[\begin{array}{c}{f}_a\\ f_b\\ f_c\end{array}\right];$

In formula:

θ _a＝ωt；

f _a: phase voltage u _aWith phase current i _aThe instantaneous value sequence;

f _b: phase voltage u _bWith phase current i _bThe instantaneous value sequence;

f _c: phase voltage u _cWith phase current i _cThe instantaneous value sequence;

Fundamental formular by quadrature Park conversion calculates:

$\left\{\begin{array}{c}{u}_d=\sqrt{\frac{2}{3}}(u_a{\cos \mathrm{\θ}}_a+u_b{\cos \mathrm{\θ}}_b+u_c{\cos \mathrm{\θ}}_c)\\ u_q=\sqrt{\frac{2}{3}}(-u_a{\sin \mathrm{\θ}}_a-u_b{\sin \mathrm{\θ}}_b-u_c{\sin \mathrm{\θ}}_c)\\ i_d=\sqrt{\frac{2}{3}}(i_a{\cos \mathrm{\θ}}_a+i_b{\cos \mathrm{\θ}}_b+i_c{\cos \mathrm{\θ}}_c)\\ i_q=\sqrt{\frac{2}{3}}(-i_a{\sin \mathrm{\θ}}_a-i_b{\sin \mathrm{\θ}}_b-i_c{\sin \mathrm{\θ}}_c)\end{array}\right.$

In formula:

u _d: d axle instantaneous voltage data sequence;

u _q: q axle instantaneous voltage data sequence;

i _d: d axle transient current data sequence;

i _q: q axle transient current data sequence;

The calculated load busbar voltage: fundamental formular is as follows:

$V=\sqrt{{u_d}^2+{u_q}^2};$

In formula:

u _d: d axle instantaneous voltage data sequence;

u _q: q axle instantaneous voltage data sequence;

V: load busbar voltage amplitude data sequence;

Calculate active power: fundamental formular is as follows:

P＝u _di _d+u _qi _q；

In formula:

u _d: d axle instantaneous voltage data sequence;

i _d: d axle transient current data sequence;

u _q: q axle instantaneous voltage data sequence;

i _q: q axle transient current data sequence;

P: actual measurement active power data sequence;

Calculate reactive power: fundamental formular is as follows:

Q＝u _qi _d-u _di _q；

In formula:

u _d: d axle instantaneous voltage data sequence;

i _d: d axle transient current data sequence;

u _q: q axle instantaneous voltage data sequence;

i _q: q axle transient current data sequence;

Q: actual measurement reactive power data sequence.

4. the power load modelling approach of distinguishing method based on total body examination as claimed in claim 3 is characterized in that:

Described step 2) obtain optimum load model parameters, comprising:

The measured value x of identification objects electric load model _i, y _iThere are two groups: (V, P), (V, Q), wherein V is actual measurement load busbar voltage, and P is actual measurement load active power, and Q is actual measurement reactive load power, y (x)=y (x; a ₀, a ₁..., a _M-1) be the functional relation of surveying independent variable and dependent variable of electric load model, wherein x _i, y _iFor measured value, a _jFor parameter to be identified; Utilize the N group data (x that has recorded _i, y _i), i=0, L, N-1, adopt levenberg-marquart algorithm to estimate unknown parameter a=(a ₀, a ₁, L, a _M-1) ^T, make target function J (a) (residual sum of squares (RSS)) minimum, namely

Minimum;

If y=is (y ₀, y ₁, L, y _N-1) ^TAnd linear y=Aa+b between a, directly solve a=(A with linear least square ^TA) ^-1A ^T(y-b);

If y is (x; a ₀, L, a _M-1) be the non-linear form of a, can only first by one group of initial solution of a, set out and carry out iterative computation, namely with nonlinear least square method, solve a;

If A (a) is function y (x)=y (x; a ₀, L, a _M-1) the Jacobian matrix, namely

$A\left(a\right)=\left[\begin{array}{cccc}\frac{{\∂y}_0}{{\∂a}_0}& \frac{{\∂y}_0}{{\∂a}_1}& L& \frac{{\∂y}_0}{{\∂a}_{M-1}}\\ \frac{{\∂y}_1}{{\∂a}_0}& \frac{{\∂y}_1}{{\∂a}_1}& L& \frac{{\∂y}_1}{{\∂a}_{M-1}}\\ M& M& & M\\ \frac{{\∂y}_{N-1}}{{\∂a}_0}& \frac{{\∂y}_{N-1}}{{\∂a}_1}& L& \frac{{\∂y}_{N-1}}{{\∂a}_{M-1}}\end{array}\right];$

Described levenberg-marquart algorithm namely solves following system of linear equations:

$(A^{T}A+\mathrm{\λI})\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ M\\ y_{N-1}-y(x_{N-1};a)\end{array}\right);$

In formula:

$\mathrm{\λI}\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ y_{N-1}-y(x_{N-1};a)\end{array}\right)$ It is the system of linear equations that solves of steepest descent method;

$\left(A^{T}A\right)\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ M\\ y_{N-1}-y(x_{N-1};a)\end{array}\right)$ It is the system of linear equations that solves of Gauss-Newton method;

Wherein, damping factor λ=10 ^-4～10 ^-2

When λ>0 is too large, may occur that program does steepest at the little place, the lowest point of descending grade and descend, { the a} convergence rate descends to make sequence;

When λ too hour, convergence domain is too small, initial approximation a ₀Restricted;

Levenberg-marquart algorithm uses steepest descent method away from minimum point the time, selects larger λ value;

Levenberg-marquart algorithm is switched to Gauss-Newton method gradually when near minimum value, reduces λ value.

5. the power load modelling approach of distinguishing method based on total body examination as claimed in claim 4 is characterized in that:

Described step 2) damped least squares iterative equation solves, and comprising:

Input solve for parameter initial value scope;

According to the motor steady-state equation, ask the initial value of induction motor state equation;

Utilize 4 rank Runge Kuttas of fixed step size to ask the electric motor state equation, i.e. the array of the differential equation;

Utilize initial parameter value to calculate active power and the reactive power of each data point, form the model output of vector form;

Solve following damped least squares iterative equation:

$(A^{T}A+\mathrm{\λI})\·\mathrm{\δa}=A^T\left(\begin{array}{c}{y}_0-y(x_0;a)\\ M\\ y_{N-1}-y(x_{N-1};\mathrm{\α})\end{array}\right),$

And judge whether to meet the condition of convergence:

If do not met convergence, the Returning utilization initial parameter value continues to calculate active power and the reactive power of each data point;

As met convergence, the result of output global convergence, and whether judgement convergence result meets the constrained condition again, namely the load model parameters of motor can not be less than 0;

If do not met constraints, return to input solve for parameter initial value scope, automatically change initial parameter value identification again;

As meet constraints, output model parameter identification value.

6. the power load modelling approach of distinguishing method based on total body examination as claimed in claim 5 is characterized in that:

Described step 3) static load model means namely the load model of functional relation of busbar voltage V of actual measurement active-power P and actual measurement reactive power Q and synchronization lower node voltage magnitude.

$\left\{\begin{array}{c}P=P_0(a_p{\left(\frac{V}{V_0}\right)}^2+b_p\left(\frac{V}{V_0}\right)+c_p)\\ Q=Q_0(a_q{\left(\frac{V}{V_0}\right)}^2+b_q\left(\frac{V}{V_0}\right)+c_q)\end{array};\right.$

a _p+b _p+c _p＝1；

a _q+b _q+c _q＝1；

In formula:

Coefficient a _p: constant-impedance (Z) component occupies the percentage of merit power;

Coefficient b _p: constant current (I) component occupies the percentage of merit power;

Coefficient c _p: permanent power (P) component occupies the percentage of merit power;

Coefficient a _q: constant-impedance (Z) component accounts for the percentage of reactive power;

Coefficient b _q: constant current (I) component accounts for the percentage of reactive power;

Coefficient c _q: permanent power (P) component accounts for the percentage of reactive power;

Described dynamic load model refers to the three order induction motor model, and its system state equation is as follows:

$\left\{\begin{array}{c}\frac{d{E}_q^{\′}}{\mathrm{dt}}=-w_{B}s{E}_d^{\′}-\frac{1}{T_0^{\′}}{E}_q^{\′}+\frac{X-X^{\′}}{T_0^{\′}}{i}_d\\ \frac{d{E}_d^{\′}}{\mathrm{dt}}=w_{B}s{E}_q^{\′}-\frac{1}{T_0^{\′}}{E}_d^{\′}-\frac{X-X^{\′}}{T_0^{\′}}{i}_q\\ \frac{{\mathrm{dw}}_m}{\mathrm{dt}}=\frac{1}{2H}(T_e-T_m)\end{array}\right.;$

In formula:

E ' _d: the transient potential of d axle;

E ' _q: the transient potential of q axle;

X: synchronous reactance;

X ': transient state reactance;

T ' ₀: the transient state open circuit time constant;

w _m: rotor speed;

w _s: synchronous speed;

i _q: the electric current in transient circuit is at the component of q axle;

i _d: the electric current in transient circuit is at the component of d axle;

T _s: be loaded into epitrochanterian electromagnetic torque;

T _m: the machine torque that rotor bears;

H: inertia time constant;

As follows through the power stage equation that derivation draws:

$\left\{\begin{array}{c}P=\frac{X_{t1}}{G}({u^2}_d-u_d{E}_d^{\′}+{u^2}_q-u_q{E}_q^{\′})+\frac{X_{t2}}{G}(u_q{E}_d^{\′}-u_d{E}_q^{\′})\\ Q=\frac{X_{t1}}{G}(u_d{E}_q^{\′}-u_q{E}_d^{\′})+\frac{X_{t2}}{G}({u^2}_d-u_d{E}_d^{\′}+{u^2}_q-u_q{E}_q^{\′})\end{array}\right.;$

In formula:

$X_{t1}=\frac{R_s(X-X^{\′})}{{R_s}^2+X^{\′2}};$

$X_{t2}=\frac{X^{\′}(X-X^{\′})}{{R_s}^2+X^{\′2}};$

G＝X-X′；

X＝X _s+X _m；

$X^{\′}=X_s+\frac{X_m{X}_r}{X_m+X_r}.$

7. the power load modelling approach of distinguishing method based on total body examination as claimed in claim 5 is characterized in that:

Described step 3) integrated load model is the comprehensive load model that comprises static load model and dynamic load model.

8. the power load modelling approach of distinguishing method based on total body examination as claimed in claim 5 is characterized in that:

The integrated load model of consideration distribution described step 3) is to have considered equivalent distribution network, reactive power compensation, and virtual bus is set between transformer and distribution network

Integrated load model, described virtual bus

With the actual load bus

Between be the equivalent impedance of transmission and distribution networks.

9. the power load modelling approach of distinguishing method based on total body examination as claimed in claim 5 is characterized in that:

The static load model of described consideration distribution described step 3) is to have considered equivalent distribution network, reactive power compensation, and virtual bus is set between transformer and distribution network

Static load model, described virtual bus

With the actual load bus

Between be the equivalent impedance of transmission and distribution networks.

10. distinguishing and it is characterized in that the power load modelling approach of method based on total body examination as any one as described in claim 1～9:

Described step 4) checking extrapolability, to set up load model with the less data of change in voltage, with the test of echoing of the larger data of change in voltage, obtain the meritorious and reactive power curve of match again, by the extrapolability of the size checking load model of error of fitting;

Described step 4) checking interpolation ability, to set up load model with the larger data of change in voltage, with the test of echoing of the less data of change in voltage, obtain the meritorious and reactive power curve of match again, by the interpolation ability of the size checking load model of error of fitting;

Described step 4) ability that checking comprehensive description different load forms is can describe the ability of the part throttle characteristics of different load composition, different change in voltage amplitudes with the module verification that the data that measure under a certain change in voltage are set up.