JP2010011620A - Apparatus, method and program for creating abridged model of electric power system - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To efficiently create a high precision abridged model of a large scale electric power system. <P>SOLUTION: A time series analysis section 142 estimates main oscillation modes before and after abridged by using an ARMA model, an eigen value analysis section 143 calculates eigen values respectively corresponding to the main oscillation modes before and after abridged by using an inverse iteration method, a mode matching section 151 specifies a pair of oscillation modes before and after abridged by using the inner product of an eigenvector, and a parameter adjusting section 152 adjusts the inertia constant, AVR constant and PSS constant in such a manner that the main oscillation modes before and after abridged match each other. <P>COPYRIGHT: (C)2010,JPO&INPIT

Description

  The present invention relates to a power system contracted model creating device, a power system contracted model creating method, and a power system contracted model creating program for creating a contracted model by contracting parts other than the main system from a power system model, In particular, the present invention relates to a power system contracted model creating apparatus, a power system contracted model creating method, and a power system contracted model creating program capable of efficiently contracting a large-scale power system with high accuracy.

  In the stability analysis of the power system, it is necessary to consider the characteristics of the control system of each generator. It takes a long time to calculate. Therefore, analysis other than the main system is performed using the reduced system.

  As a power system contraction method, a physical contraction method (see, for example, Patent Document 1, Patent Document 2, and Non-Patent Document 1) that contracts while leaving a physical image of a generator, a transmission line, and the like is used. There is a mathematical reduction method that formulates the characteristics of the target system and reduces the eigenvalue.

  FIG. 19 is an explanatory diagram for explaining a Q method (eQuivalence method) which is an example of a physical reduction method. As shown in the figure, in the Q method, when the areas other than the main system are “contracted area 1” and “contracted area 2”, “contracted area 1” is contracted to generator “G1”, The contracted area 2 is contracted to the generator “G2”. Here, the power generation capacity of “G1” is the sum of the power generation capacities in “contracted area 1”, and the control system (AVR, PSS) of “G1” is the control system of the largest generator in “contracted area 1”. The power generation capacity of “G2” is the sum of the power generation capacities in “contracted area 2”, and the control system of “G2” is the control system of the largest generator in “contracted area 2”.

  On the other hand, FIG. 20 is an explanatory diagram for explaining a mode reduction method which is an example of a mathematical reduction method. As shown in the figure, in the mode contraction method, contraction is performed by mode-decomposing the contraction target system on the mathematical formula and merging the modes obtained by the decomposition.

