CN103401238B - A kind of power load modelling approach based on Measurement-based approach - Google Patents

Abstract

Based on a power load modelling approach for Measurement-based approach, there are following steps successively: 1) by carrying out preliminary treatment to recorder data, obtain load modeling desired data; 2) the load modeling desired data obtained preliminary treatment carries out damped least squares iterative equation and solves, and obtains optimum load model parameters; 3) load model parameters is adopted to set up optimum load model; 4) carry out parameters validation to the load model parameters obtained, checking has the load model that practical value comprises the ability of good extrapolability, interpolation capability and comprehensive description different load composition, is optimum load model.Save modeling cost, load model more meets on-the-spot actual, accuracy and the confidence level of emulation are all very high, and diametrically opposite conclusion can never be drawn under critical condition, can mutually coordinate with existing generating and distribution model, make the analysis of electric power system be accurate and believable, thus significantly improve accuracy and the confidence level of electric system simulation analysis.

Description

A kind of power load modelling approach based on Measurement-based approach

Technical field

The present invention relates to Modeling for Electric Loads, particularly relate to a kind of power load modelling approach based on Measurement-based approach.

Background technology

Power system load model is the power of reflection practical power systems load ports and electric current with the math equation of its port voltage and frequency variation characteristics and corresponding parameter.Load modeling is not only to various concrete power consumption equipment element Modling model, the more important thing is the relation that the power that the overall load on research load bus absorbs changes with load busbar voltage and frequency variation, and describe the form of determining the math equation of this relation and parameter wherein.The Digital Simulation of modern power systems design, planning, decision-making and operation, all needs the load model that accuracy is higher.Because the result of the change of load model on system transient modelling, voltage stabilization and Load flow calculation has impact in various degree, under critical condition, even there is the change of matter.And load modeling not only obtains the structure and parameter of model by one group of data fitting, also must carry out by other data test of echoing, to verify its comprehensive and to noise robustness.Under change in voltage a little or less greatly, remain correct by the load model of the measurement data identification under a certain change in voltage, the i.e. extrapolation of load model and interpolation, there is the load model of the ability of good extrapolability, interpolation capability and comprehensive description different load composition, just with practical value.

Existing load model is from basic conception mostly, adopt Utopian model and invariable power, constant-impedance, constant current or three combination as integrated load model.This coarse load model with generate electricity accurately and distribution model seriously inharmonious, significantly reduce 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 of Measurement-based approach.

Summary of the invention

Technical problem to be solved by this invention is the defect making up above-mentioned prior art, provides a kind of power load modelling approach based on Measurement-based approach.

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

This power load modelling approach based on Measurement-based approach, has following steps successively:

1) by carrying out preliminary treatment to recorder data, load modeling desired data is obtained;

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, comprises fault recorder data, namely comprise the whole obvious disturbance transient process after the data of presteady state process at least 2 cycles, fault, and disturbance tends to be steady the data of process or rear steady-state process at least 5 cycles;

2) select levenberg-marquart algorithm to carry out damped least squares iterative equation as Identification of parameter to the load modeling desired data that preliminary treatment obtains to solve, obtain optimum load model parameters;

3) load model parameters is adopted to set up optimum load model, the load model of identified parameters can comprise the integrated load model of static load model, dynamic load model, integrated load model, consideration distribution, and consider the static load model of distribution;

4) carry out parameters validation to the load model parameters obtained, checking has the load model that practical value comprises the ability of good extrapolability, interpolation capability and comprehensive description different load composition, is optimum load model.

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

Described step 1) the sample frequency of fault recorder data be every cycle at least 20 point, be namely at least 1KHz, 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 is close to zero, if instantaneous value collection exists phase error, can bring very large impact to subsequent calculations.

Described step 1) preliminary treatment carried out to recorder data comprise:

Conversion perunit value: the actual value of conversion three-phase voltage instantaneous value and three-phase current instantaneous value and the ratio of 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, the 0-1-2 component of electricity and a-b-c component exist following relation:

In formula:

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

the phasor value of a phase electricity, x is voltage u or current i;

the phasor value of b phase electricity, x is voltage u or current i;

