CIS 890: Final report

Alvin Leung

December 16, 1994

Abstract

The goals of this project is to explore the internals of ETMSP 3.1 and to interface it with a parallel linear system solver. For the ETMSP 3.1 internals, we are specifically interested in the Jacobian matrix generation mechanism and the sparse linear system system solver interface. We have completed our investigation with a satisfactory result.

1 Jacobian matrix generation

To understand the Jacobian matrix generation for the EPRI ETMSP 3.1, we first explore circumstances of which the Jacobian matrix is needed. ETMSP 3.1 solves the system of differential-algebraic equations of the electrical power stability problem by decomposing the system into a ODE subsystem and a non-linear algebraic subsystem. The ODE system of equations is form by the model of the generators and the related equipment. The non-linear algebraic system is formed by the transmission network and the non-linear load. Then the program applies trapezoidal rule, a implicit integration method, to approximate the ODE system by a non-linear algebraic system. Finally, the program uses Newton method to solve the resulted algebraic systems.

The system of equations after applying the Newton method can be expressed in the form

\begin{aligned}x_{n+1}^{k+1} & = & x_{n+1}^{k} − A_{G}^{−1} [F_{n+1}^{k} + B_{G}(V_{n+1}^{k+1} − V_{n+1}^{k})] \\ (Y + Y_{G} − C_{G}A_{G}^{−1}B_{G})V_{n+1}^{k+1} & = & I_{n+1}^{k} + C_{G}A_{G}^{−1}F_{n+1}^{k} + (Y_{G} − C_{G}A_{G}^{−1}B_{G})V_{n+1}^{k} \end{aligned} ,

where

\begin{matrix}A_{G} = I − \frac{h}{2} \frac{∂f}{∂x} |_{x=x_{n+1}^{k},V =V_{n+1}^{k}} & B_{G} = − \frac{h}{2} \frac{∂f}{∂V} |_{x=x_{n+1}^{k},V =V_{n+1}^{k}} \\ C_{G} = − \frac{∂I}{∂x} |_{x=x_{n+1}^{k},V =V_{n+1}^{k}} & C_{G} = − \frac{∂I}{∂V} |_{x=x_{n+1}^{k},V =V_{n+1}^{k}} \end{matrix} .

In other words, A_{G} is the Jacobian matrix of the generator equations, B_{G} is the Jacobian matrix of the linkage between generators and the network, C_{G} is the Jacobian matrix of the linkage between the network and the generators and Y_{G} is the Jacobian matrix of the network equations.

These Jacobian matrices are generated in subroutine DGENLN as intermediate matrices. To compute the term A_{G}^{−1}B_{G}, DGENLN performs a static minimal degree ordering and a dense LU factorization on A_{G}. Then forward and back solve is used to compute U_{G}^{−1}L_{G}^{−1}B. Finally, a sparse matrix - matrix multiplication is done to obtain C_{G}A_{G}^{−1}B_{G}. In order to obtain V_{n+1}^{k+1}, a sparse LU factorization with forward and back solve is used to compute (Y + Y_{G} − C_{G}A_{G}^{−1}B_{G})^{−1}.

2 Sequential Linear Solver

The sequential linear solver for ETMSP 3.1 uses static minimum degree ordering with symbolic factorization. The subroutine ORDER performs the ordering and returns an array, IBUSO, which contains the bus number in optimal order. The linear solver package is reasonably well coded. However, the quality of maintenance is very poor. For example, computing inverse of the pivot values of the LU factorized matrix is done by a separate subroutine, instead of single line code change in the LU factorization subroutine.

The following time is obtained by running ETMSP 3.1 on RS 6000 workstation with 6000 buses data.

ETMSP 3.1 NPAC
Ordering  5.17 sec 5.15 sec
Factorization 1.80 sec 1.05 sec
Solve  0.01 sec 0.01 sec

In ETMSP 3.1 the LU linear system solver is divided into three parts. The first part handles the factorization in complex. The factorization chooses the pivot elements according to the provided optimal order. The factorization also set up the “data structure” for the next phase of the linear solver. The forward and back solve is performed within the Newton solver. As the timing data shows, the forward and back solve have a satisfactory performance.

The LU factorization, forward and back solve is performed by subroutine AFACTC, AUPDSC and ASOLSC, respectively. As shown in figure 1, the solving the sparse linear system is performed in many different phases.

Figure 1: The call graph for ETMSP 3.1 3

Figure 2: The call graph for interfacing parallel solver to ETMSP 3.1

According to the ETMSP 3.1 manual and code, the factorized matrix can be save and restore from the temporary file. This posts some inconvenience for interfacing a parallel real matrix linear solver with ETMSP 3.1.

3 Parallel Linear Solver

