JP6650333B2 - Frequency characteristic analyzer, method, and program - Google Patents

Description

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

  Instantaneous value analysis that can analyze the time change of power and current at a waveform level is effective for analysis of overvoltage and overcurrent generated in a power system, analysis of power quality, analysis of dynamic characteristics of an AC / DC converter, and the like. 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 a specific bus to the system. It is used for the analysis of transient overvoltage including harmonic instability and harmonic characteristics generated by system interaction.

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

The formulation process is a process of creating a circuit equation from information on all nodes and all circuit components in a circuit. For example, assuming that there are N _n independent nodes and N _b circuit components (branches) in the circuit, a vector u consisting of the voltages of N _n nodes in the circuit flows through N _b circuit components. 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 a spar stab blow equation, and F is referred to as a tableau matrix. F, x, and y in the above equation (1) are represented by the following equation (2).

Here, A is N _b × N _n branch pair node connection matrix that associates a connection state of N _n independent nodes and N _b number of branches in the circuit. I is a unit matrix of _{N b} × _{N b.} B _{= [B} 1 -B _2] is a branch characteristic matrix of _{N b} × 2N _b associating the branch voltage v and branch currents i of each branch. s is a vector determined by the values of the voltage and current sources in the circuit. When _{N = N n × 2N b,} x, y is the number becomes N vectors element, F is a matrix of N × N. In the AC steady-state solution calculation, all quantities of the superstab blow equation shown here are complex numbers.

  The branch characteristic matrix B is obtained from a 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.

Here, i _L and i _C are currents flowing through L and C (branch currents), and v _L and v _C are voltages across L and C (branch voltages). If L is the b-th branch, the b-th row of B and s is represented by the following equation (5).

However, the (5) of the left side of equation "|" is left of the indicia indicate the b-th row of the matrix B _1, "|" from right mark indicates b-th row of the matrix -B _2, the right-hand side b rows of s Show the eyes. Further, the (5) two arrows of Formula bottom, 1 and -1 / j.omega.L indicates that each B ₁ and b-th column of B _2, the elements omitted in the "..." is 0. When a non-linear element is arranged in the analysis circuit, the non-linear element is ignored or linearized with characteristics near the origin. Further, the voltage source in the circuit, the vector s determined by the value of the current source, when the voltage waveform voltage source produces the _{E m cos (ωt + θ)} , it is substituted for E _m e ^j.theta. The corresponding element of s. The current source for generating the _{I m cos (ωt + θ)} , substituting I m ^{e jθ.}

In the solution process, processes called ordering, LU factorization, and forward / backward substitution are sequentially executed. Ordering is an operation of rearranging matrix elements for the purpose of improving calculation accuracy. For example, the matrix elements are rearranged according to the following procedure.
(A) For the matrix F, the number NR of non-zero elements existing in each row and the number NC of non-zero elements existing in each column are obtained.
(B) NRC = (NR-1) × (NC-1) is calculated for each element, and the element having the smallest NRC is replaced with 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 (sub-matrix) excluding one row and one column is considered as F, and the procedures (a) to (c) are repeated until the last row.

  Here, the procedure (b) is a process of selecting an element having the smallest NRC (an element having a large absolute value when the NRC is the same) and making it the diagonal element. In the search, two methods are conceivable: a case where the search target is only the element of the i-th row, and a case where the search target is all elements except the row and the column for which the diagonal element has already been determined. The former is called partial pivoting, and the latter is called complete 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 more optimal diagonal elements can be selected.

The LU factorization and forward / backward substitution decompose a coefficient matrix of a circuit equation (tableau matrix F in the above equation (1)) into a lower triangular matrix L and an upper triangular matrix U, and perform operations called forward substitution and backward substitution. This is an operation for obtaining x that becomes an unknown number. Now, _let the matrix elements of the above equation (1) be a _ij (where the first suffix is the row number, the second suffix is the column number, i = 1, 2,..., N, j = 1, 2,. , N), the element l _ij of the lower triangular matrix L and the element u _ij of the upper triangular matrix U after LU decomposition of F can be calculated by the following equations (6) and (7), respectively. .

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

