JP2017194811A - Frequency characteristic analysis apparatus, method, and program - Google Patents

Abstract

A frequency characteristic analysis apparatus, method, and program capable of significantly reducing the time required for frequency characteristic calculation while maintaining calculation accuracy are provided. A frequency characteristic analyzing apparatus 1 generates a circuit equation to be analyzed according to a calculation frequency using a sparse matrix, and a circuit equation created by a circuit equation creating unit 14a. An ordering unit 14b that performs ordering by complete pivoting only on a sparse matrix of a circuit equation created in a specific case, and an LU that performs LU decomposition of the sparse matrix using the ordering result most recently performed by the ordering unit 14b A decomposition unit 14c, and a solution calculation unit 14e that performs a solution calculation using the matrix subjected to LU decomposition by the LU decomposition unit 14c. [Selection] Figure 1

Description

  The present invention relates to a frequency characteristic analysis apparatus, method, and program, and more particularly to a frequency characteristic analysis apparatus, method, and program that are useful when applied to frequency characteristic analysis of a power system.

  For analysis of overvoltage and overcurrent generated in a power system, analysis of power quality, analysis of dynamic characteristics of an AC / DC converter, instantaneous value analysis capable of analyzing temporal changes in power and current at a waveform level is effective. Many of the instantaneous value analysis programs that can perform such instantaneous value analysis have a function to calculate the frequency characteristics of impedance from the specified bus, and the AC / DC converter of the DC power transmission equipment and frequency conversion equipment and AC It is used for analysis of transient overvoltage including harmonic instability and harmonic characteristics generated by system interaction.

  FIG. 12 is a flowchart showing a calculation process of conventional frequency characteristic calculation. As shown in FIG. 12, the frequency characteristic calculation is usually a calculation that repeatedly executes a calculation for obtaining an AC steady solution of the circuit (AC steady solution calculation) while changing the frequency. The AC steady solution calculation includes a formulation process for creating a circuit equation that is a simultaneous linear equation and a solution process for solving the created circuit equation.

The formulation process is a process of creating circuit equations from information on all nodes and all circuit components in the circuit. For example, if it is assumed that N _n independent nodes and N _b circuit components (branches) exist in the circuit, the vector u and N _b circuit components composed of the voltages of N _n nodes in the circuit flow. vector i consisting of current, handling the vector v composed of the voltage across the N _b-number of circuit components as unknown variables, the circuit equation becomes the following equation (1).

The above equation (1) is referred to as the Spasterblow equation, and F is referred to as a tableau matrix. F, x, and y in the above formula (1) are represented by the following formula (2).

Here, A is an N _b × N _n branch-to-node connection matrix that associates connection states of N _n independent nodes and N _b branches in the circuit. I is a unit matrix of N _b × N _b . B = [B ₁ −B ₂ ] is a branch characteristic matrix of N _b × 2N _b that associates the branch voltage v and the branch current i of each branch. s is a vector determined by the values of the voltage source and current source in the circuit. When N = N _n × 2N _b , x and y are vectors having N elements, and F is an N × N matrix. In the AC steady solution calculation, all the quantities of the Spasterblow equation shown here are complex numbers.

  The branch characteristic matrix B is obtained from the branch characteristic equation. When the frequency is f = ω / (2π), for example, the branch characteristic equations of the inductor L and the capacitor C are expressed by the following equations (3) and (4), respectively.

However, i _L and i _C are currents (branch currents) flowing through L and C, and v _L and v _C are both-end voltages (branch voltages) of L and C, respectively. When L is the b-th branch, the b-th row of B and s is expressed by the following equation (5).

However, the left side from the “|” mark on the left side of equation (5) indicates the b-th row of the matrix B ₁ , the right side from the “|” mark indicates the b-th row of the matrix −B ₂ , and the right side is the b row of s. Showing eyes. The two arrows at the bottom of the above equation (5) indicate that 1 and −1 / jωL are the b-th columns of B ₁ and B ₂ , respectively, and the elements omitted by “...” Are 0. If a nonlinear element is arranged in the analysis circuit, the nonlinear element is ignored or handled by linearizing with characteristics near the origin. The vector s determined by the values of the voltage source and the current source in the circuit substitutes E _m e ^jθ for the corresponding element of s when the voltage waveform generated by the voltage source is E _m cos (ωt + θ). For a current source that generates I _m cos (ωt + θ), I _m e ^jθ is substituted.

In the solution process, processes called ordering, LU decomposition, and forward / backward substitution are executed in order. Ordering is an operation of rearranging matrix elements for the purpose of improving calculation accuracy. For example, the matrix elements are rearranged in the following procedure.
(A) For the matrix F, the number of non-zero elements NR present in each row and the number of non-zero elements NC present in each column are obtained.
(B) NRC = (NR−1) × (NC−1) is calculated for each element, and the element with the smallest NRC is replaced so as to be the diagonal element (pivot) in the first row. When there are a plurality of elements having the same NRC, the element having the largest absolute value is selected.
(C) The remaining matrix (submatrix) excluding 1 row and 1 column is considered as F, and the procedures (a) to (c) are repeated until the last row.

  Here, the procedure (b) is an operation for selecting an element having the smallest NRC (an element having a large absolute value when NRC is the same) and making it a diagonal element. In the search, two methods can be considered, in which the search target is only the i-row element and all elements except the row and column for which the diagonal elements have already been determined. The former is called partial pivoting, and the latter is called full pivoting. In general, partial pivoting requires a small amount of calculation because the search range is narrow, and complete pivoting requires a large amount of calculation because the search range is wide, but the calculation accuracy is high because a more optimal diagonal element can be selected.

