\documentstyle[11pt]{article}
\setlength{\oddsidemargin}{.25in}
\settowidth{\marginparwidth}{Lect 99}
\setlength{\marginparsep}{1em}
\setlength{\textwidth}{6in}
\setlength{\topmargin}{-.5in}
\setlength{\textheight}{9in}
\renewcommand{\baselinestretch}{1}
\setlength{\parindent}{0in}

\reversemarginpar

\input{psfig}


\begin{document}

\section{Differential algebraic equation (DAE)}
\label{DAE}

A simple description of DAE is ordinary differential equation with algebraic
equation.   Different ways to look at DAE can tell us how DAE can be used.
A system of DAE can be formed by a model of
a system (ODE), which we want to investigate, have some algebraic. constrains,
\[
\begin{array}{lcl}
x' & = & f(x, y, t) \\
y & = & g(x, t) \\
\end{array}
\].

This type of systems is called Constrained Variational problems.  Transient
stability problem is of this type.  DAE can also come from model reduction
and singular perturbations.  Another view of the network equation can classify
Transient stability problem as model reduction.

Model reduction is a system of ODE of which
some variables contribute a very small factor to the system.  When we are
interested in the first order approximation of its behavior, we eliminate
the some of the variables.
\[
\begin{array}{lcl}
x' & = & f(x, y, t, \epsilon) \\
\epsilon y' & = & g(x, y, t, \epsilon) \\
\end{array}
\],
where $\epsilon \rightarrow 0$.
The resulted system is a DAE system.

The difference between DAE and ODE is that DAE can be more difficult to solve.
The degree of difficulty can be indicated by the index of the DAE system.