Rikimichi Yonezawa, Taku Noda, "Development of Power System Instantaneous Value Analysis Program (Part 5)-Speed-up of XTAP by Improvement of Large-Scale Sparse Matrix LU Decomposition Algorithm-", Report of Central Research Institute of Electric Power Industry, Research Report H11009, 2012 Rikimichi Yonezawa, Taku Noda, "Development of Power System Instantaneous Value Analysis Program (Part 7)-Speeding up of XTAP by Parallelizing Internal State Update Processing of Circuit Parts-", Central Research Institute of Electric Power Industry Research Report H13005, 2014

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

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

  Here, as described above, in the ordering in the solution process, partial pivoting or complete pivoting can be used. In order to reduce the calculation time of the solution process in the frequency characteristic calculation, it is considered that partial pivoting, which requires less calculation time than complete pivoting, may be used. Although the calculation accuracy of partial pivoting is inferior to perfect pivoting in principle, the decrease in analysis accuracy is not a problem in many cases.However, there have been reports of cases where the calculation accuracy is problematic with partial pivoting, and it is not necessarily a universal method. Not really. Therefore, it is necessary to shorten the calculation time of the solution process while using complete pivoting from the viewpoint of maintaining the calculation accuracy.

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

In order to solve the above problem, a frequency characteristic analysis apparatus (1) of the present invention includes a circuit equation creation unit (14a) that creates a circuit equation to be analyzed according to a calculation frequency using a sparse matrix, Among the circuit equations created by the equation creating section, an ordering section (14b) for performing ordering by perfect pivoting only on the sparse matrix of the circuit equation created in a specific case, An LU decomposition unit (14c) for performing LU decomposition of the sparse matrix by using the ordered result, and a solution calculation unit (14e) for performing solution calculation using the matrix that has been LU-decomposed by the LU decomposition unit.
Further, in the frequency characteristic analysis apparatus according to the present invention, when the circuit equation is first created by the circuit equation creation unit, the ordering unit performs the ordering by the complete pivoting.
Further, the frequency characteristic analysis device of the present invention further includes a calculation accuracy evaluation unit (14d) for evaluating calculation accuracy using the matrix that has been LU-decomposed by the LU decomposition unit, and the ordering unit performs 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 analysis apparatus according to the present invention, the calculation accuracy evaluation unit obtains an instability factor (ρ) indicating an index of the calculation accuracy using the matrix decomposed by the LU decomposition unit, Exceeds a predetermined threshold (ρ _th ), it is evaluated that the calculation accuracy has deteriorated.
The frequency characteristic analysis method of the present invention includes a first step (S30) of creating a circuit equation to be analyzed according to a calculation frequency using a sparse matrix, and a method of specifying a circuit equation among the circuit equations created in the first step. (2) a second step (S41) of performing ordering by complete pivoting only on the sparse matrix of the circuit equation created in the case of (2), and performing LU of the sparse matrix using the ordering result most recently performed in the second step A third step (S42) of performing decomposition, a fourth step (S44) of performing a solution operation using the matrix subjected to LU decomposition in the third step, and a fourth step (S44) of changing the calculation frequency. And a fifth step (S20) of repeating.
In addition, the frequency characteristic analysis method of the present invention further includes a sixth step (S43) of evaluating calculation accuracy using the matrix that has been LU-decomposed in the third step, and the second step includes: When a circuit equation is first created in one step, or when it is evaluated that the calculation accuracy has deteriorated in the sixth step, the ordering by the complete pivoting is performed.
A frequency characteristic analysis program according to the present invention provides 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. An ordering means (14b) for performing ordering by complete pivoting only on the sparse matrix of a circuit equation created in a specific case among the equations; and LU decomposition means (14c) for performing LU decomposition of the matrix and solution calculation means (14e) for performing a solution operation using the matrix that has been LU-decomposed by the LU decomposition means function.
Further, the frequency characteristic analysis program of the present invention causes the computer to further function as calculation accuracy evaluation means (14d) for evaluating calculation accuracy using a matrix that has been LU-decomposed by the LU decomposition means. When the circuit equations are first created by the circuit equation creating means, or when the calculation accuracy evaluation means evaluates that the calculation accuracy has deteriorated, the ordering by the complete pivoting is performed.

  According to the present invention, a circuit equation to be analyzed according to a calculation frequency is created using a sparse matrix, and only the sparse matrix of the circuit equation created in a specific case is completely included in the created circuit equations. Ordering is performed by pivoting, LU decomposition of a sparse matrix is performed using the result of the most recently performed ordering, and a solution operation is performed using the matrix that has been LU-decomposed. Thereby, since ordering by perfect pivoting in cases other than the specific case can be omitted, it is possible to greatly reduce the time required for frequency characteristic calculation while maintaining calculation accuracy.