The parallel real matrix linear solver, which is currently interfaced with ETMSP 3.1, is from Toolbox by Tony Skjellum. The call graph of the parallel solver is shown in figure 2.

The initialization of the parallel computer is done only once during the first call to the parallel linear system solver. After the initialization, ETMSP 3.1 can repeatedly calling factorization, forward and back solve.

The factorized matrix data is retained in the parallel computer’s memory and will be use for the next call of the forward and back solve. As mentioned in section 2, ETMSP 3.1 can restore past factorized matrix from temporary file without notifying linear solver subroutines. As a result, a correct interface will require overhaul the ETMSP’s structure, which is beyond the scope of this project. Consequently, ETMSP with parallel linear solver will only work with limited type of test cases, such as, 9-buses, 450-buses, 1000-buses

0 1 2 3 4 5 6 7 8
0 ⊕  ⊕
1  ⊕  ⊕
2  ⊕  ⊕
3 ⊕  ⊕  ⊕  ⊕
4  ⊕  ⊕  ⊕ ⊕
5  ⊕  ⊕ ⊕  ⊖
6  ⊕  ⊕ ⊕ ⊕ ⊖
7  ⊕  ⊕ ⊕ ⊖
8  ⊕ ⊕ ⊖ ⊖ ⊖ ⊕

Figure 3: The non-zero pattern for the original matrix

and the 6000-buses system supplied by EPRI. The commonality of these cases is a simple switching pattern; hence, ETMSP can not reuse previously factorized matrix.

IN addition, the characteristics of the sequential and parallel linear solve differs significantly. The original ETMSP linear solver performs forward and back solve separately. However, the parallel forward and back solve is computed by the replacement subroutine for AUPDSC. As a result, the replacement call for ASOLSC consists of moving the result to the correct location.

In figure 3, the sparsity pattern of the matrix for a 9-bus system is shown. The reordered sparsity pattern of the matrix is shown in figure 4. In both figures, the symbol ⊕ indicates the non-zero of the matrix and the symbol ⊖ indicates the fill-ins of the matrix.

Although this example does not show any advantages on reordering, the reduction of fill-in by reordering the matrix is justified by other researchers.

In order to solve this complex matrix by using real matrix linear solver, the matrix have to be rewritten into a 2n × 2n real matrix, where n is the order of the complex matrix. The mathematical justification of the complex to real matrix is similar to the multiplication of two complex number.

\begin{aligned}x + ıy & = & (a + ıb) × (c + ıd) \\ & = & (ac − bd) + ı(bc + ad) \end{aligned}

Rewrite the multiplication in matrix form.

0 1 2 3 4 5 6 7 8
0 ⊕  ⊕ ⊕
1  ⊕  ⊕ ⊕
2  ⊕  ⊕
3  ⊕  ⊕  ⊕
4  ⊕  ⊕
5  ⊕  ⊕
6 ⊕  ⊕ ⊕  ⊕ ⊖ ⊖
7 ⊕ ⊕  ⊕ ⊕ ⊕ ⊖
8  ⊕  ⊕ ⊕  ⊖ ⊖ ⊕

Figure 4: The non-zero pattern for the ordered matrix

\begin{aligned}[ \begin{matrix}x \\ y \end{matrix} ] & = & (ac − bd) + ı(bc + ad) \\ & = & [ \begin{matrix}a & −b \\ b & a \end{matrix} ] [ \begin{matrix}c \\ d \end{matrix} ] \end{aligned}

However, using decomposition as show in the above mathematical equation will destroy the reordered sparsity pattern as shown in figure 5. As a result, large number of fill-in will generated.

The original ordering can be retained by applying the decomposition equation to each entry of the complex matrix. In other words, each entry of the original matrix is represented by a 2 × 2 matrix. In figure 6, the ⊕ symbols indicates the result of this decomposition technique on the original matrix.

An additional concert of decomposing a complex matrix into a real matrix is the destruction of diagonal dominant property. Even though the original matrix is diagonal dominant, the 2n × 2n real matrix does not have this property. In fact, the diagonal entries of the real matrix can be zero. As a result, pivoting strategy must be used to obtain best possible answer for solving the linear system. By using partial row or column pivoting, a suitable pivot can be located without destroying the ordering (figure 6).

4 Conclusion

