CN102074955A - Method based on knowledge discovery technology for stability assessment and control of electric system - Google Patents

Abstract

The invention provides a method for transient stability assessment and real-time control of an electric system. According to the method, technologies concerning knowledge discovery are introduced to realize the stability assessment and the control of the electric system; and on the basis of establishing a transient stability assessment model, the corresponding relation between the system running condition and the overall stability of the electric system is found out, the inherent running discipline of the electric system are grasped, and the fast assessment of the stable running state of the electric system is realized, therefore, the actual running management of the electric system can be guided. The method has the advantages that timely and effective stability control scheme can be established for the running state which loses stability after great disturbance, and the safe and the stable running of the system can be ensured at lowest price.

Description

Power system stability assessment and control method based on knowledge discovering technologies

Technical field

The invention belongs to power system operation and control technology field, be specifically related to a kind of power system stability assessment and control method based on knowledge discovering technologies.

Background technology

In electric power system, stable problem is put on the first place all the time, and it is determining that can electric power system normal power supply.In the power system stability accident, majority is the transient stability malicious event, and the transient stability problem of research electric power system has important practical sense.Main research method has the mixed method of time domain simulation method (or claiming step by step integration), transient energy function method (or claiming direct method) and these two kinds of methods at present.

The electric power system transient stability that utilizes intelligent method to carry out non-model is differentiated the advantages such as heuristic rule with fast, the easy generation decision-making of online computational speed, can constitute good complementation with traditional transient stability analysis method.Meanwhile, computer technology has obtained extensive use in electric power system, and electric power system can be collected a large amount of service datas by various monitoring devices, and this directly utilizes metric data analysis for the transient stability problem sufficient technical support is provided.Current outstanding problem is in the face of rich data like this, how systematically to develop strong data analysis tool, accurately and timely extracts Useful Information from initial data.

Knowledge discovering technologies is applied to electric power system, can utilize actual operating data more fully, disclose power system operation feature and rule that mass data contains behind, find the method for more reasonably dealing with problems, for the formulation of making a strategic decision provides stronger scientific basis.

Summary of the invention

The invention provides a kind of electric power system transient stability online evaluation and control method, the Knowledge Discovery correlation technique is introduced power system stability assessment and control problem, on the basis that makes up the transient stability assessment models, find out the corresponding relation between system conditions and the stability of power system, the moving law of grasp system inherence, realization is to the quick identification of power system stability running status, in order to instruct the operational management of electric power system reality; To suffering the running status of unstability after the big disturbance, provide and stablize controlling schemes timely and effectively, guarantee the safe and stable operation of system with less cost.

The present invention starts with from static and dynamic two aspects, and the input feature vector variable of the stable assessment of definition is set up stable assessment Mathematical Modeling.Generate sample by field measurement electric parameters or real-time simulation data.Use the data preconditioning technique to improve the quality of data, reduce actual assessment and calculate the used time, improve precision and the performance analyzed.Utilize SVMs to come the structural classification device, carry out the power system stability assessment.Instability status is adopted the feasible emergency control measure of mixed type particle swarm optimization algorithm search, guarantee the stable operation of system.

According to a kind of power system stability assessment of the present invention and control method, may further comprise the steps:

(1) determines the input feature vector variable of electric power system, make up the transient stability assessment models of electric power system;

(2) obtain electrical network on-line operation data and the simulation calculation of the transient stability in short-term result under given fault collection thereof in the electric power system, according to the input feature vector variable that step (1) is determined, the composing training sample space, each sample standard deviation is by input feature vector variable set x ₀Temporarily steady result of calculation y with correspondence ₀Two parts constitute;

(3) on the basis of lost data information not, initial mass data sample is carried out the data preliminary treatment, obtain new data sample set (x, y);

(4) according to training sample set (x, y), structure transient stability optimal classification discriminant function:

$f\left(x\right)=\mathrm{sgn}\{(w^{*}\·x)+b^{*}\}=\mathrm{sgn}\{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i^{*}{y}_{i}K(x_i\·x)+b^{*}\}---\left(1\right)$

Wherein sgn () is a sign function, w ^*Be the slope of once linear function, b ^*Be the intercept of once linear function, α _iBe the Lagrange multiplier corresponding, K (x with each sample _i, be Gauss's radial kernel function x), have

$K(x_i,x)=\exp (-\frac{{|x-x_i|}^2}{{\mathrm{\σ}}^2})$

Wherein σ is the width parameter of function, has controlled the radial effect scope of function.

(5) for the service data x of operation states of electric power system the unknown, after (3) passed through the data preliminary treatment set by step with it, substitution formula (1) promptly obtained its corresponding transient stability running status f (x); The value of f (x) has been represented the distance of each sample to the optimal classification hyperplane, and in the transient stability evaluation problem, the value of f (x) has been represented the distance of current running status x to the transient stability classifying face, i.e. transient stability nargin.F (x) value represent system stability, and the big more system of numerical value is stable more for just; Otherwise value is a negative indication system transient state instability, and this number big more expression of absolute value system is unstable more;

(6) step (5) is judged to be the power system operation mode of unstability, utilizes the mixed type particle swarm optimization algorithm, search for and obtain to make the emergency control measure of operation of power system recovery transient stability and cutter or cutting load amount minimum.

Wherein, determine to reflect the input feature vector x of transient characterisitics when electric power system is broken down ₀And corresponding temporarily steady result of calculation y ₀, make up the electric power system transient stability assessment models; The input feature vector variable x that step (1) adopts ₀Being divided into static nature variable and behavioral characteristics variable, also is x ₀=(static nature variable, behavioral characteristics variable) specifically comprises:

(1) the static nature variable comprises

1) generator of whole electric power system exerting oneself and loading;

2) the minimum and maximum generated output of electric power system;

3) the minimum and maximum load of electric power system;

4) the minimum and maximum value of busbar voltage;

5) the maximum meritorious through-put power of circuit, and corresponding idle through-put power;

6) key transmission cross-section and zone contact section) through-put power;

7) generated output of crucial generator;

8) each regional generator exerting oneself and loading;

(2) the behavioral characteristics variable comprises:

Below various in, δ _iBe the generator amature angle; ω _iBe generator speed,

a _iBe the generator amature angular acceleration, M _iInertia constant for generator; P _MiBe the generator mechanical output; P _EiBe generator electromagnetic power; V _KiBe generator energy; V _PiBe generator potential energy; N _GFor the generator number and, t ₀Be fault generation initial time, t _ClBe the failure removal moment;

1) fault generation initial time has the relative rotor angle of the generator of maximum work angular acceleration:

$x_1=({\mathrm{\&delta;}}_{a\max }-{\mathrm{\&delta;}}_{\mathrm{COI}}){|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

In the following formula

Be equivalent rotor angle with respect to the system inertia center, δ _AmaxIt is the rotor angle that fault generation initial time has the generator of maximum work angular acceleration;

2) failure removal has the relative rotor angle of the generator of maximum kinetic energy constantly:

$x_2=({\mathrm{\&delta;}}_{p\max }-{\mathrm{\&delta;}}_{\mathrm{COI}}){|}_{t_{\mathrm{cl}}}$

δ in the following formula _PmaxIt is the rotor angle that the fault excision has the generator of maximum kinetic energy constantly;

3) the merit angle of the machine and the machine of bringing up the rear was poor before failure removal was led constantly; Fault generation initial time, the generator with maximum rotor angular acceleration is the preceding machine of neck; Generator with minimum rotor angle acceleration is the machine of bringing up the rear;

$x_3=\max \left({\mathrm{\&delta;}}_{a\max }\right)-\min \left({\mathrm{\&delta;}}_{a\min }\right){|}_{t_{\mathrm{cl}}}$

δ in the following formula _Amax, δ _AminIt is respectively the rotor angle of the preceding machine of neck and the machine of bringing up the rear;

4) the maximum constantly generator outer corner difference of failure removal:

$x_4=\max \left({\mathrm{\&delta;}}_i\right)-\min \left({\mathrm{\&delta;}}_i\right){|}_{t_{\mathrm{cl}}}$

5) failure removal differs maximum merit angle with center of inertia COI constantly:

$x_5=\max ({\mathrm{\&delta;}}_i-{\mathrm{\&delta;}}_{\mathrm{COI}}){|}_{t_{\mathrm{cl}}}$

The inertia weighted average angle of the generator before 6) failure removal is led constantly in a group of planes and the remaining group of planes poor:

$x_6=\underset{i\&Element;s}{\mathrm{\&Sigma;}}{M}_i{\mathrm{\&delta;}}_i/\underset{i\&Element;s}{\mathrm{\&Sigma;}}{M}_i-\underset{j\&Element;A}{\mathrm{\&Sigma;}}{M}_j{\mathrm{\&delta;}}_j/\underset{j\&Element;A}{\mathrm{\&Sigma;}}{M}_j{|}_{t_{\mathrm{cl}}}$

Wherein, group of planes set before S represents to lead, A represents remaining group of planes set.To specifying each generator corner constantly to carry out descending sort, shape as:

$({\mathrm{\&delta;}}_1,{\mathrm{\&delta;}}_2,{\mathrm{\&delta;}}_3,...,{\mathrm{\&delta;}}_{N_G-2},{\mathrm{\&delta;}}_{N_G-1},{\mathrm{\&delta;}}_{N_G}),$ And ${\mathrm{\&delta;}}_1>{\mathrm{\&delta;}}_2>{\mathrm{\&delta;}}_3>...>{\mathrm{\&delta;}}_{N_G-2},{\mathrm{\&delta;}}_{N_G-1}>{\mathrm{\&delta;}}_{N_G}$

The difference of per two adjacent corners constitutes an angular clearances, shape as:

$({\mathrm{\&delta;}}_{\mathrm{1,2}},{\mathrm{\&delta;}}_{\mathrm{2,3}},...,{\mathrm{\&delta;}}_{N_G-2,N_G-1},{\mathrm{\&delta;}}_{N_G-1,N_G}),$ And

${\mathrm{\&delta;}}_{\mathrm{1,2}}={\mathrm{\&delta;}}_1-{\mathrm{\&delta;}}_2,{\mathrm{\&delta;}}_{\mathrm{2,3}}={\mathrm{\&delta;}}_2-{\mathrm{\&delta;}}_3,...,{\mathrm{\&delta;}}_{N_G-2,N_G-1}={\mathrm{\&delta;}}_{N_G-2}-{\mathrm{\&delta;}}_{N_G-1},{\mathrm{\&delta;}}_{N_G-1,N_G}={\mathrm{\&delta;}}_{N_G-1}-{\mathrm{\&delta;}}_{N_G}$