the phasor value of c phase electricity, x is voltage u or current i;

Subscript 0 represents zero-sequence component;

Subscript 1 represents positive sequence component;

Subscript 2 represents negative sequence component;

High-frequency noise in smothing filtering data: adopt balancing filter algorithm to step 2) positive sequence voltage that obtains and forward-order current carry out filtering, the high-frequency noise contained in the busbar voltage amplitude V that removal calculates, active-power P, reactive power Q, improves signal to noise ratio; Described balancing filter 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;

Orthogonal Park conversion: adopt orthogonal park transformation for mula to convert the positive sequence voltage after level and smooth and forward-order current, 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 _ainstantaneous value sequence;

F _b: phase voltage u _bwith phase current i _binstantaneous value sequence;

F _c: phase voltage u _cwith phase current i _cinstantaneous value sequence;

The fundamental formular converted by orthogonal Park 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;

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 _ihave 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 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 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, then direct linear least square solves a=(A ^ta) ^-1a ^t(y-b);

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

If A (a) is function y (x)=y (x; a ₀, L, a _m-1) 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];$

Namely described levenberg-marquart algorithm 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)$ Be steepest descent method solve system of linear equations;

$\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)$ Be Gauss-Newton method solve system of linear equations;

Wherein, damping factor λ=10 ^-4~ 10 ^-2;

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

When λ is too little, then convergence domain is too small, initial approximation a ₀restricted;

Using steepest descent method away from levenberg-marquart algorithm during minimum point, select larger λ value;

When being switched to Gauss-Newton method gradually close to levenberg-marquart algorithm during minimum value, reduce λ value.

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

Input solve for parameter initial value scope;

The initial value of induction motor state equation is asked according to motor steady-state equation;

4 rank Runge Kuttas of fixed step size are utilized to ask 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, the model forming vector form exports;

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, then Returning utilization initial parameter value continues the active power and the reactive power that calculate each data point;

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

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

As met constraints, then output model parameter identification value.

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

$\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 accounts for the percentage of active power;

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

Coefficient c _p: invariable power (P) component accounts for the percentage of active 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: invariable power (P) component accounts for the percentage of reactive power.

Described dynamic load model refers to 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 ' ₀: 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;

The power stage equation drawn through deriving is as follows:

$\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 be the comprehensive load model comprising static load model and dynamic load model.

Described step 3) the integrated load model of consideration distribution be consider equivalent distribution network, reactive power compensation, and virtual bus is set between transformer and distribution network integrated load model, described virtual bus with actual load bus between be the equivalent impedance of transmission and distribution networks.

Described step 3) the static load model of consideration distribution be consider equivalent distribution network, reactive power compensation, and virtual bus is set between transformer and distribution network static load model, described virtual bus with actual load bus between be the equivalent impedance of transmission and distribution networks.

Described step 4) checking extrapolability, set up load model by the data that change in voltage is less, carry out by the data that voltage changes greatly test of echoing again, obtain the meritorious of matching and reactive power curve, verified the extrapolability of load model by the size of error of fitting;

Described step 4) checking interpolation capability, set up load model by the data that voltage changes greatly, carry out by the data that change in voltage is less test of echoing again, obtain the meritorious of matching and reactive power curve, verified the interpolation capability of load model by the size of error of fitting;

Described step 4) the ability of checking comprehensive description different load composition, be the ability of part throttle characteristics with measuring module verification that the data that obtain set up under a certain change in voltage and can describe different load composition, different change in voltage amplitude.

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

Load modeling method of the present invention makes full use of the transient state recorder data resource of oscillograph record, without the need to installing harvester again, save load modeling cost, identification can obtain in real time load model parameters accurately, and parameter identification can be carried out, for electrical network sets up reasonable and realistic Optimal Load model to the static load model of the integrated load model of static load model, dynamic load model, integrated load model, consideration distribution, consideration distribution.The active power of the load model of this optimum and the matched curve of reactive power curve and measured curve closer to, more meet on-the-spot actual, accuracy and the confidence level of emulation are all very high, and diametrically opposite conclusion can never be drawn under critical condition, with existing electric power system generate electricity accurately and distribution model can mutually coordinate, make the analysis of electric power system be accurate and believable, thus significantly improve accuracy and the confidence level of electric system simulation analysis.