The lessons learned by exploring ETMSP 3.1

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
0 ⊕  ⊕ ⊕  ⊕  ⊕ ⊕
1  ⊕  ⊕ ⊕  ⊕  ⊕ ⊕
2  ⊕  ⊕  ⊕  ⊕
3  ⊕  ⊕  ⊕  ⊕  ⊕  ⊕
4  ⊕  ⊕  ⊕  ⊕
5  ⊕  ⊕  ⊕  ⊕
6 ⊕  ⊕ ⊕  ⊕ ⊖ ⊖ ⊕  ⊕ ⊕  ⊕ ⊖ ⊖
7 ⊕ ⊕  ⊕ ⊕ ⊕ ⊖ ⊕ ⊕ ⊖ ⊖  ⊕ ⊕ ⊕ ⊖
8  ⊕  ⊕ ⊕  ⊖ ⊖ ⊕ ⊖ ⊕ ⊖ ⊕ ⊕ ⊖ ⊖ ⊖ ⊕
9 ⊕  ⊕ ⊕ ⊖ ⊕ ⊖ ⊖ ⊖ ⊖ ⊖ ⊕ ⊕ ⊖
10  ⊕  ⊕ ⊕ ⊖ ⊕ ⊖ ⊖ ⊖ ⊖ ⊖ ⊕ ⊕
11  ⊕  ⊕ ⊖ ⊖ ⊖ ⊖ ⊕ ⊖ ⊖ ⊖ ⊕ ⊖ ⊖
12  ⊕  ⊕ ⊖ ⊕ ⊖ ⊖ ⊖ ⊕ ⊖ ⊖ ⊕ ⊖ ⊕
13  ⊕  ⊕ ⊖ ⊖ ⊖ ⊖ ⊕ ⊖ ⊖ ⊖ ⊕
14  ⊕  ⊕ ⊖ ⊖ ⊖ ⊖ ⊖ ⊖ ⊕ ⊖ ⊕ ⊖
15 ⊕  ⊕ ⊕  ⊕ ⊖ ⊖ ⊕ ⊖ ⊕ ⊕ ⊖ ⊖ ⊕ ⊖ ⊖
16 ⊕ ⊕  ⊕ ⊕ ⊕ ⊖ ⊕ ⊕ ⊖ ⊖ ⊖ ⊕ ⊕ ⊕ ⊖
17  ⊕  ⊕ ⊕  ⊖ ⊖ ⊕ ⊖ ⊕ ⊖ ⊕ ⊕ ⊖ ⊖ ⊖ ⊕

Figure 5: The non-zero pattern for the ordered matrix with block decomposition

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
0 ⊕ ⊕  ⊕ ⊕ ⊕ ⊕
1 ⊕ ⊕  ⊕ ⊕ ⊕ ⊕
2  ⊕ ⊕  ⊕ ⊕ ⊕ ⊕
3  ⊕ ⊕  ⊕ ⊕ ⊕ ⊕
4  ⊕ ⊕  ⊕ ⊕
5  ⊕ ⊕  ⊕ ⊕
6  ⊕ ⊕  ⊕ ⊕  ⊕ ⊕
7  ⊕ ⊕  ⊕ ⊕  ⊕ ⊕
8  ⊕ ⊕  ⊕ ⊕
9  ⊕ ⊕  ⊕ ⊕
10  ⊕ ⊕  ⊕ ⊕
11  ⊕ ⊕  ⊕ ⊕
12 ⊕ ⊕  ⊕ ⊕ ⊕ ⊕  ⊕ ⊕ ⊖ ⊖ ⊖ ⊖
13 ⊕ ⊕  ⊕ ⊕ ⊕ ⊕  ⊕ ⊕ ⊖ ⊖ ⊖ ⊖
14 ⊕ ⊕ ⊕ ⊕  ⊕ ⊕ ⊖ ⊖ ⊕ ⊕ ⊖ ⊖
15 ⊕ ⊕ ⊕ ⊕  ⊕ ⊕ ⊖ ⊖ ⊕ ⊕ ⊖ ⊖
16  ⊕ ⊕  ⊕ ⊕ ⊕ ⊕  ⊖ ⊖ ⊖ ⊖ ⊕ ⊕
17  ⊕ ⊕  ⊕ ⊕ ⊕ ⊕  ⊖ ⊖ ⊖ ⊖ ⊕ ⊕

Figure 6: The non-zero pattern for the ordered matrix with decoupled real and imaginary part using partial row pivoting factorization.

• parallizing ETMSP requires major overhaul of the source code. In fact, sequential ETMSP requires major overhaul to eliminate dead code segments.

• ETMSP’s operation is highly optimized to run in sequential environment. Many time saving tricks will be useless in parallel computing environment.

Currently, ETMSP uses “Very Dishonest Newton Method” to solve the non-linear algebraic systems. If parallel linear solver performs as fast as expected, then ETMSP should consider using “Dishonest Newton Method.” So that the fast convergent properties of the Newton method can be exploited. In addition, multiple initial guesses can be used with the Newton method. Hopefully, one of the initial guesses will converge early than the others. This will also increase the efficiency of the parallel forward and back solve.