JP 2006-5992 A Japanese Patent Application Laid-Open No. 10-56735 Hiraiwa, Hayashi, Uchida "Improvement of accuracy of engineering system reduction by modal analysis", Electron Theory B, Vol. 124, No. 7, 2004

  However, the physical reduction method is easy to physically understand the reduction target system, but there is no guarantee that the main oscillation mode is preserved before and after reduction, and there is a problem that analysis accuracy is poor. On the other hand, the mathematical reduction method has a problem that the mode of interest can be saved, but the reduction target system is expressed as a linear differential equation and is difficult to understand physically, and is difficult to combine with an existing simulation program.

  The present invention has been made in order to solve the above-described problems caused by the prior art, and includes a power system contraction model creation device, a power system contraction system, and a power system contraction model creation device capable of efficiently contracting a large-scale power system with high accuracy. It is an object to provide an approximately model creation method and a power system contraction model creation program.

  In order to solve the above-described problems and achieve the object, the present invention is an electric power system contracted model creating apparatus that creates a contracted model by contracting parts other than the main system from an electric power system model. A pre-contraction oscillation mode estimation means for estimating a plurality of principal oscillation modes of a previous model using an ARMA model, and a main eigenvalue is inverted from the principal oscillation modes estimated by the pre-contraction oscillation mode estimation means. Eigenvalue estimation means before reduction using an iterative method, physical reduction model generation means for creating a physical reduced model obtained by reducing the model by physical reduction, and a plurality of reduced models Post-contraction oscillation mode estimation means for estimating the main oscillation mode using an ARMA model, and estimation of main eigenvalues from the main oscillation mode estimated by the post-contraction oscillation mode estimation means using an inverse iteration method A post-reduction eigenvalue estimation means, a plurality of vibration modes estimated by the pre-contraction vibration mode estimation means, and a plurality of vibration modes estimated by the post-contraction vibration mode estimation means before and after reduction A swing mode set specifying means for specifying a plurality of sets using inner products between eigenvectors of the swing mode, and a PSS constant and an AVR constant of the generator control system of the physical contraction model created by the physical contraction model creation means And a swing mode matching means for matching each swing mode set specified by the swing mode set specifying means by adjusting the inertia constant based on the eigenvectors and eigenvalues before and after contraction.

  Further, the present invention is a power system contracted model creation method by a power system contracted model creating apparatus that creates a contracted model by contracting a part other than the main system from a model of the power system. A pre-contraction oscillation mode estimation step for estimating a plurality of principal oscillation modes of the model using the ARMA model, and a main eigenvalue from the primary oscillation mode estimated by the pre-contraction oscillation mode estimation step by performing an iterative method. A pre-contraction eigenvalue estimation step to estimate using, a physical contraction model creation step to create a physical contraction model that contracts the model by physical contraction, and multiple major fluctuations of the model after contraction A post-contraction oscillation mode estimation step that estimates a mode using an ARMA model, and a main eigenvalue is inverted from the main oscillation mode estimated by the post-contraction oscillation mode estimation step. A reduced eigenvalue estimation step estimated using a method, a plurality of oscillation modes estimated by the pre-contraction oscillation mode estimation step, and a plurality of oscillation modes estimated by the after-contraction oscillation mode estimation step A swing mode set specifying step of specifying a plurality of sets of swing modes to be matched before and after using inner products between the eigenvectors of the swing mode, and generator control of the physical contracted model created by the physical contracted model creating step And adjusting the PSS constant, AVR constant, and inertia constant of the system based on the eigenvectors and eigenvalues before and after the contraction to match each of the oscillation mode sets identified by the above-mentioned oscillation mode set identification step. It is characterized by that.

  Further, the present invention is a power system contracted model creation program for creating a contracted model by contracting a part other than the main system from a model of the power system, and comprising a plurality of main oscillation modes of the model before contraction A pre-contraction oscillation mode estimation procedure using an ARMA model, and a pre-contraction eigenvalue estimation using a reverse iteration method from the main oscillation modes estimated by the pre-contraction oscillation mode estimation procedure An estimation procedure, a physical contraction model creation procedure for creating a physical contraction model in which the model is contracted by physical contraction, and estimation of a plurality of main shaking modes of the contracted model using an ARMA model A post-contraction oscillation mode estimation procedure, a post-contraction eigenvalue estimation procedure for estimating main eigenvalues from the main oscillation modes estimated by the post-contraction oscillation mode estimation procedure using an inverse iteration method, and the reduction The inner product between the eigenvectors of the oscillation mode is set to a set of the oscillation modes that match the plurality of oscillation modes estimated by the previous oscillation mode estimation procedure and the plurality of oscillation modes estimated by the post-contraction oscillation mode estimation procedure before and after contraction. A plurality of oscillation mode set identification procedures to be used, and a PSS constant, an AVR constant and an inertia constant of the generator control system of the physical contraction model created by the physical contraction model creation procedure, and eigenvectors before and after contraction, and By adjusting based on the eigenvalue, the computer is caused to execute a swing mode matching procedure for matching each swing mode set specified by the swing mode set specifying procedure.

  According to this invention, a plurality of principal oscillation modes of the model before reduction are estimated using the ARMA model, and principal eigenvalues are estimated from the estimated principal oscillation modes using the inverse iteration method, and physical reduction is performed. Create a physical contraction model that is reduced by approximating the model, estimate multiple main vibration modes of the reduced model using the ARMA model, and repeat the main eigenvalues from the estimated main vibration modes Using the inner product between the eigenvectors of the sway mode to identify multiple sets of sway modes to be matched before and after the reduction from the estimated sway modes before and after the reduction. Because the PSS constant, AVR constant, and inertia constant of the machine control system are adjusted based on the eigenvectors and eigenvalues before and after contraction, each set of oscillation modes is matched. It is possible to save the upset mode.

  According to the present invention, since the main oscillation modes are stored before and after contraction, there is an effect that a highly accurate physical contraction model can be created.

  Exemplary embodiments of a power system contracted model creating apparatus, a power system contracted model creating method, and a power system contracted model creating program according to the present invention will be described below in detail with reference to the accompanying drawings.

  First, the outline | summary of the system | strain reduction method based on a present Example is demonstrated. FIG. 1 is an explanatory diagram for explaining the outline of the system contraction method according to the present embodiment. As shown in the figure, in the system reduction method according to the present embodiment, the main oscillation mode included between the systems is estimated by applying an ARMA (Auto-Regressive Moving Average) model to the oscillation waveform ( 1). Details of the ARMA model will be described later.

  Then, pairing of the estimated swing modes before and after contraction is performed (2). Here, the inner product of both eigenvectors is used to evaluate the degree of coincidence of the oscillation modes before and after contraction. In order to match the fluctuation mode before and after contraction, the PSS (Power System Stabilizer) constant, the AVR (Automatic Voltage Regulator) constant, and the inertia constant of the generator are set to the eigenvalue sensitivity. Use to adjust (3).

  In this way, identification of the oscillation mode using the ARMA model, evaluation based on the eigenvector inner product of the oscillation mode match before and after reduction, and adjustment of the oscillation mode before and after reduction by adjusting the PSS constant, AVR constant and inertia constant are performed. As a result, the dynamic characteristics before contraction can be automatically reproduced with high accuracy.

  Next, the configuration of the system contraction model creation device according to the present embodiment will be described. FIG. 2 is a functional block diagram illustrating the configuration of the system contraction model creation apparatus according to the present embodiment. As shown in the figure, the system contraction model creation apparatus 100 includes a data reading unit 110, a system data storage unit 120, a physical contraction unit 130, a swing mode estimation unit 140, and a swing mode insertion unit. 150 and a data output unit 160.

  The data reading unit 110 is a processing unit that reads the system configuration data defining the system model and the assumed accident data from a storage device such as a file and stores them in the system data storage unit 120.

  In addition to the system configuration data and assumed accident data read by the data reading unit 110, the system data storage unit 120 generates generator speed deviation data, interconnection line phase difference data, contracted model data, oscillation mode, It is a storage unit for storing data used for creating a contracted model such as eigenvalues.

  The physical contraction unit 130 is a processing unit that creates a contracted model by contracting a system other than the main system by a physical method using the system configuration data stored in the system data storage unit 120. The reduced model data is stored in the system data storage unit 120.

  The fluctuation mode estimation unit 140 performs transient stability calculation using the assumed accident data for each of the original system model and the contracted model, and the observed time series such as the obtained generator speed deviation or interconnection line phase difference. Is a processing unit that estimates the main oscillation mode by applying the ARMA model, and includes a transient stability analysis unit 141, a time series analysis unit 142, and an eigenvalue analysis unit 143.

  The transient stability analysis unit 141 is a processing unit that performs transient stability calculation using assumed accident data, and generates an observation time series such as a generator speed deviation and a connection line phase difference. The time series analysis unit 142 is a processing unit that estimates main oscillation modes by applying the ARMA model to the observation time series generated by the transient stability analysis unit 141.

  The eigenvalue analysis unit 143 is a processing unit that calculates main eigenvalues by applying an inverse iteration method from main oscillation modes obtained by time series analysis. Here, the inverse iterative method is one of eigenvalue analysis methods, and is a method of calculating a true eigenvalue closest to the approximate eigenvalue by giving an estimated eigenvalue as an approximate value and performing iterative calculation.

That is, when the approximate value σ for the eigenvalue λ _i of A is known, the equation (1) is established. By applying the inverse iteration method to this equation, | λ _i −σ | and | λ _k −σ | having the minimum absolute value can be calculated.