FIG. 1 is a block diagram illustrating a main configuration of a frequency characteristic analyzer according to an embodiment of the present invention. 5 is a flowchart illustrating an example of a frequency characteristic analysis method according to an embodiment of the present invention. FIG. 3 is a diagram illustrating a verification circuit used to verify a relationship between a circuit scale and a calculation speed. FIG. 4 is a diagram illustrating an equivalent circuit used for simulating the generators G2 to G10 in FIG. 3. FIG. 4 is a diagram illustrating an equivalent circuit used for simulating the transformers TR2 to TR10 in FIG. 3. FIG. 4 is a diagram illustrating calculation time required when frequency characteristic calculations are performed on “Case 1” to “Case 9” illustrated in FIG. 3. FIG. 4 is a diagram illustrating a relationship between a circuit scale and a speed improvement multiple. FIG. 4 is a diagram illustrating an instability factor at each sample point when frequency characteristic calculation is performed for “Case 9” illustrated in FIG. 3. FIG. 4 is a diagram illustrating a verification circuit used for verifying calculation accuracy and calculation speed. FIG. 10 is a diagram illustrating a calculation result of frequency characteristic calculation performed on the verification circuit of FIG. 9. FIG. 10 is a diagram illustrating an instability factor at each sample point when frequency characteristic calculation is performed on the verification circuit of FIG. 9. 9 is a flowchart showing a calculation process of a conventional frequency characteristic calculation.

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