Data prediction in load modeling method of the present invention, the symmetry of recorder data and asymmetric noisy data is utilized to carry out load modeling, by orthogonal Park conversion and space vector conversion, make the symmetry in recorder data and asymmetric noisy data all can be used for load modeling, effectively expand the available data sources of load modeling; The levenberg-marquart algorithm adopted has the advantages such as, fast convergence rate less demanding to initial value,

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 the load model parameters of its accurate true reflection part throttle characteristics can be selected for power system operation mode, simulation calculation and systems organization provide basis.

Embodiment

Below in conjunction with embodiment, the present invention will be described.

Belong to the power load modelling approach of Measurement-based approach based on recorder data, have following steps successively:

1) by carrying out preliminary treatment to recorder data, load modeling desired data is obtained;

08 point of SSH339 gathered for 25 seconds when recorder data is 2010 05 month 14 days 17, 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 the whole obvious disturbance transient process after the data of presteady state process at least 2 cycles, fault is comprised, and disturbance tends to be steady the data of process or rear steady-state process at least 5 cycles, sample frequency is every cycle 24 points, 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 is close to zero, if instantaneous value collection exists phase error, can bring very large impact to subsequent calculations;

Carry out preliminary treatment to 11DL voltage and 11DL current data to comprise:

Conversion perunit value: the fixed numbers of getting the commensurate that namely voltage reference value is selected is 57.74V, and current reference value is 1A, the actual value of scaled voltage and electric current and the ratio of 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, the 0-1-2 component of electricity and a-b-c component exist following relation:

In formula:

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

the phasor value of a phase electricity, x is voltage u or current i;

the phasor value of b phase electricity, x is voltage u or current i;

the phasor value of c phase electricity, x is voltage u or current i;

Subscript 0 represents zero-sequence component;

Subscript 1 represents positive sequence component;

Subscript 2 represents negative sequence component;

High-frequency noise in smothing filtering data: adopt 5 secondary smothing filtering algorithms to step 2) positive sequence voltage that obtains and forward-order current carry out filtering, the high-frequency noise contained in the busbar voltage amplitude V that removal calculates, active-power P, reactive power Q, improves 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;

Orthogonal park conversion: adopt orthogonal park transformation for mula to convert the positive sequence voltage after level and smooth and forward-order current, 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 _ainstantaneous value sequence;

F _b: phase voltage u _bwith phase current i _binstantaneous value sequence;

F _c: phase voltage u _cwith phase current i _cinstantaneous value sequence;

The fundamental formular converted by orthogonal Park 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;

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) select levenberg-marquart algorithm to carry out damped least squares iterative equation as Identification of parameter to the load modeling desired data that preliminary treatment obtains to solve, obtain optimum load model parameters, the measured value x of identification objects electric load model _i, y _ihave 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 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 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, then direct linear least square solves a=(A ^ta) ^-1a ^t(y-b);

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

If A (a) is function y (x)=y (x; a ₀, L, a _m-1) 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];$

Namely described levenberg-marquart algorithm 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)$ Be steepest descent method solve system of linear equations;

$\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)$ Be Gauss-Newton method solve system of linear equations;

Wherein, damping factor λ=10 ^-4~ 10 ^-2;

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

When λ is too little, then convergence domain is too small, initial approximation a ₀restricted;

Using steepest descent method away from levenberg-marquart algorithm during minimum point, select larger λ value;

When being switched to Gauss-Newton method gradually close to levenberg-marquart algorithm during minimum value, reduce λ value.

Described damped least squares iterative, comprising:

Input solve for parameter initial value scope;

The initial value of induction motor state equation is asked according to motor steady-state equation;

4 rank Runge Kuttas of fixed step size are utilized to ask 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, the model forming vector form exports;

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

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

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

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

As met constraints, then output model parameter identification value;