A true eigenvalue λ _k can be obtained by adding σ to the eigenvalue obtained in this way.

  The shaking mode combining unit 150 is a processing unit that combines the shaking mode of the original system model and the shaking mode of the contracted model using the shaking mode estimated by the shaking mode estimation unit 140, and includes a mode matching unit 151, And a parameter adjustment unit 152.

  The mode matching unit 151 is a processing unit that identifies a pair of swing modes for performing matching between the swing modes estimated from the system models before and after contraction. The mode matching unit 151 specifies a pair of rocking modes using an inner product between eigenvectors of the rocking mode.

In general, the inner product of the vector space is expressed by equation (3), and the phase angle difference expresses the degree of coincidence of the vectors. That is, if COSθ is 0, A and B are vertical, and if COSθ is 1, it means parallel.

When this is expanded as an inner product of a linear space as A = (α ₁ , α ₂ ,..., Α _n ) and B = (β ₁ , β ₂ ,..., Β _n ), the inner product of the linear space is defined by equation (4) Is done. However, α and β are complex numbers.

Unlike the vector space, the inner product of the linear space is defined as a complex number. However, when both vectors are orthogonal (A, B) = 0, the property of changing to 0 is not changed (this is called generalized orthogonality). A normalized value MAC is defined as follows. This MAC (Modal Assurance Criterion) is used in the field of mechanical vibrations for comparison between experimental results and the finite element method.

  Here, when the two vectors are orthogonal, the MAC value is 0, and when the same mode is used, the value is 1. In practice, various disturbances and estimation errors are included in the eigenvector values, so that the closer the MAC value is to 1, the more the same mode can be considered.

  FIG. 3 shows a MAC value between three eigenvectors as a sample example. It can be confirmed from the MAC value that the vectors A and B are the same, and that C is a vector different from A or B.

When applying the inner product to eigenvectors before and after contraction, the dimensions are different, so it is necessary to match the dimensions before and after contraction using a weighted average of vector elements corresponding to each contracted area as shown in equation (6). . In Equation (6), X _i represents a vector element in the non-contracted area, and therefore is applied as it is.

  Further, the mode matching unit 151 selects, as the adjustment target generator, a generator that is most sensitive to the eigenvalue before contracting, that is, a generator that has been contracted including those having a large eigenvector value.

  The parameter adjustment unit 152 adjusts the contraction model parameters, that is, the inertia constant M, the AVR constant (gain, phase compensation element), and the PSS constant (gain, phase compensation element) for each pair specified by the mode matching unit 151. It is a processing unit for adjusting the shaking mode.

The oscillation frequency in one machine infinite system is expressed by equation (7). Here, K represents the synchronizing force, and ω ₀ represents the rated angular velocity. That is, it can be seen that the influence of the inertia constant M on the oscillation frequency is large. In the Q method, the inertial constant M of the reduced generator uses a weighted average of the active power output of the generator included in the reduced system, but no linear characteristic with respect to the capacity is seen for a large-capacity inertial constant. Therefore, the inertia constant is used as an adjustment parameter for matching the oscillation frequency.