Determine maximum angular clearances, the generator on this gap constitutes the preceding group of planes S of neck, and the generator under this gap constitutes remaining group of planes A;

7) failure removal constantly with the absolute value sum of the relative rotor angle difference of fault generation initial time generator:

$x_7=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}|{\mathrm{\&theta;}}_{\mathrm{ci}}-{\mathrm{\&theta;}}_{0i}|$

In the following formula, θ _i=δ _i-δ _COI, θ _cBe the relative rotor angle of failure removal moment generator, θ ₀Be the relative rotor angle of fault generation initial time generator.

8) failure removal constantly with the maximum of the absolute value of the relative rotor angle difference of fault generation initial time generator:

$x_8=\underset{i}{\max }\left(\right|{\mathrm{\&theta;}}_{\mathrm{ci}}-{\mathrm{\&theta;}}_{0i}\left|\right)$

9) the failure removal moment and fault generation initial time center of inertia coordinate COI change in rotational speed:

$x_9={\mathrm{\&omega;}}_{\mathrm{COI}}{|}_{t_{\mathrm{cl}}}-{\mathrm{\&omega;}}_{\mathrm{COI}}{|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

In the following formula ${\mathrm{\&omega;}}_{\mathrm{COI}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i{\mathrm{\&omega;}}_i/\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i$

10) the absolute value sum of failure removal moment generator amature angular speed:

$x_{10}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{|{\mathrm{\&omega;}}_i-{\mathrm{\&omega;}}_{\mathrm{COI}}|}_{t_{\mathrm{cl}}}$

11) maximum of failure removal moment generator amature angular speed absolute value:

$x_{11}=\max \left(\right|{\mathrm{\&omega;}}_i-{\mathrm{\&omega;}}_{\mathrm{COI}}\left|\right){|}_{t_{\mathrm{cl}}}$

12) all generator amatures are in the most large and small value of fault generation initial time relative acceleration:

$x_{12-1}=(a_{\max }-a_{\mathrm{COI}}){|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

$x_{12-2}=(a_{\min }-a_{\mathrm{COI}}){|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

Wherein be a _COICenter of inertia acceleration

a _Max, a _MinIt is respectively the most large and small value of rotor angle acceleration;

13) all generator amatures are in the average of fault generation initial time acceleration:

$x_{13}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{a}_i{|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

14) all generator amatures are in the root-mean-square error of fault generation initial time acceleration:

$x_{14}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{(a_i-x_{13})}^2}{|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

15) all generator amatures are at the mean value of fault generation initial time relative acceleration:

$x_{15}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}(a_i-a_{\mathrm{COI}}){|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

16) all generator amatures are in the variance of fault generation initial time relative acceleration:

$x_{16}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{((a_i-a_{\mathrm{COI}})-x_{15})}^2}{|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

17) poor at instant of failure and fault generation initial time acceleration of machine and the machine of bringing up the rear before the neck:

$x_{17-1}=\max \left(a_i\right)-\min \left(a_i\right){|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

$x_{17-2}=\max \left(a_i\right)-\min \left(a_i\right){|}_{t_{\mathrm{cl}}}$

18) maximum of all rotor kinetic energy of the failure removal moment:

$x_{18}=\max \left(V_{\mathrm{ki}}\right){|}_{t_{\mathrm{cl}}}i=\mathrm{1,2},...,N_G$

19) failure removal has the rotor kinetic energy of the generator p of maximum rotor angle constantly:

$x_{19}=\frac{1}{2}{M}_p\&times;({{\mathrm{\&omega;}}_p}^2-1){|}_{t_{\mathrm{cl}}}$

20) total " the energy adjustment " of electric power system:

$x_{20}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}(P_{\mathrm{ai}}\&times;{\mathrm{\&delta;}}_{\mathrm{di}})$

Wherein, P _AiBe the accelerating power of i platform generator at fault generation initial time;

Relative rotor angle when being the excision of i platform generator failure;

21) mean value of all rotor kinetic energy of the failure removal moment:

$x_{21}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\frac{1}{2}{M}_i\&times;({{\mathrm{\&omega;}}_i}^2-1){|}_{t_{\mathrm{cl}}}$

22) mean value of all generator mechanical input powers before the fault:

$x_{22}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{p}_{\mathrm{mi}}$

23) variance of all generator mechanical input powers before the fault:

$x_{23}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{(P_{\mathrm{mi}}-x_{22})}^2}$

24) generator is in the average of fault generation initial time accelerating power:

$x_{24}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}){|}_{t_{\mathrm{cl}}}$

25) generator is in the variance of the relative accelerating power of fault generation initial time:

$x_{25}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{(\mathrm{\&Delta;}{P}_i-x_{24})}^2}{|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

26) generator is in the average of the relative accelerating power of fault generation initial time:

$x_{26}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&Delta;P}}_i}{P_{\mathrm{mi}}}{|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

27) generator is in the variance of the relative accelerating power of fault generation initial time:

$x_{27}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{(\frac{{\mathrm{\&Delta;P}}_i}{P_{\mathrm{mi}}}-x_{26})}^2}{|}_{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="HDA0000044564380000041.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}$

28) total kinetic energy of failure removal moment generator:

$x_{28}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\left|V_{\mathrm{KEi}}\right|{|}_{t_{\text{cl}}}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\left|\frac{1}{2}{M}_i{{\mathrm{\&omega;}}_i}^2\right|{|}_{t_{\mathrm{cl}}}$

29) maximum kinetic energy of failure removal moment generator is poor:

$x_{29-1}=(\max \left(V_{\mathrm{ki}}\right)-\min \left(V_{\mathrm{ki}}\right)){|}_{t_{\mathrm{cl}}}$

$x_{29-2}=(\max \left({\stackrel{\&CenterDot;}{V}}_{\mathrm{ki}}\right)-\min \left({\stackrel{\&CenterDot;}{V}}_{\mathrm{ki}}\right)){|}_{t_{\mathrm{cl}}}$

30) the suffered maximum active power of fault generation initial time separate unit generator is impacted:

$x_{30}=\max (P_{i1}-P_{i0}){|}_{t_0}i=\mathrm{1,2},...,N_G$

P in the formula _I1, P _I0Be respectively the meritorious output of generator i when instant of failure and before the fault;

31) the suffered minimum active power of fault generation initial time separate unit generator is impacted:

$x_{31}=\min (P_{i1}-P_{i0}){|}_{t_0}i=\mathrm{1,2},...,N_G$

32) the suffered maximum normalization active power of fault generation initial time separate unit generator is impacted:

$x_{32}=\underset{i}{\max }\frac{P_{i1}-P_{i0}}{P_{\mathrm{ie}}}{|}_{t_0}i=\mathrm{1,2},...,N_G$

P in the formula _IeFor the specified of generator i gained merit.

33) the suffered minimum normalization active power of fault generation initial time separate unit generator is impacted:

$x_{33}=\min \frac{P_{i1}-P_{i0}}{P_{\mathrm{ie}}}{|}_{t_0}i=\mathrm{1,2},...,N_G$

34) characterize the performance characteristic variable of fault to the system shock size:

$x_{34}=\left(\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{M}_i\right|P_{\mathrm{di}}\left|\right)/\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{M}_i$

In the following formula

Be the deceleration power of every generator with respect to the system inertia center; P _CiIt is fault excision moment output of a generator; P _0iIt is output of a generator before the fault;

35) generator accelerating power sum in a flash before accident is removed:

$x_{35}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{P}_{\mathrm{acc},i}^{\mathrm{BF}}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}(P_{\mathrm{mi}}^{\mathrm{BF}}-P_{\mathrm{ei}}^{\mathrm{BF}})=P_{\mathrm{COI}}^{\mathrm{BF}}$

For accident remove after the accelerating power of generator i in a flash;

For accident remove after the accelerating power in the center of inertia in a flash;

The accident that is respectively is removed the mechanical output and the electromagnetic power of back moment generator;

36) maximum of generator accelerating power absolute value in a flash after accident is removed:

$x_{36}=\underset{i}{\max }\left\{\right|P_{\mathrm{acc},i}^{\mathrm{PF}}\left|\right\}=\underset{i}{\max }\left\{\right|P_{\mathrm{mi}}^{\mathrm{PF}}-P_{\mathrm{ei}}^{\mathrm{PF}}\left|\right\}$

37) accelerating power and the generator inertia time constant ratio sum of generator in a flash before accident is removed:

$x_{37}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{P}_{\mathrm{acc},i}/M_i=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\left(\frac{P_{\mathrm{mi}}^{\mathrm{BF}}-P_{\mathrm{ei}}^{\mathrm{BF}}}{M_i}\right)$

38) accelerating power of generator and the maximum of generator inertia time constant ratio in a flash before accident is removed:

$x_{38}=\underset{i}{\max }\left\{\right|P_{\mathrm{acc},i}^{\mathrm{BF}}/M_i\left|\right\}=\underset{i}{\max }\left\{\right|\frac{P_{\mathrm{mi}}^{\mathrm{BF}}-P_{\mathrm{ei}}^{\mathrm{BF}}}{M_i}\left|\right\}$

39) accelerating power and the generator inertia time constant ratio sum of generator in a flash after accident is removed:

$x_{39}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{P}_{\mathrm{acc},i}^{\mathrm{PF}}/M_i=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\left(\frac{P_{\mathrm{mi}}^{\mathrm{PF}}-P_{\mathrm{ei}}^{\mathrm{PF}}}{M_i}\right)$

40) accelerating power of generator and the maximum of generator inertia time constant ratio in a flash after accident is removed:

$x_{40}=\underset{i}{\max }\left\{\right|P_{\mathrm{acc},i}^{\mathrm{PF}}/M_i\left|\right\}=\underset{i}{\max }\left\{\right|\frac{P_{\mathrm{mi}}^{\mathrm{PF}}-P_{\mathrm{ei}}^{\mathrm{PF}}}{M_i}\left|\right\}.$

Wherein, described mass data preprocessing process may further comprise the steps:

(1) data normalization is handled

$x=\frac{x_0-x_{0\max }}{x_{0\max }-x_{0\min }}$

X in the formula ₀, x is respectively the data sample set before and after the conversion; x _0maxBe the sample maximum; x _0minBe the sample minimum value;

(2) sample stipulations

1) according to mixing the F statistical value, determines best cluster numbers C;