LU decomposition and forward / backward substitution are operations called forward substitution and backward substitution by decomposing a coefficient matrix (tableau matrix F in the above equation (1)) into a lower triangular matrix L and an upper triangular matrix U. Is an operation for obtaining x which is an unknown. Now, the matrix element of the above equation (1) is a _ij (where the first subscript is the row number, the second subscript is the column number, i = 1, 2,..., N, j = 1, 2, .., N), the element l _ij of the lower triangular matrix L after the LU decomposition of F and the element u _ij of the upper triangular matrix U can be calculated by the following equations (6) and (7), respectively. .

Next, x is obtained by forward / backward substitution using the obtained u _ij and l _ij . A detailed description of the forward / backward substitution process following the LU decomposition will be omitted. By obtaining x, the above-mentioned node voltage, branch current, and branch voltage are obtained as an AC steady solution. The frequency characteristics of the circuit can be obtained by repeatedly performing the series of processes (formulation and solution) described above by updating the calculation frequency f. The following Non-Patent Documents 1 and 2 disclose a method for reducing the calculation time in the normal instantaneous value analysis for obtaining the time change of voltage and current, although it does not relate to the frequency characteristic calculation.

Ryomichi Yonezawa and Satoshi Noda, “Development of Instantaneous Value Analysis Program for Power System (Part 5) —Acceleration of XTAP by Improvement of Large Scale Sparse Matrix LU Decomposition Algorithm”, Research Report H11009, 2012 Ryomichi Yonezawa and Satoshi Noda, “Development of Power System Instantaneous Value Analysis Program (Part 7) —Acceleration of XTAP by Parallelizing Internal State Update Processing of Circuit Components”, Research Report H13005, 2014

  By the way, the solution based on the above-described LU decomposition is performed not only by the frequency characteristic calculation but also by the normal instantaneous value analysis for obtaining the time change of the voltage and current. Normal instantaneous value analysis is a calculation process in which the above equation (1) is solved at each calculation time step. However, when there is no time-varying element or nonlinear element in which the element value changes with time, In equation (1), only y changes for each calculation time step, and F does not change. That is, in the normal instantaneous value analysis, the ordering and LU decomposition for F need only be performed once at first, and the solution x of each calculation time step can be performed simply by performing forward substitution and backward substitution for changing y. Can be sought.

  However, in the frequency characteristic calculation, the branch characteristic equation changes by updating the calculation frequency f, so that F changes for each set calculation frequency f. For this reason, in the solution process in the frequency characteristic calculation, ordering and LU decomposition cannot be omitted as in the solution process in normal instantaneous value analysis, and it is necessary to execute ordering and LU decomposition for each calculation frequency f. As described above, in the frequency characteristic calculation, unlike the case of the instantaneous value analysis, more time is required for the solution-finding process, and the reduction of the calculation time is a problem.

  Here, as described above, partial pivoting or complete pivoting can be used for ordering in the solution process. In order to shorten the calculation time of the solution process in the frequency characteristic calculation, it may be considered to use partial pivoting which requires less calculation time than complete pivoting. Partial pivoting is theoretically inferior to perfect pivoting in principle, but in many cases, the degradation of analysis accuracy is not a problem. That's not true. Therefore, it is necessary to reduce the calculation time of the solution process while using complete pivoting from the viewpoint of maintaining calculation accuracy.

  The present invention has been made in view of the above circumstances, and an object thereof is to provide a frequency characteristic analysis apparatus, method, and program capable of significantly reducing the time required for frequency characteristic calculation while maintaining calculation accuracy. And

In order to solve the above problems, a frequency characteristic analyzing apparatus (1) of the present invention includes a circuit equation creating unit (14a) that creates a circuit equation to be analyzed according to a calculation frequency using a sparse matrix, and the circuit Of the circuit equations created by the equation creation unit, an ordering unit (14b) that performs ordering by complete pivoting only on the sparse matrix of the circuit equation created in a specific case, and the ordering unit that is most recently performed An LU decomposition unit (14c) that performs LU decomposition of the sparse matrix using the ordered result, and a solution calculation operation unit (14e) that performs solution calculation using the matrix subjected to LU decomposition by the LU decomposition unit.
In the frequency characteristic analyzing apparatus of the present invention, the ordering unit performs the ordering by the complete pivoting when the circuit equation is first created by the circuit equation creating unit.
The frequency characteristic analysis apparatus of the present invention further includes a calculation accuracy evaluation unit (14d) that evaluates calculation accuracy using the matrix subjected to LU decomposition by the LU decomposition unit, and the ordering unit includes a calculation accuracy evaluation. When it is evaluated that the calculation accuracy has deteriorated in the section, the ordering is performed by the complete pivoting.
Further, in the frequency characteristic analyzing apparatus of the present invention, the calculation accuracy evaluation unit obtains an unstable factor (ρ) indicating an index of calculation accuracy using a matrix subjected to LU decomposition by the LU decomposition unit, and the unstable factor When the value exceeds a predetermined threshold value (ρ _th ), it is evaluated that the calculation accuracy has deteriorated.
The frequency characteristic analysis method according to the present invention includes a first step (S30) in which a circuit equation to be analyzed corresponding to a calculation frequency is created using a sparse matrix, and a specification among the circuit equations created in the first step. A second step (S41) for ordering by complete pivoting only for the sparse matrix of the circuit equation created in the case of (2), and using the ordering result most recently performed in the second step, the LU of the sparse matrix A third step (S42) for performing decomposition, a fourth step (S44) for performing a solution calculation using the matrix subjected to LU decomposition in the third step, and the first to fourth steps by changing the calculation frequency. The fifth step (S20) is repeated.
The frequency characteristic analysis method of the present invention further includes a sixth step (S43) for evaluating calculation accuracy using the matrix subjected to LU decomposition in the third step, and the second step includes the first step. When the circuit equation is first created in one step, or when it is evaluated that the calculation accuracy is deteriorated in the sixth step, the ordering by the complete pivoting is performed.
The frequency characteristic analysis program of the present invention is a circuit created by a circuit equation creating means (14a) for creating a circuit equation to be analyzed according to a calculation frequency using a sparse matrix, and a circuit created by the circuit equation creating means. Among the equations, the ordering means (14b) for performing ordering by complete pivoting only on the sparse matrix of the circuit equation created in a specific case, and the sparseness using the ordering result most recently performed by the ordering means. An LU decomposition means (14c) that performs LU decomposition of a matrix and a solution calculation operation means (14e) that performs solution calculation using the matrix subjected to LU decomposition by the LU decomposition means function.
The frequency characteristic analysis program according to the present invention further causes the computer to function as a calculation accuracy evaluation unit (14d) that evaluates calculation accuracy using a matrix subjected to LU decomposition by the LU decomposition unit, and the ordering unit includes: When the circuit equation is first created by the circuit equation creating unit, or when it is evaluated that the calculation accuracy has been deteriorated by the calculation accuracy evaluating unit, the ordering by the complete pivoting is performed.

  According to the present invention, a circuit equation to be analyzed corresponding to a calculation frequency is created using a sparse matrix, and only a sparse matrix of a circuit equation created in a specific case among the created circuit equations is complete. Ordering is performed by pivoting, LU decomposition of a sparse matrix is performed using the most recent ordering result, and solution calculation is performed using the LU decomposed matrix. As a result, since ordering by complete pivoting in cases other than a specific case can be omitted, there is an effect that the time required for frequency characteristic calculation can be greatly shortened while maintaining calculation accuracy.

