CN103795071A - Wide-area damping control method based on identification of power system model - Google Patents

Abstract

The invention relates to a wide-area damping control method based on identification of a power system model. The wide-area damping control method based on identification of the power system model includes the following steps that an observing signal is collected by means of a wide-area measuring system; pre-processing is carried out on the obtained observing signal to obtain a stable time sequence and the Bayes criterion is used for carrying out order determination on an ARMA model; an identification result is obtained by means of an improved ARMA model recognition algorithm; the time-lag power system state space model is obtained and a linear matrix inequality gamma and a linear matrix inequality theta are set up for obtaining the maximum time lag upper limit h which can be tolerated by the system and a controller K; fifthly, an output signal of the controller K is attached to a generator set excitation side which is in strong correlation with a low-frequency oscillation mode to finish damping control. According to the wide-area damping control method based on identification of the power system model, the computational process is simple, the operation speed is high, the obtained maximum time lag upper limit can meet the practical requirements and engineering realization is easy to achieve.

Description

A kind of wide area damp control method based on electric power system model identification

Technical field

The invention belongs to electrical power system wide-area time-delay damping control technology field, particularly a kind of wide area damp control method based on electric power system model identification.

Background technology

In recent years, along with the continuous expansion of China's electrical network scale, the safety and stability problem of electrical network is also increasingly sophisticated, particularly the large capacity power transmission of interregional long transmission line easily causes the area oscillation problem of system, to network-wide security, stable operation has caused great threat, and has limited to a certain extent the conveying capacity of interregional circuit.And therefore the effectively head it off of the control method of electric power system based on local signal that tradition adopts is necessary to seek new technological means.The recent WAMS (Wide Area Measurement System-WAMS) in electric power system extensive use, there is the technical characterstics such as high-precise synchronization measurement, high-speed communication and fast reaction, give the safe operation of electrical network and control the developing direction that provides new.

As tackling one of measure of interregional low-frequency oscillation, electrical power system wide-area control has obtained experts and scholars' extensive concern, but also exist many technical barriers simultaneously, for example: how to solve Time Delay in wide area signal, how to tackle complicated high order system be difficult to Accurate Model and controller dimension calamity problem all needs to further investigate, this is also several problems that this patent is focused on solving.

Summary of the invention

The object of the invention is to propose a kind of wide area damp control method based on electric power system model identification, be intended to evade system is carried out to the cumbersome procedure that detailed modeling is made, and improve System Identification Accuracy by a kind of improved arma modeling identification algorithm, then adopt a kind of improvement right of freedom matrix that does not comprise nonlinear terms based on this identification model, design is a kind of to the insensitive multi-input multi-output controller of time lag in certain limit, realize the wide area control of interconnected electric power system, to suppress low-frequency oscillation between system realm.

Technical scheme of the present invention is a kind of wide area damp control method based on electric power system model identification, comprises the following steps:

The first step, collects observation signal by WAMS;

Second step, carries out preliminary treatment by the observation signal obtaining, and obtains time series stably, and adopts bayesian criterion to carry out arma modeling to determine rank, comprises that definition bayesian criterion function is as follows,

${\mathrm{\&delta;}}_{\mathrm{BIC}}\left(P\right)=N1n{\mathrm{\&sigma;}}_a^2+\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="1401241418033.png"\; state="\{\{state\}\}">P</figure-callout>}1\mathrm{nN}$

Wherein, N is observation signal length,

for the variance of observation signal, P=n+m, the number that n is auto-regressive parameter, m is the number of moving average parameter; As bayesian criterion function δ _bICthe value exponent number that hour corresponding P is arma modeling;

The 3rd step, by following sub-process, obtains identification result,

Step 1, selects new breath length, makes time sequence number t=1, establishes initial value p ₀=10 ⁶, P (0)=p ₀i,

wherein I is unit matrix, 1 _nfor ranks element is 1 n rank matrix entirely;