2) calculate the information gain of each characteristic variable, as the feature weight W in the cluster analysis to stabilization result _iNormalized sample data is carried out cluster analysis, with each sample be divided into its apart from the corresponding cluster subclass of the cluster centre of minimum, form C data sample subclass;

(3) dimension stipulations

Sample data is carried out logarithm centralization conversion and the vectorial centralization of row,,, drop to p (p＜m) by initial m with the dimension (being the characteristic variable number) of sample set by non-linear Karhunen-Loeve transformation;

1) to initial data X=(x _Ij), i=1 ..., n, j=1 ..., m (n is a sample number, and m is the input feature vector variable number) carries out logarithmic transformation, and the data set after the conversion is designated as W={w _Ij, have

w _ij＝ln?x _ij

2) go vectorial centralization, order

$y_{\mathrm{ij}}=w_{\mathrm{ij}}-\frac{1}{m}\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}$

3) compute matrix Y=(y _Ij) _{N * m}Covariance matrix and characteristic root and the standard orthogonal characteristic vector corresponding with characteristic value

4) characteristic value is carried out descending sort, λ is arranged ₁＞λ ₂＞...＞λ _m

5) according to contribution rate of accumulative total

Select contribution rate of accumulative total＞contribution rate of accumulative total threshold values preceding p (the individual characteristic value of p＜m), general, the contribution rate of accumulative total threshold values gets 0.85～0.95;

6) Karhunen-Loeve transformation forms new data sample set

F _i＝Γ _iWi＝1，2，...，p

Γ _iBe and i eigenvalue of above-mentioned covariance matrix _iCorresponding standard orthogonal characteristic vector.

Wherein, also comprise the following steps:

(1) structure POWER SYSTEM EMERGENCY CONTROL target function:

Max (f (stability margin, MVA _G, MVA _L)) (2)

When system is unstable:

F=-abs (stability margin) (3)

During system stability:

MVA _Gi＝P _gi×r _gi?i＝1，...，N _G

MVA _Lj＝P _lj×r _lj?j＝1，...，N _L

Weight is the weight index for balance and stability nargin value and controlled quentity controlled variable value in the following formula, and general value is 100; MVA _GBe the cutter amount of generator, MVA _LBe the cutting load amount, P _GiCut the meritorious of generator and exerted oneself P _LjBy the burden with power of cutting load.r _GiBe that generator cuts out force rate example, r _LjBe the cutting load ratio, satisfy:

r _Gi=0 or 1 (5)

0≤r _lj≤1 (6)

(2) define generator and load to be excised, comprise the N of faulty line sending end side _GPlatform generator and be subjected to distolateral N _LIndividual load:

(3) initialization N input particle

Each particle is (a N _G+ N _L) n dimensional vector n, preceding N _GDimension data is represented N _GThe state r of platform generator to be cut _g, in the emergency control of reality, generator has only cuts and does not cut fully two states, this N fully _GThe value right and wrong 0 of dimension data are 1; Back N _LDimension data is represented N _LThe individual excision ratio r that waits to excise load _l, the value of each dimension is a real number between [0,1]; The initial position of each particle is set in allowed limits, at random

And speed

(4) corresponding a kind of control measure of particle are utilized power system analysis software (as PSASP) to carry out transient stability and are calculated, and after judgement added the control measure of this particle correspondence, whether system was stable; According to declaring steady result, calculate the target function of each particle correspondence by formula (3) or (4);

(5) utilize mangcorn subgroup method to search optimum controlling schemes.To the k time iteration, upgrade each particle's velocity according to formula (7), to each particle the 1st～N _GThe generator dimension is upgraded corresponding dimension position according to formula (8); To each particle (N _G+ 1)～(N _G+ N _L) the load dimension, upgrade corresponding dimension position value according to formula (9)

$v_i^{(k+1)}=w\&CenterDot;v_i^{\left(k\right)}+a_1\&times;\mathrm{rand}\left(\mathrm{0,1}\right)\&times;({\mathrm{pbest}}_i^{\left(k\right)}-x_i^{\left(k\right)})+a_2\&times;\mathrm{Rand}\left(\mathrm{0,1}\right)\&times;({\mathrm{gbest}}^{\left(k\right)}-x_i^{\left(k\right)})---\left(7\right)$

$\mathrm{If}(\mathrm{rand}\left(\right)<S\left(v_i^{(k+1)}\right))\mathrm{then}{x}_i^{(k+1)}=1$

$\mathrm{else}{x}_i^{(k+1)}=0---\left(8\right)$

$x_i^{(k+1)}=x_i^{\left(k\right)}+v_i^{(k+1)}---\left(9\right)$

W is an inertia weight in the formula, a ₁And a ₂Be two study factors, rand (0,1) and Rand (0,1) they are two random numbers that are evenly distributed between (0,1), i=1, and 2 ..., m.The w definition mode is as follows:

$w=w_{\max }-\frac{w_{\max }-w_{\min }}{{\mathrm{iter}}_{\max }}\&times;\mathrm{iter}---\left(10\right)$

W wherein _Max=0.9, w _Min=0.4, iter is current iterations, and greatest iteration number is made as 100.

The invention has the beneficial effects as follows:

Utilize Knowledge Discovery theory and technology and active data preconditioning technique, can from the operation power data of magnanimity, extract the corresponding relation between operating condition and the electric power system resistance to overturning automatically automatically, find the moving law of electric power system inherence, in order to operational management and the scheduling decision that instructs electric power system reality.This method computational speed is fast, and can provide the quantizating index of power system transient stability, and the canbe used on line of the online dynamic secure estimation of electric power system is had important Practical significance.

The wherein application of data preconditioning technique when effectively reducing individual segregation prediction data space, can keep the classifying quality of initial data preferably, and this shows especially outstandingly under the huge situation of data volume.

The mixed type particle swarm optimization algorithm that the present invention proposes can be realized feasible fast emergency control measure search, can be when guaranteeing system stability, make the control cost minimum of paying, the result of calculation that obtains tallies with the actual situation, and whole computational process is finished automatically.

Description of drawings

The present invention is further described below in conjunction with accompanying drawing.

Fig. 1 is according to optimal classification face of the present invention.

Fig. 2 shows power system stability assessment and control flow.

Fig. 3 shows steady temporarily control mixed type particle swarm optimization algorithm and calculates the particle definitions example.

Fig. 4 shows the mixed type particle swarm optimization algorithm flow process of determining to stablize controlled quentity controlled variable.

Fig. 5 shows transregional electrical network regional area signal.

Fig. 6 shows the mixed type particle swarm optimization algorithm target function convergence curve that transregional electrical network is surely controlled problem.

Fig. 7 shows the element excision amount convergence curve that transregional electrical network is surely controlled problem.

Fig. 8 shows the power-angle curve that transregional electrical network does not have control measure.

Fig. 9 shows the power-angle curve that transregional electrical network has control measure.

Embodiment

Definition

The assessment of transient stability assessment electric power system transient stability belongs to the dynamic analysis problems of large-scale nonlinear dynamic system.Transient stability refers to that the synchronous generator in the interconnected electric power system is suffering serious disturbance (as the transmission line short trouble) back to keep synchronous ability.Main research method has the mixed method of time domain simulation method (or claiming step by step integration), transient energy function method (or claiming direct method) and these two kinds of methods at present.

Knowledge Discovery identifies effective, novel, potentially useful from data, and the level process of final intelligible pattern.The Knowledge Discovery process can be divided into: problem definition, data preparation and preliminary treatment, Knowledge Discovery, and result's explanation and assessment.

The preliminary treatment of data preprocessed data is exactly on the basis of lost data information not, and data are optimized processing, improves the quality of data, satisfies the particular requirement of knowledge discovery algorithm to data, improves the precision and the performance of Knowledge Discovery process thereafter.

The cluster analysis cluster is according to other process of similitude divide into several classes physics or abstract data object, i.e. " things of a kind come together, people of a mind fall into the same group ", its purpose are to make distance between the individuality that belongs to same classification as much as possible little and distance between the individuality that belongs to a different category is big as much as possible.

K-L (Karhunen-Loeve) conversion Karhunen-Loeve transformation is a kind of multivariate statistical method that original a plurality of indexs is converted into a few mutual incoherent overall target, can reach the data abbreviation, disclose relation between variable and the purpose of carrying out statistical interpretation, for character and the statistical nature of further analyzing data provides important information.

SVMs support vector machine method (Support Vector Machines, be called for short SVMs) be that the VC that is based upon Statistical Learning Theory ties up on theoretical and the structure risk minimum principle basis, it is according to limited sample information, between complexity of the model learning accuracy of specific training sample (promptly to) and learning ability (promptly discerning the ability of arbitrary sample error-free), seek optimal compromise, in the hope of obtaining best popularization ability.The basic thought of SVMs is the higher dimensional space at the training sample set place, construct the optimum hyperplane that accurately to divide Various types of data, make that the distance between hyperplane and the inhomogeneity sample set is maximum, promptly have maximum generalization ability, thus the standard of conduct classification unknown sample.The bidimensional situation explanation of concrete available Fig. 1.

In Fig. 1, it is sorting track that square dot and circular point are represented two class samples, H respectively, H ₁, H ₂Be respectively all kinds of in from the nearest sample of sorting track and be parallel to the straight line of sorting track, be called the class interval apart from r between them.So-called optimal classification line requires sorting track not only two class data correctly can be separated exactly, and makes the class interval maximum.The sorting track equation can be written as xw+b=0, and it is carried out normalization, feasible sample set (x to linear separability _i, y _i), x _i(x _i∈ R ⁿ)) be input feature value, y _i(y _i∈ R) be the class object attribute, i=1 ..., N, N are number of samples, satisfy:

y _i[(w·x _i)+b]-1≥0，i＝1，…，N (1)

W is the slope of once linear function in the following formula, and b is the intercept of once linear function.This moment class interval r=2/||w||.Class interval r maximum is equivalent to be made || w|| ²Minimum.Satisfy condition (2-3) and make

Minimum classifying face is called optimal classification face, H ₁, H ₂On the training sample point be known as support vector.

The problem of finding the solution the optimal classification face can be expressed as following constrained optimization problem:

Min(Φ(w))＝Min(||w|| ²/2)＝Min((w·w)/2) (2)

s.t.y _i(w·x _i+b)-1≥0，i＝1，...，N

Definition Lagrange function:

$L(w,b,\mathrm{\&alpha;})=(w\&CenterDot;w)/2-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i[y_i(w\&CenterDot;x_i+b)-1]---\left(3\right)$

α wherein _i＞0 is the Lagrange coefficient.Utilize the Lagrange optimization method to be converted into its dual problem (4) to above-mentioned optimal classification face problem, promptly find the solution α ₁... α _N, make

$\mathrm{Max}\left(Q\left(\mathrm{\&alpha;}\right)\right)=\mathrm{Max}(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i-\frac{1}{2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\mathrm{\&alpha;}}_j{y}_i{y}_j(x_i\&CenterDot;x_j))---\left(4\right)$

$s.t.\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{y}_i=0,{\mathrm{\&alpha;}}_i\&GreaterEqual;0,i=1,...,N$

Finding the solution the optimal classification function that obtains after the problems referred to above is

$w^{*}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{*}{y}_i{x}_i---\left(5\right)$

b ^＊＝y _k-w ^T?x _kx _k (6)

$f\left(x\right)=\mathrm{sgn}\{(w^{*}\&CenterDot;x)+b^{*}\}=\mathrm{sgn}\{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{*}{y}_i(x_i\&CenterDot;x)+b^{*}\}---\left(7\right)$

Sgn () is a sign function.b ^*Be classification thresholds, can try to achieve with any support vector (satisfying the equal sign in (1)), or get intermediate value by any a pair of support vector in two classes and try to achieve.

Top method is guaranteeing that training sample is all correctly classified, and obtains best popularization performance by the maximization class interval.Consider to exist some data can not satisfy the condition of formula (1), promptly linear inseparable, need this moment in constraints, to add a slack ξ i, problem becomes:

Min(1/2w ^Tw+CΣξi) (8)

s.t.yi(w ^Tx _i+b)≥1-ξi，i＝1，...，N

ξi≥0，i＝1，...，N

C is the constant of certain appointment, plays control divides sample punishment to mistake degree.Optimization problem herein separate with can the branch situation under separate almost completely identically, just condition becomes

0≤α _i≤C，i＝1，...，N (9)

To nonlinear problem, can the input space be transformed to a higher dimensional space by the nonlinear transformation K (x, x ') that satisfies the Mercer condition, in this new space, ask for the optimal classification face then.This nonlinear transformation is to realize by defining suitable inner product function.The majorized function of this up-to-date style (4) becomes:

$\mathrm{Max}\left(Q\left(\mathrm{\&alpha;}\right)\right)=\mathrm{Max}(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i-\frac{1}{2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{\mathrm{\&alpha;}}_j{y}_i{y}_{j}K(x_i,x_j))---\left(10\right)$

Corresponding discriminant function formula (7) also becomes

$f\left(x\right)=\mathrm{sgn}\{(w^{*}\&CenterDot;x)+b^{*}\}=\mathrm{sgn}\{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{*}{y}_{i}K(x_i\&CenterDot;x)+b^{*}\}---\left(11\right)$

Inner product kernel function K (x in the formula _i, x) adopt Gauss's radial kernel function

$K(x_i,x)=\exp (-\frac{{|x-x_i|}^2}{{\mathrm{\&sigma;}}^2})---\left(12\right)$

σ is the width parameter of function in the formula, has controlled the radial effect scope of function.

Transient Stability Control is in time carried out effective Transient Stability Control after electric power system suffers big disturbance, often can guarantee the safe and stable operation of system with less cost.Control measure commonly used mainly comprise cutter, cutting load, quick closing valve porthole, electric braking, forced compensation etc.

Particle group optimizing particle group optimizing method (Particle Swarm Optimization, be called for short PSO) essence is a kind of many agent algorithms, group that research is made up of simple individuality and the mutual-action behavior between environment and the individuality.The transition process of population is directive, use feedback principle in the search procedure and can adopt parallel computing, have higher search efficiency, can find the globally optimal solution of problem with big probability, computational efficiency is than traditional random device height, fast convergence rate.

In the PSO algorithm, represent separating of problem to be optimized with particle position, the good and bad degree of each particle performance depends on the adaptive value that problem target function to be optimized is determined, each particle determines its heading and speed size by a speed.Be located in the target search space of a d dimension, form a colony, be expressed as the position of particle i constantly the k time iteration by N particle Corresponding flying speed is expressed as

When beginning to carry out the PSO algorithm, at first need a random initializtion N particle position and speed, seek optimal solution by iteration then.With the target function of particle substitution problem definition, the quality of particle is to weigh according to the size of target function adaptive value.In iteration each time, each particle upgrades own Position And Velocity in solution space by following the tracks of two extreme values: extreme value is this self optimal solution that searches so far in iterative process of particle, is called individual extreme value, and its position is expressed as pbest _i=(pbest _I1, pbest _I2..., pbest _Id), another extreme value is the optimal solution that up to the present whole population finds, and is called global extremum, its position is expressed as gbest=(gbest ₁, gbest ₂..., gbest _d).

When the k+1 time iterative computation, particle x _iUpgrade oneself speed and position according to following rule:

(1) if the x value is continuous,

$v_i^{(k+1)}=w\&CenterDot;v_i^{\left(k\right)}+a_1\&times;\mathrm{rand}\left(\mathrm{0,1}\right)\&times;({\mathrm{pbest}}_i^{\left(k\right)}-x_i^{\left(k\right)})+a_2\&times;\mathrm{Rand}\left(\mathrm{0,1}\right)\&times;({\mathrm{gbest}}^{\left(k\right)}-x_i^{\left(k\right)})---\left(13\right)$

$x_i^{(k+1)}=x_i^{\left(k\right)}+v_I^{(k+1)}---\left(14\right)$

W is an inertia weight in the formula, a ₁And a ₂Be two study factors, rand (0,1) and Rand (0,1) they are two random numbers that are evenly distributed between (0,1), i=1, and 2 ..., m.In (14), speed v _iBe counted as correction to the position, more suitable.In addition, particle rapidity v _iTie up all by a maximal rate v at each _MaxLimit.If the acceleration of current particle causes it to surpass maximal rate on this dimension in the speed of certain one dimension, then the speed of this dimension is restricted to maximal rate.

(2) if the x value is discrete, particle position more new formula (14) is replaced by following formula (15),

$\mathrm{If}(\mathrm{rand}\left(\right)<S\left(v_i^{(k+1)}\right))\mathrm{then}{x}_i^{(k+1)}=1$

$\mathrm{else}{x}_i^{(k+1)}=0---\left(15\right)$

In binary system PSO algorithm, each dimension x of particle _IdEach dimension pbest with individual extreme value _IdBe restricted to 1 or 0, and to speed v _IdDo not do this restriction.When upgrading the position, if v with speed _IdHigher, particle position x _IdMore likely select 1; Otherwise, if v _IdHang down a bit, then x _IdSelect 0.The function S of threshold value between [0,1] in the formula (15) (v) be the sigmoid function, S (v)=1/ (1+e ^-v).

Fig. 2 is stable assessment of the present invention and control flow chart.Stable assessment of the present invention and control method comprise the steps:

Step 1 definition input feature vector variable is set up the electric power system transient stability assessment models.Electric power system transient stability assessment models feature is divided into steady state characteristic and behavioral characteristics, and behavioral characteristics is divided into merit corner characteristics, angular speed feature, angular acceleration feature and power, energy feature.The characteristic variable of the present invention's definition is as follows:

The step 2 data acquisition, sample generates.Among the present invention, utilize time-domain-simulation that the actual electric network system model is provided with various forecast accidents, carry out a large amount of simulation calculation, utilize initial steady state flow data and fault that 1～2 second interior multidate information of finite time in back takes place, produce simulation sample, therefrom extract the characteristic quantity data that need then.

The preliminary treatment of step 3 data.The operation of power networks characteristic data dimension that collects is inconsistent, and each characteristic value size differences is big.The feature that value is little in the computational process is easily flooded by the big data of value.In addition, measured data is accumulated in time, and data volume is big.By the data preconditioning technique, can effectively improve the quality of data, reduce the subsequent calculations required time.The present invention mainly adopts following data preconditioning technique:

1. data normalization is handled the characteristic bi-directional scaling, makes it to fall into (as 0.0 to 1.0) between a little given zone.The feature that standardization processing can effectively prevent to have big initial codomain is compared the weight problems of too of appearance with the feature with less initial codomain, can avoid the influence to result of calculation of the dimension of each characteristic variable and the order of magnitude.To raw sample data, the present invention carries out standardization processing by following formula:

$x=\frac{x_0-x_{0\max }}{x_{0\max }-x_{0\min }}---\left(16\right)$

X in the formula ₀, x is respectively the data sample set before and after the conversion; x _0maxBe the sample maximum; x _0minBe the sample minimum value

2. the data reduction is compressed data set, obtains much smaller data representation, but can produce same analysis result.Usually, for the data that exist with the bivariate table form, successively carry out stipulations from horizontal sample dimension and vertical characteristic dimension respectively.

(1) the present invention of sample stipulations utilizes the cluster analysis technology that service data is carried out the sample stipulations.Can guarantee under the constant situation of analysis result original large data sets to be divided into some less data sets, newly-generated small data set, each self-contained sample is similar each other, and is different each other between each small data set.But each small data set independent parallel is handled, to improve computational efficiency.Clustering problem can be expressed as an iteration optimizing problem.

If the sample space X={x after the standardization that constitutes by n sample, a m characteristic variable _Ij, x _IjBe the value of sample j on characteristic variable i, i=1,2 ..., m, j=1,2 ..., n.If the fuzzy clustering center matrix is S={s _Hi} _{C * m}, s _HiBe the cluster centre of i variable in the h class sample, h=1,2 ..., c.

If sample space X is gathered into the C class.Best cluster value C determines according to mixing the F statistical value.

$\mathrm{Mixed}\_F=\frac{m}{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}1/F\left(i\right)}---\left(17\right)$

$F\left(i\right)=\frac{\underset{h=1}{\overset{c}{\mathrm{\&Sigma;}}}n{(s_{\mathrm{hi}}-{\stackrel{\&OverBar;}{s}}_i)}^2(n-c)}{\underset{h=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n_h}{\mathrm{\&Sigma;}}}{(x_{\mathrm{ijh}}-s_{\mathrm{hi}})}^2(c-1)}(i=\mathrm{1,2},...,m)---\left(18\right)$

N in the formula _hIt is h class sample number;

Be the mean value of i variable cluster centre; x _IjhBe i variate-value of j sample in the h cluster.The Mixed_F value is big more, and the cluster effect is good more.

Ambiguity in definition cluster matrix U={ u _Hj, u _HjFor sample j belongs to the relative degree of membership of classification h, satisfy

0≤u _hj≤1

$\underset{h=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{hj}}=1---\left(19\right)$

$\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{hj}}>0$

If the fuzzy clustering center matrix is S={s _Ih, s _IhBe the standardization value of cluster centre h on characteristic variable i, i=1,2 ..., m, h=1,2 ..., c.Consider the difference of each characteristic variable, establish weight vectors W=(w the cluster effect ₁..., w _m) ^T, and satisfy

$\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{w}_i=1---\left(20\right)$

The fuzzy clustering problem can be expressed as an iteration optimizing problem:

$\min \{F\left(u_{\mathrm{hj}}\right)=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{h=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left[u_{\mathrm{hj}}\right|\left|w_i(x_{\mathrm{ij}}-s_{\mathrm{ih}})\right|\left|\right]}^2\}---\left(21\right)$

And satisfy constraints (19) and (20).Separating of the problems referred to above can obtain the iteration expression formula of U battle array and S battle array by the structure Lagrangian:

$u_{\mathrm{hj}}=\frac{1}{\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left[w_i(x_{\mathrm{ij}}-s_{\mathrm{ih}})\right]}^2}{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left[w_i(x_{\mathrm{ij}}-s_{\mathrm{ik}})\right]}^2}}---\left(22\right)$