It is a block diagram which shows the principal part structure of the frequency characteristic analyzer by one Embodiment of this invention. It is a flowchart which shows an example of the frequency characteristic analysis method by one Embodiment of this invention. It is a figure which shows the circuit for a verification used in order to verify the relationship between a circuit scale and calculation speed. It is a figure which shows the equivalent circuit used for simulation of the generators G2-G10 in FIG. It is a figure which shows the equivalent circuit used for simulation of transformer TR2-TR10 in FIG. It is a figure which shows the calculation time required when the frequency characteristic calculation was performed about "Case 1"-"Case 9" shown in FIG. It is a figure which shows the relationship between a circuit scale and speed improvement multiple. It is a figure which shows the instability factor in each sample point at the time of performing frequency characteristic calculation about "Case 9" shown in FIG. It is a figure which shows the circuit for verification used in order to verify calculation accuracy and calculation speed. It is a figure which shows the calculation result of the frequency characteristic calculation performed with respect to the circuit for verification of FIG. It is a figure which shows the instability factor in each sample point at the time of performing frequency characteristic calculation about the verification circuit of FIG. It is a flowchart which shows the calculation process of the conventional frequency characteristic calculation.

  Hereinafter, a frequency characteristic analyzing apparatus, method, and program according to an embodiment of the present invention will be described in detail with reference to the drawings. In the following, the case where the analysis target is a power system will be described as an example, but the present invention is not limited to the power system, and any electric circuit can be the analysis target. .

<Frequency characteristics analyzer>
FIG. 1 is a block diagram showing a main configuration of a frequency characteristic analyzing apparatus according to an embodiment of the present invention. As shown in FIG. 1, the frequency characteristic analysis apparatus 1 of this embodiment includes an input device 11, an output device 12, an auxiliary storage device 13, and an arithmetic processing device 14, and calculates frequency characteristics of a power system as an analysis target. I do. Such a frequency characteristic analyzing apparatus 1 is realized by, for example, a desktop personal computer.

  The input device 11 includes, for example, a keyboard and a pointing device, and outputs an instruction according to an operation of a user who uses the frequency characteristic analysis device 1 to the arithmetic processing device 14. The output device 12 includes, for example, a liquid crystal display device, a printer, and the like, and outputs various information output from the arithmetic processing device 14. The auxiliary storage device 13 includes, for example, an HDD (hard disk drive), an SSD (solid state drive), and the like, and various types of information (ordering result OR and evaluation parameter PR: details will be described later) used in analyzing the power system. ) And a frequency characteristic analysis program PG and the like.

  The arithmetic processing unit 14 uses the various information stored in the auxiliary storage device 13 according to the frequency characteristic analysis program PG stored in the auxiliary storage device 13 based on the operation instruction input from the input device 11, and the frequency characteristics of the power system. Perform the calculation. The arithmetic processing unit 14 includes a circuit equation creation unit 14a (circuit equation creation unit), an ordering unit 14b (ordering unit), an LU decomposition unit 14c (LU decomposition unit), a calculation accuracy evaluation unit 14d (calculation accuracy evaluation unit), and The solution calculation unit 14e (solution calculation means) is provided.