However, although it is the inertia constant of the reduced generator, it is not appropriate to freely adjust the constant expressing the physical characteristics. Therefore, adjustment of the inertia constant is limited to the extent that parameter adjustment by AVR is made efficient.

  AVR plays a role of maintaining the synchronizing power of the generator by keeping the generator terminal voltage constant. That is, it can be generally said that the influence on the oscillation frequency is strong as shown in the equation (7).

  When contracting by the Q method, the excitation system provided in the maximum capacity generator in the reduction system or the most installed type excitation system is applied as the excitation system of the equivalent reduction generator. This problem is not limited to the Q method, but is a problem that all reduction methods have, but it cannot be said that the replaced excitation system is necessarily appropriate.

  Therefore, the thyristor type super-fast response excitation shown in FIG. 4 is provided for the adjustment generator after contraction, and the parameters are adjusted to match the oscillation frequency. The three adjustment parameters are phase lag, lead element and gain.

  The parameter adjustment unit 152 adjusts the PSS constant after adjusting the inertia constant and AVR. Since the PSS varies the terminal voltage as an auxiliary input of the AVR and increases the equivalent damping torque in the positive direction to stabilize the power fluctuation, the influence on the attenuation rate of the eigenvalue is large. The adjustment effect on the frequency is small. In consideration of the auxiliary signal for AVR, the adjustment order is adjusted separately to the inertia constant, AVR, and PSS.

  In order to adjust the PSS constant, a PSS is provided to the generator for adjustment in the same manner as the AVR. FIG. 5 shows a control block diagram of the PSS. The control block employs a single-stage ΔP type PSS and corresponds to a high gain control block as shown in the figure. Therefore, when switching from the excitation system already installed in the adjustment generator to FIG. 5, the waveform and the dominant eigenvalue (main eigenvalue) may change greatly. As adjustment parameters, the phase delay, the leading element, and the gain are targeted as in the case of AVR.

  FIG. 6 is a diagram illustrating a parameter adjustment processing procedure by the parameter adjustment unit 152. As shown in the figure, the parameter adjustment unit 152 adjusts the inertia constant, AVR, and PSS constant for each generator independently in this order using a predetermined parameter adjustment logic.

  In the parameter adjustment logic, parameters are adjusted using the eigenvalue sensitivity. It is assumed that the system matrix A changes by ΔA with respect to a minute change in the system parameter α. Then, the eigenvalue also slightly changes by Δλ. As a result, the system changes as Δλ / Δα. This is the eigenvalue sensitivity. Using this, the evaluation function is minimized so that the target eigenvalue (before reduction) matches the adjusted eigenvalue (after reduction) and these eigenvectors match. The evaluation function is shown in equation (9). The parameter adjustment unit 152 calculates a minimum j among the evaluation functions related to the parameter j and K at that time, and repeats until the evaluation function satisfies the convergence condition.

  Returning to FIG. 2, the data output unit 160 is a processing unit that reads out and outputs the data of the reduced model after parameter adjustment from the system data storage unit 120, and the output data of the reduced model is stability analysis or the like Used in

  Next, a processing procedure of the system contraction model creation device 100 according to the present embodiment will be described. FIG. 7 is a flowchart illustrating the processing procedure of the system contraction model creation device 100 according to the present embodiment. As shown in the figure, in the system contraction model creation device 100, the data reading unit 110 reads the system configuration data and the assumed accident data and stores them in the system data storage unit 120 (step S1).

  Then, the transient stability analysis unit 141 performs system stability analysis with reference to the system data storage unit 120, calculates the generator speed deviation and interconnection line phase difference before contraction, and stores them in the system data storage unit 120. Store (step S2).

  Then, the time series analysis unit 142 applies the ARMA model to the generator speed deviation or interconnection line phase difference before contraction, performs time series analysis, and estimates the main oscillation mode before contraction. That is, the oscillation frequency, attenuation rate, and eigenvector of the main oscillation mode before contraction are calculated (step S3). Then, the eigenvalue analysis unit 143 performs eigenvalue analysis using the inverse iteration method, and calculates eigenvalues corresponding to the main oscillation modes before contraction (step S4).

  Then, the physical contraction unit 130 reads the system configuration data from the system data storage unit 120 to perform physical contraction, and stores the contracted model data in the system data storage unit 120 (step S5).

  Then, the transient stability analysis unit 141 performs system stability analysis with reference to the system data storage unit 120, and calculates the reduced generator speed deviation and interconnection line phase difference (step S6). Then, the time series analysis unit 142 applies the ARMA model to the generator speed deviation or the interconnection phase difference after contraction to perform time series analysis, and estimates the main oscillation mode after contraction. That is, the oscillation frequency, attenuation rate, and eigenvector of the main oscillation mode after contraction are calculated (step S7). Then, the eigenvalue analysis unit 143 performs eigenvalue analysis using the inverse iteration method, and calculates eigenvalues corresponding to the main shake modes after contraction (step S8).

  Then, the mode matching unit 151 performs the swing mode pairing before and after the contraction based on the inner product of the eigenvectors before and after the contraction, and selects the adjustment generator (step S9). Then, the parameter adjustment unit 152 adjusts the inertia constant, the AVR constant, and the PSS constant using the eigenvalue sensitivity, and updates the contracted model of the system data storage unit 120 (step S10). When the parameter adjustment by the parameter adjustment unit 152 is completed, the data output unit 160 reads the reduced model data from the system data storage unit 120 and outputs it (step S11).

  As described above, by adjusting the parameters so that the eigenvectors and eigenvalues of the main oscillation modes before and after the reduction match, the accuracy of the reduction model can be improved.

  Next, the parameter adjustment logic shown in FIG. 6 will be described. FIG. 8 is a flowchart showing the parameter adjustment logic taking the adjustment of the PSS constant as an example. As shown in FIG. 8, the parameter adjustment unit 152 first performs eigenvalue calculation as a base value using an inverse iteration method (step S21).

  Then, eigenvalue calculation is performed on values obtained by changing each PSS constant (delay time constant, advance time constant, gain constant) by a predetermined value, and eigenvalue sensitivity relating to each PSS constant is calculated (steps S22 to S24). ).

  When the eigenvalue sensitivity for each PSS constant is calculated (No in step S25), the evaluation function is minimized using the eigenvalue sensitivity (step S26), and the PSS constant that minimizes the evaluation function and its constant multiple K are obtained. (Step S27).

  Then, it is determined whether or not the evaluation function is larger than a predetermined threshold (step S28). If the evaluation function is larger than the predetermined threshold, the PSS constant obtained in step S26 is multiplied by K (step S29), and the process returns to step S21. . On the other hand, if it is not greater than the predetermined threshold, the adjustment of the PSS constant is terminated.

  In this way, the PSS constant can be adjusted while obtaining the PSS constant that minimizes the evaluation function using the eigenvalue sensitivity.

Next, time series analysis using the ARMA model by the time series analysis unit 142 will be described. For the observation time series x _t , the output at a certain point in time is a linear combination of past outputs (x _t , x _t−1 , x _t−2 ...) And white noise (e _t , e _t−1 , e _If the system can be expressed by a linear combination of _t-2 ...), the system is called an ARMA model and is expressed by the following equation.

と定義すると(１２)式は(１５)式のように表わされ、このシステムの伝達関数Ｈ(ｚ^-1)は(１６)式で表現される（図９参照)。 Here, using the delay operator z ^-1

(12) is expressed as (15), and the transfer function H (z ⁻¹ ) of this system is expressed as (16) (see FIG. 9).

At this time, the output response of the system is determined by the arrangement of characteristic roots of the transfer function. I.e. p _k is a vibration waveform when the complex conjugate, p _k will indicate the monotonous response if real roots. Since the application to the fluctuation waveform of the power system is considered here, the former vibration waveform is mainly targeted.

If the ARMA model is high order, characteristic roots corresponding to the order of the model are calculated as in equation (17). Therefore, it is necessary to specify a shaking mode having a large influence from the estimated shaking modes. In the case of the ARMA model, the influence degree of each mode can be expressed by the power spectral density shown in the equation (18). This represents the influence degree of the vibration component included in the single observation waveform, and it corresponds to the multi-point data. Correlation cannot be expressed. Therefore, a cross spectrum is applied to obtain the correlation between the two points.

The ARMA model at time t is defined by the following equation, where x (t) is a response time series (reference waveform) at a specific point in the power system and y (t) is a response time series at an arbitrary point.

Here, a (i) represents the AR coefficient, and b _x (j) and b _y (j) represent the MA coefficients of the reference waveform and the observed waveform y (t), respectively. Moreover, e _x (t) and e _y (t) are white noises in the respective observed waveforms and satisfy the following conditions.

Σ _exy ^{2 in the} equation (21) means that there is a correlation when e _x and e _y are the same time. In the ARMA model defined by equation (19), the same AR coefficient a (i) is applied to two observed waveforms. This is because the AR coefficient is a vibration characteristic (natural frequency) included in the power system. It is assumed that the AR coefficients are equal regardless of all observed waveforms.

Here, a cross-correlation function R _xy (τ) is applied to obtain the correlation between x (t) and y (t). As described above, when a transfer function is obtained from each time series, it is expressed as follows.

