%make two column with no page numbering, default is 10 point
\documentstyle[twocolumn]{article}
\pagestyle{empty}

%set dimensions of columns, gap between columns, and space between paragraphs
\setlength{\textheight}{8.75in}
\setlength{\columnsep}{2.0pc}
\setlength{\textwidth}{6.8in}
\setlength{\footheight}{0.0in}
\setlength{\topmargin}{0.25in}
\setlength{\headheight}{0.0in}
\setlength{\headsep}{0.0in}
\setlength{\oddsidemargin}{-.19in}
\setlength{\parindent}{1pc}

%I copied stuff out of art10.sty and modified them to conform to IEEE format

\makeatletter
%as Latex considers descenders in its calculation of interline spacing,
%to get 12 point spacing for normalsize text, must set it to 10 points
\def\@normalsize{\@setsize\normalsize{12pt}\xpt\@xpt
\abovedisplayskip 10pt plus2pt minus5pt\belowdisplayskip \abovedisplayskip
\abovedisplayshortskip \z@ plus3pt\belowdisplayshortskip 6pt plus3pt
minus3pt\let\@listi\@listI} 

%need an 11 pt font size for subsection and abstract headings
\def\subsize{\@setsize\subsize{12pt}\xipt\@xipt}

%make section titles bold and 12 point, 2 blank lines before, 1 after
\def\section{\@startsection {section}{1}{\z@}{24pt plus 2pt minus 2pt}
{12pt plus 2pt minus 2pt}{\large\bf}}

%make subsection titles bold and 11 point, 1 blank line before, 1 after
\def\subsection{\@startsection {subsection}{2}{\z@}{12pt plus 2pt minus 2pt}
{12pt plus 2pt minus 2pt}{\subsize\bf}}
\makeatother

\begin{document}
%don't want date printed
\date{}

%make title bold and 14 pt font (Latex default is non-bold, 16 pt)
%\title{\Large\bf My Wonderful Article in IEEE Format}
\title{\Large\bf Parallel Differential-Algebraic Equation Solvers for\\
Power System Transient Stability Analysis{\thanks{This work has been
supported in part by Niagara Mohawk Power Corporation, the New York State 
Science and Technology Foundation, the NSF under co-operative agreement No.
CCR-9120008, and ARPA under contract No. DABT63-91-K-0005.}} } 

%for single author (just remove % characters)
%\author{I. M. Author \\
%  My Department \\
%  My Institute \\
%  My City, ST, zip}
 
%for two authors (this is what is printed)
 \author{\begin{tabular}[t]{c@{\extracolsep{2em}}c}
Kamala Anupindi~~~~~ Anthony Skjellum** & Paul Coddington~~~~~ Geoffrey Fox\\
 \\
    $kamala@npac.syr.edu$ & $tony@aurora.cs.msstate.edu$\\
    Northeast Parallel Architectures Center  & ** NSF/Engineering Research
 Center \\
   Syracuse University & Mississippi State University  \\
  Syracuse, NY 13244 & Mississippi State, MS 39762\\
 \end{tabular}}

\maketitle

%I don't know why I have to reset thispagesyle, but otherwise get page numbers
\thispagestyle{empty}

\subsection*{\centering Abstract}
%IEEE allows italicized abstract
\noindent
{\em Real-time or faster-than-real-time power system transient stability 
simulations will have significant impact on the future design and
operations of both individual electrical utility companies and large
interconnected power systems. The analysis involves solution of extremely
large systems of differential and algebraic equations. Differential-Algebraic
Equation (DAE) solvers have been used to solve problems similar in 
nature to the transient stability analysis (TSA) problem. This paper 
discusses the possibility of the use of the existing DAE solvers
to solve the transient stability analysis application. We also discuss our
research in developing a scalable, parallel DAE solver for use by the
power system community and in related applications \cite{Skj90d}. }

\section{Introduction}
\noindent
In engineering and science, many problems give rise to smooth and/or
discontinuous systems of differential equations coupled with nonlinear
algebraic equations. These coupled systems of differential and algebraic
equations are commonly known as differential/algebraic equations (DAEs).
Research in the field of DAEs started about a decade ago and is still 
in its early stages. A number of engineering and science problems have 
been solved with the help of DAE solvers developed by \cite{Brenan,
Deuf, Hairer, Alvin, Skj90d}. But to this date, there has been no 
investigation of the application of DAE solvers to power systems
problems such as the Transient Stability Analysis (TSA) \cite{Anderson}.\\

Transient Stability Analysis is a compute-intensive problem requiring
the solution of nonlinear, discontinuous systems of DAEs for the 
purpose of simulating such electrical power systems phenomena as
network component and transmission line failures. Existing DAE solvers 
often use higher order implicit Runge-Kutta integration techniques or 
Backward Differentiation Formulas (BDF), techniques that yield more 
accurate solutions or permit larger time-steps.  DAE solvers have 
been especially developed to handle discontinuities in functions, 
a condition encountered in transient stability calculations.\\