The circuit equation creation unit 14a creates a circuit equation of the power system to be analyzed using a sparse matrix. Specifically, the circuit equation creating unit 14a creates the sputter blow equation represented by the above-described equation (1) according to the calculation frequency f. The above sparse matrix is a tableau matrix F shown in the above-described equation (2). When the configuration of the power system to be analyzed does not change and the calculation frequency f changes, the branch characteristic matrix B = [B ₁ -B ₂ ] in the tableau matrix F shown in the above-described equation (2) changes.

  The ordering unit 14b performs ordering by complete pivoting on the tableau matrix F of the circuit equation created by the circuit equation creating unit 14a. Here, the ordering unit 14b performs ordering by complete pivoting only on the tableau matrix F of the circuit equation created in a specific case among the circuit equations created by the circuit equation creation unit 14a. This is to reduce the time required for frequency characteristic analysis.

  Specifically, the ordering unit 14b performs ordering by complete pivoting when the circuit equation is first created by the circuit equation creation unit 14a or when the calculation accuracy evaluation unit 14d evaluates that the calculation accuracy has deteriorated. Do. Note that the result of ordering by complete pivoting performed by the ordering unit 14b is stored in the auxiliary storage device 13 as the ordering result OR.

  The LU decomposition unit 14 c performs LU decomposition of the tableau matrix F using the ordering result OR stored in the auxiliary storage device 13. Here, the ordering result OR stored in the auxiliary storage device 13 is updated every time the ordering unit 14b performs ordering by complete pivoting. For this reason, it can be said that the LU decomposition unit 14c performs LU decomposition of the tableau matrix F using the ordering result most recently performed by the ordering unit 14b.

The calculation accuracy evaluation unit 14d evaluates the calculation accuracy using the matrix subjected to LU decomposition by the LU decomposition unit 14c. Specifically, the calculation accuracy evaluation unit 14d obtains an unstable factor ρ indicating an index of calculation accuracy shown in the following equation (8), and when the unstable factor ρ becomes equal to or greater than a predetermined threshold ρ _th. In addition, it is evaluated that the calculation accuracy has deteriorated. In the following equation (8), the element l _ij of the lower triangular matrix L and the element u _ij of the upper triangular matrix U are collectively expressed as lu _ij . The threshold value ρ _th is stored in the auxiliary storage device 13 as the evaluation parameter PR.

As apparent from the above equation (8), the instability factor ρ is an increase ratio of the infinite norm of the matrix elements before and after the LU decomposition. Normally, the denominator of the above equation (8) becomes a value near 1 by a preprocessing called scaling (or normalization, equalization), so the instability factor ρ is approximately equal to the value of the maximum element of lu _ij. Become. In ordering, when an element that is not appropriate as a diagonal element, such as a value having a small absolute value, is selected as the diagonal element, the denominator of the above-described equation (7) becomes a small value, so that l _ij , u _ij increases and the instability factor ρ increases. Therefore, it is possible to determine whether or not the ordering is appropriate by using the unstable factor ρ.

The solution calculation unit 14e performs a solution calculation using the matrix (lower triangular matrix L, upper triangular matrix U) subjected to LU decomposition by the LU decomposition unit 14c. Specifically, x shown in the above equations (1) and (2) is obtained by forward / backward substitution using the element u _ij of the upper triangular matrix U and the element l _ij of the lower triangular matrix L.

  The frequency characteristic analysis program PG is recorded on a computer-readable recording medium such as a CD-ROM or DVD (registered trademark) -ROM, or downloaded via a network such as the Internet. 1 is stored in the auxiliary storage device 13. By reading and executing the frequency characteristic analysis program PG stored in the auxiliary storage device 13, the function of each block of the frequency characteristic analysis device 1 (for example, the circuit equation creation unit 14a, the ordering unit 14b, the LU decomposition unit) 14c, the calculation accuracy evaluation unit 14d, and the solution calculation unit 14e) are realized by software. That is, these functions are realized by cooperation of software and hardware resources.

<Frequency characteristics analysis method>
Next, a frequency characteristic analysis method using the above-described frequency characteristic analysis apparatus 1 will be described. When analyzing the frequency characteristic, first, the user inputs the power system to be analyzed. Specifically, the user operates the input device 11 of the frequency characteristic analysis apparatus 1 and performs a process of inputting information on all nodes and circuit components constituting the power system to the arithmetic processing device 14 in accordance with a user instruction. Done.

Next, an analysis condition is input by the user. For example, a user operating the input device 11 of the frequency characteristic analyzer 1 calculates the calculation start frequency f _begin , the calculation end frequency f _end , the number of sample points K, and the calculation frequency (equal intervals or equal intervals on the logarithmic axis). Etc. are input.

When the above input is completed and an instruction to start analysis is given by the user, calculation of frequency characteristics of the power system is started according to the flowchart shown in FIG. FIG. 2 is a flowchart showing an example of a frequency characteristic analysis method according to an embodiment of the present invention. When the frequency characteristic calculation is started, first, various initial processes are performed by the arithmetic processing unit 14 shown in FIG. 1 (step S10). For example, a variable k for specifying a sample point is initialized (a value is set to 1), and a calculation frequency f _k at each sample point is obtained.

Next, processing for updating the calculation frequency f is performed by the arithmetic processing unit 14 (step S20). Here, since the value of the variable k in the process at the step S10 is set to 1, calculating the frequency f is updated to calculate the frequency f ₁ (set). The calculation frequency f ₁ is, for example, the calculation start frequency f _begin input by the user. Thereby, the first sample point is specified.

  Next, a formulation process (creation of power system circuit equations) is performed (step S30: first step). Specifically, using the information of all nodes and circuit components constituting the power system input by the user, the Superstar blow equation corresponding to the calculated frequency f updated (set) in step S20 ((1) described above) , (2), the circuit equation creating unit 14a performs the process of creating the equation (2). When the formulation process is completed, a solution finding process is performed (step S40).

  When the solution process starts, first, the ordering unit 14b performs ordering by complete pivoting on the tableau matrix F of the circuit equation created by the circuit equation creation unit 14a (step S41: second step). Note that the result of ordering by complete pivoting performed by the ordering unit 14b is stored in the auxiliary storage device 13 as the ordering result OR. Next, the LU decomposition unit 14c performs LU decomposition of the tableau matrix F using the ordering result OR stored in the auxiliary storage device 13 (step S42: third step). That is, the first processing (calculation) of the LU decomposition unit 14c is performed after the ordering by the ordering unit 14b is performed.