<Frequency characteristic analyzer>
FIG. 1 is a block diagram showing a main configuration of a frequency characteristic analyzer according to an embodiment of the present invention. As shown in FIG. 1, the frequency characteristic analyzer 1 of the present embodiment includes an input device 11, an output device 12, an auxiliary storage device 13, and an arithmetic processing device 14, and calculates a frequency characteristic of a power system to be analyzed. I do. Such a frequency characteristic analyzer 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 corresponding to an operation of a user using the frequency characteristic analysis device 1 to the arithmetic processing device 14. The output device 12 includes, for example, a liquid crystal display device and a printer, and outputs various information output from the arithmetic processing device 14. The auxiliary storage device 13 includes, for example, a hard disk drive (HDD), a solid state drive (SSD), and the like, and various types of information (an ordering result OR and an evaluation parameter PR, which will be described later) used in analyzing a power system. ) And a frequency characteristic analysis program PG are stored.

  The arithmetic processing device 14 uses the various information stored in the auxiliary storage device 13 in accordance with 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 Perform calculations. The arithmetic processing device 14 includes a circuit equation creating unit 14a (circuit equation creating 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 It includes a solution calculation unit 14e (solution calculation means).

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 creation unit 14a creates the sparse stave blow equation shown in the above-described equation (1) according to the calculation frequency f. Note that the sparse matrix is a tableau matrix F shown in the above-described equation (2). In the case where the configuration of the power system to be analyzed does not change, when 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 orders the tableau matrix F of the circuit equation created by the circuit equation creating unit 14a by complete pivoting. 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 creating 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 creation unit 14a first creates a circuit equation or when the calculation accuracy evaluation unit 14d evaluates that calculation accuracy has deteriorated. Do. The result of the ordering performed by the 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 14c performs the 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 the 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 that has been LU-decomposed by the LU decomposition unit 14c. Specifically, the calculation accuracy evaluation unit 14d obtains an instability factor ρ indicating an index of the calculation accuracy shown in the following equation (8), and when the instability factor ρ is equal to or _larger than a predetermined threshold ρth. Then, 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 described as lu _ij . Further, the above-mentioned threshold value ρ _th is stored in the auxiliary storage device 13 as the evaluation parameter PR.

As is clear from the above equation (8), the instability factor ρ is the increase ratio of the infinite norm of the matrix element before and after LU decomposition. Normally, the scaling (or normalization, equalization) by pretreatment called, the denominator of the equation (8) to become a value close to 1, the unstable factor [rho, roughly equal to the value of the maximum element of the lu _ij Become. In the 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 equation (7) becomes a small value, so 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 instability factor ρ.

The solution calculation unit 14e performs a solution calculation using the matrix (lower triangular matrix L, upper triangular matrix U) that has been LU-decomposed 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, for example, a program recorded on a computer-readable recording medium such as a CD-ROM or a DVD (registered trademark) -ROM, or a program 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 creating unit 14a, the ordering unit 14b, the LU decomposition unit 14c, a calculation accuracy evaluation unit 14d, and a solution calculation unit 14e) are implemented as software. That is, these functions are realized by cooperation between software and hardware resources.

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

Next, the user inputs analysis conditions. For example, the user who operates the input device 11 of the frequency characteristic analysis apparatus 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 a logarithmic axis). Is input.

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

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). Note that this calculation frequency _{f 1} is the calculated starting frequency _{f the begin,} for example entered by a user. As a result, the first sample point is specified.

  Next, a formulation process (creation of a circuit equation of the power system) 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 super stab blow equation according to the calculation frequency f updated (set) in step S20 (see (1) above) , (2) is performed by the circuit equation creating unit 14a. Upon completion of the formulation process, a solution process is performed (step S40).

  When the solution process is started, first, the tableau matrix F of the circuit equation created by the circuit equation creating unit 14a is ordered by the complete pivoting by the ordering unit 14b (step S41: second step). The result of the ordering performed by the complete pivoting performed by the ordering unit 14b is stored in the auxiliary storage device 13 as the ordering result OR. Next, a process of performing LU decomposition of the tableau matrix F using the ordering result OR stored in the auxiliary storage device 13 is performed by the LU decomposition unit 14c (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 the calculation accuracy using the matrix that has been LU-decomposed by the LU decomposition unit 14c is performed by the calculation accuracy evaluation unit 14d (step S43: sixth step). Specifically, by using a matrix which is LU decomposed by LU decomposition unit 14c, obtains the [rho unstable factors indicating an index of a calculation accuracy shown in the equation (8) above, the instability factor [rho is the threshold [rho _th The calculation accuracy evaluation unit 14d performs a process of determining whether the value is also small. The above-mentioned threshold value ρ _th is stored in the auxiliary storage device 13 as the evaluation parameter PR.

If it is determined that the instability factor ρ is smaller than the threshold value ρ _th (if the determination result in step S43 is “yes”), a solution calculation using the matrix that has been LU-decomposed by the LU decomposition unit 14c Is performed by the solution calculation unit 14e (step S44: fourth step). Specifically, by performing 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 calculation for finding x shown in the above-described equations (1) and (2) is a solution calculation. 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 number of sample points K (step S50). When it is determined that the value of the variable k is equal to or smaller than the number of sample points K (when the determination result in step S50 is “no”), a 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, a process of creating a circuit equation according to the updated calculation frequency f is performed by the circuit equation creation unit 14a (Step S30). After that, the LU decomposition of the tableau matrix F is performed by the LU decomposition unit 14c using the ordering result OR stored in the auxiliary storage device 13 without performing the ordering by the ordering unit 14b (step S42). In other words, since the processing (calculation) performed by the LU decomposition unit 14c corresponds to the second and subsequent calculations, the LU decomposition of the tableau matrix F is performed using the latest ordering result 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). Then, when it is determined that the instability factor ρ is smaller than the threshold value ρ _th (when the determination result of step S43 is “yes”), the solution calculation (step S44) and the value of the variable k are performed. The comparison with the number of sample points K (step S50) is performed in order. If it is determined that the value of the variable k is equal to or smaller than the number K of sample points (if the determination result in step S50 is “no”), the process of updating the calculation frequency f after incrementing the value of the variable k is an arithmetic process. It is performed again by the device 14 (step S20).

As described above, the instability factor ρ is smaller than the threshold value ρ _th (the determination result in step S43 is “yes”), and the value of the variable k is equal to or smaller than the number of sample points K (the determination result in step S50 is If “no”, the ordering unit 14b does not perform the ordering (step S41), updates the calculation frequency f (step S20), creates a circuit equation according to the calculation frequency f (step S30), The LU decomposition (step S42) and the solution calculation (step S44) are repeatedly performed.

If the instability factor ρ obtained in step S43 becomes _{equal to} or _larger than the threshold value ρth while the above processing is being performed, the determination result in step S43 becomes “no”. Then, ordering by complete pivoting is performed by the ordering unit 14b (step S41), and the result of the ordering 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 _larger than the threshold value ρth, the ordering by the ordering unit 14b is performed, and LU decomposition is performed using the newly obtained ordering result OR. .

  When the value of the variable k is equal to or smaller 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 at each sample point are sequentially changed. Desired. Note that the series of processes illustrated in FIG. 2 ends when the value of the variable k is determined to be equal to or smaller than the number K of sample points (when the determination result in step S50 is “no”).

  As described above, in the present embodiment, when the difference between the frequency at which the calculation is completed and the frequency to be calculated next is small, the change in the circuit equation obtained at each frequency is also small, and the ordering unit 14 b The ordering is not performed for each sample point, and the result of the ordering performed in the past is reused. For this reason, the time required for frequency characteristic calculation can be significantly reduced.

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 needed. Specifically, when it is evaluated that the calculation accuracy has deteriorated (when the instability factor ρ becomes equal to or _larger than the threshold ρth), the ordering unit 14b performs the ordering. Therefore, the calculation accuracy can be maintained.

<Verification of calculation speed and calculation accuracy>
The inventor of the present application actually performed the frequency characteristic analysis of the power system using the above-described frequency characteristic analysis device 1, and verified the calculation speed and the calculation accuracy. In the frequency characteristic calculation of the present embodiment, when the calculation accuracy does not matter, the ordering is omitted to improve the efficiency of the calculation. In general, however, the calculation amount of the matrix operation such as the ordering affects the matrix order N. Receive. Therefore, in the following, the relationship between the matrix order N and the calculation speed will be verified first, and then the calculation accuracy and calculation speed will be verified while comparing with the conventional example. Incidentally, 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) Therefore, no problem be considered as circuit scale. Therefore, hereinafter, N is referred to as a circuit scale.

(A) Relationship between circuit scale and calculation speed Since the amount of calculation in ordering by perfect pivoting is considered to be affected by circuit scale N, when circuit scale N changes, the speed improvement factor (computation speed compared to the conventional method) The index indicating how many times the number has improved) is also considered to change. Therefore, for the purpose of confirming the relationship between the circuit scale N and the speed improvement rate, the speed improvement multiple when the circuit scale N was changed was calculated.

  FIG. 3 is a diagram illustrating a verification circuit used to verify the relationship between the circuit scale and the calculation speed. The circuit for verification 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 respectively connected to the transmission network via the transformers TR2 to TR10. This is a circuit connected to a (transmission network composed of transmission lines TL1 to TL16). The circuit scale of this verification circuit was changed in nine cases, “Case 1” to “Case 9” shown in the figure, and the frequency characteristics were calculated for each case.

  For example, “Case 1” is a circuit of a scale including 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 of a scale composed of electric wires TL1, TL2, TL3, and TL11. “Case 9” is a circuit of a scale including all of the generators G2 to G10, the transformers TR2 to TR10, and the transmission lines TL1 to TL16 shown in FIG. Table 1 below shows the number of nodes and the number of branches in each of "Case 1" to "Case 9".

  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 the simulation of the generators G2 to G10 is a circuit in which an inductance L corresponding to the next transient reactance is connected in series to an 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 is considered that an equivalent circuit shown in FIG. 4B can be simulated. Here, since the main purpose is to evaluate the calculation speed, the next transient impedance is assumed to have no loss component, and is simply simulated by the equivalent circuit shown in FIG. The inductance L of the generators G2 to G7 and G9 was 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. As with the generator, the leakage impedance of the transformer has a characteristic that the loss increases due to the skin effect. For example, it is considered that the equivalent circuit shown in FIG. 5B can be simulated. 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 transformer ratio was 500 [kV]: 22 [kV], and the stray capacitance to the ground of the transformer winding was simulated with a capacitance of 2000 [pF] for each transformer. The transmission lines TL1 to TL16 should originally be simulated realistically using a frequency-dependent distributed constant line model. However, in this study, since the main purpose is to evaluate the calculation speed, the transmission lines TL1 to TL16 are simplified. A π-type equivalent circuit model (1 stage) having no interphase coupling was used.

  In all cases of “Case 1” to “Case 9” shown in FIG. 3, the frequency characteristics of the impedance were calculated in view of the system from the bus A in the figure. Here, the wider the calculation frequency is, the larger the difference of the circuit equation between the frequency samples becomes, which is a disadvantageous condition for the frequency characteristic calculation in the present embodiment. Therefore, in the verification here, the effect of the frequency characteristic calculation in the present embodiment was confirmed by setting the frequency range from 0.1 [Hz] to 1 [MHz], which is wider than the frequency range usually performed. The frequency sampling points were 2,000 at equal intervals on the logarithmic axis. For the calculation, a desktop personal computer that is currently widely used was used.

Here, in numerical calculations using a general computer, it is standard to use double precision arithmetic, and the number of significant digits is about 16 digits. In other words, when the above-mentioned instability factor ρ exceeds 16 digits, information drop occurs due to 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 reduced. Therefore, in this calculation, the threshold ρ _{th is set} to 10 ^{11 in} consideration of a margin of five digits. Note that “Case 1” to “Case 9” used here are cases created for verification, and that the frequency characteristics themselves obtained as calculation results have no special meaning.

  FIG. 6 is a diagram illustrating the calculation time required when the frequency characteristics are calculated for “Case 1” to “Case 9” illustrated 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 characteristic calculation is performed. 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 illustrating the relationship between the circuit scale and the speed improvement multiple. FIG. 7 is a plot of the ratio of the calculation time (the value obtained by dividing the calculation time of the conventional method by the calculation time of the proposed method) for each case shown in FIG. 6 as a speed improvement multiple. The circuit scale in each case shown in FIG. 6 is obtained from Table 1 described above. Referring to FIG. 7, it can be seen that the speed-up factor increases as the circuit size increases. In particular, in the case of “Case 9” shown in FIG. 6 (when the circuit scale is 1353), the speed improvement multiple is nearly 60 times, and it can be seen that the effect of the speed improvement is extremely large. Therefore, it is considered that the larger the circuit, the more the calculation time becomes a problem, the greater the effect of speeding up and the more practical the circuit.

FIG. 8 is a diagram showing an instability factor at each sample point when frequency characteristic calculation is performed for “Case 9” shown in FIG. Referring to FIG. 8, the instability factors of the conventional method and the proposed method have substantially the same value until the calculation frequency exceeds about 5000 [Hz]. It can be seen that increases rapidly. Referring to FIG. 8, it can be confirmed that the instability factor of the proposed method rapidly decreases around 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 value ρ _th (10 ¹¹ ) when the calculation frequency is about 500 [kHz], and the ordering is performed again. As is clear from FIG. 8, the number of times of ordering including the first calculation in the entire calculation process is two, and the ordering of 1998 times out of the calculation of all 2,000 points (the number of frequency sample points) is omitted. .

(B) Verification of Calculation Accuracy and Calculation Speed FIG. 9 is a diagram illustrating a verification circuit used for verifying calculation accuracy and calculation speed. The verification circuit illustrated in FIG. 9 includes generators G1 to G3, transformers TR1 to TR3, a frequency converter FC, and transmission lines TL1 to TL3, and the generators G1 to G3 respectively operate the transformers TR1 to TR3. This is a circuit of a 275 [kV] system that is connected to a power transmission network (a power transmission network including transmission lines TL1 to TL3) via a power supply and a frequency conversion station FC is connected to the power transmission network. This verification circuit is a circuit assuming calculation of frequency characteristics of an AC system, which is performed in analysis of a harmonic instability phenomenon and a harmonic expansion phenomenon accompanying an AC / DC converter. In FIG. 9, for simplicity of illustration, it is drawn in a single-line connection diagram, but all are actually configured by three-phase models.

  For the frequency converter FC, the circuit configuration of a separately-excited frequency converter (a circuit configuration including a converter, an AC filter, and a phase adjustment facility) was used. Table 2 below shows the parameters of the frequency converter FC.

Each of the generators G1 to G3 was 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 was represented by a circuit including the inductances L ₁ and L ₂ and the resistor R ₂ shown in FIG. 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 the analysis of harmonics.

As shown in FIG. 5B, the transformers TR1 to TR3 were simulated by an equivalent circuit in which a leakage impedance was connected in series to an ideal transformer connected in Δ-Y. Leakage impedance of the transformer TR1~TR3 was expressed by a circuit consisting of the inductance _L 1, _{L 2} and resistor _{R 2} shown in Figure 5 (b). 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 transformer ratio was 275 [kV]: 22 [kV], and the stray capacitance to ground of the transformer winding was simulated with a capacitance of 2000 [pF] for each transformer.

  The transmission lines TL1 to TL3 were all overhead lines, and the columns were assumed to be 275 [kV] standard columns. The ground resistivity was set to 100 [Ωm] and simulated by a frequency-dependent line model. The transmission lines TL1 to TL3 have a length of 5 [km], 15 [km], and 100 [km], respectively. The sample of the line constant calculation frequency at the time of creating the frequency-dependent line model is 0.1 [Hz] to 400 points were set in the range of 10 [MHz].

With respect to the above-described verification circuit, the frequency characteristics of the impedance of the AC system viewed from the point P in FIG. 9 were calculated from 1 [Hz] to 2000 [Hz] in increments of 1 [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 showing 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 is reduced at the frequency of the resonance order of each frequency filter. Further, in FIG. 10, both the calculation result of the conventional method and the calculation result of the proposed method are drawn. However, both calculation results match well, 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 has been confirmed that the proposed method can calculate at a higher speed than the conventional method.

FIG. 11 is a diagram illustrating an instability factor at each sample point when frequency characteristic calculation is performed on the verification circuit of FIG. Referring to FIG. 11, there was almost no increase in the instability factor at each sample point, and the transition also coincided 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 was clarified that the result of the ordering can be reused, thereby increasing the calculation speed.

  The frequency characteristic analysis device, method, and program according to one 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 case, for example, the analysis target and the analysis conditions are transmitted to the server from the terminal device operated by the user, and the analysis result of the analysis performed by the server is returned to the terminal device.

  Further, in the above-described embodiment, the frequency characteristic analysis device 1 has been described as a dedicated device for performing the frequency characteristic calculation. However, the frequency characteristic analysis device 1 has 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 is applied to analysis in other fields requiring sparse matrix operation (eg, structural analysis, fluid analysis, and the like). It is also possible to apply to. INDUSTRIAL APPLICABILITY The present invention can be generally applied to applications for speeding up a sparse matrix operation for a sparse matrix whose contents are changed when the calculation frequency is changed.

DESCRIPTION OF SYMBOLS 1 Frequency characteristic analyzer 14a Circuit equation creation unit 14b Ordering unit 14c LU decomposition unit 14d Calculation accuracy evaluation unit 14e Solution calculation unit ρ Instability factor ρ _th threshold

Claims (8)

A circuit equation creation unit that creates a circuit equation to be analyzed according to the 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 on the sparse matrix of the circuit equation created in a specific case,
An LU factorization unit that performs an LU factorization of the sparse matrix using an ordering result most recently performed by the ordering unit;
A solution calculation operation unit that performs a solution operation using the matrix that has been LU-decomposed by the LU decomposition unit;
A frequency characteristic analysis device comprising:   The frequency characteristic analysis apparatus according to claim 1, wherein the ordering unit performs the ordering by the complete pivoting when the circuit equation is first created by the circuit equation creating unit. The image processing apparatus further includes a calculation accuracy evaluation unit that evaluates calculation accuracy using the matrix that has been LU-decomposed by the LU decomposition unit,
The ordering unit, if the calculation accuracy evaluation unit in deterioration of the calculation accuracy has been evaluated to have occurred, the full pivoting performs ordering by, the frequency characteristic analyzer according to claim 1 or claim 2, wherein.   The calculation accuracy evaluation unit obtains an instability factor indicating an index of calculation accuracy using the matrix that has been LU-decomposed by the LU decomposition unit, and calculates a calculation accuracy when the instability factor exceeds a predetermined threshold. The frequency characteristic analysis device according to claim 3, wherein the frequency characteristic analysis device evaluates that deterioration of the frequency has occurred. A first step in which a circuit equation creation unit creates a circuit equation to be analyzed according to the calculation frequency using a sparse matrix;
A second step in which the ordering unit performs ordering by complete pivoting only on the sparse matrix of the circuit equation created in a specific case among the circuit equations created by the circuit equation creating unit ;
A third step in which an LU factorization unit performs an LU factorization of the sparse matrix using the ordering result most recently performed by the ordering unit ;
A fourth step of performing a solution operation using the matrix that has been LU-decomposed by the LU factorization unit ,
A fifth step in which an arithmetic processing unit repeats the first to fourth steps while changing the calculation frequency;
A frequency characteristic analysis method having the following. The calculation accuracy evaluation unit further includes a sixth step of evaluating the calculation accuracy using the matrix that has been LU-decomposed by the LU decomposition unit ,
The ordering unit, when a circuit equation is first created in the circuit equation creation unit , or when it is evaluated that the calculation accuracy has deteriorated in the calculation accuracy evaluation unit , performs the ordering by the complete pivoting, claim. Item 6. The frequency characteristic analysis method according to Item 5. Computer
Circuit equation creating means for creating a circuit equation to be analyzed according to the calculation frequency using a sparse matrix,
Among the circuit equations created by the circuit equation creating means, ordering means for performing ordering by complete pivoting only on the sparse matrix of the circuit equation created in a specific case,
LU decomposition means for performing LU decomposition of the sparse matrix using an ordering result most recently performed by the ordering means;
Solution calculation means for performing a solution operation using the matrix that has been LU-decomposed by the LU decomposition means;
Frequency analysis program to function as a function. The computer further functions as calculation accuracy evaluation means for evaluating calculation accuracy using the matrix that has been LU-decomposed by the LU decomposition means,
The ordering means, when the circuit equation is first created by the circuit equation creation means, or when it is evaluated that the calculation accuracy has deteriorated by the calculation accuracy evaluation means, performs the ordering by the complete pivoting. Item 7. The frequency characteristic analysis program according to Item 7.

Images