However, a (0) = 1, and each transfer characteristic is as follows.

From the Wiener-Khintchine theorem, the Fourier transform of the autocorrelation function is a power spectrum, whereas the cross-correlation function R _xy (τ) is defined as a cross spectrum S _xy . For discrete signals x (t) and y (t), the cross spectrum S _xy is defined as follows:

By expanding the equation (27) using the impulse response sequence of the observed waveform, the cross spectrum of the ARMA model is expressed using the transfer characteristic of the system as shown in the equation (28).

  Here, the physical meaning of the spectrum estimated from the waveform between two points is as follows. When the original signal is shifted by a time interval p (referred to as lag in time series analysis) to check whether there is periodicity for a certain observation signal, p = T, 2T, 3T, etc. If a large correlation is seen, it can be seen that this signal contains a large component of period T. That is, the periodicity of the waveform can be estimated from the correlation with the waveform shifted by the lag p. This is called an autocorrelation function, and a power spectrum is expressed by an energy density function in a frequency band.

On the other hand, when it is desired to know the degree of relationship between two signals, the correlation between the two signals is confirmed by a cross-correlation function. However, the cross-correlation function does not always have a maximum value at the lag p = 0, and is generally not symmetric with respect to the lag p = 0. The cross-correlation function is expressed by an energy density function in the frequency band, which is a cross spectrum, and the expression (27) expressed using the delay operator z ₋₁ is the frequency spectrum X of x (t) and y (t). When expressed as a function of frequency ω using (ω) and Y (ω), it is expressed by equation (29).

である。よってその位相角は

となる。 Here, X (ω) is a complex function, so

It is. So the phase angle is

It becomes.

  FIG. 10 shows a diagram of the phase angle of the cross spectrum. In order to deal with physical phenomena, the real axis component actually makes sense. That is, the phase difference between the two greatly affects the absolute value of the cross spectrum. By calculating the cross spectrum in this way, the phase relationship between the frequency components can be grasped. Here, by using x (t) as a reference waveform and y (t) as various observation data, the phase relationship with respect to a specific frequency obtained from the characteristic root and its spectral density are grasped from the cross spectrum, and from each observation point. By comparing the obtained results on the same plane, the fluctuation relation of each observation data, that is, the eigenvector can be estimated.

  In order to grasp the transfer characteristic of the ARMA model expressed by the equation (19), it is necessary to estimate the AR coefficient and the MA coefficient. Here, the Burg method, which is one of MEM (maximum entropy method), is applied to determine the AR coefficient. The Burg method is known to have a high resolution for observation time series having a short data length. When the transient phenomenon of the power system handled here is considered to be about 10 to 20 seconds, the sampling frequency is determined to be about 100 Hz and the analysis method of the Y method. Therefore, the number of data points is about several hundred to 2,000, and the usefulness of the Burg method can be expected when the sampling frequency of the observed waveform is small. The Burg method is described in J. P. Burg, “Maximum Entropy Spectral Analysis”, Proc 37th Meeting Sco. Of Explaration Geophysists Oklahoma, 1967.

By applying the algorithm of Levinson-Durbin (hereinafter LD) to the AR model (Nakamizo “Signal Analysis and System Identification”, Corona, 2000), the m-th order AR coefficient is sequentially determined from the m-1th order coefficient. Can be calculated. Burg method (see (34)) before addition to the location prediction error e _t and minimizes the backward prediction error r _t with respect LD method is a method of obtaining the reflection coefficients k _m. (34) by substituting e _t and r _t in equation (33) of the type, which can be obtained k _m that minimizes Differentiating, determining the m-th order AR coefficient by repeating this sequentially Can do.

The cross spectrum is effective for specifying the main oscillation mode from the observed waveform. In order to obtain the cross spectrum as shown by the equation (28), in addition to the AR coefficient, the MA coefficient and the mean square of the residual σ _exy ² must be estimated. The AR coefficient can be obtained by the Burg method, but in general, the MA coefficient is difficult to obtain analytically, and can be obtained approximately by using a numerical calculation method such as a maximum likelihood method. However, if spectral estimation is the purpose, the spectrum can be estimated without the need for MA coefficients.

When variable conversion of i = 1 and j = 1 + k is performed from the expressions (25) and (26), the following expression is obtained.

On the other hand, the ARMA model is divided into variables r (t) and s (t) as follows.

Here, the cross-correlation function of r (t) and s (t) is expressed as follows using the AR term.

R _xy is a cross-correlation function of x (t) and y (t) and can be calculated from the observation time series. Similarly, using the MA term in equation (36) and equation (21),

(38) compared to those obtained by multiplying the z _-l to (35), defines a cross-spectrum without the MA terms as those (28) by substituting the equation (39) below that the equation Can do.

The vibration characteristics of the system are estimated based on the AR coefficient obtained from the observation time series and the cross spectrum. By clarifying the relationship between the ARMA model cross-spectrum and the vibration system as shown in the literature (Kanazawa “Structural damage detection method based on microtremor measurement”, Electromagnetic Research Institute report, T00026, 2001/4) Can be identified. In Equation (39), if the pole of A (z ⁻¹ ) = 0 is z = z ₁ , z ₂ ,..., Z _p , the pole of A (z) = 0 is z = z ₁ ⁻¹ , z ₂ ⁻¹ ,..., Z _p ⁻¹ , and when expressed as a (0) = 1, it is expressed by equation (40).

と表わされ、(４１)式をラプラス平面に変換すると次式で表現することができる。ここでｓはラプラス演算子である。

Furthermore, if partial fractional decomposition of equation (40)

When the equation (41) is converted into a Laplace plane, it can be expressed by the following equation. Here, s is a Laplace operator.

β _j and γ _j are residues for the poles −s _j and s _j and are given by the following equations.