Next, processing for evaluating calculation accuracy using the matrix subjected to LU decomposition by the LU decomposition unit 14c is performed by the calculation accuracy evaluation unit 14d (step S43: sixth step). Specifically, the unstable factor ρ indicating the calculation accuracy index shown in the above-described equation (8) is obtained using the matrix subjected to the LU decomposition by the LU decomposition unit 14c, and the unstable factor ρ is determined from the threshold value ρ _th . The calculation accuracy evaluation unit 14d performs a process for determining whether the value is smaller. The above threshold value ρ _th is stored in the auxiliary storage device 13 as the evaluation parameter PR.

When it is determined that the instability factor ρ is smaller than the threshold ρ _th (when the determination result in step S43 is “yes”), a solution calculation using the matrix decomposed by the LU decomposition unit 14c. Is performed by the solution calculation unit 14e (step S44: fourth step). Specifically, the operation for obtaining x shown in the above-described equations (1) and (2) by forward / backward substitution using the element u _ij of the upper triangular matrix U and the element l _ij of the lower triangular matrix L is a solution operation. This is performed by the unit 14e.

When the above processing is completed, the arithmetic processing unit 14 determines whether or not the value of the variable k has become larger than the sample score K (step S50). When it is determined that the value of the variable k is equal to or less than the number K of sample points (when the determination result of step S50 is “no”), the process of updating the calculation frequency f after incrementing the value of the variable k is calculated. This is performed by the processing device 14 (step S20: fifth step). Here, the calculation frequency f is updated to calculate the frequency f _2. Thereby, the next sample point is specified.

  When the calculation frequency f is updated, the circuit equation creation unit 14a performs a process of creating a circuit equation corresponding to the updated calculation frequency f (step S30). Thereafter, the LU decomposition unit 14c performs the LU decomposition of the tableau matrix F using the ordering result OR stored in the auxiliary storage device 13 without performing the ordering by the ordering unit 14b (step S42). That is, since the processing (calculation) performed by the LU decomposition unit 14c corresponds to the second and subsequent calculations, LU decomposition of the tableau matrix F is performed using the ordering result most recently performed by the ordering unit 14b.

Subsequently, a process for evaluating the calculation accuracy is performed by the calculation accuracy evaluation unit 14d (step S43). When it is determined that the unstable factor ρ is smaller than the threshold ρ _th (when the determination result in step S43 is “yes”), the solution calculation (step S44), the value of the variable k, and Comparison with the number K of sample points (step S50) is performed in order. When it is determined that the value of the variable k is less than or equal to the number K of sample points (when the determination result of step S50 is “no”), the process of updating the calculation frequency f after incrementing the value of the variable k It is performed again by the device 14 (step S20).

As described above, the instability factor ρ is smaller than the threshold ρ _th (the determination result in step S43 is “yes”), and the value of the variable k is equal to or less than the sample score K (the determination result in step S50 is In the case of “no”), the ordering unit 14b does not perform the ordering (step S41), the calculation frequency f is updated (step S20), and the circuit equation corresponding to the calculation frequency f is created (step S30). LU decomposition (step S42) and solution calculation (step S44) are repeatedly performed.

If the instability factor ρ obtained in step S43 is _{equal to} or greater than the threshold ρ _th during the above processing, the determination result in step S43 is “no”. Then, ordering by complete pivoting is performed by the ordering unit 14b (step S41), and the ordering result is stored in the auxiliary storage device 13 as the ordering result OR. Thereafter, LU decomposition using the ordering result OR stored in the auxiliary storage device 13 is performed by the LU decomposition unit 14c (step S42). As described above, when the instability factor ρ obtained in step S43 is equal to or greater than the threshold ρ _th , the ordering unit 14b performs ordering, and LU decomposition is performed using the newly obtained ordering result OR. .

  When the value of the variable k is equal to or less than the number K of sample points (when the determination result in step S50 is “no”), the series of processes described above is repeated, and the frequency characteristics of the power system for each sample point are sequentially determined. Desired. The series of processes shown in FIG. 2 ends when it is determined that the value of the variable k is equal to or less than the number K of sample points (when the determination result in step S50 is “no”).

  As described above, in this embodiment, when the difference between the frequency at which the calculation is completed and the frequency to be calculated next is small, attention is paid to the fact that the change in the circuit equation obtained at each frequency is small. The ordering result performed in the past is reused without performing the ordering for each sample point. For this reason, it is possible to significantly reduce the time required for frequency characteristic calculation.

In the present embodiment, the ordering by the ordering unit 14b is not always omitted, but the ordering by the ordering unit 14b is performed as necessary. Specifically, when it is evaluated that the calculation accuracy has deteriorated (when the unstable factor ρ becomes equal to or greater than the threshold ρ _th ), the ordering unit 14b performs ordering. For this reason, calculation accuracy can also be maintained.

<Verification of calculation speed and calculation accuracy>
The inventor of the present application actually performed frequency characteristic analysis of the power system using the frequency characteristic analysis apparatus 1 described above, and verified the calculation speed and calculation accuracy. In the frequency characteristic calculation of the present embodiment, the calculation efficiency is improved by omitting the ordering when the calculation accuracy is not a problem. In general, the calculation amount of the matrix operation such as the ordering affects the matrix order N. Receive. Therefore, in the following, first, the relationship between the matrix order N and the calculation speed will be verified, and then the calculation accuracy and calculation speed will be verified in comparison with the conventional example. The matrix order N is determined by the number of nodes N _n and the number of branches N _b (N = N _n × 2N _b ), so it can be considered as a circuit scale. Therefore, N is hereinafter referred to as a circuit scale.

