\documentstyle[11pt,psfig]{article}

\title{CIS 890: Final report}
\author{Alvin Leung}

\begin{document}
\maketitle

\begin{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.

\end{abstract}

\section{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{math}
\begin{array}{lcl}
x_{n + 1}^{k + 1} & = & x_{n + 1}^{k} - A_{G}^{-1} \left[ F^{k}_{n + 1}
+ B_{G}(V_{n + 1}^{k + 1} - V_{n + 1}^{k}) \right] \\
(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{array}
\end{math},

where

\begin{math}
\begin{array}{cc}
A_{G} = I - \left.\frac{h}{2} \frac{\partial f}{\partial x}
\right|_{x=x_{n + 1}^{k}, V=V_{n + 1}^{k}} & 
B_{G} = - \left.\frac{h}{2} \frac{\partial f}{\partial V}
\right|_{x=x_{n + 1}^{k}, V=V_{n + 1}^{k}} \\
C_{G} = - \left. \frac{\partial I}{\partial x} \right|_{x=x_{n + 1}^{k},
V=V_{n + 1}^{k}} &
C_{G} = - \left. \frac{\partial I}{\partial V} \right|_{x=x_{n + 1}^{k},
V=V_{n + 1}^{k}} \\
\end{array}
\end{math}.

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 {\it DGENLN} as intermediate
matrices.  To compute the term $A_{G}^{-1}B_{G}$, {\it 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}$.

\section{Sequential Linear Solver}
\label{seq_linear}

  The sequential linear solver for ETMSP 3.1 uses static minimum
  degree ordering with symbolic factorization.  The subroutine {\it ORDER}
  performs the ordering and returns an array, {\it 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.

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
               & ETMSP 3.1 & NPAC \\ \hline
 Ordering      & 5.17 sec  & 5.15 sec \\ \hline
 Factorization & 1.80 sec  & 1.05 sec \\ \hline
 Solve         & 0.01 sec  & 0.01 sec \\ \hline
\end{tabular}
\end{center}

  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.

\begin{figure}
\begin{center}
\begin{picture}(350,250)
\psfig{figure=etmsp_call_graph.ps,width=350pt,height=227pt,angle=270}
\end{picture}
\end{center}
\caption{The call graph for ETMSP 3.1 \label{ETMSP_call_graph}}
\end{figure}

  The LU factorization, forward and back solve is performed by subroutine
  {\it AFACTC}, {\it AUPDSC} and {\it ASOLSC}, respectively.  As shown in
  figure~\ref{ETMSP_call_graph}, the solving the sparse linear system is
  performed in many different phases.  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.
   
\section{Parallel Linear Solver}

  The parallel real matrix linear solver, which is currently interfaced
  with ETMSP 3.1, is from {\bf Toolbox} by Tony Skjellum.  The call graph
  of the parallel solver is shown in figure~\ref{ETMSP_par_solver}.

\begin{figure}
\begin{center}
\begin{picture}(217,198)
\psfig{figure=etmsp_pallel_linear_call_graph.eps,width=217pt,height=198pt}
\end{picture}
\end{center}
\caption{The call graph for interfacing parallel solver to ETMSP 3.1
\label{ETMSP_par_solver}}
\end{figure}

   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~\ref{seq_linear},
   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 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 {\it ASOLSC}
   consists of moving the result to the correct location.

   In figure~\ref{Original_A}, the sparsity pattern of the matrix for
   a 9-bus system is shown.  The reordered sparsity pattern
   of the matrix is shown in figure~\ref{Ordered_complex_A}.  In both
   figures, the symbol $\oplus$ indicates the non-zero of the matrix
   and the symbol $\ominus$ indicates the fill-ins of the matrix.

\begin{figure}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|} \hline
  &  0&  1&  2&  3&  4&  5&  6&  7&  8\\ \hline
 0&  $\oplus$&  &  &  $\oplus$&  &  &  &  &  \\ \hline
 1&  &  $\oplus$&  &  &  $\oplus$&  &  &  &  \\ \hline
 2&  &  &  $\oplus$&  &  &  &  $\oplus$&  &  \\ \hline
 3&  $\oplus$&  &  &  $\oplus$&  &  $\oplus$&  &  &  $\oplus$ \\ \hline
 4&  &  $\oplus$&  &  &  $\oplus$&  &  &  $\oplus$&  $\oplus$ \\ \hline
 5&  &  &  &  $\oplus$&  &  $\oplus$&  $\oplus$&  &  $\ominus$ \\ \hline
 6&  &  &  $\oplus$&  &  &  $\oplus$&  $\oplus$&  $\oplus$& $\ominus$ \\ \hline
 7&  &  &  &  &  $\oplus$&  &  $\oplus$&  $\oplus$& $\ominus$  \\ \hline
 8&  &  &  &  $\oplus$&  $\oplus$& $\ominus$ & $\ominus$ & $\ominus$ &  $\oplus$\\ \hline
\end{tabular}
\end{center}
\caption{The non-zero pattern for the original matrix
\label{Original_A}}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|} \hline
  &  0&  1&  2&  3&  4&  5&  6&  7&  8\\ \hline
 0&  $\oplus$&  &  &  &  &  &  $\oplus$&  $\oplus$&  \\ \hline
 1&  &  $\oplus$&  &  &  &  &  &  $\oplus$&  $\oplus$\\ \hline
 2&  &  &  $\oplus$&  &  &  &  $\oplus$&  &  \\ \hline
 3&  &  &  &  $\oplus$&  &  &  $\oplus$&  &  $\oplus$\\ \hline
 4&  &  &  &  &  $\oplus$&  &  &  &  $\oplus$\\ \hline
 5&  &  &  &  &  &  $\oplus$&  &  $\oplus$&  \\ \hline
 6&  $\oplus$&  &  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\ominus$&
$\ominus$\\ \hline
 7&  $\oplus$&  $\oplus$&  &  &  &  $\oplus$&  $\oplus$&  $\oplus$&
$\ominus$ \\ \hline
 8&  &  $\oplus$&  &  $\oplus$&  $\oplus$&  &  $\ominus$&  $\ominus$&
$\oplus$\\ \hline
\end{tabular}
\end{center}
\caption{The non-zero pattern for the ordered matrix
\label{Ordered_complex_A}}
\end{figure}

   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 \times 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{math}
\begin{array}{lcl}
x + \imath y & = & (a + \imath b) \times (c + \imath d) \\
& = & (ac - bd) + \imath (bc + ad) \\
\end{array}
\end{math}

   Rewrite the multiplication in matrix form.

\begin{math}
\begin{array}{lcl}
\left[
\begin{array}{c}
x \\
y \\
\end{array}
\right]
& = & (ac - bd) + \imath (bc + ad) \\
& = & \left[
\begin{array}{cc}
a & -b \\
b & a \\
\end{array}
\right]
\left[
\begin{array}{c}
c \\
d \\
\end{array}
\right]
\end{array}
\end{math}

   However, using decomposition as show in the above mathematical equation
   will destroy the reordered sparsity pattern as shown in
   figure~\ref{Ordered_real_A_block}. 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 \times 2$ matrix.  In
   figure~\ref{Ordered_A_real_LUed}, the $\oplus$ 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 \times 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~\ref{Ordered_A_real_LUed}).


\begin{figure}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline
  &    0     &  1&  2&  3&  4&  5&  6&  7&  8&  9& 10& 11& 12& 13& 14& 15& 16& 17\\ \hline
 0& $\oplus$ &  &  &  &  &  &  $\oplus$&  $\oplus$& & $\oplus$&  &  &  &  &  &
  $\oplus$&  $\oplus$& \\ \hline
 1&  &  $\oplus$&  &  &  &  &  &  $\oplus$&  $\oplus$ & & $\oplus$ &  &  &  &
  &  &  $\oplus$&  $\oplus$ \\ \hline
 2&  &  &  $\oplus$&  &  &  &  $\oplus$&  & & &  &  $\oplus$&  &  &  &
  $\oplus$& & \\ \hline
 3&  &  &  &  $\oplus$&  &  &  $\oplus$&  &  $\oplus$ & &  &  &  $\oplus$&  &
  &  $\oplus$&  &  $\oplus$ \\ \hline
 4&  &  &  &  &  $\oplus$&  &  &  &  $\oplus$ & &  &  &  &  $\oplus$&  &  &  &
  $\oplus$ \\ \hline
 5&  &  &  &  &  &  $\oplus$&  &  $\oplus$& &  & &  &  &  &  $\oplus$&  &
  $\oplus$&  \\ \hline
 6&  $\oplus$&  &  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\ominus$& $\ominus$
 & $\oplus$& &  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\ominus$& $\ominus$
 \\ \hline
 7&  $\oplus$&  $\oplus$&  &  &  &  $\oplus$&  $\oplus$&  $\oplus$&
$\ominus$ & $\oplus$&  $\oplus$& $\ominus$ & $\ominus$ &  & $\oplus$ &
$\oplus$&  $\oplus$& $\ominus$ \\ \hline
 8&  & $\oplus$ &  & $\oplus$ & $\oplus$&  & $\ominus$&  $\ominus$&
$\oplus$ & $\ominus$ &  $\oplus$& $\ominus$ &  $\oplus$&  $\oplus$& $\ominus$
&  $\ominus$&  $\ominus$& $\oplus$ \\ \hline
 9& $\oplus$ &  &  &  &  &  &  $\oplus$& $\oplus$& $\ominus$ & $\oplus$ &
$\ominus$ & $\ominus$ & $\ominus$ & $\ominus$ & $\ominus$ & $\oplus$& $\oplus$
& $\ominus$ \\ \hline
 10&  &  $\oplus$&  &  &  &  &  &  $\oplus$ & $\oplus$ & $\ominus$ &
$\oplus$ & $\ominus$ & $\ominus$ & $\ominus$ & $\ominus$ & $\ominus$ &
$\oplus$ & $\oplus$ \\ \hline
 11&  &  & $\oplus$ &  &  &  & $\oplus$ & $\ominus$ & $\ominus$ & $\ominus$ &
 $\ominus$ & $\oplus$& $\ominus$ & $\ominus$ & $\ominus$ &
  $\oplus$& $\ominus$ & $\ominus$ \\ \hline
 12&  &  &  &  $\oplus$&  &  &  $\oplus$& $\ominus$ & $\oplus$ &
 $\ominus$ & $\ominus$ & $\ominus$ & $\oplus$ & $\ominus$ & $\ominus$ &
 $\oplus$ & $\ominus$ &  $\oplus$ \\ \hline
 13&  &  &  &  &  $\oplus$&  &  &  &  $\oplus$ & $\ominus$ & $\ominus$ &
 $\ominus$ & $\ominus$ & $\oplus$ & $\ominus$  & $\ominus$ & $\ominus$ &
 $\oplus$ \\ \hline
 14&  &  &  &  &  &  $\oplus$&  &  $\oplus$& $\ominus$ & $\ominus$ &
 $\ominus$ & $\ominus$ & $\ominus$ & $\ominus$ & $\oplus$ & $\ominus$ &
 $\oplus$ & $\ominus$  \\ \hline
 15&  $\oplus$&  &  $\oplus$&  $\oplus$&  &  & $\oplus$ & $\ominus$ &
 $\ominus$ & $\oplus$ & $\ominus$ & $\oplus$ & $\oplus$ & $\ominus$ &
 $\ominus$ & $\oplus$ & $\ominus$ & $\ominus$ \\ \hline
 16&  $\oplus$&  $\oplus$&  &  &  &  $\oplus$&  $\oplus$&  $\oplus$&
$\ominus$ & $\oplus$&  $\oplus$& $\ominus$ & $\ominus$ & $\ominus$ &
$\oplus$&  $\oplus$&  $\oplus$& $\ominus$ \\ \hline
 17&  &  $\oplus$ & & $\oplus$ & $\oplus$ &  &  $\ominus$ &  $\ominus$ &
$\oplus$ & $\ominus$ & $\oplus$ & $\ominus$ & $\oplus$ & $\oplus$& $\ominus$
&  $\ominus$&  $\ominus$& $\oplus$ \\ \hline
\end{tabular}
\caption{The non-zero pattern for the ordered matrix with block decomposition
\label{Ordered_real_A_block}}
\end{figure}

\begin{figure}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline
  &  0&  1&  2&  3&  4&  5&  6&  7&  8&  9& 10& 11& 12& 13& 14& 15& 16& 17\\ \hline
 0&  $\oplus$&  $\oplus$&  &  &  &  &  &  &  &  &  &  &  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$&  &  \\ \hline
 1&  $\oplus$&  $\oplus$&  &  &  &  &  &  &  &  &  &  &  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$&  &  \\ \hline
 2&  &  &  $\oplus$&  $\oplus$&  &  &  &  &  &  &  &  &  &  &  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$\\ \hline
 3&  &  &  $\oplus$&  $\oplus$&  &  &  &  &  &  &  &  &  &  &  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$\\ \hline
 4&  &  &  &  &  $\oplus$&  $\oplus$&  &  &  &  &  &  &  $\oplus$&  $\oplus$&  &  &  &  \\ \hline
 5&  &  &  &  &  $\oplus$&  $\oplus$&  &  &  &  &  &  &  $\oplus$&  $\oplus$&  &  &  &  \\ \hline
 6&  &  &  &  &  &  &  $\oplus$&  $\oplus$&  &  &  &  &  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\oplus$\\ \hline
 7&  &  &  &  &  &  &  $\oplus$&  $\oplus$&  &  &  &  &  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\oplus$\\ \hline
 8&  &  &  &  &  &  &  &  &  $\oplus$&  $\oplus$&  &  &  &  &  &  &  $\oplus$&  $\oplus$\\ \hline
 9&  &  &  &  &  &  &  &  &  $\oplus$&  $\oplus$&  &  &  &  &  &  &  $\oplus$&  $\oplus$\\ \hline
10&  &  &  &  &  &  &  &  &  &  &  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\oplus$&  &  \\ \hline
11&  &  &  &  &  &  &  &  &  &  &  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\oplus$&  &  \\ \hline
12&  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$&  &  &  &  &  $\oplus$&  $\oplus$&  $\ominus$&  $\ominus$&  $\ominus$&  $\ominus$\\ \hline
13&  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$&  &  &  &  &  $\oplus$&  $\oplus$&  $\ominus$&  $\ominus$&  $\ominus$&  $\ominus$\\ \hline
14&  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$&  &  &  &  &  &  &  $\oplus$&  $\oplus$&  $\ominus$&  $\ominus$&  $\oplus$&  $\oplus$&  $\ominus$&  $\ominus$\\ \hline
15&  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$&  &  &  &  &  &  &  $\oplus$&  $\oplus$&  $\ominus$&  $\ominus$&  $\oplus$&  $\oplus$&  $\ominus$&  $\ominus$\\ \hline
16&  &  &  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$&  &  &  $\ominus$&  $\ominus$&  $\ominus$&  $\ominus$&  $\oplus$&  $\oplus$\\ \hline
17&  &  &  $\oplus$&  $\oplus$&  &  &  $\oplus$&  $\oplus$&  $\oplus$&  $\oplus$&  &  &  $\ominus$&  $\ominus$&  $\ominus$&  $\ominus$&  $\oplus$&  $\oplus$\\ \hline
\end{tabular}
\caption{The non-zero pattern for the ordered matrix with decoupled real and
imaginary part using partial row pivoting factorization.
\label{Ordered_A_real_LUed}}
\end{figure}

\section{Conclusion}

   The lessons learned by exploring ETMSP 3.1

\begin{itemize}

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

\item ETMSP's operation is highly optimized to run in sequential environment.
      Many time saving tricks will be useless in parallel computing
      environment.
\end{itemize}

   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.

\end{document}