A system of ODE is classify as index 0.
Next, we have index 1 DAEs.  These DAE can subdivided into two classes,
semi-explicit
\[
\begin{array}{lcll}
x' & = & f(x, y, t) \\
0 & = & g(x, y, t) & \exists \frac{\partial g}{\partial y} \\
\end{array},
\]
and full-implicit system,
\[
\begin{array}{lcll}
0 & = & f(x, x', y, t) \\
0 & = & g(x, y, t) & \exists \frac{\partial g}{\partial y} \\
\end{array}
\].  There are higher index problem; however, we will not introduce them in
this report.

The underlining model of electrical transient analysis is based on the
simple RLC circuit model.  The transmittion line, transformer and
capator is modeled by
a perfect components, such as resistor, capticler and inductor.
The infinite voltage and current sources is replace by
the classical generator model and the 2-axis generator model.
A classical generator model is a perfect voltage source with bookkeeping
for the rotor angle added.  The 2-axies generator model is more detail model
of a generator.  As a result, the rotor angle keep track by a 2-axies
generator
model is more accurate.  The two type of generators contribute ODEs to the
model's equations.  The network is contribute the algebraic. system of equation.
The model's equations are
\[
\begin{array}{lcl}
x' & = & f(x, V, I) \\
I & = & Y(x) V \\
\end{array}
\].  As a result, electrical transient problem consists of solving
semi-explicit index one differential algebraic equations.

\subsection{Integration methods}

In PTI software, modified Euler method is used.  There are two types of
modified Euler methods, explicit and implicit.  Explicit modified Euler
method is not stiff accurate.
In other words, the explicit modified Euler method requires small step size
to maintain the desire accuracy when solving stiff systems~\cite{DE}.
A stiff system
is generally referred to a system, which some of the variables varies
slowly and some variables varies very fast.  In our case, the generator's
rotor angle varies slowly (governed by the ordinary differential equations),
while the
voltage of the network varies instantaneously (governed by the algebraic equations).
Since the explicit modified Euler method is not stiff accurate, the numerical
solution of the DAE system do not automatically satisfy the algebraic part of
the DAE, ie, $I_{n} - Y(x_{n}) V \neq 0$.  The formulation
of implicit modified Euler method is the same as the trapezoidal rule.
Implicit modified Euler method, on the other hand,
have a very different properties.

In ETMSP, electrical transient analysis problem is currently solved by
applying trapezoidal rule to the system~\cite{ETMSPPM}.
The system of equation becomes,
\begin{eqnarray}
\left[
\begin{array}{c}
F \\
G \\
\end{array}
\right]
& = &
\left[
\begin{array}{c}
x_{n + 1} - x_{n} - \frac{h}{2}(f(x_{n + 1}, V_{n + 1}) + f(x_{n}, V_{n})) \\
YV_{n + 1} - I(x_{n + 1}, V_{n + 1}) \\
\end{array}
\right]. \label{Trap} \\
\end{eqnarray}

Trapezoidal rule is a stuff accurate, order 2 method~\cite{SODEII}.  A stuff
accurate method implies that the step size used on solving the problem have
reduced impact on the accuracy of the numerical result.  The order 2 accurate
means the order convergence is proportional to $h^{2}$.  In addition, the
trapezoidal
rule is A-stable, so that the behavior of the system model will not be affected
by the analytical behavior of the numerical method.
Since trapezoidal rule is stuff accurate, the algebraic part of the DAE system
is satisfied~\cite{Num-Sol}.  Consequently, the numerical result of this method
is more
accurate and the method allows to use large step size.  Furthermore, new
numerical model can be developed without any concern on the analytical
behavior of the numerical DAE solver's method.  As a result, ETMSP is a good
research tools on developing electrical component's models.

The implicit non-linear algebraic equations~\ref{Trap} are solved by
Newton method.
After applying Newton method to the system, the system becomes
\begin{eqnarray}
\left[
\begin{array}{c}
F \\
G \\
\end{array}
\right]
& = &
\left[
\begin{array}{cc}
A & B \\
C & D \\
\end{array}
\right]
\left[
\begin{array}{c}
\triangle x \\
\triangle  V \\
\end{array}
\right], \label{Trap-Newton} \\
\end{eqnarray}
where
$A$ is the Jacobian matrix for the generator equations,
$B$ is the Jacobian matrix for the coupling between generators,
$C$ is the Jacobian matrix for the coupling between generators and
$D$ is the Jacobian matrix for the network equations.
The system then further decoupled into two systems,
\begin{eqnarray}
x_{n + 1}^{k + 1} & = & x_{n + 1}^{k} - A^{-1}(F_{n + 1}^{k}
+ B(V_{n + 1}^{k + 1} - V_{n + 1}^{k}) \\
(Y + D - CA^{-1}B)(V_{n + 1}^{k + 1} - V_{n + 1}^{k}) & = & -G_{n + 1}^{k} +
CA^{-1}F_{n + 1}^{k} \\
\end{eqnarray}
The first system consists of all the generator equations.  The sparcity pattern
of $A$ is block diagonal.  Each diagonal block represents ODE of the
generator's model. As a result, the  generator equation is mutually independent
of other generator.  The nonzero values of the second matrix contributed mostly
by the network equations, $Y$; however, additional nonzero values may be added
by $D - CA^{-1}B$, which is introduced by the decoupling process.

The advantages of this decoupling reduce the over all size of the matrices
that need to be factorized.  This reduction can be show by the following
calculation.  Let $a$ be the dimension of the Jacobian matrix of the generator
equations, $A$ and $d$ be the dimension of the Jacobina matrix of the network
equations, $D$. 

\begin{math}
\begin{array}{lcl}
(a + d)^{3} & = & (a^{2} + 2ad + d^{2})(a + d) \\
& = & a^{3} + 3a^{2}d + 3ad^{2} + d^{3} \\
\end{array}
\end{math}
This computation assumes that the matrices, $A$ and $D$, are dense.
The reduction is discounted by using sparse matrix algorithms.  However,
this reduction can still be significant for sequential processing.
The trade off for this processing method is the destruction of the
sparcity pattern offered by the system.


\subsection{DAE solver}

As we see from the previous section, the trapezoidal rule is a excellent tool
for developing models for new components.  However, once the component is
analyzed and determined stable.  Trapezoidal rule becomes the limiting
factor of the computation speed.  In this section, we will introduce a new
DAE solver package to replace trapezoidal rule.  
DASSL~\cite{DASSLSW}~\cite{DASPKSW} is specified designed for solving index
one DAEs.
The features of this packages are
variable convergence order, variable step size and clean interface to
linear sparse solver.

DASSL~\cite{Num-Sol}, the DAE solver package, uses implicit backward
difference formula.  BDF provides
formula for convergence from order 1 to 5.  The analytical stability of
order 1 and 2 BDF formula are A-stable.  In fact, order 1 BDF is Euler
method and
order 2 BDF is modify Euler method.  However, the gain in
computational speed are expected form using the higher order of convergence
formula.  The solver package have automatic control for shifting from the
lower formula to higher order formula, when the function exhibit stable
behavior for higher order formula.  At higher convergence order, the error
of the result is reduced.  Hence, larger step
size can be used to compute the result with the same accuracy.  However,
the trade-off is that the higher order formula is not A-stable.  The stability
of these higher formula is A($\alpha$)-stable~\cite{SODEII}.
\begin{tabular}{|c|c|c|c|c|c|}
\hline
order    &  1 &  2 &  3    &   4   &   5  \\ \hline
$\alpha$ & 90 & 90 & 86.03 & 73.35 & 51.84 \\ \hline
\end{tabular}
As a result, using DASSL for model development requires in-depth knowledge of
DASSL and the working knowledge of the model.  However, if the model's behavior
is known, then DASSL can be used to provide fast simulations.

DASSL~\cite{Num-Sol} package also have a automatic variable step size control.
Although
ETMSP only allows fixed step size, its users usually use different step
size for different phase of the simulation.  For example, the user uses
larger step size before fault, smaller step size just after the fault and
larger step size after the disturbance is expected to be settled fault.
In DASSL, the selection of the step size is based on the acceptable error,
which is defined by the user, and the behavior of the function.  If the
history of the function indicates the function is ``slow'' varying, then the
step size is increase until either the error is larger than specified or
unstable result is computed.  As a result, we expected further reduction
of the runtime of a simulation.

Furthermore, the DASSL package have a well defined interface to linear system
solver~\cite{SLCDASSL}.  This interface allows DASSL to call any parallel
sparse matrix solver with minimum code changes.  Currently, DASSL can use two
major types of linear solver, direct solver and semi-iterative solver.

\section{Sparse Solver}

\subsection{ETMSP}

   ETMSP uses a inhouse developed sparse matrix factorization
package.
In this section, we briefly introduce facilities of this sparse matrix
package.  The sparse package provides five major functions~\cite{ETMSPSW},

   \begin{itemize}
   \item Matrix reordering
   \item Symbolic factorization
   \item LU factorization
   \item Forward and Back solve
   \item Interface and Supporting functions
   \end{itemize}

The matrix reordering strategy consists of static minimum degree ordering.
The matrix is reordered once before the first factorization.  The
reordering is redone only, when there is significant changes to the matrix
structure.  Minimum degree ordering~\cite{DM_SMat} is referred to placing the
bus with
minimum number of off diagonal elements before the bus other buses.
This step is repeated for the $n - 1 \times n - 1$ submatrix,
where $n$ is the dimension of the original matrix.

The symbolic factorization can be used to discover the fill-ins patterns of
the matrix during the factorization with specific ordering.  Symbolic
factorization is computed
in a similar way to a factorizing a matrix; however, no numerical computation
is done.  As a result, symbolic factorization is much cheaper than a full
factorization.  The symbolic factorization function in the ETMSP sparse
solver packages is not used.

The LU factorization function of this package offers three functionalities.
First, the LU factorization perform a full LU factorization on a symmetric
structured matrix with ordering given by the matrix reordering function.
Second, the function can partially un-factorization the matrix.  This
function is intend to provide a way to reduce the cost of factorization, when
the changes of the matrix occur at the later half of the matrix.
Third, partial factorize the previously partial un-factorization matrix.
However, ETMSP only use the function for full LU factorization.

Forward and Back solve functions is used in the very dishonest newton function.
The interface function provides interface services for ETMSP to access the
sparse package.  Additional function to store and to restore the solver's
state is provided. 

\section{Semi-iterative solver}

   In the course of solving DAEs, DASSL~\cite{DASSLSW}~\cite{DASPKSW} uses
   the Newton's method to
   solve systems of non-linear equations.  Each Newton iteration requires to
   solve a system of linear equations, $Ax = b$.  In the
   original DASSL~\cite{DASSLSW}, LU factorization is employed to solve the
   system of linear
   equations.   However, the recent advances on Krylov subspace iterative
   linear solvers, for example, CGS, QMR and GMRES, has prompted the
   DASSL's authors
   to incorporate GMRES iterative solver into DASSL~\cite{DASPKSW}.
   In~\cite{Peter1},
   ~\cite{Peter2} and ~\cite{Peter3}, the authors show DASSL and Krylov
   subspace iterative solver combination, the sequential software package
   is renamed to DASPK, offers
   reduced memory storage for large problems.  Experience shows that
   DASSL with semi-iterative solver computes the
   solution faster than DASSL with LU factorization.  The test
   program is
   solving heat equations, so the timing of the test program has no
   bearing on predicting the DASSL's
   performance on solving the electrical power problem.  However, the test
   does demonstrate a possible way for exploration of fast on-line transient
   stability program. Currently, Concurrent DASPK is capable to use three
   iterative solvers, CGS, TFQMR and GMRES.

   When computing a Newton iteration, DASPK does not compute the Jacobian
   matrix ($J$). Instead, an approximation vector of $Jy$ is generated by
   using a difference quotient~\cite{Peter3},

\[ Jy = F'(x)y \approx \frac{F(x + \sigma) - F(x)}{\sigma}, \sigma \in \Re. \]
   Using this approximation, DASSL has shown to converge theoretically and
   numerically,
   when DASSL is applied to a system of ODEs.  Although no difficulties
   are expected
   from using DASPK on the electrical power problem, a closer study of the
   interaction of DASPK with the electrical power DAEs is needed.

   In order to further speed up the convergence of the semi-iterative solving
   $Ax = b$, a preconditioner is used.
The goal of preconditioner is to generate a
precondition matrix ($C$) as close to the inverse of Jacobian matrix as possible
with acceptable amount of computations.

The basic function of the precondition matrix can be illustrate by the equation,

\[ CJx = Cb. \]
If $C = J^{-1}$, then the system of linear equations is solved.  However,
$J^{-1}$ is very difficult to generate, so $C$ is a good
approximation of $J^{-1}$ at best.  In this situation, the precondition matrix
serves as a mechanism to cluster the eigenvalues of $J$ as shown in
figure~\ref{ex1}. 
Although $C \approx J^{-1}$ does not solve the system, it maps the system of
linear equations to a better behaved system in the point of view of the Krylov
subspace methods.

\begin{figure}
\begin{minipage}{2in}
\begin{math}
\begin{array}{lcl}
A & = &
\left[
\begin{array}{cccc}
3 & 1 & 0 & 0 \\
1 & 3 & 1 & 0 \\
0 & 1 & 3 & 1 \\
0 & 0 & 1 & 3 \\
\end{array}
\right]
\end{array}
\end{math}
\end{minipage}
\hspace{1em}
\begin{minipage}{2in}
\begin{math}
\begin{array}{lcl}
C & = & 
\left[
\begin{array}{cccc}
\frac{1}{3} & 0 & 0 & 0 \\
0 & \frac{1}{3} & 0 & 0 \\
0 & 0 & \frac{1}{3} & 0 \\
0 & 0 & 0 & \frac{1}{3} \\
\end{array}
\right]
\end{array}
\end{math}
\end{minipage}


\begin{math}
\begin{array}{lcl}
\Lambda_{A} & = & 
\left[
\begin{array}{cccc}
3.6180 & 0 & 0 & 0 \\
0 & 4.6180 & 0 & 0 \\
0 & 0 & 2.3820 & 0 \\
0 & 0 & 0 & 1.3820 \\
\end{array}
\right]
\end{array}
\end{math}

\begin{math}
\begin{array}{lcl}
\Lambda_{CA} & = & 
\left[
\begin{array}{cccc}
1.2060 & 0 & 0 & 0 \\
0 & 1.5393 & 0 & 0 \\
0 & 0 & 0.7940 & 0 \\
0 & 0 & 0 & 0.4607 \\
\end{array}
\right]
\end{array}
\end{math}


\caption{The eigenvalues of $A$ is clustered after the application of
the precondition matrix C. \label{ex1}}

\end{figure}

The underline process of the Krylov subspace methods is to
find a function ($f$) such that $f(0) = 1$ and the roots of $f$ are the
eigenvalues of the matrix $J$.  The function is found by
increasing the index of the current function by one for each iteration as
shown in figure~\ref{ex2}.
Another important property of iterative methods is that the large eigenvalues
are found with relative ease.  When a iterative solver is terminated at the
$i$th
iteration and $i$ is less than $\mbox{dim }J$, the user accepts the function,
$f_{i}$, as a ``good'' approximation to the function, $f_{n}$.
Consequently, the clustered eigenvalues allows
a lower order function to approximate the higher order function better,
which translates into fast convergence of the linear solver.

\begin{figure}

\begin{minipage}{2.4in}
\begin{picture}(173,121)
%\BoundingBox{173}{121}
\psfig{figure=min_poly.ps,height=121pt,width=173pt}
\end{picture}
\end{minipage}
\begin{minipage}{2.4in}
\begin{picture}(173,121)
%\BoundingBox{173}{121}
\psfig{figure=poly_0.ps,height=121pt,width=173pt}
\end{picture}
\end{minipage}

\begin{minipage}{2.4in}
\begin{picture}(173,121)
%\BoundingBox{173}{121}
\psfig{figure=poly_1.ps,height=121pt,width=173pt}
\end{picture}
\end{minipage}
\begin{minipage}{2.4in}
\begin{picture}(173,121)
%\BoundingBox{173}{121}
\psfig{figure=poly_2.ps,height=121pt,width=173pt}
\end{picture}
\end{minipage}

\begin{minipage}{2.4in}
\begin{picture}(173,121)
%\BoundingBox{173}{121}
\psfig{figure=poly_3.ps,height=121pt,width=173pt}
\end{picture}
\end{minipage}
\begin{minipage}{2.4in}
\begin{picture}(173,121)
%\BoundingBox{173}{121}
\psfig{figure=poly_4.ps,height=121pt,width=173pt}
\end{picture}
\end{minipage}

\caption{The minimal polynomial plot shows a polynomial with the eigenvalues as
roots.  The degree $i$ polynomial shows the shape of the current polynomial at
$i$ iteration of the linear solver.  \label{ex2}}
\end{figure}

\newpage

\bibliography{report}
\bibliographystyle{plain}

\end{document}
\endinput