(A) Relationship between circuit scale and calculation speed Since the calculation amount of ordering by complete pivoting is considered to be affected by the circuit scale N, when the circuit scale N changes, the speed improvement multiple (computation speed compared with the conventional method) It is considered that the index that shows how many times improvement of Therefore, for the purpose of confirming the relationship between the circuit scale N and the speed improvement rate, a speed improvement multiple when the circuit scale N is changed was calculated.

  FIG. 3 is a diagram showing a verification circuit used to verify the relationship between the circuit scale and the calculation speed. The verification circuit shown in FIG. 3 includes generators G2 to G10, transformers TR2 to TR10, and transmission lines TL1 to TL16, and the generators G2 to G10 are connected to the power transmission network via the transformers TR2 to TR10, respectively. It is a circuit connected to (a power transmission network composed of power transmission lines TL1 to TL16). The circuit scale of this verification circuit was changed to nine cases from “Case 1” to “Case 9” shown in the figure, and frequency characteristics were calculated for each case.

  For example, “Case 1” is a circuit having a scale composed of a generator G2, a transformer TR2, and transmission lines TL1 and TL2, and “Case 2” is a generator G2, G3, transformers TR2, TR3, and a transmission line. This is a circuit having a scale composed of the electric wires TL1, TL2, TL3, and TL11. “Case 9” is a circuit having a scale including all of the generators G2 to G10, the transformers TR2 to TR10, and the power transmission lines TL1 to TL16 illustrated in FIG. The number of nodes and the number of branches in each of “Case 1” to “Case 9” are shown in Table 1 below.

  FIG. 4 is a diagram showing an equivalent circuit used for simulating the generators G2 to G10 in FIG. As shown in FIG. 4A, the equivalent circuit used for simulating the generators G2 to G10 is a circuit in which an inductance L corresponding to the next transient reactance is connected in series to the electromotive force E of the generator. The next transient impedance of the generator has a characteristic that the loss increases due to the skin effect as the frequency increases. For example, it can be simulated by the equivalent circuit shown in FIG. Here, since the main purpose is to evaluate the calculation speed, it is assumed that the next transient impedance has no loss component and is simply simulated by the equivalent circuit shown in FIG. The inductances L of the generators G2 to G7 and G9 were 32.1 [μH], the inductance L of the generator G8 was 64.2 [μH], and the inductance L of the generator G10 was 10.7 [μH]. The power supply voltage E was short-circuited.

FIG. 5 is a diagram showing an equivalent circuit used for simulating the transformers TR2 to TR10 in FIG. As shown in FIG. 5 (a), an equivalent circuit used to simulate the transformer TR2~TR10 is ideal transformer delta-Y connection, a circuit leakage inductance L ₁ connected in series. Similarly to the generator, the transformer's leakage impedance also has a characteristic that loss increases due to the skin effect, and it is considered that it can be simulated by, for example, the equivalent circuit shown in FIG. Here, for the same reason as the simulation of the generator, it is assumed that the leakage impedance has no loss component and is simply simulated by the equivalent circuit shown in FIG.

Transformer TR2～TR7, leakage inductance _{L 1} of TR9 is a 9.28 [mH], leakage inductance _{L 1} of the transformer TR8 is a 18.6 [mH], leakage inductance _{L 1} of the transformer TR10 3.10 [MH]. The transformation ratio was 500 [kV]: 22 [kV], and the ground floating capacitance of the transformer winding was simulated with a capacitance of 2000 [pF] for all transformers. The transmission lines TL1 to TL16 should be realistically simulated using a frequency-dependent distributed constant line model. However, since the main purpose is to evaluate the calculation speed in this study, the transmission lines TL1 to TL16 are simplified. A π-type equivalent circuit model (one stage) having no coupling between phases was used.

  In any case of “Case 1” to “Case 9” shown in FIG. 3, impedance frequency characteristics were calculated from the bus A in FIG. Here, the wider the calculation frequency, the greater the difference in circuit equations between frequency samples, which is a disadvantageous condition for frequency characteristic calculation in the present embodiment. Therefore, in the verification here, the frequency characteristic calculation effect in this embodiment is confirmed by setting the frequency range from 0.1 [Hz] to 1 [MHz], which is wider than the frequency range that is normally performed. The frequency sample points were 2000 points at equal intervals on the logarithmic axis. For the calculation, a desktop computer, which is currently widely used, was used.

Here, in a numerical calculation using a general computer, it is standard to use a double precision operation, and the number of significant digits is about 16 digits. That is, when the above-mentioned unstable factor ρ exceeds 16 digits, information loss occurs in the difference calculation in the calculation executed during the LU decomposition (the calculation shown in the above-described equations (6) and (7)). However, there is a possibility that the calculation accuracy is lowered. For this reason, in this calculation, the threshold ρ _{th is set} to 10 ¹¹ with an allowance of 5 digits. It should be noted that “case 1” to “case 9” used here are cases created for verification, and the frequency characteristics obtained as a calculation result have no special meaning.

  FIG. 6 is a diagram showing calculation time required when frequency characteristics are calculated for “Case 1” to “Case 9” shown in FIG. The legend “conventional method” in FIG. 6 is obtained when the frequency characteristics are calculated by the conventional method shown in FIG. 12, and the legend “proposed method” is the method according to the present embodiment shown in FIG. This is when frequency characteristics are calculated. Referring to FIG. 6, in the conventional method, the calculation time increases in the order of “Case 1” to “Case 9” (that is, as the circuit size increases), whereas in the proposed method, the circuit size increases. However, it can be seen that the calculation time hardly changes.

  FIG. 7 is a diagram showing the relationship between the circuit scale and the speed improvement multiple. FIG. 7 is a plot of the ratio of calculation time (value obtained by dividing the calculation time of the conventional method by the calculation time of the proposed method) as a speed improvement multiple for each case shown in FIG. The circuit scale of each case shown in FIG. 6 is obtained from Table 1 described above. Referring to FIG. 7, it can be seen that the speed enhancement factor increases as the circuit scale increases. In particular, in the case of “Case 9” shown in FIG. 6 (when the circuit scale is 1353), the speed improvement multiple is close to 60 times, and it can be seen that the speed improvement effect is extremely large. Therefore, it can be considered that the larger the large-scale circuit in which the calculation time is a problem, the greater the effect of speeding up and the more practical.