This paper describes research on the use of existing DAE solvers to 
solve the power system TSA problem. In section 2, a detailed description 
of the power system problem is given. In section 3, we discuss various 
existing DAE solvers and their potential for application to the TSA 
problem. Our ongoing research in the benchmarking of various sequential 
solvers and the development of parallel solvers is discussed in section 
4. Lastly, in section 5, our conclusions on the application of DAE 
solvers for TSA problem are discussed.

\section{Power System Transient Stability Analysis Problem}
\noindent
Transient Stability Analysis (TSA) examines the dynamic behavior of a
power system for as much as several seconds following a power-flow
disturbance. The general nature of the transient stability problem is
explained in \cite{Mack} as follows: Consider a simple mechanical 
analog where a number of masses, representing the generators in an 
electric power system, are suspended from a ``network'' consisting of 
elastic strings representing the electric transmission lines. Now one 
of the strings is suddenly cut, corresponding to a sudden loss of a 
transmission line. Thereafter the masses will experience transient 
coupled motion, and the forces in the strings will fluctuate. This 
sudden disturbance may cause one of the following two effects: 

\begin{enumerate}
\item The system will settle down to a new equilibrium state, 
characterized by a new set of string forces corresponding to the power 
voltages in the electricity case.

\item  Because of the transient forces, one or more additional strings 
will break, causing a weaker network, resulting in a chain reaction of 
broken strings and total collapse of the system. If the system has the 
strength to survive the disturbance and settle into a new steady state, 
it is referred to as ``transient stable for the fault in question.''
\end{enumerate}

\noindent
The transient stability analysis problem can be described by a system of
DAEs where the generators are represented by differential equations and 
the power system network interconnecting the generators is represented 
by nonlinear algebraic equations.  The reaction of a single generator to 
variations in its load can be modeled using the following Ordinary 
Differential Equations (ODEs):