The constant term S ₀ is given by the following equation.

  Here, considering the poles of the ARMA model, since these are obtained as complex numbers, naturally, the conjugate complex number is also obtained as a characteristic root. However, since the imaginary part of the pole represents the frequency of the vibration mode in the real function, the imaginary part cannot be negative in the physical phenomenon. Therefore, these correspond to so-called fake modes observed on the analysis model. Further, as shown in the equation (42), a mode in which the signs of the real parts are opposite to each other is also estimated, but this also represents the attenuation rate of the waveform (usually takes a negative value and becomes an attenuation stable waveform). Therefore, it is unlikely that both modes will appear in the observed waveform. That is, these modes are also observed as fake modes. Thus, the mode observed in the actual waveform with respect to the order of the model is a part thereof.

The attenuation rate and oscillation frequency of the main mode obtained from s _j are expressed as follows.

At this time, since the value obtained from the reference waveform x (t) is used for the AR coefficient, s _j , that is, the mode attenuation rate and the oscillation frequency are common (actually, the AR coefficient is different for each observed waveform). Because they are different, the estimated attenuation rate and oscillation frequency are not the same value).

  Therefore, by treating the residual value γi for each mode of each waveform as a vector, it is possible to express the magnitude and phase relationship of fluctuations between points with respect to the frequency. As a result, the eigenvector can be estimated from only the observed waveform, and the vibration characteristics of the main oscillation modes can be grasped (see FIG. 11).

  As described above, if the order is increased, so many modes are calculated, but when applied to the power system, the vibration components that can be confirmed as the system oscillation phenomenon are limited. , Filtering mode.

  The target of estimation is the fluctuation phenomenon in the electric power system, especially the inter-system fluctuation mode that contributes to the stability in a wide area. Since these oscillation frequencies are usually about 1 Hz, the upper and lower limits of the oscillation frequency are set to 0.0 to 2.0 Hz.

  For 10-second data, if the sampling frequency is 10 Hz, the number of data points is 100, and even if the data has a maximum frequency (Nyquist frequency) of 5 Hz and a minimum frequency of 0.1 Hz, it is sufficient system characteristics. Can be grasped.

  A vibration component with a large damping rate is very stable and often cannot be observed as a waveform. Generally, in an electric power system, if the attenuation rate is greater than -1.0 [1 / sec], the waveform is likely to be detected, so the upper limit is set to -3.0 [1 / sec] (comparison using only the attenuation rate) Rather than doing this, evaluation is often made with a damping ratio expressed by an attenuation rate corresponding to the oscillation frequency, but this applies because the S method is based on the attenuation rate).

In order to rank vibration components included in the reference waveform in descending order, the cross spectrum S _xy is used. That is, the spectrum absolute value | S _xy (f0) | with respect to the oscillation frequency f0 specified from the AR coefficient is used as a criterion. However, it is conceivable that a false mode that does not actually contribute to fluctuation is calculated as the order of the model increases. Therefore, in order to avoid such erroneous determination of the mode, as shown in FIG. 12, it is determined that the cross spectrum of the designated frequency substantially shows a peak value as the main oscillation mode. Here, it is almost assumed that f0 is a peak value obtained from the power spectrum of the reference waveform, and the peak value of the cross spectrum of each waveform includes an estimation error of the time series model, and these peak values do not necessarily match. Because.

  In addition, when considering the application of time series analysis to fluctuations in the power system, various observation time series can be selected as the reference waveform, but in order to increase the model identification accuracy, It is desirable to apply. In system analysis, assumed accident failures correspond to these disturbances, but voltage, current, reactive power, etc. are likely to be affected by the accident and are suitable for phenomenon analysis at the time of the accident. In the case of analysis, it is not suitable because it is easily affected by disturbance. On the other hand, since the phase angle and the phase difference have small steep fluctuations compared to the above-mentioned various amounts, it can be expected that the identification accuracy can be maintained while being influenced by the disturbance.

  FIG. 13 shows a flowchart of the eigenvector calculation procedure by the time series analysis unit 142. FIG. 13 shows a case where the interconnection phase difference and the generator speed deviation are used as the observation time series. As illustrated in FIG. 13, the time series analysis unit 142 estimates the AR coefficient by applying the Burg method to the interconnection line phase difference (step S31).

  And the power spectrum regarding a connection line phase difference is computed (step S32), and the frequency and attenuation factor of a system oscillation mode are estimated (step S33). Then, the residue γ is obtained from the speed deviation of each generator (steps S34 to S36). That is, a cross spectrum is calculated for each generator speed deviation (step S34), and the residue γ of each oscillation mode is estimated (step S35). Then, an eigenvector is estimated from each γ (step S37).

  Next, the effect of contraction by the system contraction model creation apparatus 100 according to the present embodiment will be described. Here, the effect of contraction by the system contraction model creation apparatus 100 according to the present embodiment was confirmed for the case where the area 1 to 5 is contracted by three machines with FIG. 14 as the original system. Three units are reduced as one equivalent generator in areas 1 and 2, and similarly, one unit for areas 4 and 5 and one generator for area 3; It is a model that expresses.

  15 (a) and 15 (c) show the results of the conventional contraction, and it can be seen that the interconnection effective tidal current and the internal phase difference are greatly deviated from the results of the original system. On the other hand, FIGS. 15B and 15D show reduction results by the system reduction model creation apparatus 100 according to the present embodiment. From the waveform, it can be seen that the degree of coincidence with the original system is improved compared to the conventional contraction. That is, it can be seen that the dynamic characteristics of the original system are reproduced even after reduction.

  FIG. 16 shows the transition of the shaking mode at this time. In the figure, ◯ represents the eigenvalue before reduction and before adjustment, and x represents the dominant eigenvalue before reduction. As can be seen from the figure, the change in the vertical axis direction, that is, the fluctuation frequency is adjusted by adjusting the inertia constant and the AVR constant, and the attenuation rate of each mode is adjusted by adjusting the PSS constant. In this example, the AVR constant is not adjusted. This is because the oscillation frequency is sufficiently improved by the adjustment based on the inertia constant, and it is not necessary to adjust the AVR. From such characteristics, when adjusting the inertia constant and the AVR constant, only the oscillation frequency is considered in the evaluation function.

  The dominant eigenvalue before contraction and the eigenvalue after adjustment are as shown in FIG. In this example, two dominant eigenvalues were confirmed. However, since it is known from the eigenvector that the influence of areas 4 and 5 is very strong for both modes, the reduced generator 3140 is adjusted in both modes. It is said.

  As described above, in the present embodiment, the time series analysis unit 142 estimates the main swing modes before and after contraction using the ARMA model, and the eigenvalue analysis unit 143 corresponds to the main swing modes before and after contraction, respectively. The eigenvalues to be calculated are calculated using the inverse iteration method, the mode matching unit 151 identifies a pair of main oscillation modes before and after contraction using the inner product of the eigenvectors, and the parameter adjustment unit 152 performs the main oscillation mode before and after the contraction. Since the inertia constant, the AVR constant, and the PSS constant of the contracted model are adjusted so as to match, a highly accurate contracted model can be created efficiently.

  In the present embodiment, the system contraction model creation apparatus has been described. However, by realizing the configuration of the system contraction model creation apparatus with software, a system contraction model creation program having the same function can be obtained. it can. Therefore, a computer that executes this system contraction model creation program will be described.

  FIG. 18 is a functional block diagram illustrating the configuration of a computer that executes a system contraction model creation program according to the present embodiment. As shown in the figure, the computer 200 includes a RAM 210, a CPU 220, an HDD 230, a LAN interface 240, an input / output interface 250, and a DVD drive 260.

  The RAM 210 is a memory that stores a program and a program execution result, and the CPU 220 is a central processing unit that reads the program from the RAM 210 and executes the program. The HDD 230 is a disk device that stores programs and data, and the LAN interface 240 is an interface for connecting the computer 200 to other computers via the LAN. The input / output interface 250 is an interface for connecting an input device such as a mouse or a keyboard and a display device, and the DVD drive 260 is a device for reading / writing a DVD.

  The system contraction model creation program 211 executed in the computer 200 is stored in the DVD, read from the DVD by the DVD drive 260, and installed in the computer 200. Alternatively, the system contraction model creation program 211 is stored in a database or the like of another computer system connected via the LAN interface 240, read from these databases, and installed in the computer 200. The installed system reduction model creation program 211 is stored in the HDD 230, read into the RAM 210, and executed by the CPU 220.

  In this embodiment, the system contraction model creation apparatus has been described. However, the present invention is not limited to this, and the system model after contraction by inputting the system model before and after physical contraction is the main. The present invention can be similarly applied to a reduced model adjusting apparatus that adjusts a reduced model so as to preserve the mode. That is, it is possible to remove the physical contraction unit 130 from the system contraction model creation apparatus 100 and input data of the physical contraction model.

  As described above, the power system contracted model creation device, the power system contracted model creation method, and the power system contracted model creation program according to the present invention are useful for power system stability analysis, etc. This is suitable for analyzing a complex power system.