Step 2, input observation data { x _t;

Step 3, establishes X (p, t) for the new corresponding accumulation output vector that ceases the data that length is p, V (p, t) is the new corresponding accumulation error vector that ceases the data that length is p, Φ (p, t) be the new corresponding information matrix that ceases the data that length is p

for piling up the estimated value of error V (p, t),

for the transposed matrix of the estimated value of information matrix Φ (p, t); With formula one and formula two construct X (p, t) and

formula one

$\widehat{\mathrm{\&Phi;}}(p,t)=[X(p,t-1),X(p,t-2),\&CenterDot;\&CenterDot;\&CenterDot;X(p,t-n),\hat{V}(p,t-1),\hat{V}(p,t-2),\&CenterDot;\&CenterDot;\&CenterDot;\hat{V}(p,t-m)]$ Formula two

Wherein,

for the dimension vector space that is p;

Step 4, establishes multi-input multi-output system and is decomposed into the subsystems of the single output of the many inputs of l, and arbitrary subsystem is designated as f, f=1, and 2 ..., l, calculates covariance matrix P (t) by formula three, refreshes the estimated value of parameter vector θ in the t moment by formula four

$P_f^{-1}\left(t\right)=P_f^{-1}(t-1)+{\widehat{\mathrm{\&Phi;}}}_f(p,t){\widehat{\mathrm{\&Phi;}}}_f^T(p,t)$ Formula three

${\widehat{\mathrm{\&theta;}}}_f\left(t\right)={\widehat{\mathrm{\&theta;}}}_{f-1}(t-1)+P_f\left(t\right){\widehat{\mathrm{\&Phi;}}}_f^T(p,t)[X_f(p,t)-{\widehat{\mathrm{\&Phi;}}}_f^T(p,t){\widehat{\mathrm{\&theta;}}}_{f-1}\left(t\right)]$ Formula four

Wherein, P _f(t) be the covariance matrix of subsystem f,

for the parameter vector θ of subsystem f is in the estimated value in t moment

x _f(p, t) is the corresponding accumulation output vector of the data that in subsystem f, new breath length is p,

for information matrix Φ in subsystem f _fthe transposed matrix of the estimated value of (p, t);

Step 5, calculates by formula five

return to step 2 iteration and carry out, until

while being less than corresponding predetermined threshold value or t > L, finish sub-process,

${\hat{V}}_f(p,t)=X_f(p,t)-{\widehat{\mathrm{\&Phi;}}}_f^T(p,t){\widehat{\mathrm{\&theta;}}}_f\left(t\right)$ Formula five

The 4th step, according to above-mentioned steps gained

with

with arma modeling X (p, t)=Φ ^t(p, t) θ+V (p, t), conversion obtains suc as formula the time-lag power system state-space model shown in six;

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=\mathrm{Ax}+\mathrm{BKu}(t-d\left(t\right))\\ y=\mathrm{Cx}\end{array}\right.$ Formula six

Wherein,

for the first derivative of x, d (t) in wide area signal time become time lag amount, K is state feedback controller, A, B, C is the sytem matrix that identification obtains, x is state variable, u (t) is control inputs amount, y is output variable;

If there is matrix U=U ^t> 0, Z=Z ^t> 0,

and the matrix M of any appropriate dimension ₁, M ₂and V, by set up suc as formula seven and formula eight shown in LMI Γ, Θ, and utilize Matlab generalized eigenvalue solver to solve, the patient maximum time lag upper limit h of the system that obtains and controller K;

$\mathrm{\&Gamma;}=\left[\begin{array}{ccc}J+J^T+h{Y}_{11}& A_{\mathrm{<figure-callout\; id="2"\; label="d\; M"\; filenames="1401241418034.png,1401241418035.png"\; state="\{\{state\}\}">d</figure-callout>}}{M}_{}{}_2+U-M_1^T+h{Y}_{12}& h{J}^T\\ *& -M_2-M_2^T+h{Y}_{22}& h{M}_2^T{A}_d^T\\ *& *& -\mathrm{hZ}\end{array}\right]<0$ Formula seven

$\mathrm{\&Theta;}=\left[\begin{array}{cc}{Y}_{11}& Y_{12}\\ *& Y_{22}-Z\end{array}\right]\&GreaterEqual;0$ Formula eight

In formula seven, matrix J=AU+A _dm ₁+ BV; A _d=BK is time lag state matrix;

The 5th step, is attached to controller K output signal with the generating unit excitation side of low frequency oscillation mode strong correlation and completes damping control.

The wide area damp control method for designing computational process that the present invention provides is simple, and the speed of service is very fast, and the maximum time lag upper limit obtaining can meet actual requirement, is easy to Project Realization.Compared with prior art, the present invention has following beneficial effect:

1, the present invention, by adopting identification technique, has avoided the cumbersome procedure to large complicated power network modeling, and compared with prior art, the present invention has improved the precision of identification by improved arma modeling identification algorithm, has reduced iterations.Obtain system transter by carrying out identification, avoided, in conventional method, the detailed modeling of system is carried out to linearisation model reduction again and the complex process that causes.By adopting many new breath thought to improve arma modeling identification algorithm, the precision of identification the iterations reducing are improved, for next step wide area damping control design is had laid a good foundation.

The method of the improvement right of freedom matrix that 2, the present invention adopts, design one the insensitive wide area damping control of time lag in certain limit, the method, containing nonlinear terms, has not been avoided the complicated iterative process causing because comprising nonlinear terms in existing right of freedom matrix method.By solving LMI (Linear Maxtrix Inequality, LMI), can obtain to time become the insensitive state feedback controller of time lag, can realize the effective inhibition to the low-frequency oscillation of interconnected electric power system region, improve the damping of electric power system.

Accompanying drawing explanation

Fig. 1 is the multi-input multi-output system submodel STRUCTURE DECOMPOSITION figure that the embodiment of the present invention adopts.

Fig. 2 is the EPRI-8 machine 36 node test system line charts that the embodiment of the present invention adopts.

Fig. 3 selects different new breath length the error comparison diagrams with traditional arma modeling discrimination method in the discrimination method that adopts of the present invention.

Fig. 4 is the merit angular difference dynamic response comparison diagram of lower No. 1 and No. 8 generator of fault 1 situation in the embodiment of the present invention.

Fig. 5 is the meritorious dynamic response figure of interconnection between fault 1 situation lower node 22 and node 20 in the embodiment of the present invention.

Fig. 6 is the merit angular difference dynamic response comparison diagram of lower No. 1 and No. 8 generator of fault 2 situations in the embodiment of the present invention

Fig. 7 is the meritorious dynamic response figure of interconnection between fault 2 situation lower nodes 22 and node 20 in the embodiment of the present invention.

Embodiment

Describe technical solution of the present invention in detail below in conjunction with drawings and Examples.

When technical solution of the present invention is specifically implemented, can adopt software engineering to realize automatic operational process by those skilled in the art.The flow process of embodiment comprises the following steps:

The first step: by WAMS systematic collection observation signal sequence, conventionally choose important interconnection meritorious or with low-frequency oscillation strong correlation generator group merit angular difference etc. as observation signal.

Second step: the observation signal obtaining is gone to the preliminary treatment such as trend, zero-mean, obtain time series stably, and adopt bayesian criterion (Bayesian Information Criterion, BIC) to carry out arma modeling to determine rank, definition bayesian criterion function is:

${\mathrm{\&delta;}}_{\mathrm{BIC}}\left(P\right)=N1n{\mathrm{\&sigma;}}_a^2+\mathrm{<figure-callout\; id="1"\; label="P"\; filenames="1401241418033.png"\; state="\{\{state\}\}">P</figure-callout>}1\mathrm{nN}$

Wherein, N is observation signal length, for the variance of observation signal; P=n+m is the exponent number of model.As bayesian criterion function δ _bICvalue hour corresponding P is the exponent number of arma modeling.

Arma modeling (Auto-Regressive and Moving Average Model) is the important method of search time sequence, is basis " mixing " formation by autoregression model (being called for short AR model) and moving average model (being called for short MA model).

The 3rd step: by following sub-process, obtain identification result.

Step 1. is selected new breath length, makes t=1, establishes initial value p ₀=10 ⁶, P (0)=p ₀i, here, I is unit matrix, 1 _nfor ranks element is 1 n rank matrix entirely.

Step 2. is collected observation data { x _t;

Formula for step 3. (7) and formula (15) structure X (p, t) and

Formula for step 4. (21) is calculated P (t), refreshes by formula (20)

Formula for step 5. (22) is calculated

return to step 2 iteration and carry out, until Identification Errors is enough little (

be less than predetermined threshold value) or the complete data window node of iteration L(right _t=Lexecute iteration, after t=t+1, be greater than t > L), stop sub-process, enter the 4th step.

If iterate to data window node, through above flow process, to { x _t, t=1,2 ..., L has carried out processing successively.

The 4th step: by above-mentioned steps (first step～three step), can obtain suc as formula the time-lag power system state-space model shown in (23).By setting up suc as formula the LMI shown in (28) and formula (29), and utilize Matlab generalized eigenvalue solver (GEVP) to solve, the patient maximum time lag upper limit h of the system that obtains and controller K.

The 5th step: controller K output signal is attached to the generating unit excitation side of low frequency oscillation mode strong correlation and completes damping control.

The wide area damp control method based on electric power system model identification is applied to EPRI-8 machine 36 node test systems by the present embodiment, as shown in Figure 2, wherein there is synchronous machine No. 1, No. 2, No. 3, No. 4, No. 5, No. 6, No. 7, No. 8, and bus 1-9,11-14,16,18-31,33,34,50-52, in EPRI-8 machine 36 node test systems, number discontinuous.Designed wide area control signal append to synchronous generator No. 1, No. 3, No. 5, No. 7 excitation system and on.And controller input signal, wide area measurement system signal is chosen respectively the meritorious P between bus 9 and 24 _9-24, the meritorious P between bus 19 and 30 _19-30, the merit angular difference δ between No. 1 synchronous machine and No. 3 synchronous machines _1-3, and merit angular difference δ between No. 1 synchronous machine and No. 7 synchronous machines _1-7.

Obtain needed wide area controller and maximum time lag upper limit 435.2ms through above-mentioned design cycle, can meet actual requirement of engineering.

For the identification effect of check this patent identification algorithm, in the improvement arma modeling algorithm that this patent is adopted, choose different new breath length p=3,5,8 iteration error and contrast with traditional arma modeling algorithm p=1, as shown in Figure 3, abscissa is t, unit is second (s), ordinate is iteration error delta, and unit is percentage %.

In order to check the final effect of this patent control method, test macro is carried out respectively to the l-G simulation test of two class faults, system responses is respectively as Fig. 4,5 and Fig. 6,7, and abscissa is t (s), and ordinate is respectively: merit angular difference, unit is degree deg; Meritorious, unit is the p.u. of standard unit.

Fault 1:2 generator terminal voltage raises 10%.

Fault 2: three phase short circuit fault, fault recovery after 100ms occur at node 16 places.

Below the specific implementation of committed step in embodiment is elaborated:

1, the 3rd step has been used improvement arma modeling identification algorithm

In the middle of electric power system actual moving process, due to the random variation at any time such as load switching, Unit Commitment, in system, each parameter exists the small size disturbance that is similar to noise, and such noise-like signal can be described as steadily, the time series { x of zero-mean _t, t=1,2 ..., L, t is time sequence number, L is data window node, its Auto Regressive Moving Average(ARMA) model can be expressed as:

$x_t-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}_i{x}_{t-i}=v_t-\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{d}_j{v}_{t-j}---\left(1\right)$

Wherein, a _i(i=1,2 ..., n) be auto-regressive parameter (Auto Regressive, AR), the number that n is auto-regressive parameter; d _j(j=1,2 ..., m) being moving average parameter (Moving Average, MA), m is the number of moving average parameter; { v _tbe residual sequence, formula (1) is denoted as ARMA (n, m) model.

Introduce the backward shift operator z of unit ^-1, formula (1) can be expressed as:

$(1-a_1{z}^{-1}-a_2{z}^{-2}-\&CenterDot;\&CenterDot;\&CenterDot;-a_n{z}^{-n})x_t=(1-d_1{z}^{-1}-d_2{z}^{-2}\&CenterDot;\&CenterDot;\&CenterDot;-d_n{z}^{-n})v_t---\left(2\right)$

Define parameter vector θ to be identified and information vector

as follows:

Formula (1) can be written as following form:

In above formula, x (t) is signal to be identified, and v (t) is noise signal.Adopt least square method or random gradient identification algorithm parameter vector θ to be estimated to have following form to formula (4):

$\widehat{\mathrm{\&theta;}}\left(t\right)=\widehat{\mathrm{\&theta;}}(t-1)+L\left(t\right)e\left(t\right)---\left(5\right)$

Wherein,

for parameter vector θ is in the estimated value in t moment,

for gain vector,

for scalar newly ceases, i.e. single new breath, wherein representation vector space,

for the dimension vector space that is n.And many new breath identification rules are to revise on identification idea basis at single new breath, scalar is newly ceased

be generalized to vector information

for the dimension vector space that is p, and gain vector

transform into gain matrix

wherein p is new breath length,

for the dimension vector space that is p × n,, formula (5) can be converted into following many new breath forms:

$\widehat{\mathrm{\&theta;}}\left(t\right)=\widehat{\mathrm{\&theta;}}(t-1)+\mathrm{\&Gamma;}(p,t)E(p,t)---\left(6\right)$

This patent is selected recurrence method identification method, and in conjunction with many new breath identification thought, obtain many new breath RLSs (Multi-Innovation Recursive Least Square Algorithm, MIRLS algorithm), its concrete discrimination method is as follows:

The data that are p for new breath length, define corresponding accumulation output vector X (t), pile up error vector V (t), information matrix Φ (t) and be:

Wherein X (p, t) be the new corresponding accumulation output vector that ceases the data that length is p, V (p, t) is the corresponding accumulation error vector of the data that newly breath length is p, Φ (p, t) is the new corresponding information matrix that ceases the data that length is p.

Arma modeling can be expressed as:

X(p,t)＝Φ ^T(p,t)θ+V(p,t) (10)

Owing to having comprised the accumulation error vector V (p, t) that can not survey noise in the middle of information matrix Φ (p, t), therefore need by its estimated value replace:

$\hat{V}(p,t)=X(p,t)-{\widehat{\mathrm{\&Phi;}}}^T(p,t)\widehat{\mathrm{\&theta;}}\left(t\right)---\left(11\right)$

In above formula,

for the transposed matrix of the estimated value of information matrix Φ (p, t),

for parameter vector θ is in the estimated value in t moment, now get criterion function J (θ):

$J\left(\mathrm{\&theta;}\right)=\underset{t=1}{\overset{L}{\mathrm{\&Sigma;}}}{\left|\right|X(p,t)-{\mathrm{\&Phi;}}^T(p,t)\mathrm{\&theta;}\left|\right|}^2---\left(12\right)$

In above formula, X (p, t) and Φ ^t(p, t) is respectively the t row element of the transposed matrix of piling up output matrix and information matrix.Try to achieve the estimated value of parameter vector θ by minimization J (θ) with information vector

get final product parameter vector θ in the estimated value in t moment

$\widehat{\mathrm{\&theta;}}\left(t\right)=\widehat{\mathrm{\&theta;}}(t-1)+P\left(t\right){\widehat{\mathrm{\&Phi;}}}^T(p,t)[X(p,t)-{\widehat{\mathrm{\&Phi;}}}^T(p,t)\widehat{\mathrm{\&theta;}}(t-1)]---\left(13\right)$

$P^{-1}\left(t\right)=P^{-1}(t-1)+\widehat{\mathrm{\&Phi;}}(p,t){\widehat{\mathrm{\&Phi;}}}^T(p,t)---\left(14\right)$

$\widehat{\mathrm{\&Phi;}}(p,t)=[X(p,t-1),X(p,t-2),\&CenterDot;\&CenterDot;\&CenterDot;X(p,t-n),\hat{V}(p,t-1),\hat{V}(p,t-2),\&CenterDot;\&CenterDot;\&CenterDot;\hat{V}(p,t-m)]---\left(15\right)$

In above formula, P (t) represents covariance matrix.Then above-mentioned algorithm is generalized to polynary form, consideration dimension is l, the polynary identification model that data window length is r, shown in (16):

The subsystem that above multi-input multi-output system can be decomposed into l many single outputs of input (Multi-Input Single Output, MISO) conventionally carries out identification, as shown in Figure 1: l input v ₁(t) ... v _l(t) obtain l output x through ssystem transfer function ₁(t) ... x _l(t), can be decomposed into l input v ₁(t) ... v _l(t) obtain exporting x through subsystem 1 ₁(t) ... l input v ₁(t) ... v _l(t) obtain exporting x through subsystem l _l(t).

If wherein arbitrary subsystem f is input as v _f(t), be output as x _f(t), f=1,2 ..., l, formula (4) can be converted into:

Wherein, the information vector of subsystem f

with parameter vector θ _f:

θ _f＝[a _f1,a _f2,...a _fn；d _f11,d _f12,...d _f1m；d _f21,d _f22,...d _f2m；...d _fl1,d _fl2,...d _flm] ^T (19)

Wherein, a _fi(i=1,2 ..., n) be the auto-regressive parameter of subsystem f, the number that n is auto-regressive parameter; d _fqj(q=1,2 ..., l, j=1,2 ..., m) be the moving average parameter of subsystem f.

In above formula, f ∈ (1, l) represent f subsystem, derive according to formula (13)～formula (15), have:

${\widehat{\mathrm{\&theta;}}}_f\left(t\right)={\widehat{\mathrm{\&theta;}}}_{f-1}(t-1)+P_f\left(t\right){\widehat{\mathrm{\&Phi;}}}_f^T(p,t)[X_f(p,t)-{\widehat{\mathrm{\&Phi;}}}_f^T(p,t){\widehat{\mathrm{\&theta;}}}_{f-1}\left(t\right)]---\left(20\right)$

$P_f^{-1}\left(t\right)=P_f^{-1}(t-1)+{\widehat{\mathrm{\&Phi;}}}_f(p,t){\widehat{\mathrm{\&Phi;}}}_f^T(p,t)---\left(21\right)$

${\hat{V}}_f(p,t)=X_f(p,t)-{\widehat{\mathrm{\&Phi;}}}_f^T(p,t){\widehat{\mathrm{\&theta;}}}_f\left(t\right)---\left(22\right)$

Wherein, P _f(t) be the covariance matrix of subsystem f,

for the parameter vector θ of subsystem f is in the estimated value in t moment

for the corresponding accumulation output vector of the data that in subsystem f, new breath length is p,

for piling up error V in subsystem f _fthe estimated value of (p, t), for information matrix Φ in subsystem f _fthe transposed matrix of the estimated value of (p, t).

Therefore, this patent, by the circulation of following steps, can obtain identification result.

Step 1. is selected new breath length, makes t=1, establishes initial value p ₀=10 ⁶, P (0)=p ₀i,

here, I is unit matrix, 1 _nfor ranks element is 1 n rank matrix entirely.

Step 2. is inputted observation data { x _t;

Formula for step 3. (7) and formula (15) structure X (p, t) and

Step 4. is established multi-input multi-output system and is decomposed into the subsystems of the single output of the many inputs of l, arbitrary subsystem is designated as f, f=1, and 2 ..., l, calculates each subsystem by formula (21) obtain P (t), calculate each subsystem by formula (20)

refresh

Formula for step 5. (22) is calculated

comprise and use respectively formula (22) to calculate to each subsystem

return to step 2 iteration and carry out, until

enough little (is each subsystem

all be less than corresponding predetermined threshold value) or t > L, step 2 no longer returned to.When concrete enforcement, can be designed to, in step 5, first use formula (22) to calculate

judging whether to meet termination condition, be to finish this sub-process, otherwise t=t+1 returns to step 2 again.

2, the 4th step is based on improving right of freedom matrix design wide area damping control

First,, by above-mentioned discrimination method, obtain

with

now can obtain suc as formula the system arma modeling shown in (10), by model transferring, can obtain suc as formula the state space form shown in (23):

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=\mathrm{Ax}+\mathrm{BKu}(t-d\left(t\right))\\ y=\mathrm{Cx}\end{array}\right.---\left(23\right)$

This model is the electric power system universal model containing time lag.Wherein,

for the first derivative of x, d (t) in wide area signal time become time lag amount, K is state feedback controller, A, B, C is the sytem matrix that identification obtains, x is state variable, u (t) is control inputs amount, y is output variable.

Then, pass through Newton Leibniz formula

$x(t-d\left(t\right))=x\left(t\right)-{\&Integral;}_{t-d\left(t\right)}^t\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}---\left(24\right)$

Wherein, s is integral sign.

Can obtain:

$2[x^T\left(t\right){\mathrm{<figure-callout\; id="1"\; label="N"\; filenames="1401241418033.png"\; state="\{\{state\}\}">N</figure-callout>}}_1+x^T(t-d\left(t\right))N_2]\&times;[x\left(t\right)-{\&Integral;}_{t-d\left(t\right)}^t\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}-x(t-d\left(t\right))]=0---\left(25\right)$

Wherein, N ₁and N ₂for the matrix of any appropriate dimension, be right of freedom matrix, it has overcome the fixing conservative that weight matrix brings of use in conventional method, and its optimal value can be passed through LMI (Linear Maxtrix Inequality, LMI) solution obtains, and in formula (25)

can replace by through type (24).

, if there is P=P in given scalar h > 0 ^t> 0, Q=Q ^t> 0,

and the matrix N of any appropriate dimension ₁, N ₂, following LMI Φ, Ψ are set up:

$\mathrm{\&Phi;}=\left[\begin{array}{ccc}{\mathrm{\&Phi;}}_{11}& {\mathrm{\&Phi;}}_{12}& h{A}^T\\ *& {\mathrm{\&Phi;}}_{22}& h{A}_d^T\\ *& *& -h{Q}^{-1}\end{array}\right]<0---\left(26\right)$

$\mathrm{\&Psi;}=\left[\begin{array}{ccc}{X}_{11}& X_{12}& {\mathrm{<figure-callout\; id="1"\; label="N"\; filenames="1401241418033.png"\; state="\{\{state\}\}">N</figure-callout>}}_1\\ *& X_{22}& {\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="1401241418034.png,1401241418035.png"\; state="\{\{state\}\}">N</figure-callout>}}_2\\ *& *& Q\end{array}\right]\&GreaterEqual;0---\left(27\right)$

Wherein, P, Q, X are the process matrix in LMI computational process, and * represents upper transposed matrix corresponding to diagonal matrix in matrix.If be limited to h in maximum time lag, the formula (23) that meets time lag constraint 0≤d (t)≤h is progressive stable, wherein

$\begin{array}{c}{\mathrm{\&Phi;}}_{11}=\mathrm{PA}+A^{T}P+{\mathrm{<figure-callout\; id="1"\; label="N"\; filenames="1401241418033.png"\; state="\{\{state\}\}">N</figure-callout>}}_1+N_1^T+\mathrm{<figure-callout\; id="11"\; label="h\; X"\; filenames="1401241418032.png"\; state="\{\{state\}\}">h</figure-callout>}{X}_{}{}_{11}\\ {\mathrm{\&Phi;}}_{12}=P{A}_d-{\mathrm{<figure-callout\; id="1"\; label="N"\; filenames="1401241418033.png"\; state="\{\{state\}\}">N</figure-callout>}}_1+N_2^T+h{X}_{12}\\ {\mathrm{\&Phi;}}_{22}=-N_2-N_2^T+\mathrm{<figure-callout\; id="22"\; label="h\; X"\; filenames="1401241418032.png"\; state="\{\{state\}\}">h</figure-callout>}{X}_{}{}_{22}\end{array}$

Inference: given scalar h > 0, if there is matrix U=U ^t> 0, Z=Z ^t> 0,

and the matrix M of any appropriate dimension ₁, M ₂and V, following LMI Γ, Θ are set up:

$\mathrm{\&Gamma;}=\left[\begin{array}{ccc}J+J^T+h{Y}_{11}& {\mathrm{<figure-callout\; id="2"\; label="A\; d\; M"\; filenames="1401241418034.png,1401241418035.png"\; state="\{\{state\}\}">A</figure-callout>}}_d{M}_{}{}_2+U-M_1^T+h{Y}_{12}& h{J}^T\\ *& -M_2-M_2^T+h{Y}_{22}& h{M}_2^T{A}_d^T\\ *& *& -\mathrm{hZ}\end{array}\right]<0---\left(28\right)$

$\mathrm{\&Theta;}=\left[\begin{array}{cc}{Y}_{11}& Y_{12}\\ *& Y_{22}-Z\end{array}\right]\&GreaterEqual;0---\left(29\right)$

In formula (28), matrix J=AU+A _dm ₁+ BV; A _d=BK is time lag state matrix.

Existence feedback controller K, the formula (23) that makes to meet time lag constraint 0≤d (t)≤h is progressive stable, and K=VL ^-1, this state feedback controller K is wide area damping control to be designed.

Prove:

If matrix

$H=\left[\begin{array}{cc}P& 0\\ N_1^T& N_2^T\end{array}\right],A_K=\left[\begin{array}{cc}A& \mathrm{BK}\end{array}\right],I_2=\left[\begin{array}{cc}I& -I\end{array}\right],\stackrel{-}{A}=\left[\begin{array}{c}{A}_{\mathrm{<figure-callout\; id="2"\; label="K\; I"\; filenames="1401241418034.png,1401241418035.png"\; state="\{\{state\}\}">K</figure-callout>}}& \\ I_{}{}_2\end{array}\right]$

By A _kbring formula (26) into,

$\mathrm{\&Xi;}=\left[\begin{array}{cc}{H}^T\stackrel{-}{A}+{\stackrel{-}{A}}^T\stackrel{}{H+\mathrm{hX}}& h{A}_K^T\\ h{A}_K& -h{Z}^{-1}\end{array}\right]<0$

If

$H^{-1}={\left[\begin{array}{cc}P& 0\\ N_1^T& N_2^T\end{array}\right]}^{-1}=\left[\begin{array}{cc}L& 0\\ M_1& {\mathrm{<figure-callout\; id="2"\; label="M"\; filenames="1401241418034.png,1401241418035.png"\; state="\{\{state\}\}">M</figure-callout>}}_2\end{array}\right]$

Ξ represents LMI, uses respectively diag{H ^-T, I} and diag{H ^-1, I} premultiplication and right multiplier (28), wherein diag{} represents diagonal matrix, and makes following variable and replace, and makes variable Y, R, V be respectively:

Y＝H ^-TXH ^-1,R＝Z ^-1,V＝KL

Use respectively diag{H ^-T, I} and diag{H ^-1, I} premultiplication and right multiplier (27), and utilize taper complement fixed reason to get final product to obtain formula (29), card is finished.

Taper (Schur) complement fixed reason: to given symmetrical matrix

wherein S ₁₁∈ R ^{r × r}, R ^{r × r}represent that ranks are the square formation of r.Three conditions are of equal value below:

(1)S＜0；

(2) ${\mathrm{<figure-callout\; id="11"\; label="S"\; filenames="1401241418032.png"\; state="\{\{state\}\}">S</figure-callout>}}_{11}<0,S_{22}-S_{12}^T{S}_{11}^{-1}{S}_{12}<0;$

(3) ${\mathrm{<figure-callout\; id="22"\; label="S"\; filenames="1401241418032.png"\; state="\{\{state\}\}">S</figure-callout>}}_{22}<0,S_{11}-S_{12}{S}_{22}^{-1}{S}_{12}^T<0$

Owing to not comprising any nonlinear terms in formula (28), (29), therefore do not need by transforming, the non-linear minimization problem that can be converted into based on LMI is carried out iterative, do not need to carry out parameter adjustment yet, can sum up and become a LMI generalized eigenvalue problem (Generalized Eeigenvalue Problem, GEVP), finally solve controlled device K.

Specific embodiment described herein is only to the explanation for example of the present invention's spirit.Those skilled in the art can make various modifications or supplement or adopt similar mode to substitute described specific embodiment, but can't depart from spirit of the present invention or surmount the defined scope of appended claims.

Claims (1)

1. the wide area damp control method based on electric power system model identification, is characterized in that, comprises the following steps:

The first step, collects observation signal by WAMS;

Second step, carries out preliminary treatment by the observation signal obtaining, and obtains time series stably, and adopts bayesian criterion to carry out arma modeling to determine rank, comprises that definition bayesian criterion function is as follows,

${\mathrm{\&delta;}}_{\mathrm{BIC}}\left(P\right)=N1n{\mathrm{\&sigma;}}_a^2+P1\mathrm{nN}$

Wherein, N is observation signal length,

for the variance of observation signal, P=n+m, the number that n is auto-regressive parameter, m is the number of moving average parameter; As bayesian criterion function δ _bICthe value exponent number that hour corresponding P is arma modeling;

The 3rd step, by following sub-process, obtains identification result,

Step 1, selects new breath length, makes time sequence number t=1, establishes initial value p ₀=10 ⁶, P (0)=p ₀i,

wherein I is unit matrix, 1 _nfor ranks element is 1 n rank matrix entirely;

Step 2, input observation data { x _t;

Step 3, establishes X (p, t) for the new corresponding accumulation output vector that ceases the data that length is p, V (p, t) is the new corresponding accumulation error vector that ceases the data that length is p, Φ (p, t) be the new corresponding information matrix that ceases the data that length is p for piling up the estimated value of error V (p, t),

for the transposed matrix of the estimated value of information matrix Φ (p, t); With formula one and formula two construct X (p, t) and

formula one

$\widehat{\mathrm{\&Phi;}}(p,t)=[X(p,t-1),X(p,t-2),\&CenterDot;\&CenterDot;\&CenterDot;X(p,t-n),\hat{V}(p,t-1),\hat{V}(p,t-2),\&CenterDot;\&CenterDot;\&CenterDot;\hat{V}(p,t-m)]$ Formula two

Wherein,

for the dimension vector space that is p;

Step 4, establishes multi-input multi-output system and is decomposed into the subsystems of the single output of the many inputs of l, and arbitrary subsystem is designated as f, f=1, and 2 ..., l, calculates covariance matrix P (t) by formula three, refreshes the estimated value of parameter vector θ in the t moment by formula four

$P_f^{-1}\left(t\right)=P_f^{-1}(t-1)+{\widehat{\mathrm{\&Phi;}}}_f(p,t){\widehat{\mathrm{\&Phi;}}}_f^T(p,t)$ Formula three

${\widehat{\mathrm{\&theta;}}}_f\left(t\right)={\widehat{\mathrm{\&theta;}}}_{f-1}(t-1)+P_f\left(t\right){\widehat{\mathrm{\&Phi;}}}_f^T(p,t)[X_f(p,t)-{\widehat{\mathrm{\&Phi;}}}_f^T(p,t){\widehat{\mathrm{\&theta;}}}_{f-1}\left(t\right)]$ Formula four

Wherein, P _f(t) be the covariance matrix of subsystem f,

for the parameter vector θ of subsystem f is in the estimated value in t moment

x _f(p, t) is the corresponding accumulation output vector of the data that in subsystem f, new breath length is p,

for information matrix Φ in subsystem f _fthe transposed matrix of the estimated value of (p, t);

Step 5, calculates by formula five

return to step 2 iteration and carry out, until

while being less than corresponding predetermined threshold value or t > L, finish sub-process,

${\hat{V}}_f(p,t)=X_f(p,t)-{\widehat{\mathrm{\&Phi;}}}_f^T(p,t){\widehat{\mathrm{\&theta;}}}_f\left(t\right)$ Formula five

The 4th step, according to above-mentioned steps gained

with with arma modeling X (p, t)=Φ ^t(p, t) θ+V (p, t), conversion obtains suc as formula the time-lag power system state-space model shown in six;

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=\mathrm{Ax}+\mathrm{BKu}(t-d\left(t\right))\\ y=\mathrm{Cx}\end{array}\right.$ Formula six

Wherein,

for the first derivative of x, d (t) in wide area signal time become time lag amount, K is state feedback controller, A, B, C is the sytem matrix that identification obtains, x is state variable, u (t) is control inputs amount, y is output variable;

If there is matrix U=U ^t> 0, Z=Z ^t> 0,

and the matrix M of any appropriate dimension ₁, M ₂and V, by set up suc as formula seven and formula eight shown in LMI Γ, Θ, and utilize Matlab generalized eigenvalue solver to solve, the patient maximum time lag upper limit h of the system that obtains and controller K;

$\mathrm{\&Gamma;}=\left[\begin{array}{ccc}J+J^T+h{Y}_{11}& A_d{M}_2+U-M_1^T+h{Y}_{12}& h{J}^T\\ *& -M_2-M_2^T+h{Y}_{22}& h{M}_2^T{A}_d^T\\ *& *& -\mathrm{hZ}\end{array}\right]<0$ Formula seven

$\mathrm{\&Theta;}=\left[\begin{array}{cc}{Y}_{11}& Y_{12}\\ *& Y_{22}-Z\end{array}\right]\&GreaterEqual;0$ Formula eight

In formula seven, matrix J=AU+A _dm ₁+ BV; A _d=BK is time lag state matrix;

The 5th step, is attached to controller K output signal with the generating unit excitation side of low frequency oscillation mode strong correlation and completes damping control.

Images