\begin{center}
{\em
E +$x_{d}I_{d} = E^{'}_{q} + x^{'}_{d}$\\
\medskip
$E_d + x_{q}I_{q} = E^{'}_d + x^{''}_{q}I_{q}$\\
\medskip
${\tau}^{'}_{q0}{\dot{E}}^{'}_{d} = -E^{'}_{d} - (x_{q} - x^{'}_{q} )I_{q}$
\\\medskip
$\dot{E}^{'}_{q} = {\frac{1}{T^{'}_{d0}}}(E_{FD} - E)$ \\
\medskip
$T_{e} = E^{'}_{d}I_{d} + E^{'}_{q}I_{q} - (L^{'}_{q} -
L^{'}_{d})I_{d}I_{q}$\\
\medskip
${\tau}_{j}{\dot\omega} = T_{m} - D{\omega}-[E^{'}_{d}I_{d}+E^{'}_{q}I_{q} -
(L^{'}_{q} - L^{'}_{d})I_{d}I_{q}]$\\
\medskip
$\dot{\delta} = \omega - 1 $
\medskip
}
\end{center}

\noindent
The equations described above roughly describe the ``structure'' of the
differential equations involved in the transient stability analysis DAEs.
The notation used in the above equations is as follows:
\begin{itemize}
\item $E$ is the magnitude of source voltage of the generator,
\vspace{-0.09in}
\item $I$ is the field current of the generator,
\vspace{-0.09in}
\item $\delta$ is the rotor angle, 
\vspace{-0.09in}
\item $T_{e}$ and $T_{m}$ are the electrical and mechanical torques 
 respectively, 
\vspace{-0.09in}
\item $\omega$ is the shaft angular velocity,
\vspace{-0.09in}
\item $x$ is the reactance of the generator, 
\vspace{-0.09in}
\item $L$  is inductance involved, 
\vspace{-0.09in}
\item ${\tau}_{j}$ and ${\tau}^{'}_{q0}$ are the
time constants involved in the derivation of these equations.
\end{itemize}
\noindent
Subscripts $d, q, F$, and $D$ represent the respective $d$-axis, $q$-axis,
$F$-axis and $D$-axis components of the elements described.
Detailed explanation on the notation involved in the 
described equations can be found in \cite{Anderson}.\\

There are multiple generators supplying power to the network and the 
generators are coupled through an electrical power network. The power
network can be modeled by nonlinear algebraic equations as follows:
\begin{center}
{\em I = V Y }
\end{center}
\noindent
where $I$ is the network current, $V$ is the network voltage and $Y$ is 
the network admittance.  Apart from these two sets of equations, TSA 
involves equations for loads (which can be differential or algebraic) 
and equations for the generator control (that are differential). All 
these equations together form a system of DAEs that must be solved 
simultaneously. 

\paragraph*{ Structure of the TSA DAEs:}
The transient stability analysis DAEs are
non-symmetric in nature. They are of bordered-block-diagonal form wherein
blocks of generator equations along the diagonal are coupled with the 
power system distribution network. The admittance matrix involved is 
extremely large, complex and sparse. Research is being conducted 
\cite{David} to reorder this matrix into block-diagonal-bordered form 
in order to exploit the structure for parallelism.

\paragraph*{ Order of the TSA DAEs:}
The number of equations involved in the transient stability analysis 
solution is extremely large. Consider an example of 20 differential 
equations to describe each generator and two equations to describe 
the complex voltage/current at each network bus. Then an interconnected 
power system with 2000 buses and 300 generators could generate a sparse,
unsymmetric system of 10,000 nonlinear algebraic equations that must be 
solved simultaneously. For regional power systems, the number could be 
five times higher. 

\paragraph*{ Index of the TSA DAEs:}
The differential index of a system is the minimal number of analytical 
differentiations needed such that the differential equations of the DAE 
can be transformed by algebraic manipulations into an explicit ODE 
system. The transient stability problem leads to an index-one or 
higher-index DAE system with the possibility of the ODEs involved being 
stiff in nature.

\paragraph*{ Algebraic Equations of TSA:}
The network equations of the TSA are
complex in nature (the voltages and currents involved are both complex).
These equations are also nonlinear. Research in developing a solver that 
can be used directly to solve a system of complex, nonlinear algebraic 
equations has been limited and leaves the field wide open for further
research. 

\paragraph*{ Scalability:}
Since a regional power network can consist of many individual power 
systems, it is highly desirable that the DAE solver used for the 
transient stability analysis be scalable. Also, the power networks can 
vary in size depending on the region under consideration.
 
\section{Differential-Algebraic Equation Solvers}
\noindent
The power system community has been solving the transient stability 
problem approximately as a decoupled system of differential and 
algebraic equations. As explained in \cite{Brenan}, this method of 
reduction of large DAE system to a system of explicit ODEs not only 
destroys the sparsity of the system but also prevents the exploitation 
of the system structure (symmetry of the Jacobian). Several sequential 
DAE solvers exist that use various numerical techniques for their 
solution. Prominent methods among the existing solvers are as follows:
\begin{itemize}
\item Extrapolation methods \cite{Deuf,Hairer},
\item Implicit Runge-Kutta methods \cite{Hairer},
\item Backward Differentiation Formulas (BDF) methods \cite{Brenan,nmg-Hi82,nmg-h83,nmg-Hi85b},
\item Rosenbrock Methods \cite{Hairer}.
\end{itemize}
\noindent
Our research is concentrated on analyzing the following DAE solvers:

\begin{itemize}
\item {\bf LIMEX} \cite{Deuf} uses an extrapolation technique. As 
discussed in \cite{Brenan}, one-step methods such as extrapolation 
techniques have an inherent advantage over BDFs methods when applied to 
DAEs with frequent discontinuities. LIMEX is a code developed for 
semi-explicit index one DAEs. It implements an extrapolation of the 
semi-implicit Euler method and can handle stiff ODEs. This code has 
been successfully used in the solution of problems arising in combustion 
modeling where there are frequent discontinuities in time. Because of 
its ability to handle stiff differential equations and frequent 
discontinuities, LIMEX is an excellent choice for the power system TSA 
solution.

\item {\bf SODEX} \cite{Hairer} and {\bf SEULEX} \cite{Hairer} also use 
extrapolation technique for the solution of DAEs.

\item {\bf RADAU5} \cite{Hairer} uses an Implicit Runge-Kutta (IRK) 
method of order five. Research on solution of stiff ODEs has focused a 
great deal on Implicit Runge-Kutta methods. RADAU5 is designed to solve 
index one, two, and three systems.  IRK methods have a definite advantage
over multistep methods such as BDFs methods when the DAEs exhibit 
frequent discontinuities. Because of their one-step nature, IRK methods 
can be started at a higher order after every discontinuity, whereas 
multistep methods must be restarted, usually at a low order, after 
every discontinuity. This fact helps the IRK methods to be more 
efficient. Also, with IRK methods, it is possible to construct high 
order A-stable or nearly A-stable IRK formulas, which is important 
while solving stiff ODEs with eigenvalues lying close to the imaginary 
axis. The generator ODEs in the transient stability analysis may be of 
this type. Since the TSA DAEs can be of index higher than one, RADAU5, 
with its ability to handle higher index systems makes an excellent 
choice for the transient stability analysis problem.

\item {\bf DASSL} \cite{Brenan} uses a multistep, backward differentiation 
technique. It uses a variable step-size, variable-order 
fixed-leading-coefficient implementation of BDFs formulas to advance 
solution from one time step to the next. DASSL can solve DAEs of index 
zero and one. BDFs methods suffer no order reduction for index one 
systems. Also, BDFs methods achieve the same order of convergence for 
this class of DAEs as they do for ODEs. DASSL has a robust order 
selection strategy for DAEs with eigenvalues close to the imaginary 
axis in the complex plane. DASSL can also be used to solve semi-explicit
index-two systems of DAEs. DASSL by itself cannot handle discontinuities,
but a variation of DASSL, DASRT, has root finding capabilities to locate
discontinuities when they are sufficiently large that DASSL cannot 
integrate through without intervention. Research is being conducted at 
the Northeast Parallel Architectures Center (NPAC) \cite{Alvin} to solve 
the TSA problem using DASSL. 

\item {\bf LSODI} \cite{nmg-Hi82,nmg-h83,nmg-Hi85b} uses BDFs. It uses 
a fixed-coefficient implementation of BDFs formulas to solve linearly 
implicit system of DAEs. The fixed-coefficient methods can be implemented
efficiently for smooth problems, but suffer from inefficiency and 
possible instability for problems requiring frequent step-size 
adjustments. LSODI solves DAEs of index zero but if the coefficient 
matrix is singular, then it solves the resulting index one DAE system. 
It can solve systems of stiff and non-stiff ODEs as well. LSODI has been 
successfully used to solve systems arising from PDE problems. Another 
version of LSODI, LSODA, switches automatically  between stiff (BDF) and
non-stiff (Adams) methods.  Because of its ability to handle stiff ODEs 
and frequent discontinuities, LSODI presents a viable alternative to 
other solvers discussed for the power system TSA problem.

\item {\bf RODAS} \cite{Hairer} uses an embedded Rosenbrock technique of 
order four. Rosenbrock methods have a big advantage over IRK methods in 
that they completely avoid non-linear systems of equations, at the same 
time providing the advantages of accurate solutions with stiff 
differential equations. RODAS can solve DAEs that are expressible in 
semi-explicit form and that are of index one. Because of the lack of 
data about the application of RODAS to any significant problem, and 
because of the varied advantages of other solvers discussed above, our 
research will not involve RODAS, even though it has the qualities of
being a good contender.
\end{itemize}

\section{Research Status}
\noindent
Currently at NPAC, we are involved in the following two issues:
Transient Stability Analysis Benchmarks, and Parallel DAE Solvers.

\subsection{Transient Stability Analysis Benchmarks}
\noindent
We are setting up transient stability analysis benchmarks with 
interfaces to different sequential solvers like LIMEX, LSODI, RADAU5 
and DASSL. This is to test our empirical conclusion that these solvers 
are best suited for the TSA problem. We will also present a comparative 
review of their speed and performance. The test data being used is the 
standard IEEE 118-bus data but we will also benchmark larger systems at 
a later stage.

\subsection{Parallel DAE Solvers}
\noindent
CDASSL (Concurrent DASSL) \cite{Skj90d} is the first attempt that has 
been made at parallelizing these DAE solvers. More recently, attempts 
have also been made to develop a data-parallel and message-passing 
versions of DASPK \cite{Linda} on the CM5. Our research is examining 
the scalability of the various DAE solution techniques. We are 
concentrating on selecting the most suitable solver for the power system 
application, and creating a scalable library for an application based on 
it. The parallel version of this DAE solver is being developed for the 
use by the power system community, and in related applications such as 
the chemical plant simulations ({\em e.g., }\cite{Skj90d}) and the 
electrical circuit simulations \cite{Brenan}. We are using the {\em 
Multicomputer Toolbox} \cite{Skj93c} for the development of the parallel 
solver. This software is currently based on the portable message-passing 
system {\sl Zipcode\/} \cite{Skj93b}; the {\em Toolbox } will port to 
{\em MPI } \cite{MPI} next year and hence our programs will run on all 
reasonable message-passing multicomputers.

\section{Conclusions}
\noindent
We have discussed the power system transient stability analysis problem 
(TSA). We also presented a review of some of the existing DAE solvers 
that we consider best suited for this problem. Benchmarking of sequential
DAE solvers and the development of a parallel, scalable DAE solver are in 
progress. This research is being closely coordinated with parallel sparse 
matrix solver \cite{David} research being performed at the Northeast 
Parallel Architectures Center (NPAC) at Syracuse University.

%this is how to do an unnumbered subsection
\subsection*{Acknowledgments}
\noindent
We wish to thank David Koester and Alvin Leung, of NPAC, for their help 
with this research.
\noindent
\bibliographystyle{plain}
\bibliography{final}

\vspace{1.5in}

\end{document}