It is explanatory drawing for demonstrating the system | strain reduction method of the system | strain reduction model creation apparatus which concerns on a present Example. It is a functional block diagram which shows the structure of the system | strain reduction model production apparatus which concerns on a present Example. It is a figure which shows the sample example of a MAC value. It is an AVR control block diagram. It is a PSS control block diagram. It is a figure which shows the process sequence of the parameter adjustment by a parameter adjustment part. It is a flowchart which shows the process sequence of the system | strain reduction model creation apparatus which concerns on a present Example. It is a flowchart which shows the parameter adjustment logic which makes adjustment of a PSS constant an example. It is a figure which shows the transfer function of an ARMA model. It is a figure which shows phase angle (theta) _xy of a cross spectrum. It is a figure which shows the eigenvector calculation method from each waveform. It is a figure which shows the determination method of main spectrum components. It is a flowchart which shows the eigenvector calculation procedure by a time series analysis part. It is a figure which shows an IEEJ WEST30 machine standard system model. It is a figure which shows the effect of reduction by the system | strain reduction model creation apparatus which concerns on a present Example. It is a figure which shows transition of the shaking mode by parameter adjustment. It is a figure which shows the adjustment result of a parameter. It is a functional block diagram which shows the structure of the computer which executes the system | strain reduction model creation program which concerns on a present Example. It is explanatory drawing for demonstrating Q method which is an example of the physical reduction method. It is explanatory drawing for demonstrating the mode reduction method which is an example of the mathematical reduction method.

Explanation of symbols

100 System Reduction Model Creation Device 110 Data Reading Unit 120 System Data Storage Unit 130 Physical Reduction Unit 140 Oscillation Mode Estimation Unit 141 Transient Stability Analysis Unit 142 Time Series Analysis Unit 143 Eigenvalue Analysis Unit 150 Oscillation Mode Integration Unit 151 Mode Matching unit 152 Parameter adjusting unit 160 Data output unit 200 Computer 210 RAM
211 System Reduction Model Creation Program 220 CPU
230 HDD
240 LAN interface 250 I / O interface 260 DVD drive

Claims (8)