$s_{\mathrm{ih}}=\frac{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{hj}}^2{w}_i^2{x}_{\mathrm{ij}}}{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{hj}}^2{w}_i^2}---\left(23\right)$

For the transient stability assessment, be one two (or many) pattern classifications problem, each characteristic variable has inherent law to the effect degree of cluster analysis.In information theory, feature with the compression of the highest information gain or maximum entropy can make the required amount of information minimum of classification, contribution maximum to classification problem, the present invention is with based on the information gain of the entropy tolerance as feature weight, and can avoid weights W is the problem that the artificial given result of calculation that causes is influenced by subjective factor.

If X is the set that comprises n data sample, in temporarily steady evaluation problem, all samples can be divided into stable and unstable two classes, are designated as C _i(i=1,2).n _iBe to belong to class C among the X _iSample number, characteristic variable A has v different value { a ₁, a ₂..., a _v, can X be divided into v subclass { X with attribute A ₁, X ₂..., X _v, wherein, X _jComprise sample such among the X, their values on attribute A are a _jThe information gain of characteristic variable A is expressed as:

$W_A=\mathrm{Gain}\left(A\right)$

$=-\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}{p}_i{\log }_2\left(p_i\right)-\underset{j=1}{\overset{v}{\mathrm{\&Sigma;}}}\frac{n_{1,j}+{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="HDA0000044564380000011.png"\; state="\{\{state\}\}">n</figure-callout>}}_{2,j}}{n}(-\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}{\log }_2\left(p_{\mathrm{ij}}\right))---\left(24\right)$

P in the formula _iBe that arbitrary sample belongs to C _iProbability, p _IjIt is subset X _jIn sample belong to class C _iProbability,

|| the sample number of expression set, n _IjIt is subset X _jIn belong to class C _iSample number.

Find the solution above-mentioned optimization problem, obtain fuzzy clustering center matrix S.For each operation sample x _i(i=1 ..., n),, be divided in the cluster subclass of the minimum cluster centre correspondence of distance with it according to the distance of itself and each cluster centre vector.

(2) the dimension stipulations are deleted or are reduced incoherent characteristic attribute, reach the purpose that improves computational efficiency.The present invention utilizes improved non-linear Karhunen-Loeve transformation method to carry out the dimension stipulations of transient stability assessment.Karhunen-Loeve transformation is a kind of original a plurality of indexs to be converted into the multivariate statistical method of a few mutual incoherent overall target, can reach the data abbreviation, disclose relation between variable and the purpose of carrying out statistical interpretation.

In practical problem, often be non-linear relation between the characteristic variable and between overall target and the initial data, therefore be necessary " linearisation " in traditional Karhunen-Loeve transformation improved.For by n sample, m the data set X=(x that feature constitutes _Ij) _{N * m}, the present invention adopts the non-linear Karhunen-Loeve transformation method of " logarithm centralization ", and main way is as follows:

(1) to initial data x _IjCarry out logarithmic transformation, the data set after the conversion is designated as W={w _Ij, have

w _ij＝lnx _ij (25)

(2) go vectorial centralization, order

$y_{\mathrm{ij}}=w_{\mathrm{ij}}-\frac{1}{m}\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}---\left(26\right)$

(3) compute matrix Y=(y _Ij) _{N * m}Covariance matrix V, find the solution the characteristic value of V and corresponding standard orthogonal characteristic vector, according to contribution rate, determine to satisfy contribution rate of accumulative total

Preceding p (the individual overall target of p＜m);

(4) set up the overall target equation, calculate comprehensive index value, form new data set

F _i＝Γ _iW?i＝1，2，...，p (27)

In the following formula, Γ _iBe and i eigenvalue of above-mentioned covariance matrix _iCorresponding standard orthogonal characteristic vector.

The input feature vector variable that utilizes Karhunen-Loeve transformation to choose can drop to p (p＜m), help to improve the computational efficiency of transient stability assessment by m with the dimension of data set.

Step 4 transient stability assessment the present invention adopts the support vector machine technology to carry out towards the Knowledge Discovery of transient stability assessment.The purpose of transient stability assessment is according to system's service data of collecting, determines that corresponding system running state is stable or unsettled.To the history data of collecting, at first utilize algorithm of support vector machine, determine to distinguish optimal classification line (or face) mathematic(al) representation (11) of steady and unstable operating point in the sample space.

For the service data x of running status the unknown,, can obtain its running status with its substitution formula (11).The output of SVMs has represented test data to the distance that hyperplane is arranged most, exports the big more hyperplane that shows apart from far away more.Make that SVMs is that system stability is represented in timing, then should the big more systems of number stable more, vice versa.

The step 5 Transient Stability Control is from the angle of mathematics, and Transient Stability Control belongs to discrete control category; From the angle of optimization subject, can be summed up as a some extreme value problem, promptly, how to provide best controlled quentity controlled variable to improving this controlled target of power system transient stability.

(1) Mathematical Modeling of steady control problem solved by the invention can be described as:

1) target function

Max (f (stability margin, MVA _G, MVA _L)) (28)

Target function f is system stability nargin, cutter amount MVA _GWith cutting load amount MVA _LFunction, be according to steady control amount and implement stablize control measure after system mode definite.Obtain system running state according to above transient stability appraisal procedure, when system's instability, target function lays particular emphasis on the variation of stability margin; When system mode return to stable after, target is when guaranteeing system stability, to find the control mode of excision amount minimum, so take into account the factor of cutter cutting load amount in the target function at this moment.Target function defines suc as formula (29) and formula (30):

When system is unstable:

F=-abs (stability margin) (29)

During system stability:

MVA _Gi＝P _gi×r _gi，i＝1，...，N _G

MVA _Lj＝P _lj×r _lj，j＝1，...，N _L

Weight is the weight index for balance and stability nargin value and controlled quentity controlled variable value in the following formula, and general value is 100; MVA _GBe the cutter amount of generator, MVA _LBe the cutting load amount, P _GiCut the meritorious of generator and exerted oneself P _LjBy the burden with power of cutting load.r _GiBe that generator subtracts the ratio of exerting oneself, r _LjIt is the cutting load ratio.

2) constraints

r _Gi=0 or 1 (31)

0≤r _lj≤1 (32)

Before optimizing, need definition to wait to ask separating of problem, i.e. particle in the PSO algorithm.For temporarily steady controlled quentity controlled variable problem identificatioin, it is the vector of one (NG+NL) dimension that each particle is expressed as.The 1st～N _GDimension is the generator dimension, N _G+ 1～N _G+ N _LDimension is the load dimension.Because in the emergency control of reality, generator has only to be cut fully and not to cut two states fully, and so the problem of what not being cut is this N _GThe value right and wrong 0 of dimension data are 1.Back N _LDimension data is represented N _LIndividual excision ratio of waiting to excise load, the value of each dimension are real numbers between [0,1].

(2) Transient Stability Control process

1) definition generator and load to be excised.In principle, all generators in the system and load all can be used as controlled device, but do like this, do not meet actual conditions on the one hand, also can cause huge computation burden on the other hand.Must clearly excise the scope of element at the position of fault point for this reason.Be to excise the sending end side generator of faulty line and be subjected to distolateral load in principle.

2) initialization N input particle

Each particle is (a N _G+ N _L) n dimensional vector n.The initial position of each particle is set in allowed limits, at random

And speed

3) corresponding a kind of control measure of particle, the disturbance input file (as the ST.S12 in the PSASP software) in the corresponding temporarily steady calculation procedure.According to the content of each particle, write the disturbance input file, the temporarily steady program of operation, after judgement added the control measure of this particle correspondence, whether system was stable.

4) according to declaring steady result, calculate the target function of each particle correspondence by formula (29) or (30).

5) when k=1, i.e. iteration for the first time, each particle

Coordinate is set to its current location,

Calculate its corresponding individual extreme value (being the target function value of individual extreme point); When k＞1, then calculate the target function value of particle, if be better than the current individual extreme value of this particle, Then pbest is set to this particle position, promptly

And upgrade current individual extreme value; To each particle, when k=1, i.e. iteration for the first time, global extremum (being the target function value of global extremum point) is best in all individual extreme values, gbest is set to the current location of this best particle; When k＞1, if preferably in the individual extreme value of all particles be better than current global extremum, then gbest is set to this particle position,

And upgrade current global extremum;

6) upgrade each particle's velocity according to formula (13), upgrade each particle position according to formula (14) and (15); If

Then

If

Then

7) whether the cycle-index of interpretation problems has surpassed maximum number of times restriction, if surpass, then stops iteration, otherwise turns to step 3);

8) the overall utmost point figure of merit of output, this value be exactly each machine and each load should the excision amount.

Below be one embodiment of the present of invention, carry out emulation experiment with certain transregional electrical network and be embodiment, further specify as follows:

(1) sample generates