FIG. 8 is a diagram showing instability factors at each sample point when frequency characteristics are calculated for “Case 9” shown in FIG. 3. Referring to FIG. 8, the instability factors of the conventional method and the proposed method are almost the same values until the calculation frequency exceeds about 5000 [Hz], but if the calculation frequency becomes higher than that, the instability factor in the proposed method. It can be seen that increases rapidly. Referring to FIG. 8, it can be confirmed that the instability factor of the proposed method decreases rapidly around the frequency of about 500 [kHz] and returns to the same value as the instability factor of the conventional method. This is because the instability factor of the proposed method exceeds the threshold ρ _th (10 ¹¹ ) around the calculation frequency of about 500 [kHz], and the ordering is performed again. As is apparent from FIG. 8, the number of calculations for ordering including the first calculation in all the calculation processes is two, and the ordering for 1998 is omitted from the calculation of all 2000 points (number of frequency sample points). .

(B) Verification of Calculation Accuracy and Calculation Speed FIG. 9 is a diagram showing a verification circuit used for verifying the calculation accuracy and calculation speed. The verification circuit shown in FIG. 9 has generators G1 to G3, transformers TR1 to TR3, a frequency converter station FC, and power transmission lines TL1 to TL3, and the generators G1 to G3 have transformers TR1 to TR3, respectively. And a frequency conversion station FC is a circuit of a 275 [kV] system connected to the power transmission network (a power transmission network including the transmission lines TL1 to TL3). This verification circuit is a circuit that assumes the calculation of the frequency characteristics of the AC system, which is performed in the analysis of the harmonic instability phenomenon and the harmonic expansion phenomenon associated with the AC / DC converter. In FIG. 9, for simplicity of illustration, a single-line diagram is drawn, but in actuality, all are configured by a three-phase model.

  For the frequency conversion station FC, the circuit configuration of a separately excited frequency conversion station (circuit configuration including a converter, an AC filter, and phase adjusting equipment) was used. The parameters of this frequency converter station FC are shown in Table 2 below.

The generators G1 to G3 were simulated by an equivalent circuit in which an initial transient impedance was connected in series to a voltage source representing an electromotive force of the generator, as shown in FIG. The initial transient impedance is expressed by a circuit composed of inductances L ₁ and L ₂ and a resistor R ₂ shown in FIG. 4B in consideration of the skin effect. Specifically, for the generator G1, the inductance _{L 1} and 0.283 [mH], the inductance _{L 2} and 0.566 [mH], and the resistor _{R 2} and 0.577 [Omega]. Also, the generator G2, the inductance _{L 1} and 0.566 [mH], the inductance _{L 2} and 1.13 [mH], and the resistor _{R 2} and 1.15 [Omega]. Also, the generator G3, the inductance _{L 1} and 54.5 [.mu.H], the inductance _{L 2} and 0.109 [mH], and the resistor _{R 2} and 0.111 [Omega]. The power supply voltage E was short-circuited because of harmonic analysis.

As shown in FIG. 5B, the transformers TR1 to TR3 are simulated by an equivalent circuit in which a leakage impedance is connected in series to an ideal transformer of Δ-Y connection. The leakage impedance of the transformers TR1 to TR3 is expressed by a circuit including the inductances L ₁ and L ₂ and the resistor R ₂ shown in FIG. Specifically, for the transformer TR1, the inductance _{L 1} and 24.2 [mH], the inductance _{L 2} and 4.78 [mH], and the resistor _{R 2} and 12.7 [Omega]. Also, the transformer TR2, an inductance _{L 1} and 46.5 [mH], the inductance _{L 2} and 11.5 [mH], and the resistor _{R 2} and 29.6 [Omega]. Also, the transformer TR3, the inductance _{L 1} and 4.84 [mH], the inductance _{L 2} and 0.955 [mH], and the resistor _{R 2} and 2.55 [Omega]. The transformation ratio was 275 [kV]: 22 [kV], and the floating capacitance to ground of the transformer winding was simulated with a capacitance of 2000 [pF] for all transformers.

  The transmission lines TL1 to TL3 are all overhead lines, and the mounting poles are assumed to be 275 [kV] standard mounting poles, and the ground resistivity is set to 100 [Ωm], and simulated by a frequency-dependent line model. The lengths of the transmission lines TL1 to TL3 are 5 [km], 15 [km], and 100 [km], respectively, and the sample of the line constant calculation frequency when creating the frequency-dependent line model is 0.1 [Hz] to 400 points were set within a range of 10 [MHz].