3) adopting load model parameters to set up optimum load model is the integrated load model considering distribution, namely considers equivalent distribution network, reactive power compensation, and between transformer and distribution network, arranges virtual bus integrated load model, virtual bus with actual load bus between be the equivalent impedance of transmission and distribution networks, this embodiment is set up integrated load model and the invariable power of the consideration distribution obtained based on recorder data on-line identification, constant-impedance, the combination of constant current or three compares as the power matched curve of the idealized load model of integrated load model, the active power of this embodiment and the matched curve of reactive power curve and measured curve closer to, more meet on-the-spot actual, accuracy and the confidence level of emulation are all very high, and diametrically opposite conclusion can never be drawn under critical condition, with existing electric power system generate electricity accurately and distribution model can mutually coordinate, the analysis of electric power system is made to be accurate and believable, thus significantly improve accuracy and the confidence level of electric system simulation analysis,

4) carry out parameters validation to the load model parameters obtained, checking has the load model that practical value comprises the ability of good extrapolability, interpolation capability and comprehensive description different load composition, is optimum load model.

Checking extrapolability: the data less by change in voltage set up load model, then carry out by the data that voltage changes greatly test of echoing, obtains the meritorious of matching and reactive power curve, is verified the extrapolability of load model by the size of error of fitting;

Checking interpolation capability: set up load model by the data that voltage changes greatly, then carry out by the data that change in voltage is less test of echoing, obtains the meritorious of matching and reactive power curve, is verified the interpolation capability of load model by the size of error of fitting;

The ability of checking comprehensive description different load composition: by the ability of part throttle characteristics measuring module verification that the data that obtain set up under a certain change in voltage and can describe different load composition, different change in voltage amplitude.

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 general technical staff of the technical field of the invention; make some equivalent alternative or obvious modification without departing from the inventive concept of the premise; and performance or purposes identical, all should be considered as belonging to the scope of patent protection that the present invention is determined by submitted to claims.

Claims (7)

1., based on a power load modelling approach for Measurement-based approach, it is characterized in that:

There are following steps successively:

1) by carrying out preliminary treatment to recorder data, load modeling desired data is obtained;

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, comprises fault recorder data, namely comprise the whole obvious disturbance transient process after the data of presteady state process at least 2 cycles, fault, and disturbance tends to be steady the data of process and rear steady-state process at least 5 cycles;

The sample frequency of described fault recorder data is every cycle at least 20 point, and be namely at least 1KHz, 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;

Describedly preliminary treatment carried out to recorder data comprise:

Conversion perunit value: the actual value of conversion three-phase voltage instantaneous value and three-phase current instantaneous value and the ratio of 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, the 0-1-2 component of electricity and a-b-c component exist following relation:

$\left[\begin{array}{c}{\stackrel{\·}{x}}_a^0\\ {\stackrel{\·}{x}}_a^1\\ {\stackrel{\·}{x}}_a^2\end{array}\right]=\frac{1}{3}\left[\begin{array}{ccc}1& 1& 1\\ 1& \mathrm{\α}& {\mathrm{\α}}^2\\ 1& {\mathrm{\α}}^2& \mathrm{\α}\end{array}\right]\left[\begin{array}{c}{\stackrel{\·}{x}}_a\\ {\stackrel{\·}{x}}_b\\ {\stackrel{\·}{x}}_c\end{array}\right];$

In formula:

Operator α=e ^{j120 °};

the phasor value of a phase electricity, x is voltage u or current i;

the phasor value of b phase electricity, x is voltage u or current i;

the phasor value of c phase electricity, x is voltage u or current i;

Subscript 0 represents zero-sequence component;

Subscript 1 represents positive sequence component;

Subscript 2 represents negative sequence component;

High-frequency noise in smothing filtering data: adopt balancing filter algorithm to step 2) positive sequence voltage that obtains and forward-order current carry out filtering, the high-frequency noise contained in the busbar voltage amplitude V that removal calculates, active-power P, reactive power Q, improves signal to noise ratio; Described balancing filter 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;