Adopting certain transregional electrical network annual running mode to calculate and using real data, scale of power is 4936 nodes, 660 generator bus, 1646 load buses, 4290 alternating current circuits and 3 DC line.Under a certain operational mode, calculate the through-put power of each bar circuit, preceding 10 circuits of getting the through-put power maximum are as critical circuits; Preceding 8 generators of the maximum of taking-up power similarly, are as crucial generator.Add up to and defined 86 static nature variablees and 42 temporarily steady characteristic variables.The scene provides 15 kinds of trend operational modes, has comprised 39 groups of three phase short circuit fault in the given error listing, by excision faulty line excision fault.Get 14 kinds of trend operational modes at random, calculate 80 seconds whole process simulation of 14 * 39=546 temporarily steady operation of its correspondence, extract the quiet/behavioral characteristics variable of fault front and back and corresponding temporarily steady result of calculation (stable state), the formation training data.Utilize one group of remaining running mode data, calculate 80 seconds whole process simulation of 39 temporarily steady operations of its correspondence, obtain corresponding characteristic vector value, constitute test data.What deserves to be mentioned is that test data does not participate in the training process of SVMs, belong to the service data of new the unknown.

(2) data preliminary treatment

1), the sample data of collecting is carried out standardization processing according to formula (16);

2) the sample stipulations are established cluster numbers C and are changed to 10 from 2, calculate Mixed_F, the results are shown in Table 1.

The Mixed_F value of the different cluster numbers correspondences of table 1

As can be seen from Table 1, when cluster numbers was 3, the value maximum of Mixed_F correspondence was so can gather into the training set sample 3 classes.After the cluster, the sample number in each class is respectively 147,187 and 212, and sample space is respectively 73.08%, 65.75% and 6117% of a former data set.For each the operation sample in the test set, distance according to itself and each cluster centre vector, be divided into it in the cluster subclass of the minimum cluster centre correspondence of distance, the sample number in each group test set is respectively 7,16 and 16, has so just formed 3 groups of training-test subclass.

3) the dimension stipulations are to then carrying out the dimension stipulations through the data of sample stipulations.In nonlinear method, the number of features in each data subset is reduced to 25,31 and 28 by 128.

(3) transient stability assessment

3 groups of training subclass (comprising 147,187 and 212 training samples respectively) are trained the temporarily steady evaluation prediction model of each self-structuring.3 groups of training subclass 12,100,14 stable support vectors of model difference and 6,11,7 unstable support vectors are formed.Support vector is the data that really work when prediction, and its quantity is much smaller than the quantity of training sample.Can executed in parallel to the training-test of three groups of subclass, get computing time wherein the most consuming time one group as final computing time.

Than the time-domain-simulation method, based on the transient stability assessment of knowledge discovering technologies, maximum advantage is that its computational speed is fast.Utilize two kinds of methods, successively state's adjusting data is carried out simulation calculation, the result is shown in table 2 and table 3:

Table 2 carries out the pretreated stability forecast result of data

Table 3 Knowledge Discovery method and time-domain-simulation method compare computing time

Whether stable from table 2 and table 3 as can be seen, obtain system's operation conclusion, time-domain-simulation calculates needs 40 seconds dynamic processes of emulation, calculates 116.64 seconds consuming time; If the working knowledge discover method only needs 13.812 seconds.Forecasting accuracy sees Table 4.

Table 4 predicts the outcome based on the transient stability of SVMs model

The SVMs value is stable for just representing current temporarily steady operation, be the negative indication instability, and the absolute value of SVMs value is big more, and the temporarily steady operation of expression correspondence is big more apart from the distance of SVMs classification hyperplane; It is stable to stablize classification 1 expression, and 2 expressions are unstable.In temporarily steady evaluation problem, mistake is divided and is referred to be judged to be instability with stablizing sample, leaks to divide to refer to be judged to be unstable sample stable.As shown in table 4: as to have 1 operation to be leaked branch.Leak to divide the SVMs value of operation to equal 0.5008, near hyperplane, i.e. y=0.0.Those are in the sample on neutrality border, and available picking out utilizes the time-domain-simulation method to carry out detailed transient stability analysis, to guarantee all data final analysis results' reliability.

(4) Transient Stability Control

Fig. 5 is transregional electrical network regional area schematic diagram, and the circle among the figure is represented each power generation region or electricity consumption zone, and directed line segment is represented interregional interconnection and direction of tide thereof.Wherein connect by AC line 3600 between the zone 1 and 2.Be located on the interconnection 54 and 58 between the zone 1 and 10, near regional 10 places, three-phase shortcircuit ground connection permanent fault take place, fault took place back 0.2 second, open-circuit line:

Interconnection 53 and 57 between zone 10 and the zone 11;

Interconnection 56 between zone 1 and the zone 7;

Circuit 59 in the zone 1.

Simulation calculation shows, cut-offs above-mentioned circuit and can not make system restoration stable.Now adopt and mix the control measure that the PSO algorithm search can guarantee system stability.

1) control area is determined

For the big system of actual motion, if with all generators in the system and load all as the excision object, promptly each generator and the load in each dimension data of particle correspondence system in the PSO algorithm will reduce algorithm efficiency.Before carrying out the PSO search, the scope of necessary clear and definite search volume, and this scope is as much as possible little.In the present invention, the search volume is made up of cut element (cause the unbalanced sending end generator of power and be subjected to end load).In this example, determine that finally the generator in power sending end zone (11,33,34,37) and the load bus that power is subjected to end regions 2 constitute the initial ranging space of PSO algorithm.

2) PSO algorithm data initialization

Initially consider not cutting load,, then consider the excision load again if only depend on the excision generator can't make system restoration stable.Therefore, calculate in the particle, the corresponding position of generator initially puts 1, and the corresponding position of load initially puts 0

The Control Parameter of mixing in the PSO algorithm is chosen as: population N=20, study factor a ₁=a ₂=2.0, velocity constraint v _Max=1.0, v _Min=-1.0, iterations=100.Inertia weight w adopts the definition mode of formula (33):

$w=w_{\max }-\frac{w_{\max }-w_{\min }}{{\mathrm{iter}}_{\max }}\&times;\mathrm{iter}---\left(33\right)$

W in the formula _Max=0.9, w _Min=0.4, iter is current iterations, and greatest iteration number is made as 100.

3) stablizing control strategy determines

Carry out institute's determining step, obtain making stable final cutter of power system restoration and cutting load amount scheme.Set the action after fault takes place 0.22 second of all steady control means.Table 5 has provided the final steady control result that operation PSO algorithm obtains.

Table 5 makes transregional network system recover stable cutter and cutting load amount

Fig. 6 represents target function iteration convergence curve, and Fig. 7 represents element excision amount convergence curve.From Fig. 6 and Fig. 7, as seen, utilize mixing PSO algorithm to determine the emergency control regulation scheme soon, in the present embodiment, just can converge to stationary value 10 times transregional network system iteration.Fig. 8 and Fig. 9 represent to add before and after the control measure, the time-domain response result of generator's power and angle.From Fig. 9 as seen, add steady control means after, it is stable that the place system can recover in 2 seconds, and merit angle unstable phenomenon as shown in Figure 8 can not take place.In PSO algorithm of the present invention, each particle in the iteration is represented a kind of steady control means each time, its fitness function value, and just steady control means all will calculate by transient stability the influence of the stability of a system.Resemble the so big system of transregional electrical network for one, obtain final optimal and separate, need repeatedly call transient stability computation program.If all carry out 40 seconds overall process time-domain-simulations at every turn, the required time of so whole calculating will be very considerable.In order to improve computational efficiency, in the steady control emulation to transregional electrical network, the present invention adopts the SVMs method to replace time-domain-simulation, carries out the quick judgement of running status.The SVMs value that obtains as stability margin substitution formula (29) or (30), is calculated the fitness target function, as the foundation of further iteration.Each particle, each iteration if declare surely with the time domain simulation method, needs with 34 seconds; If declare surely, only need time less than 2 seconds with the SVMs method.If the population size is 15, carry out iteration 10 times, obtain finally to separate, operation time-domain-simulation 34 seconds * 15 * 10=5100 consuming time second=85 minute, and operation SVMs method needs 2 seconds * 15 * 10=300 second=5 minute.

Claims (4)

1. a power system stability is assessed and control method, it is characterized in that may further comprise the steps:

(1) determines the input feature vector variable of electric power system, make up the transient stability assessment models of electric power system;

(2) obtain electrical network on-line operation data and the simulation calculation of the transient stability in short-term result under given fault collection thereof in the electric power system, according to the input feature vector variable that step (1) is determined, the composing training sample space, each sample standard deviation is by input feature vector variable set x ₀Temporarily steady result of calculation y with correspondence ₀Two parts constitute;

(3) on the basis of lost data information not, initial mass data sample is carried out the data preliminary treatment, obtain new data sample set (x, y);

(4) according to training sample set (x, y), structure transient stability optimal classification discriminant function:

$f\left(x\right)=\mathrm{sgn}\{(w^{*}\&CenterDot;x)+b^{*}\}=\mathrm{sgn}\{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i^{*}{y}_{i}K(x_i\&CenterDot;x)+b^{*}\}---\left(1\right)$

Wherein sgn () is a sign function, w ^*Be the slope of once linear function, b ^*Be the intercept of once linear function, α _iBe the Lagrange multiplier corresponding, K (x with each sample _i, be Gauss's radial kernel function x), have

$K(x_i,x)=\exp (-\frac{{|x-x_i|}^2}{{\mathrm{\&sigma;}}^2})$

Wherein σ is the width parameter of function, has controlled the radial effect scope of function;

(5) for the service data x of operation states of electric power system the unknown, after (3) passed through the data preliminary treatment set by step with it, substitution formula (1) promptly obtained its corresponding transient stability running status f (x); The value of f (x) has been represented the distance of each sample to the optimal classification hyperplane, and in the transient stability evaluation problem, the value of f (x) has been represented the distance of current running status x to the transient stability classifying face, i.e. transient stability nargin.F (x) value represent system stability, and the big more system of numerical value is stable more for just; Otherwise value is a negative indication system transient state instability, and this number big more expression of absolute value system is unstable more;

(6) step (5) is judged to be the power system operation mode of unstability, utilizes the mixed type particle swarm optimization algorithm, search for and obtain to make the emergency control measure of operation of power system recovery transient stability and cutter or cutting load amount minimum.

2. the method for claim 1 is characterized in that definite input feature vector x that can reflect transient characterisitics when electric power system is broken down ₀And corresponding temporarily steady result of calculation y ₀, make up the electric power system transient stability assessment models; The input feature vector variable x that step (1) adopts ₀Being divided into static nature variable and behavioral characteristics variable, also is x ₀=(static nature variable, behavioral characteristics variable) specifically comprises:

(1) the static nature variable comprises

1) generator of whole electric power system exerting oneself and loading;