With respect to the above verification circuit, the frequency characteristic of the impedance of the AC system estimated from the point P in FIG. 9 was calculated in increments of 1 [Hz] from 1 [Hz] to 2000 [Hz]. For the calculation, the desktop personal computer used in the above-described verification (verification for the verification circuit shown in FIG. 3) was used, and the threshold ρ _th was set to 10 ¹¹ .

  FIG. 10 is a diagram illustrating a calculation result of frequency characteristic calculation performed on the verification circuit of FIG. Referring to FIG. 10, it can be confirmed that the impedance decreases at the frequency of the resonance order of each frequency filter. In FIG. 10, both the calculation result of the conventional method and the calculation result of the proposed method are drawn, but the two calculation results are in good agreement, and the difference cannot be determined from FIG. 10. The calculation time was about 156 seconds for the conventional method and about 8 seconds for the proposed method. This is about 20 times faster than the conventional method, and it was confirmed that the proposed method can calculate faster than the conventional method.

FIG. 11 is a diagram showing instability factors at each sample point when frequency characteristics are calculated for the verification circuit of FIG. Referring to FIG. 11, there was almost no increase in instability factors at each sample point, and the transition was also consistent at the waveform level. The size of the instability factor is 10 approximately ⁵ at the maximum, ordering is not repeated in the proposed method, it became the most computational small case. Thus, in this verification, it became clear that the result of ordering can be reused, and thereby the calculation speed can be increased.

  The frequency characteristic analyzing apparatus, method, and program according to an embodiment of the present invention have been described above. However, the present invention is not limited to the above-described embodiment, and can be freely changed within the scope of the present invention. . For example, in the above-described embodiment, an example in which the frequency characteristic analysis device 1 is realized by a desktop personal computer has been described. However, the frequency characteristic analysis device may be realized in the form of a server connected to a network. In such a form, for example, the analysis target and the analysis condition are transmitted to the server from the terminal device operated by the user, and the analysis result of the analysis performed on the server is returned to the terminal device.

  In the above-described embodiment, the frequency characteristic analysis apparatus 1 has been described as a dedicated apparatus for performing frequency characteristic calculation. However, the frequency characteristic analysis apparatus 1 includes functions other than the frequency characteristic calculation function (for example, instantaneous value analysis). Function) may be provided. Further, in the above-described embodiment, the case where the frequency characteristic calculation of the power system is performed has been described as an example. However, the present invention analyzes in other fields that require sparse matrix calculation (for example, structural analysis, fluid analysis, etc.) It is also possible to apply to. The present invention can be applied to general applications for speeding up sparse matrix operations on sparse matrices whose contents are changed when the calculation frequency is changed.

DESCRIPTION OF SYMBOLS 1 Frequency characteristic analyzer 14a Circuit equation preparation part 14b Ordering part 14c LU decomposition | disassembly part 14d Calculation precision evaluation part 14e Solving calculation part ρ Instability factor ρ _th threshold

Claims (8)

A circuit equation creation unit for creating a circuit equation to be analyzed according to a calculation frequency using a sparse matrix;
Among the circuit equations created by the circuit equation creation unit, an ordering unit that performs ordering by complete pivoting only for the sparse matrix of the circuit equation created in a specific case,
An LU decomposition unit that performs LU decomposition of the sparse matrix using an ordering result most recently performed by the ordering unit;
A solution calculation unit that performs a solution calculation using the matrix that has been LU decomposed by the LU decomposition unit;
A frequency characteristic analyzer comprising:   The frequency characteristic analysis apparatus according to claim 1, wherein the ordering unit performs ordering by the complete pivoting when a circuit equation is first created by the circuit equation creation unit. A calculation accuracy evaluation unit for evaluating calculation accuracy using the matrix decomposed by the LU decomposition unit;
The frequency characteristic analysis apparatus according to claim 1, wherein the ordering unit performs the ordering by the complete pivoting when it is evaluated that the calculation accuracy is deteriorated by the calculation accuracy evaluation unit.   The calculation accuracy evaluation unit obtains an unstable factor indicating an index of calculation accuracy using the matrix decomposed by the LU decomposition unit, and when the unstable factor exceeds a predetermined threshold, the calculation accuracy The frequency characteristic analysis apparatus according to claim 3, wherein the frequency characteristic analysis apparatus evaluates that the deterioration has occurred. A first step of creating a circuit equation to be analyzed according to a calculation frequency using a sparse matrix;
A second step of performing ordering by complete pivoting only on the sparse matrix of the circuit equation created in a specific case among the circuit equations created in the first step;
A third step of performing LU decomposition of the sparse matrix using the ordering result most recently performed in the second step;
A fourth step of performing a solution calculation using the matrix subjected to LU decomposition in the third step;
A fifth step of repeating the first step to the fourth step by changing the calculation frequency;
A frequency characteristic analysis method comprising: A sixth step of evaluating the calculation accuracy using the matrix subjected to LU decomposition in the third step;
The second step performs the ordering by the complete pivoting when the circuit equation is first created in the first step or when it is evaluated that the calculation accuracy is deteriorated in the sixth step. 5. The frequency characteristic analysis method according to 5. Computer
Circuit equation creation means for creating a circuit equation to be analyzed according to the calculation frequency using a sparse matrix;
Ordering means for performing ordering by complete pivoting only for the sparse matrix of the circuit equation created in a specific case among the circuit equations created by the circuit equation creation means;
LU decomposition means for performing LU decomposition of the sparse matrix using the ordering result most recently performed by the ordering means;
Solution calculation means for performing solution calculation using the matrix subjected to LU decomposition by the LU decomposition means;
Frequency characteristic analysis program to function. Further causing the computer to function as calculation accuracy evaluation means for evaluating calculation accuracy using the matrix subjected to LU decomposition by the LU decomposition means;
The ordering unit performs the ordering by the complete pivoting when the circuit equation is first created by the circuit equation creation unit or when it is evaluated that the calculation accuracy is deteriorated by the calculation accuracy evaluation unit. Item 8. The frequency characteristic analysis program according to Item 7.

Images