Orthogonal Park conversion: adopt orthogonal park transformation for mula to convert the positive sequence voltage after level and smooth and forward-order current, 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}{\mathrm{cos\θ}}_a& {\mathrm{cos\θ}}_b& {\mathrm{cos\θ}}_c\\ -{\mathrm{sin\θ}}_a& -{\mathrm{sin\θ}}_b& -{\mathrm{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 _ainstantaneous value sequence;

F _b: phase voltage u _bwith phase current i _binstantaneous value sequence;

F _c: phase voltage u _cwith phase current i _cinstantaneous value sequence;

The fundamental formular converted by orthogonal Park calculates:

$\left\{\begin{array}{c}{u}_d=\sqrt{\frac{2}{3}}(u_a{\mathrm{cos\θ}}_a+u_b{\mathrm{cos\θ}}_b+u_c{\mathrm{cos\θ}}_c)\\ u_q=\sqrt{\frac{2}{3}}(-u_a{\mathrm{sin\θ}}_a-u_b{\mathrm{sin\θ}}_b-u_c{\mathrm{sin\θ}}_c)\\ i_d=\sqrt{\frac{2}{3}}(i_a{\mathrm{cos\θ}}_a+i_b{\mathrm{cos\θ}}_b+i_c{\mathrm{cos\θ}}_c)\\ i_q=\sqrt{\frac{2}{3}}(-i_a{\mathrm{sin\θ}}_a-i_b{\mathrm{sin\θ}}_b-i_c{\mathrm{sin\θ}}_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;

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) select levenberg-marquart algorithm to carry out damped least squares iterative equation as Identification of parameter to the load modeling desired data that preliminary treatment obtains to solve, obtain optimum load model parameters;

3) load model parameters is adopted to set up optimum load model, the load model of identified parameters can comprise the integrated load model of static load model, dynamic load model, integrated load model, consideration distribution, and consider the static load model of distribution;

4) carry out parameters validation to the load model parameters obtained, checking has the load model that practical value comprises the ability of good extrapolability, interpolation capability and comprehensive description different load composition, is optimum load model.

2., as claimed in claim 1 based on the power load modelling approach of Measurement-based approach, it 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 _ihave 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 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 recorded _i, y _i), i=0 ..., N-1, adopts levenberg-marquart algorithm to estimate unknown parameter a=(a ₀, a ₁..., a _m-1) ^t, make target function J (a) minimum, namely minimum value.

3., as claimed in claim 2 based on the power load modelling approach of Measurement-based approach, it is characterized in that:

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

Input solve for parameter initial value scope;

The initial value of induction motor state equation is asked according to motor steady-state equation;

4 rank Runge Kuttas of fixed step size are utilized to ask 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, the model forming vector form exports.

4., as claimed in claim 3 based on the power load modelling approach of Measurement-based approach, it is characterized in that:

Described step 3) integrated load model be the comprehensive load model comprising static load model and dynamic load model.

5., as claimed in claim 1 based on the power load modelling approach of Measurement-based approach, it is characterized in that:

Described step 3) the integrated load model of described consideration distribution be consider equivalent distribution network, reactive power compensation, and virtual bus is set between transformer and distribution network integrated load model, described virtual bus with actual load bus between be the equivalent impedance of transmission and distribution networks.

6., as claimed in claim 1 based on the power load modelling approach of Measurement-based approach, it is characterized in that:

Described step 3) the static load model of described consideration distribution be consider equivalent distribution network, reactive power compensation, and virtual bus is set between transformer and distribution network static load model, described virtual bus with actual load bus between be the equivalent impedance of transmission and distribution networks.

7. the power load modelling approach based on Measurement-based approach of any one as described in claim 1 ~ 6, is characterized in that:

Described step 4) checking extrapolability, set up load model by the data that change in voltage is less, carry out by the data that voltage changes greatly test of echoing again, obtain the meritorious of matching and reactive power curve, verified the extrapolability of load model by the size of error of fitting;

Described step 4) checking interpolation capability, set up load model by the data that voltage changes greatly, carry out by the data that change in voltage is less test of echoing again, obtain the meritorious of matching and reactive power curve, verified the interpolation capability of load model by the size of error of fitting;

Described step 4) the ability of checking comprehensive description different load composition, be the ability of part throttle characteristics with measuring module verification that the data that obtain set up under a certain change in voltage and can describe different load composition, different change in voltage amplitude.