2) the minimum and maximum generated output of electric power system;

3) the minimum and maximum load of electric power system;

4) the minimum and maximum value of busbar voltage;

5) the maximum meritorious through-put power of circuit, and corresponding idle through-put power;

6) key transmission cross-section and zone contact section) through-put power;

7) generated output of crucial generator;

8) each regional generator exerting oneself and loading;

(2) the behavioral characteristics variable comprises:

Below various in, δ _iBe the generator amature angle; ω _iBe generator speed, α _iBe the generator amature angular acceleration,

M _iInertia constant for generator; P _MiBe the generator mechanical output; P _EiBe generator electromagnetic power; V _KiBe generator energy; V _PiBe generator potential energy; N _GBe generator number, t ₀Be fault generation initial time, t _ClBe the failure removal moment;

1) fault generation initial time has the relative rotor angle of the generator of maximum work angular acceleration:

$x_1=({\mathrm{\&delta;}}_{a\max }-{\mathrm{\&delta;}}_{\mathrm{COI}}){|}_{t_0}$

In the following formula

Be equivalent rotor angle with respect to the system inertia center, δ _AmaxIt is the rotor angle that fault generation initial time has the generator of maximum work angular acceleration;

2) failure removal has the relative rotor angle of the generator of maximum kinetic energy constantly:

$x_2=({\mathrm{\&delta;}}_{p\max }-{\mathrm{\&delta;}}_{\mathrm{COI}}){|}_{t_{\mathrm{cl}}}$

δ in the following formula _PmaxIt is the rotor angle that the fault excision has the generator of maximum kinetic energy constantly;

3) the merit angle of the machine and the machine of bringing up the rear was poor before failure removal was led constantly; Fault generation initial time, the generator with maximum rotor angular acceleration is the preceding machine of neck; Generator with minimum rotor angle acceleration is the machine of bringing up the rear;

$x_3=\max \left({\mathrm{\&delta;}}_{a\max }\right)-\min \left({\mathrm{\&delta;}}_{a\min }\right){|}_{t_{\mathrm{cl}}}$

δ in the following formula _Amax, δ _AminIt is respectively the rotor angle of the preceding machine of neck and the machine of bringing up the rear;

4) the maximum constantly generator outer corner difference of failure removal:

$x_4=\max \left({\mathrm{\&delta;}}_i\right)-\min \left({\mathrm{\&delta;}}_i\right){|}_{t_{\mathrm{cl}}}$

5) failure removal differs maximum merit angle with center of inertia COI constantly:

$x_5=\max ({\mathrm{\&delta;}}_i-{\mathrm{\&delta;}}_{\mathrm{COI}}){|}_{t_{\mathrm{cl}}}$

The inertia weighted average angle of the generator before 6) failure removal is led constantly in a group of planes and the remaining group of planes poor:

$x_6=\underset{i\&Element;s}{\mathrm{\&Sigma;}}{M}_i{\mathrm{\&delta;}}_i/\underset{i\&Element;s}{\mathrm{\&Sigma;}}{M}_i-\underset{j\&Element;A}{\mathrm{\&Sigma;}}{M}_j{\mathrm{\&delta;}}_j/\underset{j\&Element;A}{\mathrm{\&Sigma;}}{M}_j{|}_{t_{\mathrm{cl}}}$

Wherein, group of planes set before S represents to lead, A represents remaining group of planes set.To specifying each generator corner constantly to carry out descending sort, shape as:

$({\mathrm{\&delta;}}_1,{\mathrm{\&delta;}}_2,{\mathrm{\&delta;}}_3,...,{\mathrm{\&delta;}}_{N_G-2},{\mathrm{\&delta;}}_{N_G-1},{\mathrm{\&delta;}}_{N_G}),$ And ${\mathrm{\&delta;}}_1>{\mathrm{\&delta;}}_2>{\mathrm{\&delta;}}_3>...>{\mathrm{\&delta;}}_{N_G-2},{\mathrm{\&delta;}}_{N_G-1}>{\mathrm{\&delta;}}_{N_G}$

The difference of per two adjacent corners constitutes an angular clearances, shape as:

$({\mathrm{\&delta;}}_{\mathrm{1,2}},{\mathrm{\&delta;}}_{\mathrm{2,3}},...,{\mathrm{\&delta;}}_{N_G-2,N_G-1},{\mathrm{\&delta;}}_{N_G-1,N_G}),$ And

${\mathrm{\&delta;}}_{\mathrm{1,2}}={\mathrm{\&delta;}}_1-{\mathrm{\&delta;}}_2,{\mathrm{\&delta;}}_{\mathrm{2,3}}={\mathrm{\&delta;}}_2-{\mathrm{\&delta;}}_3,...,{\mathrm{\&delta;}}_{N_G-2,N_G-1}={\mathrm{\&delta;}}_{N_G-2}-{\mathrm{\&delta;}}_{N_G-1},{\mathrm{\&delta;}}_{N_G-1,N_G}={\mathrm{\&delta;}}_{N_G-1}-{\mathrm{\&delta;}}_{N_G}$

Determine maximum angular clearances, the generator on this gap constitutes the preceding group of planes S of neck, and the generator under this gap constitutes remaining group of planes A;

7) failure removal constantly with the absolute value sum of the relative rotor angle difference of fault generation initial time generator:

$x_7=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}|{\mathrm{\&theta;}}_{\mathrm{ci}}-{\mathrm{\&theta;}}_{0i}|$

In the following formula, θ _i=δ _i-δ _COI, θ _cBe the relative rotor angle of failure removal moment generator, θ ₀Be the relative rotor angle of fault generation initial time generator.

8) failure removal constantly with the maximum of the absolute value of the relative rotor angle difference of fault generation initial time generator:

$x_8=\underset{i}{\max }\left(\right|{\mathrm{\&theta;}}_{\mathrm{ci}}-{\mathrm{\&theta;}}_{0i}\left|\right)$

9) the failure removal moment and fault generation initial time center of inertia coordinate COI change in rotational speed:

$x_9={\mathrm{\&omega;}}_{\mathrm{COI}}{|}_{t_{\mathrm{cl}}}-{\mathrm{\&omega;}}_{\mathrm{COI}}{|}_{t_0}$

In the following formula ${\mathrm{\&omega;}}_{\mathrm{COI}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i{\mathrm{\&omega;}}_i/\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i$

10) the absolute value sum of failure removal moment generator amature angular speed:

$x_{10}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{|{\mathrm{\&omega;}}_i-{\mathrm{\&omega;}}_{\mathrm{COI}}|}_{t_{\mathrm{cl}}}$

11) maximum of failure removal moment generator amature angular speed absolute value:

$x_{11}=\max \left(\right|{\mathrm{\&omega;}}_i-{\mathrm{\&omega;}}_{\mathrm{COI}}\left|\right){|}_{t_{\mathrm{cl}}}$

12) all generator amatures are in the most large and small value of fault generation initial time relative acceleration:

$x_{12-1}=(a_{\max }-a_{\mathrm{COI}}){|}_{t_0}$

$x_{12-2}=(a_{\min }-a_{\mathrm{COI}}){|}_{t_0}$

Wherein be a _COICenter of inertia acceleration

a _Max, a _MinIt is respectively the most large and small value of rotor angle acceleration;

13) all generator amatures are in the average of fault generation initial time acceleration:

$x_{13}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{a}_i{|}_{t_0}$

14) all generator amatures are in the root-mean-square error of fault generation initial time acceleration:

$x_{14}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{(a_i-x_{13})}^2}{|}_{t_0}$

15) all generator amatures are at the mean value of fault generation initial time relative acceleration:

$x_{15}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}(a_i-a_{\mathrm{COI}}){|}_{t_0}$

16) all generator amatures are in the variance of fault generation initial time relative acceleration:

$x_{16}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{((a_i-a_{\mathrm{COI}})-x_{15})}^2}{|}_{t_0}$

17) poor at instant of failure and fault generation initial time acceleration of machine and the machine of bringing up the rear before the neck:

$x_{17-1}=\max \left(a_i\right)-\min \left(a_i\right){|}_{t_0}$

$x_{17-2}=\max \left(a_i\right)-\min \left(a_i\right){|}_{t_{\mathrm{cl}}}$

18) maximum of all rotor kinetic energy of the failure removal moment:

$x_{18}=\max \left(V_{\mathrm{ki}}\right){|}_{t_{\mathrm{cl}}}i=\mathrm{1,2},...,N_G$

19) failure removal has the rotor kinetic energy of the generator p of maximum rotor angle constantly:

$x_{19}=\frac{1}{2}{M}_p\&times;({{\mathrm{\&omega;}}_p}^2-1){|}_{t_{\mathrm{cl}}}$

20) total " the energy adjustment " of electric power system:

$x_{20}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}(P_{\mathrm{ai}}\&times;{\mathrm{\&delta;}}_{\mathrm{di}})$

Wherein, P _AiBe the accelerating power of i platform generator at fault generation initial time; Relative rotor angle when being the excision of i platform generator failure;

21) mean value of all rotor kinetic energy of the failure removal moment:

$x_{21}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\frac{1}{2}{M}_i\&times;({{\mathrm{\&omega;}}_i}^2-1){|}_{t_{\mathrm{cl}}}$

22) mean value of all generator mechanical input powers before the fault:

$x_{22}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{p}_{\mathrm{mi}}$

23) variance of all generator mechanical input powers before the fault:

$x_{23}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{(P_{\mathrm{mi}}-x_{22})}^2}$

24) generator is in the average of fault generation initial time accelerating power:

$x_{24}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}){|}_{t_{\mathrm{cl}}}$

25) generator is in the variance of the relative accelerating power of fault generation initial time:

$x_{25}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{(\mathrm{\&Delta;}{P}_i-x_{24})}^2}{|}_{t_0}$

26) generator is in the average of the relative accelerating power of fault generation initial time:

$x_{26}=\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&Delta;P}}_i}{P_{\mathrm{mi}}}{|}_{t_0}$

27) generator is in the variance of the relative accelerating power of fault generation initial time:

$x_{27}=\sqrt{\frac{1}{N_G}\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{(\frac{{\mathrm{\&Delta;P}}_i}{P_{\mathrm{mi}}}-x_{26})}^2}{|}_{t_0}$

28) total kinetic energy of failure removal moment generator:

$x_{28}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\left|V_{\mathrm{KEi}}\right|{|}_{t_{\text{cl}}}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\left|\frac{1}{2}{M}_i{{\mathrm{\&omega;}}_i}^2\right|{|}_{t_{\mathrm{cl}}}$

29) maximum kinetic energy of failure removal moment generator is poor:

$x_{29-1}=(\max \left(V_{\mathrm{ki}}\right)-\min \left(V_{\mathrm{ki}}\right)){|}_{t_{\mathrm{cl}}}$

$x_{29-2}=(\max \left({\stackrel{\&CenterDot;}{V}}_{\mathrm{ki}}\right)-\min \left({\stackrel{\&CenterDot;}{V}}_{\mathrm{ki}}\right)){|}_{t_{\mathrm{cl}}}$

30) the suffered maximum active power of fault generation initial time separate unit generator is impacted:

$x_{30}=\max (P_{i1}-P_{i0}){|}_{t_0}i=\mathrm{1,2},...,N_G$

P in the formula _I1, P _I0Be respectively the meritorious output of generator i when instant of failure and before the fault;

31) the suffered minimum active power of fault generation initial time separate unit generator is impacted:

$x_{31}=\min (P_{i1}-P_{i0}){|}_{t_0}i=\mathrm{1,2},...,N_G$

32) the suffered maximum normalization active power of fault generation initial time separate unit generator is impacted:

$x_{32}=\underset{i}{\max }\frac{P_{i1}-P_{i0}}{P_{\mathrm{ie}}}{|}_{t_0}i=\mathrm{1,2},...,N_G$

P in the formula _IeFor the specified of generator i gained merit.

33) the suffered minimum normalization active power of fault generation initial time separate unit generator is impacted:

$x_{33}=\min \frac{P_{i1}-P_{i0}}{P_{\mathrm{ie}}}{|}_{t_0}i=\mathrm{1,2},...,N_G$

34) characterize the performance characteristic variable of fault to the system shock size:

$x_{34}=\left(\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{M}_i\right|P_{\mathrm{di}}\left|\right)/\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{M}_i$

In the following formula

Be the deceleration power of every generator with respect to the system inertia center; P _CiIt is fault excision moment output of a generator; P _0iIt is output of a generator before the fault;

35) generator accelerating power sum in a flash before accident is removed:

$x_{35}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{P}_{\mathrm{acc},i}^{\mathrm{BF}}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}(P_{\mathrm{mi}}^{\mathrm{BF}}-P_{\mathrm{ei}}^{\mathrm{BF}})=P_{\mathrm{COI}}^{\mathrm{BF}}$

For accident remove after the accelerating power of generator i in a flash;

For accident remove after the accelerating power in the center of inertia in a flash;

The accident that is respectively is removed the mechanical output and the electromagnetic power of back moment generator;

36) maximum of generator accelerating power absolute value in a flash after accident is removed:

$x_{36}=\underset{i}{\max }\left\{\right|P_{\mathrm{acc},i}^{\mathrm{PF}}\left|\right\}=\underset{i}{\max }\left\{\right|P_{\mathrm{mi}}^{\mathrm{PF}}-P_{\mathrm{ei}}^{\mathrm{PF}}\left|\right\}$

37) accelerating power and the generator inertia time constant ratio sum of generator in a flash before accident is removed:

$x_{37}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{P}_{\mathrm{acc},i}/M_i=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\left(\frac{P_{\mathrm{mi}}^{\mathrm{BF}}-P_{\mathrm{ei}}^{\mathrm{BF}}}{M_i}\right)$

38) accelerating power of generator and the maximum of generator inertia time constant ratio in a flash before accident is removed:

$x_{38}=\underset{i}{\max }\left\{\right|P_{\mathrm{acc},i}^{\mathrm{BF}}/M_i\left|\right\}=\underset{i}{\max }\left\{\right|\frac{P_{\mathrm{mi}}^{\mathrm{BF}}-P_{\mathrm{ei}}^{\mathrm{BF}}}{M_i}\left|\right\}$

39) accelerating power and the generator inertia time constant ratio sum of generator in a flash after accident is removed:

$x_{39}=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}{P}_{\mathrm{acc},i}^{\mathrm{PF}}/M_i=\underset{i=1}{\overset{N_G}{\mathrm{\&Sigma;}}}\left(\frac{P_{\mathrm{mi}}^{\mathrm{PF}}-P_{\mathrm{ei}}^{\mathrm{PF}}}{M_i}\right)$

40) accelerating power of generator and the maximum of generator inertia time constant ratio in a flash after accident is removed:

$x_{40}=\underset{i}{\max }\left\{\right|P_{\mathrm{acc},i}^{\mathrm{PF}}/M_i\left|\right\}=\underset{i}{\max }\left\{\right|\frac{P_{\mathrm{mi}}^{\mathrm{PF}}-P_{\mathrm{ei}}^{\mathrm{PF}}}{M_i}\left|\right\}.$

3. the method for claim 1 is characterized in that wherein said mass data preprocessing process may further comprise the steps:

(1) data normalization is handled

$x=\frac{x_0-x_{0\max }}{x_{0\max }-x_{0\min }}$

X in the formula ₀, x is respectively the data sample set before and after the conversion; x _0maxBe the sample maximum; x _0minBe the sample minimum value;

(2) sample stipulations

1) according to mixing the F statistical value, determines best cluster numbers C;

2) calculate the information gain of each characteristic variable, as the feature weight W in the cluster analysis to stabilization result _iNormalized sample data is carried out cluster analysis, with each sample be divided into its apart from the corresponding cluster subclass of the cluster centre of minimum, form C data sample subclass;

(3) dimension stipulations

Sample data is carried out logarithm centralization conversion and the vectorial centralization of row,,, drop to p (p＜m) by initial m with the dimension (being the characteristic variable number) of sample set by non-linear Karhunen-Loeve transformation;

1) to initial data X=(x _Ij), i=1 ..., n, j=1 ..., m (n is a sample number, and m is the input feature vector variable number) carries out logarithmic transformation, and the data set after the conversion is designated as W={w _Ij, have

w _ij＝ln?x _ij

2) go vectorial centralization, order

$y_{\mathrm{ij}}=w_{\mathrm{ij}}-\frac{1}{m}\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}$

3) compute matrix Y=(y _Ij) _{N * m}Covariance matrix and characteristic root and the standard orthogonal characteristic vector corresponding with characteristic value

4) characteristic value is carried out descending sort, λ is arranged ₁＞λ ₂＞...＞λ _m

5) according to contribution rate of accumulative total

(the individual characteristic value of p＜m), general, the contribution rate of accumulative total threshold values gets 0.85～0.95 to select the preceding p of contribution rate of accumulative total＞contribution rate of accumulative total threshold values

6) Karhunen-Loeve transformation forms new data sample set

F _i＝Γ _iWi＝1，2，...，p

Γ _iBe and i eigenvalue of above-mentioned covariance matrix _iCorresponding standard orthogonal characteristic vector.

4. the method for claim 1 is characterized in that also comprising the following steps:

(1) structure POWER SYSTEM EMERGENCY CONTROL target function:

Max (f (stability margin, MVA _G, MVA _L)) (2)

When system is unstable:

F=-abs (stability margin) (3)

During system stability:

MVA _Gi＝P _gi×r _gi?i＝1，...，N _G

MVA _Lj＝P _lj×r _lj?j＝1，...，N _L

Weight is the weight index for balance and stability nargin value and controlled quentity controlled variable value in the following formula, and general value is 100; MVA _GBe the cutter amount of generator, MVA _LBe the cutting load amount, P _GiCut the meritorious of generator and exerted oneself P _LjBy the burden with power of cutting load.r _GiBe that generator cuts out force rate example, r _LjBe the cutting load ratio, satisfy:

r _Gi=0 or 1 (5)

0≤r _lj≤1 (6)

(2) define generator and load to be excised, comprise the N of faulty line sending end side _GPlatform generator and be subjected to distolateral N _LIndividual load:

(3) initialization N input particle

Each particle is (a N _G+ N _L) n dimensional vector n, preceding N _GDimension data is represented N _GThe state r of platform generator to be cut _g, in the emergency control of reality, generator has only cuts and does not cut fully two states, this N fully _GThe value right and wrong 0 of dimension data are 1; Back N _LDimension data is represented N _LThe individual excision ratio r that waits to excise load _l, the value of each dimension is a real number between [0,1]; The initial position of each particle is set in allowed limits, at random

And speed

(4) corresponding a kind of control measure of particle are utilized power system analysis software (as PSASP) to carry out transient stability and are calculated, and after judgement added the control measure of this particle correspondence, whether system was stable; According to declaring steady result, calculate the target function of each particle correspondence by formula (3) or (4);

(5) utilize mangcorn subgroup method to search optimum controlling schemes.To the k time iteration, upgrade each particle's velocity according to formula (7), to each particle the 1st～N _GThe generator dimension is upgraded corresponding dimension position according to formula (8); To each particle (N _G+ 1)～(N _G+ N _L) the load dimension, upgrade corresponding dimension position value according to formula (9)

$v_i^{(k+1)}=w\&CenterDot;v_i^{\left(k\right)}+a_1\&times;\mathrm{rand}\left(\mathrm{0,1}\right)\&times;({\mathrm{pbest}}_i^{\left(k\right)}-x_i^{\left(k\right)})+a_2\&times;\mathrm{Rand}\left(\mathrm{0,1}\right)\&times;({\mathrm{gbest}}^{\left(k\right)}-x_i^{\left(k\right)})---\left(7\right)$

$\mathrm{If}(\mathrm{rand}\left(\right)<S\left(v_i^{(k+1)}\right))\mathrm{then}{x}_i^{(k+1)}=1$

$\mathrm{else}{x}_i^{(k+1)}=0---\left(8\right)$

$x_i^{(k+1)}=x_i^{\left(k\right)}+v_i^{(k+1)}---\left(9\right)$

W is an inertia weight in the formula, a ₁And a ₂Be two study factors, rand (0,1) and Rand (0,1) they are two random numbers that are evenly distributed between (0,1), i=1, and 2 ..., m.The w definition mode is as follows:

$w=w_{\max }-\frac{w_{\max }-w_{\min }}{{\mathrm{iter}}_{\max }}\&times;\mathrm{iter}---\left(10\right)$

W wherein _Max=0.9, w _Min=0.4, iter is current iterations, and greatest iteration number is made as 100.

Images