A power system contracted model creation device that creates a contracted model by contracting parts other than the main system from a power system model,
A pre-contraction oscillation mode estimation means for estimating a plurality of principal oscillation modes of the model before reduction using an ARMA model;
A pre-contraction eigenvalue estimation means for estimating a main eigenvalue from the main vibration modes estimated by the pre-contraction vibration mode estimation means;
Physical contraction model creation means for creating a physical contraction model in which the model is contracted by physical contraction;
A post-contraction oscillation mode estimation means for estimating a plurality of principal oscillation modes of the contracted model using the ARMA model;
A post-contraction eigenvalue estimation means for estimating a main eigenvalue from the main oscillation modes estimated by the post-contraction shake mode estimation means using an inverse iteration method;
A swing mode that specifies a plurality of sets of swing modes that are matched before and after contraction from a plurality of swing modes estimated by the pre-contraction swing mode estimation means and a plurality of swing modes estimated by the post-contraction swing mode estimation means A pair identification means;
Each of the oscillation mode set identification means specified by adjusting the parameters of the generator control system of the physical reduction model created by the physical reduction model creation means based on eigenvectors and eigenvalues before and after reduction. A power system contraction model creating apparatus comprising: a shaking mode matching means for matching a shaking mode group.   2. The power system contraction model creation device according to claim 1, wherein the fluctuation mode set identification unit identifies a fluctuation mode set using an inner product between eigenvectors of the fluctuation mode.   The swing mode matching means adjusts the PSS constant, AVR constant and inertia constant of the generator control system of the physical contraction model created by the physical contraction model creation means by adjusting the swing mode set specifying means. 3. The power system contraction model creation device according to claim 1 or 2, wherein each specified oscillation mode set is matched. A power system contracted model creation device that creates a contracted model by contracting parts other than the main system from a power system model,
Pre-contraction oscillation mode estimation means for estimating a plurality of main oscillation modes of the model before reduction using time series analysis;
A pre-contraction eigenvalue estimation means for estimating a main eigenvalue from the main vibration modes estimated by the pre-contraction vibration mode estimation means;
Physical contraction model creation means for creating a physical contraction model in which the model is contracted by physical contraction;
A post-contraction oscillation mode estimation means for estimating a plurality of main oscillation modes of the contracted model using time series analysis;
A post-contraction eigenvalue estimation means for estimating a main eigenvalue from the main oscillation modes estimated by the post-contraction shake mode estimation means using an inverse iteration method;
Between the eigenvectors of the oscillation mode, a set of oscillation modes that match the plurality of oscillation modes estimated by the oscillation mode estimation means before contraction and the plurality of oscillation modes estimated by the oscillation mode estimation means after contraction before and after contraction. Oscillating mode set specifying means for specifying a plurality using inner products of
Each of the oscillation mode set identification means specified by adjusting the parameters of the generator control system of the physical reduction model created by the physical reduction model creation means based on eigenvectors and eigenvalues before and after reduction. A power system contraction model creating apparatus comprising: a shaking mode matching means for matching a shaking mode group.
  The swing mode matching means adjusts the PSS constant, AVR constant and inertia constant of the generator control system of the physical contraction model created by the physical contraction model creation means by adjusting the swing mode set specifying means. 5. The power system contraction model creation device according to claim 4, wherein each specified oscillation mode set is matched. A power system contracted model creation device that creates a contracted model by contracting parts other than the main system from a power system model,
Pre-contraction oscillation mode estimation means for estimating a plurality of main oscillation modes of the model before reduction using time series analysis;
A pre-contraction eigenvalue estimation means for estimating a main eigenvalue from the main vibration modes estimated by the pre-contraction vibration mode estimation means;
Physical contraction model creation means for creating a physical contraction model in which the model is contracted by physical contraction;
A post-contraction oscillation mode estimation means for estimating a plurality of main oscillation modes of the contracted model using time series analysis;
A post-contraction eigenvalue estimation means for estimating a main eigenvalue from the main oscillation modes estimated by the post-contraction shake mode estimation means using an inverse iteration method;
A swing mode that specifies a plurality of sets of swing modes that are matched before and after contraction from a plurality of swing modes estimated by the pre-contraction swing mode estimation means and a plurality of swing modes estimated by the post-contraction swing mode estimation means A pair identification means;
By adjusting the PSS constant, AVR constant and inertia constant of the generator control system of the physical contraction model created by the physical contraction model creation means based on the eigenvectors and eigenvalues before and after the contraction, A power system contraction model creating apparatus comprising: a swing mode matching means for matching each swing mode set specified by the specifying means. A power system contracted model creation method by a power system contracted model creating apparatus that creates a contracted model by contracting parts other than the main system from a power system model,
A pre-contraction oscillation mode estimation step for estimating a plurality of principal oscillation modes of the model before reduction using an ARMA model;
A pre-contraction eigenvalue estimation step for estimating a main eigenvalue from the main oscillation modes estimated by the pre-contraction oscillation mode estimation step using an inverse iteration method;
A physical contraction model creation step for creating a physical contraction model in which the model is contracted by physical contraction;
A post-contraction oscillation mode estimation step for estimating a plurality of main oscillation modes of the reduced model using the ARMA model;
A post-contraction eigenvalue estimation step for estimating a main eigenvalue from the main oscillation mode estimated by the post-contraction oscillation mode estimation step using an inverse iteration method;
Between the eigenvectors of the oscillation mode, a set of oscillation modes that match the plurality of oscillation modes estimated by the pre-contraction oscillation mode estimation step and the plurality of oscillation modes estimated by the after-contraction oscillation mode estimation step before and after the contraction A swing mode set specifying step for specifying a plurality of times using an inner product of
By adjusting the PSS constant, AVR constant and inertia constant of the generator control system of the physical contraction model created by the physical contraction model creation step based on the eigenvectors and eigenvalues before and after contraction, A method for creating a power system contraction model, characterized by comprising: a swing mode matching step for matching each swing mode set specified by the specific step. A power system contracted model creation program for creating a contracted model by contracting parts other than the main system from a power system model,
A pre-contraction oscillation mode estimation procedure for estimating a plurality of main oscillation modes of the model before reduction using an ARMA model;
A pre-contraction eigenvalue estimation procedure for estimating a main eigenvalue from the main oscillation mode estimated by the pre-contraction oscillation mode estimation procedure using an inverse iteration method;
A physical reduced model creation procedure for creating a physical reduced model obtained by reducing the model by physical reduction;
A post-contraction oscillation mode estimation procedure for estimating a plurality of main oscillation modes of the reduced model using the ARMA model;
A post-contraction eigenvalue estimation procedure for estimating a main eigenvalue from the main oscillation mode estimated by the post-contraction oscillation mode estimation procedure using an inverse iteration method;
Between the eigenvectors of the oscillation mode, a set of the oscillation modes that are matched by the plurality of oscillation modes estimated by the pre-contraction oscillation mode estimation procedure and the plurality of oscillation modes estimated by the after-contraction oscillation mode estimation procedure before and after reduction. Oscillating mode group identification procedure for identifying a plurality using inner product of
By adjusting the PSS constant, AVR constant and inertia constant of the generator control system of the physical contraction model created by the physical contraction model creation procedure based on the eigenvectors and eigenvalues before and after contraction, A power system contraction model creation program that causes a computer to execute a swing mode matching procedure for matching each set of swing modes specified by a specific procedure.

Images