
% =====

% Executive Summary
% Parallel DAE solvers
% Sparse matrix
% Survey
% Appendix on DAEs and their solution

\documentstyle[11pt]{article}

\input epsf

\textwidth=6.0in
\textheight=8.5in
\topmargin=-0.5in
\oddsidemargin=0.25in
 
\begin{document}

\begin{center}

{\LARGE \bf Final Report\\} 
\vspace*{1.0cm}
{\LARGE \bf
High Performance Computing and Parallel \\
Processing in the Utility Industry \\}
\vspace*{0.8cm}
{\Large \bf NMPC Project 07-9199} \\
\vspace*{1.0cm}
{\Large \bf Phase 1 -- Feasibility Study} \\
\vspace*{1.5cm}
{\Large Northeast Parallel Architectures Center \\
\vspace*{0.2cm}
Syracuse University} \\
\vspace*{1.5cm}
{\Large January 16, 1993} \\

% ABSTRACT??

\end{center}


\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgements}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The people at Syracuse University who have worked on this project are: 
\begin{itemize}
\item {\bf Northeast Parallel Architectures Center} 
 \vspace*{0.1cm} \\ 
 Prof. Geoffrey Fox (Principal Investigator) \\
 Dr. Paul Coddington (Project Leader) \\
 Prof. Sanjay Ranka \\
 Prof. Salim Hariri \\
 Prof. Alok Choudhary \\
 Dr. James Allwright \\
 Kamala Anupindi \\
 David Koester \\
 Alvin Leung \\
 Siva Narendra 
\item {\bf School of Information Studies} 
 \vspace*{0.1cm} \\
 Prof. Mike Nilan \\
 Richard Southwick  \\
 Caren Obad 
\end{itemize}

\noindent
We have received invaluable help and guidance from the following people
at Niagara Mohawk Power Corporation: 
\begin{description}
 \item ~~\vspace*{-0.6cm}\\
 Sid Lipton (Project Co-ordinator) \\
 Larry Eng    \\
 Bill Merritt \\
 Wes Yeomans  \\
 Jack Brennan \\
 Danielle Dylan  
\end{description}

\noindent
The following people gave their time for interviews to develop the knowledge
base for the transient stability analysis expert system:
\begin{description}
 \item ~~\vspace*{-0.6cm}\\
 Larry Eng (NMPC), 
 Ricardo Gonzalez, Jonathon Mayo, Robert Waldele (NYPP),
 Michael Henderson, Larry Hochverg, Paul Roman, George Dunn (NYPA). 
\end{description}

\noindent
We have also benefited greatly from discussions with the following people: 
\begin{description}
 \item ~~\vspace*{-0.6cm}\\
 Tony Skjellum (Mississippi State University)  \\
 Joel Saltz (University of Maryland)  \\
 Chris Pottle (Cornell University)  \\
 David Wong (Ontario Hydro)  \\
 Nelson Simons (General Electric) \\ 
 Rana Mukerji (General Electric)  
\end{description}


\newpage

\tableofcontents

\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{List of Acronyms}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{tabular}{ll}
{\bf AI}      & Artificial intelligence \\
{\bf BDF}     & Backward Differentiation Formula \\
{\bf CRPC}    & Center for Research in Parallel Computing \\
{\bf DAE}     & Differential Algebraic Equation \\
{\bf DEC}     & Digital Equipment Corporation \\
{\bf DSA}     & Dynamic Security Assessment \\
{\bf EMS}     & Energy Management System \\
{\bf EPRI}    & Electric Power Research Institute \\
{\bf ESEERCO} & Empire State Electrical Energy Research Corporation \\
{\bf ETMSP}   & Extended Transient Mid-term Stability Program \\
{\bf GE}      & General Electric \\
{\bf HPF}     & High Performance Fortran \\
{\bf IBM}     & International Business Machines \\
{\bf IST}     & School of Information Studies \\
{\bf MAPS}    & Multi-Area Production System \\
{\bf MW}      & Megawatts \\
{\bf NMPC}    & Niagara Mohawk Power Corporation \\
{\bf NPAC}    & Northeast Parallel Architectures Center \\
{\bf NUG}     & Non-Utility Generator \\
{\bf NYPP}    & New York Power Pool \\
{\bf NYPA}    & New York Power Authority \\
{\bf ODE}     & Ordinary Differential Equation \\
{\bf OLTSP}   & On-Line Transient Stability Program \\
{\bf PTI}     & Power Technology Inc. \\
{\bf SGI}     & Silicon Graphics Inc. \\
{\bf STAG}    & ??? \\
{\bf SVC}     & ??? \\
{\bf TASC}    & The Analytic Sciences Corporation \\
{\bf VAR}     & ??? \\
\end{tabular}



\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\addcontentsline{toc}{section}{Executive Summary}
\section*{Executive Summary}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MOTIVATION

OVERVIEW 

WHAT WE HAVE ACHIEVED 

WHAT WE ARE DOING / WILL DO BY THE END OF PHASE 2 OF THE PROJECT

ADD DIAGRAMS, PICTURES FIGURES, TABLES TO MAKE THINGS MORE UNDERSTANDABLE??


\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:Intro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsection{Project Objectives}
\label{sec:Objectives} 

The goal of this project is to evaluate the possible use of high-performance
computing, and in particular parallel processing, in the utility industry.
This is being done in two ways: 
\begin{enumerate}
\item A general survey of utility applications (in particular those used
      at NMPC) to identify those for which parallel processing could be
      effectively and beneficially utilized.
\item A detailed study of two particular applications, unit commitment and 
      transient stability analysis, and the development of parallel algorithms
      for these two applications.
\end{enumerate}

The main focus of this project 
is the development of a parallel transient stability program and associated
expert system which can be used on-line in a control room 
Energy Management System (EMS) to set 
dynamic transient stability limits.
This is perceived as potentially being of substantial benefit to NMPC,
in terms of cost savings and improved network reliability.
The techniques and algorithms used to develop such a system should also
be applicable to other areas of power systems simulation.

\begin{description}
\item [] Phase 1 of this project was a study of the feasibility of using
      parallel processing and expert systems for these applications. 
%
\item [] Phase 2 of the project is to develop and benchmark prototype parallel 
      programs and expert systems software.
%
\item [] Phase 3 involves developing, documenting and testing production-grade
      software which can be used at NMPC.
\end{description}


\subsection{Overview of the Report}
\label{sec:Overview} 

In this report we present the results of Phase 1 of this project, 
the feasibility study, and also discuss some of the work planned for
the future.

In section~\ref{sec:UnitComm} we study two different optimization methods for 
solving the unit commitment problem, and the use of parallel processing
to find better solutions in a shorter time.

In the following sections we examine the use of parallel processing to 
improve the performance of the security model for the on-line EMS used
at NMPC.

Section~\ref{sec:ContEval} investigates the use of parallel processing to 
improve contingency evaluation simulations.

Section~\ref{sec:TSA} examines methods for speeding-up transient
stability analysis, including improved sequential algorithms and the
requirements for developing a parallel transient stability program.

Section~\ref{sec:Sparse} presents our research on the development of
parallel algorithms for solving sparse matrix equations, 
which are a key element of transient stability analysis 
and many other power systems applications.

Section~\ref{sec:ExpertSystem} discusses the development of an expert system
for analyzing transient stability data and setting stability limits,
for both on-line simulations and off-line studies.

Section~\ref{sec:OLTS} integrates the results from the previous four sections
and discusses the feasibility of an on-line system for calculating
transient stability limits for use within the EMS security model.
A similar EPRI project to develop an on-line program for dynamic security 
assessment is also discussed.

Section~\ref{sec:Parallel} gives an overview of the relevant issues involved
in using parallel computers, and in particular the identification of
the best parallel hardware, languages, software, and algorithms for
power utility applications.

Section~\ref{sec:Survey} gives preliminary results on a survey of 
utility applications which would benefit from parallel processing.
This survey concentrates on the main algorithms used in these 
applications, and identifies whether efficient parallel algorithms 
already exist for these applications, or would need to be developed.



\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Unit Commitment}
\label{sec:UnitComm} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Utilities have available many generators with different generation
and transmission capacities and costs. 
Generators can be shut down for scheduled 
maintenance or unscheduled problems, or when there is not enough demand
for them to run. The demand (or load) varies substantially 
through daily and yearly cycles. However generators often take hours
to shut down or bring on-line, so it is necessary to start up generators
well in advance of when they are needed. The start-up and shut-down costs
must also be considered. Unit commitment is the problem of scheduling 
which generators to commit to power production in order to minimize the 
cost of meeting the forecast load for a specified time period, which 
could be a day, a week, or a year.

The cost is calculated by first doing a simplified economic dispatch
at each time interval, which determines the cost of distributing the 
load among the committed generators in the most economic way.
Generator start-up and shut-down costs, as well as any other costs, are
also added. It would also be possible to factor in other associated
``costs'' which may not be purely economic, for example environmental
effects, thus including environmental dispatch.

Unit commitment is a difficult combinatorial optimization problem. 
Suppose we wish to solve the unit commitment problem for one day of 
operation of a utility with 100 generators, given hourly load forecasts. 
The solution must specify the state of all 100 generators for each of 
the 24 one-hour time intervals, so the number of variables in the problem
is 2400. For each variable we must specify whether the generator is on or off,
so the number of possible solutions is $2^{2400}$, an unimaginably large 
number. It is clearly not possible to find the solution by an exhaustive 
search of all possible states.

The same is true for all hard optimization problems. Many different methods
(generally referred to as {\it heuristics}) have therefore been devised
to solve these problems by greatly reducing the number of possible
solutions which are checked (the {\it search space}).
Since these methods do not search the whole solution space, it is usually
not guaranteed that the result produced by the heuristic optimization
algorithm is the best possible (i.e. minimum cost) solution.
However for most optimization problems it is not necessary to find the 
optimal solution -- a good near-optimal solution is quite sufficient. 
This is certainly the case for the unit commitment problem, since both the
the load forecast and the pricing models used will only be approximate.
Our goal is therefore to identify a heuristic algorithm which will give 
a good near-optimal solution to the problem in a specified amount of 
computer time.
A standard approach to solving the unit commitment problem is the use
of prioritized lists, dynamic programming technique, or a combination
of the twos~\cite{WoodWoll}.

Parallel algorithms for unit commitment have previously been studied at
Syracuse University by Professor Yehuda Wallach and David Koester.
They took a simple unit commitment program and developed a parallel 
program for a shared memory multiprocessor machine by implementing a 
parallel dynamic programming algorithm~\cite{Koester90}.
The aim of this section of the project was to improve on this work by porting
their program to a distributed memory parallel computer, which is a more
widely used computer architecture (see section~\ref{sec:ParComp}),
and to investigate the performance of an alternative optimization technique, 
simulated annealing, for this problem.


\subsection{Dynamic Programming Approach}
\label{sec:DynProg} 

A standard method of approaching the unit commitment problem is the use
of dynamic programming optimization techniques~\cite{WoodWoll}.
Dynamic programming involves a search at each time interval,
which is done over a list of the best solutions up to the previous time 
interval.  %???
This dramatically restricts the search space to a
manageable subset of the total parameter space.

This method can also incorporate some basic knowledge of the system,
such as the use of a prioritized list of generators, with a higher priority 
given to the generators which produce the cheapest power.
This greatly reduces the search space, and makes it feasible to 
tackle even quite large problems. 
Simple prioritized-list methods are often used by utilities for unit
commitment. This is what is currently done at NMPC 
(see section~\ref{sec:NimoUnitComm}). 
However combining these with dynamic programming techniques will generally
give better solutions for large problems.
Other techniques, such as mixed integer-linear programming, could also
be used for this problem.

% DIAGRAM AS PER WOOD & WOLL ??


\subsection{Parallel Dynamic Programming}
\label{sec:ParDynProg} 

The search space used in the dynamic programming algorithm can easily be
distributed over processors, so the search can be done in parallel. 
Substantial communication is required, since each processor requires 
data from every other processor in order to search the entire space.
However this is {\it block communication}, that is, the data can be sent all
at once in a single large block, rather than one piece at a time in a lot of
small blocks, which is often required in parallel algorithms for other 
problems. Block communication is much more efficient for a distributed
memory parallel computer, since it minimizes the {\it latency}, or the 
overhead required to initiate the sending of a message. A simple analogy
would be posting a long letter -- it is much more efficient to post the 
whole letter in a single envelope, than to post each page individually in
separate envelopes. 

Results for the parallel dynamic programming algorithm on shared memory
parallel computers have shown good speed-ups for small test 
problems~\cite{Koester90}, and the efficiency of this method should 
improve with larger problem sizes and more complex cost functions,
since this would increase the ratio of work (calculation) to overhead
(communication).

We are in the process of porting this program to a distributed memory 
parallel computer. Since the algorithm can use block communication,
it should also give good speed-ups on a distributed memory architecture.



\subsection{Simulated Annealing Approach}
\label{sec:SimAnneal} 

Simulated annealing is very general optimization technique which has
been applied to many different problems since it was first introduced
10 years ago~\cite{SA}. The main advantage of this technique is that
it is very simple to apply to any problem. It generally gives quite
good results, which can be improved by using more computer time.

Simulated annealing is based on the simple optimization technique of
{\it hill climbing}. To find the highest peak in a mountain range, 
the simplest approach would be to always walk in a direction with an
upward slope, which will eventually lead you to the top of a mountain. 
However it may not be the {\it highest} mountain in the range. You may
have found only a {\it local} maximum, not a {\it global} maximum.
A simple hill climbing algorithm will rarely find the optimal value.
Avoiding getting trapped in a local maximum (or minimum, if we are
looking for the smallest value rather than the largest) is a major 
difficulty with hard optimization problems. Simulated annealing 
avoids this problem by occasionally allowing you to go in a direction with
a downward slope, so that you are free to explore other peaks. 

Simulated annealing works in an analogous way to the annealing of a 
metal. When a hot metal is cooled very slowly ({\it annealed}), 
the atoms form into a stable crystalline lattice which is very strong,
since it is a global minimum of the energy. However
if the hot metal is cooled very quickly ({\it quenched}), the atoms only 
have time to form a semi-stable lattice with many defects,
which is usually brittle and not as strong, since it is only a local 
minimum of the energy. Quenching is thus the equivalent of hill climbing.
Simulated annealing works by regarding the cost function of the problem
as being like the energy in real annealing, 
and introducing a parameter which is like temperature.
At high temperatures there is a large probability of making changes that 
increase the energy (cost function), whereas at low temperature only 
changes that reduce the cost are allowed. Simulated annealing starts 
with a high temperature value, and then slowly reduces it to zero,
at which point the cost function should be nearly optimal.

There has been quite a lot of work at NPAC over the past year to 
develop an improved simulated annealing method, and in particular
to develop good parallel algorithms. 
This work is aimed at developing a general optimization program which 
can be applied to any hard optimization problem.
We are currently working on applying this program to unit commitment, 
and it may also be applicable to other utility optimization problems.

We have also been investigating other work in this area, in particular
public domain general-purpose simulated annealing programs. 
We are currently working with the Adaptive Simulated Annealing (ASA)
program of Ingber~\cite{Ingber89}. 
Our goal is to use this well established program to get initial results,
and then use the sequential program we have developed to see if it can 
do better. Finally, we plan to use the parallel simulated annealing
algorithms to try to improve the speed of the program and/or the quality of
the results for unit commitment.

Simulated annealing requires the specification of 3 things:
%
\begin{enumerate}
\item How to vary the parameters (the {\it update})
\item How to change the temperature (the {\it cooling schedule})
\item The function to be minimized (the {\it cost function})
\end{enumerate}
%
For the unit commitment problem, the parameters are the state (on or off)
of every generator at every time interval. We will call a particular set of
parameters the {\it configuration} of the system. We can easily specify
a simple method for changing the parameters, by just reversing
the state of a given generator at a given time period. The particular 
generator and time period can just be chosen at random.
More sophisticated updates are possible and would likely improve
performance, but for simplicity we will start with this approach. 
The variation of the parameters is also affected by {\it constraints}. 
For example, some generators will be out of service and must therefore be
fixed to be off,
whereas some generators may have {\it must run} status, meaning that for 
some reason (often economic) they must always be on. Of course there is 
also the constraint that the total capacity of all the generators which are
on must exceed the forecast load, and that the load must be satisfied 
without exceeding the transmission capacity of the network.

Our goal is to rapidly find a configuration with a near-optimal cost
function. Simulated annealing does this in the following way. 
An initial configuration is specified, which can have any parameter
values, as long as the constraints are satisfied. For example, we could
set all generators in service to be on. An initial (large) temperature
value is specified. The temperature is then slowly decreased to zero,
following the specified cooling schedule, and one or more updates are
done for each new temperature value. 
In simulated annealing, a randomly chosen update (change of the 
parameters) is {\it proposed}. If the update leads to a reduction in 
the cost function, it is {\it accepted} (i.e. the change is made).
If the cost function is increased, it is only accepted with a certain
probability, which depends on the temperature. As the temperature goes
to zero, this probability decreases, so that eventually only changes which
decrease the cost are accepted. 

The only variation between different simulated annealing algorithms is
in the specification of the cooling schedule. The update and the cost
function are specified by the user for the particular problem to be
solved. In many simulated annealing programs, the user also has to 
specify the cooling schedule, which has to be tuned by trial and error
for the particular problem. Cooling too fast will give poor results, 
while cooling very slowly will give very good results but take far too
much computer time. Choosing a good cooling schedule can be very difficult,
but is crucial to obtain good performance from the annealing algorithm.

More effective simulated annealing programs use a dynamic or adaptive
cooling schedule, which uses measured data from the program as it 
is running to work out the next temperature change to be made.
The algorithm we have developed uses an improved method for finding the
adaptive schedule, that should give a cooling schedule which is close
to optimal for the particular problem being considered~\cite{PaulAnneal}. 
This is expected to greatly improve the running time and the final result.
We have implemented a preliminary version of our method, however there
are a number of improvements to this program which could still be made.
We are therefore also using a more mature public domain simulated 
annealing program (Ingber's ASA program) for our unit commitment studies.

We are currently creating subroutines for the update, the cost function,
and the constraints, in order to interface with the simulated annealing 
programs. These are based on our unit commitment program that uses
dynamic programming.

% RESULTS?? or estimate of when we will get them?


\subsection{Parallel Simulated Annealing}
\label{sec:ParSimAnneal} 

We have also developed a parallel version of our general-purpose simulated 
annealing algorithm which should further improve the cooling schedule and 
the quality of the result.
The parallelism is simply obtained by running a different annealing 
on each processor. Since simulated annealing is a probabilistic 
procedure, each processor will produce a different result, and the 
best result from all processors can be selected. Running multiple
simulations also helps in determining the adaptive schedule. 
The biggest problem with adaptive annealing algorithms is that there
is not enough data to accurately determine a good cooling schedule.
Running multiple simulations gives much more data, which can improve
this calculation. 

This method should give little speed-up over the sequential algorithm, 
but improve the quality of the results. Another goal would be to use 
parallelism to get the same quality result in a shorter time.
Some speed-up should be achievable by using this method with a
faster cooling schedule. However large speed-ups could only be 
achieved by parallelizing the computation of the cost function, 
which takes up most of the time in an annealing run. 
This is a more difficult task that is under investigation.

The parallel simulated annealing program, and the distributed memory
parallel program using dynamic programming, are both written in C 
using Express, a portable message passing language \cite{Express}.



\subsection{Unit Commitment at NMPC}
\label{sec:NimoUnitComm} 

Danielle Dylan's group in Operations does day-ahead unit commitment
for NMPC. They construct a 24 hour load forecast with hourly time periods,
based on current conditions (weather forecasts, etc.) and historical data.
They must follow NYPP procedures which require them to show that the 
committed units can meet peak demand plus reserve.
The amount of reserve is determined by NYPP, based on excess 
capacity required to cope with the most severe contingencies 
(failures of generators and transmission lines) in the whole region,
taking into account transmission stability limits.
The additional reserve requirement is usually 300-500 MW.

Costs for unit commitment are figured by using tables of average
generation costs, plus incremental costs (additional costs for each
additional MW), for each generator. 
The generators are prioritized in order of increasing cost.
Start-up and shut-down costs are also considered. 

Most NMPC generators are not ``cyclable'', that is, it takes a lot
of time and money to commit or decommit them. They are therefore
kept on-line for as long as possible, and usually only decommitted
on weekends, when demand is lower. Some units are taken off at night
and brought on in the morning. There are only two ``marginal''
generating units, Albany and Oswego, which will often need to be
committed on a daily basis. The commission of these units will also
depend on the ability to buy or sell excess power at favorable rates
to neighboring utilities or NUGs. This is factored in to the unit
commitment calculation.

The current practice of unit commitment at NMPC is to use a ``lambda
sheet'', basically a spreadsheet consisting of a prioritized list of
generators. The few most expensive generators are committed or
decommitted depending on the load forecast and the buying or selling 
of power. This is currently done manually.
Since most units are not cyclable, they are constrained to ``must-run''
status, which greatly reduces the number of possibilities to be considered,
and allows the task to be handled manually. However the process is clearly
amenable to automation.

EPRI have developed a unit commitment program called DYNAMICS, and
NMPC is a test utility for this project. This program handles combined
unit commitment and economic dispatch, using the same dynamic programming
method as we have used for unit commitment, and Lagrangian relaxation for 
economic dispatch~\cite{WoodWoll}.
It also handles emissions (environmental dispatch) and transactions
(buying and selling of power). This program runs the NMPC day-ahead 
unit commitment problem in about 1 minute on an IBM RS/6000 workstation,
and can handle week-ahead problems for the largest utilities in the
country (a couple of hundred generating units) in less than 5 minutes.
The program results have been studied at NMPC by ``after-the-fact''
analysis, and the unit commitment group are very happy with its performance. 
The main problem in testing DYNAMICS at NMPC was in integrating the 
data from the Energy Management System.

Since the new EPRI unit commitment program is so fast and effective, 
it seems that implementing a parallel unit commitment program would not
be of benefit to NMPC. 
However studies are being made on the the possible use of pool-wide unit
commitment, which is expected to lead to overall savings. Unit commitment 
and economic dispatch for the entire power pool is a much more complicated
problem, and one which would be done on-line as well as day-ahead
(economic dispatch is revised by NYPP on timescales of the order of a minute).
A faster unit commitment program may be of benefit to NYPP in this case.
We will investigate this possibility in Phase 2.
The parallel algorithms developed in this project should be applicable 
to any unit commitment program used by NYPP. 

The focus of the work in this section is on parallel algorithms for 
simulated annealing and dynamic programming, which are general
optimization techniques. It is probable that these parallel algorithms
(or variations of them) could therefore be of use in other utility 
applications, for example the MAPS and ENGINEOUS programs from General
Electric, which use optimization techniques to try to find the most
cost-effective network configurations for planning studies.
These are discussed in section~\ref{sec:ENGINEOUS}.




\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Contingency Evaluation}
\label{sec:ContEval} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

An important aspect of the operation of an electrical power network
is the security model, which simulates the response of the network to 
disturbances such as component failure. 
These simulations allow control room operators to identify
and remedy possible sources of cascading failures,
which could cause a disruption of power generation or transmission. 

We have talked at some length with William Merritt of the Operations 
group at NMPC, in order to clarify the use of the Energy Management 
System (EMS) security model software used by the Operations center, 
and to identify the areas where Operations would like to improve upon the 
current programs and practices.

The security model is run every few ?? minutes using data on the
state of the network which is provided every few ?? seconds by
remote telemetry units positioned throughout the network. 
This raw data gives an over-determination of the variables required to
solve the electrical network equations. This is necessary since not all
telemetry units will be accurate, or even functioning. 
A {\it state estimator} program is first run to make a statistical 
analysis of the raw data to determine the most likely actual state of 
the network, and this data is used for all subsequent simulations.

The security model simulates the effect on network stability of a 
number of different {\it contingencies}, or possible faults.
The NMPC Operations group analyzes a list of around 200 contingencies.
Most of these contingencies would cause very little disruption to the network, 
so to reduce computation time a {\it contingency screening} is performed. 
This is a fast, approximate calculation which allows the contingencies 
to be ranked in order of increasing severity, and to dismiss those
which are very stable.
Around 30 of the most severe cases are evaluated using a much more 
accurate model. The results are compared to the system stability limits.
If this contingency evaluation shows a contingency to be
so severe that it could cause problems,
some preventive action is taken by the operator.

The Operations group at NMPC would like to be able to use a more accurate 
screening procedure; study extra contingencies, including the occurrence
of two contingencies at one time; and to run a contingency evaluation more
frequently, on the order of once a minute.
The use of parallel processing would allow these extra capabilities
to be achieved.


\subsection{Parallel Contingency Evaluation}
\label{sec:ParContEval} 

Contingency evaluation is easily handled in parallel, 
since all the contingencies are independent.
This type of problem is often described as being
{\it embarrassingly parallel}. The independent cases can be sent to
different processors and computed in parallel by just using the sequential
algorithm on each processor. 
The fast screening of all the contingencies can be handled in this way
with almost perfect speed-up, i.e. a speed-up of $N$ on $N$ processors.
The results can be sent to one of the processors for ranking, which is
so fast that it need not be done in parallel. The full contingency 
evaluation on the most severe contingencies will also have almost 
perfect speed-up, as long as the number of cases is greater than or 
equal to the number of processors, which should be the case.

The only communication required in this embarrassingly parallel method
is access to the data produced by the state estimator. This would need
to be broadcast to each processor once at the beginning of the calculations,
which should not be a large overhead, since it would then be used for 
screening and evaluation of many contingencies, as there will be many
more contingencies than processors. This method should thus work well 
not just on a parallel computer, where communication speed between 
processors is very fast, but even when the communications channels are 
fairly slow, for example on a network of workstations.

The EMS program used by NMPC was developed by STAG.  % Acronym??
It is currently run on IBM 3090 or ES/9000 mainframes.
Parts of the program are written in non-standard languages
(IBM assembler and JCL), so the current program would be extremely 
difficult to port to a parallel computer. 
However STAG are in the process of porting the EMS to an IBM RS/6000
workstation environment. This version of the program would be very simple
to port to a parallel computer, especially to the IBM SP-1, which is 
basically just an integrated collection of RS/6000 workstations
(see section~\ref{sec:ParComp}).

% BENEFITS OF WORKSTATION-BASED PROGRAM vs MAINFRAME.

Networks of high-performance workstations are rapidly taking over the
role previously held by mainframes, since they offer improved 
performance at a lower price, are cheaper to maintain, are based on
portable industry standards rather than proprietary software, and are 
generally more reliable, which is a critical concern for a utility
Operations control center. 
The NMPC Operations mainframe is up 99.6\% of the time, 
however Operations would like to improve upon this figure.
A workstation-based EMS would be much easier to interface with newly
developed application programs, such as the EPRI DYNAMICS program for 
unit commitment and economic dispatch (see section~\ref{sec:NimoUnitComm}),
and programs being developed for on-line transient stability analysis 
(see section \ref{sec:OLTS}).

Much of the STAG EMS software has already been ported to the IBM
workstation platform. 
A pre-release version of the contingency evaluation program 
is currently undergoing testing, and
is expected to be available in the second quarter of this year. 
% ???


\subsection{State Estimation}
\label{sec:StateEst} 

Before the contingency evaluation can be performed, the state
estimation program must be run. 
This makes a statistical analysis of the data from the telemetry units
to determine the most likely state of the network.
If the sequential program is too time-consuming, it may be necessary
to parallelize it in order to avoid a bottleneck in what might otherwise
be a near-real time security analysis program.
Currently the state estimation program runs on an IBM 3090
in about 2 minutes, and on an IBM ES/9000 in well under a minute.
On a high-performance IBM RS/6000 workstation this will probably only
require a few seconds, so parallelizing the state estimator will
probably not be necessary.

However a more accurate, complex and numerically intensive 
state estimation program could perhaps be implemented if parallel
processing were used to make it run fast enough.
For example, a dynamic state estimation program could be used,
which looks at correlations in data over time, and is able to extract
information on the degradation or malfunction of telemetry units
\cite{TASC,WoodWoll}.

% J.J. Allemong et al, A fast and reliable state estimation program
% for AEP's new control center, IEEE Transactions, Vol. PAS-101, 1982, p933.
% Sci Tech TK1.P77

% also Simoes-Costa, IEEE Transactions, Vol. PAS-100, 1981, p691.

State estimation, like many other power network computations, involves 
the solution of large sparse matrix equations\cite{WoodWoll}. 
The sparse matrix is representative of the power grid,
which is very sparsely connected.
We expect that the parallel sparse matrix solver that we have developed
(see section~\ref{sec:Sparse})
could also be used as the basis for a parallel state estimation algorithm.

We are planning to investigate this problem in more detail, however we 
have so far been unable to obtain access to documentation on the 
algorithms used in the STAG state estimation program.

A simpler use of parallelism to speed up the overall run-time for the 
security model would be to overlap state estimation with contingency
evaluation, using a ``leapfrog'' technique. One processor (or workstation) 
could be performing a state estimation on the latest telemetry data, 
while a separate group of processors could be independently running a 
contingency evaluation using the results of the previous state estimation.



\subsection{Dynamic Stability Limits}
\label{sec:DynStabLim} 

Contingency evaluation programs do not simulate the network dynamics 
(the transient effects caused by a fault) which can give rise to instability.
They give only the static, steady-state behavior of the network, 
by using approximations ({\it distribution factors}~\cite{WoodWoll})
to model the response of the network.

A more accurate contingency evaluation would utilize a full {\it transient 
stability analysis}, which calculates the dynamic response of the network
to a fault. This is currently not done, due to the large computation times
required for transient stability analysis. EPRI has begun a project to 
develop an on-line program for {\it dynamic security assessment} (DSA),
which will incorporate transient stability analysis into the security
model contingency evaluation. This is discussed in section~\ref{sec:DSA}.

Transient stability analysis is currently done off-line at NMPC and NYPP,
by planning engineers. It is used to set the component and interface
stability limits used by NMPC Operations. These are set by very infrequent
studies (every few months), using a ``worst-case'' scenario.
An on-line program to provide stability limits on a near real-time basis
using the current 
operating state of the power network would be of great value to Operations,
since it would reduce the amount of excess generation (spinning reserve)
required to satisfy security constraints, and increase safe transmission
capability to allow greater transfer of power from areas with lower 
generation costs.
William Merritt has estimated that NMPC could save as much as \$3.7M % ???
per year in this way,
and EPRI have given similar arguments for large savings as 
motivation for developing their DSA program~\cite{DSA}.

The Operations group would also benefit from the ability to use an
on-line transient stability program in study mode, providing quick answers
to their occasional requests for certain stability limits which currently 
take days or weeks to be processed by the Planning department.
The feasibility of an on-line program to set transient stability limits
is discussed in section~\ref{sec:OLTS}.


\subsection{Further Improvements}
\label{sec:ContEvalImp} 

As well as allowing improvement of current simulations by the Operations
group, high performance computing could also enable the development of
new tools to help control room operators.
One possible area for further study is the development of methods and
software to help an operator respond to a situation where the
occurrence of a contingency would violate stability limits.
The security model used by Operations flags two types of security problems:
% ???
%
\begin{itemize}
\item {\bf A Minor Limit Violation}, for which the operator would have
approximately 15 minutes to respond to the occurrence of the contingency
before overloaded components could fail, causing the system to become
unstable. The operator must devise a strategy to quickly bring the
system back to stable operation should the contingency occur.
\item {\bf A Major Limit Violation}, for which the system could become quickly
unstable should the contingency occur. In this case the operator must
take immediate action to alter the operating conditions to avoid this
possibility.
\end{itemize}
%
This process should be amenable to automation.
A program (probably including an expert system) could be used to 
simulate many possible scenarios in order to 
identify a workable strategy to deal with a stability problem. 
A program of this kind could also be used to deal with the
more critical case of an unforseen disturbance driving the system
rapidly toward instability, giving the operator very little time to 
react to control the problem.

If the security model were able to integrate with a
pricing model incorporating unit commitment and economic dispatch, an even
better approach would be to find not just a viable strategy for keeping
the system stable, but one that would minimizing incurred costs.
The EPRI DYNAMICS program discussed in section~\ref{sec:NimoUnitComm}
should be fast enough to be used on-line in this way. Multiple scenarios
could be identified for bringing the system into a stable operating mode,
and the cost of implementing each of these strategies could be calculated
in parallel, with a different strategy being analyzed on each processor.
The most cost-effective strategy could then be given to the operator,
or implemented automatically.



\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Transient Stability Analysis}
\label{sec:TSA}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Transient stability analysis involves the study of the electrical 
network and generator dynamics for several seconds following a disturbance. 
The dynamics of each generator in the system is described by multiple 
differential equations, while the coupling of the generators to the 
electrical network is described by algebraic equations. 
Together these lead to a system of differential-algebraic equations (DAEs) 
that must be solved simultaneously. Since a large system of DAEs needs to 
be solved at small time increments to minimize the errors in numerical 
integration of the differential equations, the transient stability analysis
is a highly compute intensive problem.

One of the main goals of this project is to develop an on-line 
program for transient stability analysis. 
This is viewed as a very desirable tool by power utilities, 
and is the subject of much current research, particularly by EPRI
(the EPRI project is discussed in section~\ref{sec:DSA}).

There are two problems to be overcome to create an on-line 
transient stability analysis program. 
The first problem is to substantially speed up the program to make 
results available within a reasonable time period. 
We believe that this goal can be achieved by the use of parallel
processing and possibly by developing improved algorithms for 
solving the transient stability equations.
We are pursuing both these avenues of research.
The second problem is to develop an on-line program to analyze the output 
of the transient stability program and classify the results as stable or
unstable. Our work in this area is addressed in section~\ref{sec:ExpertSystem}.

There are a number of commercial programs available for transient stability
analysis. Two of the most widely used programs are the EPRI 
Extended Transient Mid-term Stability Program (ETMSP), and the PTI PSS/E
program, which is used by NMPC and NYPP.

Transient stability calculations at NMPC currently take around 
20 minutes % ?? 
on the IBM mainframe % ??? 
for each contingency.
On an IBM RS/6000 this is reduced to around ??? minutes for the 
PTI PSS/E program.
The latest version (3.1) of ETMSP
runs a 2000 bus problem on an IBM RS/6000 model 580 in 
52 seconds. 


\subsection{Simulation of Power Network Dynamics}
\label{sec:Equations} 

An electrical power network is modeled by a set of ordinary differential
equations (ODEs) that describe the internal dynamics of the ``dynamic'' 
devices (like generators, SVCs, DC links, governors, etc.), and a set of 
non-linear algebraic equations that describe the transmission network 
through which the dynamic devices and the system loads are connected.
Transient stability analysis therefore involves the solution
of a set of combined differential-algebraic equations (DAEs). 
The detailed specification of these equations can be found in 
Refs. \cite{ETMSP}. % others??
Here we will give a simplified overview of the problem of formulating
and solving the DAEs.

A simple description of a DAE is ordinary differential equations combined
with algebraic equations. However there are many different ways of specifying
a DAE, which can give us information on how the DAE can be used and solved.
A system of DAEs can be formed by a model consisting of a system of ODEs 
which have some algebraic constraints:
\[
\begin{array}{lcl}
x' & = & f(x, y, t) \\
y & = & g(x, t) \\
\end{array}
\]
This type of system is called a {\it constrained variational problem}. 
The transient stability problem is of this type, since the equations 
describing power system dynamics can be written in the form
\[
\begin{array}{lcl}
x' & = & f(x, V, I) \\
I & = & Y(x) V \\
\end{array}
\] 
where $I$ is the current, $V$ the voltage, $Y$ is the admittance matrix
describing the network, and $x$ is a vector containing the variables 
describing the dynamics of the generators.


The main practical difference between a DAE and ODE is that a DAE can be
more difficult to solve. The degree of difficulty can be inferred from the
{\it index} of the DAE system.
A system of ODEs is classified as index 0.
Next we have index 1 DAEs, which can be subdivided into two classes,
depending on whether they are specified as a semi-explicit
\[
\begin{array}{lcll}
x' & = & f(x, y, t) \\
0 & = & g(x, y, t) \\
\end{array}
\]
or fully-implicit system
\[
\begin{array}{lcll}
0 & = & f(x, x', y, t) \\
0 & = & g(x, y, t) \\
\end{array}
\]  
The electrical transient stability problem thus consists of solving
semi-explicit index-one differential-algebraic equations.
There are higher index problems; however, we will not introduce them here,
since they are not believed to be applicable to the transient stability 
problem (although this is the subject of some debate).

Another factor which makes these equations difficult to solve is that they
are infinitely stiff. A {\it stiff} ODE describes a system in which
some of the variables change quite slowly, while others vary rapidly,
making the solution more difficult to compute numerically~\cite{DE}.
In the case of power system DAEs, the relative generator rotor angles 
vary comparatively slowly (governed by the ODEs) while the voltages vary 
instantaneously (governed by the algebraic equations). 

Commercial transient stability programs such as ETMSP and PSS/E
use the partitioned or decoupled method 
to obtain the time response of the system for the 
specified fault and switching condition. 
The basic idea is to alternate between numerical integration 
to solve the differential equations describing the generators, 
and an iterative Newton method to solve the non-linear 
algebraic equations describing the network \cite{ETMSP}.

Fig.~\ref{fig:TSA} gives an outline of the method used in ETMSP 
for solving the transient stability equations. 
The generator equations are integrated using either an explicit Runge-Kutta 
method or a predictor-corrector technique \cite{DE},
depending on the circumstances.
Both are standard techniques for numerical integration of differential
equations.
The network equations are solved at the end of the integration
(for next step), and also during the integration in the corrector steps.
In ETMSP, an explicit Euler method is used for the predictor, and the 
trapezoidal method is used as a corrector method.
The system is then solved using the Newton, 
or Very Dishonest Newton (VDHN), method
(this is explained in detail in Refs. \cite{ETMSP}). % others ??
Note that the network equations are solved multiple times for each 
timestep in the Runge-Kutta integration, and also once after computing
the initial conditions.


\begin{figure}
\label{fig:TSA}
\hspace*{0.25in} DO ( Loop for each network condition or contingency )\\
\hspace*{0.5in}  Update the admittance matrix to reflect the system operation \\
\hspace*{0.5in}  DO ( Loop for each integration timestep in the time interval studied ) \\
\hspace*{0.75in} Check whether appropriate to change to predictor-corrector\\
\hspace*{0.75in} IF ( not appropriate to change ) \\
\hspace*{1.0in}    Use Runge-Kutta as starter method \\
\hspace*{0.75in}  ELSE \\
\hspace*{1in}      {\it Use predictor-corrector method} \\
\hspace*{1in}      Calculate derivatives of state variables \\
\hspace*{1in}      Use the predictor portion of method first \\
\hspace*{1.25in}   WHILE ( integration error $>$ tolerance limit ) \\
\hspace*{1.5in}      Get network solution \\
\hspace*{1.5in}      Calculate derivatives of state variables  \\
\hspace*{1.5in}      Use corrector portion of numerical integration method \\
\hspace*{1.25in}   END WHILE\\
\hspace*{0.75in}  END IF \\
\hspace*{0.75in}  Get network solution in preparation for next step\\
\hspace*{0.75in}  Save solution for next step \\
\hspace*{0.75in}  Write any intermediate output \\
\hspace*{0.75in}  Increment the time step\\
\hspace*{0.75in}  IF ( there are any changes to the system topology ) THEN \\
\hspace*{1in}   Update the admittance matrix to reflect the system operation \\
\hspace*{0.5in}   END DO\\
\hspace*{0.25in} END DO \\
%
\caption{An outline of the method of solution of the network 
dynamics by the EPRI ETMSP transient stability program.}
\end{figure}


The trapezoidal rule is a stiff accurate, order 2 method~\cite{SODEII}.  
Stiff accurate implies that the method will work for stiff systems
(such as transient stability problems),
since the size of the integration timestep $h$ used in solving the problem 
has a reduced impact on the accuracy of the numerical result. 
Order 2 means that the order of convergence (or accuracy) is proportional
to $h^{2}$. The explicit Euler method is order 1, and is not stiff accurate.
The higher the order of the method, the larger the timestep that can be used.
In addition, the Euler method and trapezoidal rule are analytically stable
(or {\it A-stable}),
which means the analytical behavior of the numerical method will not
adversely affect the accuracy of the solution.
There is also the problem of {\it numerical stability}, which concerns the
possibility that numerical round-off error due to finite-precision 
calculations on the computer can become amplified and lead to large errors 
in the solution. This problem does not appear to have been carefully 
addressed in current solution methods. 

The PTI PSS/E transient stability program uses an explicit modified Euler
solution method, which is order 2 but is not stiff accurate. 
In other words, it requires a very small step size to maintain the desired 
accuracy when solving stiff problems such as transient stability~\cite{DE}.

We have been exploring the possibility of using improved ODE solvers to
solve the generator equations.
Prominent among these are the solvers present in ODEPACK \cite{hindmarsh},
such as LSODI (implicit method ODE solver), LSODE 
(explicit method ODE solver), and LSODES (ODE solver with sparse matrix 
solving capabilities).
LSODE has been used in the latest version of the EPRI program DIRECT,
which uses direct energy methods rather than time-domain methods to
study transient stability.
An implicit method such as LSODI would be much better for solving the
time-domain transient stability DAEs (this is discussed further in
section~\ref{sec:DAE}).

We are also investigating some improved methods for decoupling DAEs
into ODEs and algebraic equations~\cite{Brenan}, to see if they could 
feasibly be applied to the transient stability problem.




\subsection{DAE Solvers}
\label{sec:DAE} 

The power systems community has traditionally solved the transient stability 
problem numerically by approximating it as a decoupled system of differential
and algebraic equations, as described in the previous section. 
This method of reduction of a large DAE system to a system of explicit ODEs 
reduces the accuracy of the solution and
prevents the full exploitation of the system structure in the solution
technique \cite{Brenan}. 
In the last few years many methods have been developed which solve the
full coupled system of DAEs simultaneously, which is believed to be a 
more efficient and accurate approach. 

Several sequential public-domain DAE solvers exist that use various
numerical methods for their solution. Prominent among these techniques are:
\begin{itemize}
\item Extrapolation methods \cite{Deuf,Hairer}
\item Implicit Runge-Kutta methods \cite{Hairer}
\item Backward Differentiation Formula (BDF) methods 
      \cite{Brenan,nmg-Hi82,nmg-h83,nmg-Hi85b}
\item Rosenbrock methods \cite{Hairer}
\end{itemize}
%
A brief description of some available DAE solvers 
and their applicability to the transient stability problem is given
below (a more detailed discussion can be found in \cite{KamalaSLC}): 
%
\begin{itemize}
\item {\bf LIMEX} \cite{Deuf} uses an extrapolation technique. 
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. Because of its ability to handle stiff differential equations and 
frequent discontinuities, LIMEX is an excellent choice for power 
system DAEs.
%
\item {\bf SODEX} \cite{Hairer} and {\bf SEULEX} \cite{Hairer} also use 
extrapolation techniques for the solution of DAEs.
%
\item {\bf RADAU5} \cite{Hairer} uses an implicit Runge-Kutta 
method of order five. RADAU5 is designed to solve index one, two, 
and three systems. It is possible that more sophisticated modeling of
the transient stability problem would result in DAEs of index higher than
one, in which case RADAU5 would be a good choice.
%
\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 to advance 
solution from one time step to the next. 
DASSL has the ability to do frequent step-size adjustments which
could be very useful for power systems problems.
DASSL can solve DAEs of index zero and one, and semi-explicit index two 
systems. 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 NPAC to solve the transient stability problem 
using a parallel implementation of DASSL (see section~\ref{sec:ParDAE}). 
%
\item {\bf LSODI} \cite{nmg-Hi82,nmg-h83,nmg-Hi85b} uses BDFs. It uses 
a fixed-coefficient implementation of BDFs (which can be implemented
very efficiently) to solve linearly implicit systems of DAEs.
LSODI solves stiff and non-stiff ODEs (DAEs of index zero), but if the 
coefficient matrix is singular, it can also solve the resulting index 
one DAE system. IS THIS THE CASE FOR THE TSA PROBLEM??? 
Because of its ability to handle stiff ODEs and frequent discontinuities, 
LSODI presents a viable alternative to the other solvers discussed for
power system problems.
%
\item {\bf RODAS} \cite{Hairer} uses an embedded Rosenbrock technique of 
order four. 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}
%
LIMEX, DASSL (or DASRT), and possibly LSODI appear to be the DAE solvers
which are best suited to the transient stability problem.
These have been used extensively for other problems, and have proven to
be very efficient. We are initially focussing on using DASSL, however we
hope to compare this with other methods in order to find the best solver
for this problem.


We believe the use of these DAE solvers is a promising new approach to solving
the transient stability problem, which offers the possibility of improved
accuracy and reduced solution time.
More efficient and accurate DAE solvers would be particularly effective
for long-term stability studies, where the time interval studied is on the 
order of minutes rather than seconds.
Currently the time taken to run these simulations is extremely large.
Good DAE solvers would allow longer integration timesteps without sacrificing
numerical accuracy, which is very important for long-term stability 
simulations.

However, since the approach of using a DAE solver is so different from
traditional methods used for solving power system equations, 
it has proved to be much more difficult than expected to integrate a DAE
solver into a real production code such as ETMSP.
In ETMSP, setting up the equations to be solved, in particular for the 
many different types of generator models used, is ``hard-wired'' to the
current solution technique used, in order to optimize the program for speed.
To use alternate methods such as a DAE solver would require substantial
additional effort in setting up the system equations in a general form
which could be used by any solution method. 
This may require rewriting sections of the program for each of the hundreds
of combinations of generators, exciters, etc. which are handled by the 
program. However we are hoping that we can implement a simpler workaround
to this problem, by utilizing the input format of the ``user defined model''
provided for in ETMSP~\cite{ETMSP}.
We are currently investigating this approach, and hope to implement 
a prototype DAE solver during Phase 2 of the project.



\subsection{Parallel Processing for Transient Stability Analysis}
\label{sec:ParallelTSA} 

There has been quite a lot of previous work attempting to use
parallel processing to speed up transient stability analysis, % \cite{??}
including research by PTI and 
EPRI~\cite{Bose93}. % \cite{previous EPRI work} ??
This work has so far not produced good results.
We will give a short overview of this work, briefly discuss some of the 
reasons for the poor results, and outline our plan for parallelizing the
EPRI ETMSP program (the comments on parallelization would mostly apply 
to the PTI PSS/E program also).

Initial attempts to utilize parallel processing consisted of using 
a parallelizing compiler to create a parallel program from the 
sequential Fortran program (either PSS/E or ETMSP). 
Because of the structure of these programs, this approach is doomed
to failure, and indeed the parallel programs constructed in this way
were invariably no faster (and sometimes slower) than the sequential
programs.
This is because parallelizing compilers are quite rudimentary, and 
work mainly by identifying loops which have little or no data 
dependencies, so that every instance of the loop is independent 
(or almost independent) of any other. These loops can therefore be
parallelized, with different instances of the loop being assigned to
different processors.
This approach will only work for problems which are very regular,
have regular data structures, and for which the programs are written
in a structured way with simple loops which make the regularity
apparent to the compiler.
Commercial transient stability programs do not satisfy these 
requirements.


The decoupling methods used in ETMSP do allow for some parallelism.
The simplest is in the solution of the generator equations. 
These can be solved independently for each generator, making the solution
embarrassingly parallel, i.e. the sequential ODE solver can be used for a
different generator (or set of generators) on each processor.
As long as the number of generators is larger than the number of processors,
which would certainly be the case for any realistic problem,
this method would give almost perfect speed-up for solving the generator 
equations, that is, the speed-up would be almost equal to the number of
processors used.

The only problem with this scenario arises when there are data
dependencies among the generators. This occurs only in the following
two cases:
%
\begin{enumerate}
\item When controllers with feedback from other devices are added to 
      the system. This will give a coupling between the data for the 
      controller and the feedback device.
\item When the system consists of both low pressure and high pressure
      governors, there is a coupling between the two.
\end{enumerate}
%
Since case 1 is not currently implemented in ETMSP (although there are
plans to add it in a future release), it does not pose a problem to our
initial parallel implementation. 
%
Case 2 is handled by the ETMSP code, however it
can be worked around by having the data for low pressure and
high pressure governors on the same processor, so that they can easily
exchange data, in the same way as in the sequential program. 
If there are a large number of these devices, they may have to be
distributed over different processors, in which case data would need to be 
communicated between processors. This would still allow parallelism, 
but the efficiency of the program would be slightly reduced due to the 
communications overhead.
Implementation of case 1 could be worked around in a similar fashion.

Since the generator equations are completely independent, 
they should be very simple to parallelize. 
However as noted in the report of EPRI research project 3103-1 on parallel 
methods for transient stability analysis~\cite{Bose93}, the EPRI ETMSP 
program is written in such a way that this independence is obscured.
There is not a ``clean'' loop over generators, or even a call
to a separate subroutine to solve the generator equations.
The solution of the generators is entwined into the overall solution
in a rather messy way.
Restructuring this part of the program would greatly simplify the task
of parallelizing the program, or indeed improving it in any way,
for example by using improved algorithms for solving the ODEs.
One of the main recommendations of EPRI project 3103-1 was that future 
versions of EPRI programs should be written in a simpler, more structured
and modular fashion, with clean loops. 

Other possible approach would be to use of a parallel ODE solver,
which may also allow for a reduction in the error introduced in 
decoupling the equations~\cite{Bose93}.

In ETMSP, about $75\%$ of the computational time goes into solving 
the generator equations, which can be efficiently parallelized, 
and about $25\%$ into solving the network equations. 
So even if massive parallelism could be used to reduce the time to
solve the generator equations to almost zero, the total time would only
be reduced by a factor of 4, since the network equations still have to 
be solved, which originally took $1 \over 4$ of the time.
Only if these can also be solved efficiently in parallel can we get a
speed-up of more than 4 (this maximum possible speed-up will be dependent 
on the transient stability program used and the network to be solved, 
but the idea is the same in each case, and the number will be of this order).

The network equations are solved in the standard way, by using an iterative
Newton method. At each iteration, a linearized version of the
equations if formed, which are solved using a sparse linear solver.
A parallel network solution thus requires a parallel sparse matrix solver,
which is a very difficult problem, especially for power systems applications
where the matrices are very sparse.
Previous work on parallel sparse matrix solvers for power systems 
has given very poor results, with speed-ups typically in the range 2-5. 
This has been the bottleneck in using parallel computers to speed up
power systems applications.
To our knowledge, parallel transient stability programs have never achieved 
speed-ups of greater than 10, no matter how many processors were used.

In section~\ref{sec:Sparse} we describe the development of an improved
parallel sparse matrix solver which exploits the structure of power systems
matrices.
Preliminary results indicate that this solver should provide better 
performance than previous methods, with speed-ups of around 10 on 16 
processors appearing possible. % ???.
If we can get perfect speed-up in solving the generator equations, 
this would produce a total speed-up for the whole ETMSP
program of 14 on 16 processors, % ???
which is very close to optimal.

This speed-up would be useful not just for an on-line program,
but also for off-line transient stability studies, and for integrated 
planning tools which use transient stability analysis, for example
General Electric's project to use ENGINEOUS for power systems planning 
(see section~\ref{sec:ENGINEOUS}).



\subsection{Implementing a Parallel Transient Stability Program}
\label{sec:ParallelETMSP} 

We are currently working on developing a prototype parallel implementation
of ETMSP. So far we have studied the parallelization of the generator
equations, and the use of parallel sparse matrix routines to solve the
network equations.
This has required us to spend a substantial amount of time familiarizing
ourselves with the ETMSP code, in order to understand exactly how the 
currently used ODE and matrix solvers interface with the rest of the program, 
how they are initialized and called, and what changes they make to the state
of the program.
As noted previously, the program is not written in a simple modular
fashion, so understanding the detailed structure of the program is not
an easy task.

We are currently writing software to interface parallel sparse matrix 
programs to ETMSP. We are planning to do this for a number of different
parallel sparse matrix solvers, in order to compare their performance 
on a real power systems application.
These will including the solver we have developed, 
some general-purpose parallel sparse solvers \cite{CMSSL,Toolbox91},
and the programs being developed at the University of Maryland 
for an EPRI research project on parallel sparse matrix methods.

We have obtained some sample 1000, 2000 and 6000 bus data for the 
north-east US power network from EPRI, and will initially be benchmarking 
the parallel programs using this real data.
Larry Eng of NMPC has agreed to work on the NMPC generator data,
which is currently in PTI format, so that it will interface with ETMSP.

We will initially benchmark the parallel programs on the 32 processor
Thinking Machines CM-5 at NPAC, since the sparse matrix solver currently
runs on this machine. We are in the process of porting this to
Zipcode~\cite{Zipcode},
a portable parallel language which will run on any distributed memory
parallel computer. We can then also benchmark the performance of the 
parallel programs on other parallel computers such as the IBM SP-1.

We will begin the parallelization of the generator equations 
as soon as we have received the latest version (3.1) of ETMSP, 
which is now being shipped. 
The first task will be to restructure the current sequential code for 
solving the generator equations, so that it is in a modular form with a 
clean loop over generators.
Substantial effort will be involved to restructure the ETMSP program in 
this way. This restructuring has been done in the case of a research-level
transient stability code developed at Arizona State University, and also 
for an in-house production-level code used by Hydro Quebec.
This allowed parallelization of the generator equations in these programs, 
resulting in almost perfect speed-up for this part of the 
programs~\cite{Bose93}.
However this has not been attempted for PSS/E or ETMSP, which are much 
larger and more complex programs (ETMSP has approximately 200,000 lines
of Fortran).
David Wong, who wrote most of the ETMSP program, has estimated that the
restructuring of ETMSP could be done in about 6 months by someone with
a good understanding of the problem and the program.
We believe this is a reasonable timescale for this undertaking.

It would be beneficial to integrate the parallelization of ETMSP with
the overall development of this program by EPRI, i.e. to get some
official backing (and possibly co-funding) from EPRI for this project.



\subsection{Parallel DAE Solvers}
\label{sec:ParDAE} 

Sequential DAE solvers offer the promise of reduced simulation time and
improved accuracy over current solution methods for transient stability
simulations. However another good reason for studying the use of 
DAE solvers is that parallel implementations of some DAE solvers
already exist, such as Cdassl,
a concurrent version of DASSL~\cite{Toolbox93,SLCdassl}.
This has been used very effectively in simulating the operation of
a petrochemical plant~\cite{TonyThesis}.
DASSL~\cite{Num-Sol,DASSLSW,DASPKSW} is a software package specifically 
designed for solving index one DAEs. The features of this package are
variable convergence order, variable step size and a clean interface to
the linear sparse solver.

DASSL uses an implicit backward difference formula, and provides formulae
for convergence from order 1 to 5.  
The order 1 BDF is the Euler method and order 2 BDF is the implicit 
modified Euler method, both of which A-stable.  
(Unlike the {\it explicit} modified Euler method, the {\it implicit} 
modified Euler method is stiff accurate, since it has the same formulation 
as the trapezoidal rule.)
We expect a gain in computational speed from using formulae with even higher
orders of convergence, since a larger step size can be used to compute the 
result with the same accuracy.  The trade-off is that the formulae with
order greater than 2 are not A-stable \cite{SODEII}. 
DASSL has automatic control for shifting from a lower order formula to a
higher order formula, when the function exhibits stable behavior for 
higher order formula. However the stability problem needs to be considered 
very carefully, so using DASSL for model development requires an in-depth
knowledge of DASSL and a good working knowledge of the model. 
However, if the behavior of the model is well understood, then DASSL 
can be used to provide fast simulations.

DASSL also has an automatic variable step size control. Although
ETMSP only allows fixed step size, its users often use different step
sizes for different phases of the simulation.  For example, a larger step 
size can be used before the fault, a smaller step size just after the fault,
and a larger step size after the disturbance is expected to be settled.
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 time
history of the function indicates the function is slowly varying, then the
step size is increase until either the error is larger than specified or an
unstable result is computed.  Variable step size should provide a further
reduction in the run-time of a simulation.

In solving the DAEs, DASSL uses Newton's method to solve systems of
non-linear equations. Each Newton iteration requires the solution of
a system of linear equations, $Ax = b$. For the transient stability 
application, this requires calling a parallel sparse matrix solver.
Cdassl includes a general purpose parallel sparse matrix solver.
During the course of this project we have done a lot of work on improving
Cdassl and making it easier to interface with other programs, in particular
providing a well-defined interface to the linear system solver~\cite{SLCdassl}. 
This interface allows DASSL to call any parallel sparse matrix solver with 
minimum code changes. In particular, are planning to interface to the solver
we have developed for power systems applications 
(see section~\ref{sec:Sparse}).
It would also be possible to interface the parallel program to a DAE solver 
other than DASSL, if another method proved to be more effective for power
systems problems.

In the original DASSL~\cite{DASSLSW}, direct methods using LU factorization 
were employed to solve the system of linear equations  
(see section~\ref{sec:Sparse} for a discussion of direct methods). 
However, the recent advances in semi-iterative linear solvers such as 
CGS, TFQMR and GMRES \cite{Templates}, 
have prompted DASSL's authors to incorporate the GMRES 
semi-iterative solver into DASSL. The resulting package, known as
DASPK~\cite{DASPKSW}, has proved to work faster than DASSL for certain
problems.
Cdassl can interface to both the direct and semi-iterative solvers.
We are planning to study the use of semi-iterative methods as a possible
improved solution technique for the transient stability problem.

We believe Cdassl promises to be the best method for solving the transient
stability problem accurately in parallel, however interfacing Cdassl into
a real production code such as ETMSP is quite problematic, as noted in 
section~\ref{sec:DAE}.




\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Parallel Sparse Matrix Algorithms}
\label{sec:Sparse} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Alvin's work on the ETMSP solvers. ??

\subsection{Future Work}
\label{sec:SparseNext} 

iterative and semi-iterative linear solvers




\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Expert System and User Interface}
\label{sec:ExpertSystem} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Currently transient stability analysis is done off-line by planning
engineers. There are two reasons for this. The first is the large
amount of computer time required for a transient stability simulation.
We have attempted to address this problem in the previous sections.
The second reason is that analyzing the output for a certain
contingency as being either transient stable or unstable is not a simple
task, and requires a trained engineer to interpret the simulation results.

In order to provide an on-line transient stability program, the second
problem must also be addressed, which requires that the procedure for
analyzing the output of the stability program be automated.
There are various ways of approaching this problem, and a number of
different methods are being studied in related EPRI projects for
on-line dynamic security assessment~\cite{DSA}. 
The approach we are taking is to develop an expert system which 
approaches the problem in the same way a human engineer would.

This expert system should be useful not only for on-line applications,
but also in helping the engineer in traditional off-line transient 
stability studies. Programs for transient stability analysis produce
tremendous amounts of graphical output which must be analyzed.
An expert system would be able to prioritize or filter this data,
for example by flagging results which show clear instability,
filtering out results which are clearly stable, and presenting the
engineer those results which are most useful for making a decision.
An improved user interface may also allow these results to be presented
in a way that would aid the engineer in interpretating the data.

Note that for both on-line and off-line use, the expert system 
must be conservative. 
It is not acceptable to classify as stable any results which may
actually be unstable.
This would not adversely affect an on-line system, since
conservative dynamic limits would generally 
still be lower than the current static limits.


Analysis of the output data is only one part of setting transient stability 
limits. The engineer must also choose appropriate contingencies, and 
adjust parameter values, in order to find the limiting value for a certain
component or interface, above which the network becomes unstable. 
This is another (more difficult) process which needs to be automated in
order to achieve our goal of developing an on-line system for setting
dynamic stability limits.
We are also beginning to develop an expert system for this process. 



\subsection{Expert Systems in the Electric Power Industry}

Artificial intelligence (AI) technology -- using expert system and neural 
network design -- has generated both a tremendous amount of interest 
as well as a proliferation of computer-based applications within the 
electric power industry.  Moreover, it would appear that this trend is 
likely to continue, at an accelerated pace, well into the future.  This fact 
was made readily apparent by a review of the current literature in the 
power industry and was further corroborated in conversations and 
interviews with industry experts during the course of Phase 1 research.  
While system security will remain as the primary criteria for power 
system operation and a significant hurdle in the implementation of 
automated decision systems (i.e. artificial intelligence technologies), 
new technical, organizational, and economic realities within the power 
utility industry are creating new challenges which must be addressed.  
It is likely that AI technology will play an important role in providing 
solutions to these challenges.

The perceived importance of expert system technology within the 
field was convincingly demonstrated at an EPRI-sponsored 
international conference and exhibition 
on Expert System Applications for the Electric Power Industry,
in Phoenix on December 8-10, 1993.  
Numerous software applications using AI technologies were presented 
during this conference.  Some of the existing and proposed expert 
system and neural network based applications included systems for:
\begin{itemize}
\item	fault location and detection
\item	equipment diagnostics (e.g. transformers, relays, etc.)
\item	communication system analysis
\item	dynamic security assessment
\item	transient stability analysis
\item	post-fault system restoration
\item	alarm analysis
\item	theft of service detection
\item	voltage/VAR control
\item	EMS software integration          
\end{itemize}

There are a variety of reasons for the enthusiasm generated by 
these real and potential applications of AI technology.  
They may be summarized as follows:
\begin{enumerate}
%
\item	The electric power industry has become more conscious of the need 
to maximize productivity in the face of budgetary constraints.
%
\item	The expertise of personnel has become an increasingly important 
commodity in the face of downsizing and retirement.  Expert systems 
are viewed as a means of capturing such expertise, especially for the 
training of new or less experienced employees.
%
\item	The Energy Policy Act of 1992 has raised concerns about the 
ramifications of a deregulated utility industry.  There appears to be 
an increasing awareness that competition will be a fact of life and 
key to survival in the future.  Expert systems are viewed as one 
component of an overall need to incorporate computerized 
technology in order to remain competitive.
%
\item	Expert systems are also conceived of as potentially more efficient 
in dealing with more stringent environmental standards.
%
\item	Planners and, especially, operators are confronted with an 
increasingly complex and voluminous amount of data in today's 
energy management systems.  Expert systems and improved graphical 
user interfaces may be able to aid planners and operators in making 
more accurate, timely and cost-efficient decisions.  
%
\end{enumerate}

EPRI are developing an expert system to integrate all the 
separate programs involved in DSA, and to analyze the results and 
provide a simple user interface to the operator.
A similar expert system and user interface program is being developed 
to integrate EPRI PSAPAC (Power Systems Analysis Package) 
programs for off-line planning.

There has been some very interesting work done by ABB Systems Control 
(EPRI Research Project 3103-2) on methods for determining stability
which use numerical and neural network techniques
rather than expert system rules. 
These rely on around a dozen ``critical indices'' (such as maximum voltage
deviation) which are obtained from the ETMSP program. 
They have identified correlations between these indices and experts'
stability classifications, and trained a neural network to take these
indices as input and output a stability classification.
This would be a useful complement to expert system rules.
They have also shown interesting results for using internal program 
variables to predict longer term stability problems, where the system
appears to be heading toward a stable state, but then suddenly goes 
unstable at some large time.
This type of instability can often not be detected from standard outputs
of ETMSP or PSS/E mid-term stability simulations.


\subsection{The Knowledge Base for the Decision Model}
\label{sec:KnowBase} 

The research team from the Syracuse University School of Information 
Studies (IST) has been designated to conduct research toward the  
development of a prototype expert system.  The expert system would 
assist planning engineers in the task of off-line transient stability 
analysis.  During the initial phase (Phase 1) of the project, we have 
concentrated our efforts on the acquisition of the expert knowledge to 
be incorporated into the analytical software of the prototype expert 
system.  In order to accomplish this objective, we have arranged 
interviews with planning engineers from Niagara Mohawk Power 
Corporation (NMPC), the New York Power Pool (NYPP), and the New 
York Power Authority (NYPA).  

In the conduct of these interviews we have applied methods for 
eliciting the knowledge required for the prototype expert system.  This 
knowledge is defined by the processes and parameters used by the 
engineers as they regularly approach the task of analyzing transient 
stability in their work.  We were specifically  interested in the 
determination of stability limits (interface transfer limits) through 
examination of the graphical outputs produced by the time domain 
module of the PTI PSS/E program.  This is the primary means used by 
NMPC planning engineers for off-line evaluation of transient stability 
for simulated cases.   

As an outcome of the interviews, we were able to develop a decision 
model based on the cognitive processes of the engineers.  This model 
was incorporated into design of the prototype expert system along with 
the associated decision parameters which formed the content of the rule 
base.  The preliminary results of the initial implementation of the 
prototype expert system would appear to indicate that we were 
successful in ``capturing'' the knowledge structure of the engineers in 
our decision model -- for a sample data set, the selection of transient 
stability limits were consistent for both the planning engineers and the 
expert system.             

As a result of preliminary conversations with NMPC engineers, it was 
decided that a system designed to assist planning engineers in the task 
of transient stability analysis could be a useful and realizable objective 
for implementing a prototype expert system.  The process involved with 
conducting transient stability analysis places considerable time 
demands on the NMPC planning engineers.  Using the PTI PSS/E 
software application, an off-line simulation is typically run for selected 
{\it base cases} in order to determine appropriate transfer limits for 
system interfaces.  The input parameters selected for initialization of 
the base case are based upon the existing interface limit plus 10\%, or 
200 MW -- in essence, a worst-case scenario.  The base case is tested 
using a set of contingencies (system failures) rank ordered from worst 
case (most detrimental effect) to best (least detrimental effect).  The 
contingencies are normally input by the planning engineers in this 
order (i.e. worst to best).  The outputs produced for each contingency 
are in a graphical form (for the time domain analysis module of the 
PSS/E application in the form of various types of post-fault oscillations 
for system units).  These outputs must be carefully examined and 
evaluated by the engineers for patterns which may predict system 
instability for specified interfaces.  The run-time involved with a typical 
analysis of transient stability is quite long and the paper output 
produced is considerable in volume.  The end result involves a 
significant expenditure in terms of time for the engineer, both in 
waiting for the computer runs and in physically examining the output.

The end result of the process has important ramifications for the 
Operations sector of NMPC.  The interface stability limits determined by 
the planning engineer's transient stability analysis are input into the 
EMS software used in Operations for real-time (or nearly real-time) 
security analysis for the entire NMPC power system.  System operators 
are constrained by these limits in that they are required to take timely 
corrective actions whenever contingencies create circumstances in 
which the limits are violated.  In fact, NYPP actively monitors the 
activities within the power grid and has the authority to order specific 
actions (e.g. shutting down power generation) under such 
circumstances, although such action is rarely implemented.  In sum, the 
transfer limits determined by the planning engineers may have a 
significant impact over the long run in terms of the cost of operations 
and overall efficiency of the power system as operators attempt to work 
within the constraints of the limits.  This would be especially true if the 
transfer limits were overly conservative, as they are likely to be given 
the concern for system security.          

The research team has operated on the hypothesis that an expert 
system may be able to significantly reduce the time involved in the 
current method of performing transient stability analysis (i.e. time 
domain analysis; a common off-line method in the industry), and, 
possibly, improve the accuracy of the interpretations.  Ideally, this may 
be achieved by 
\begin{enumerate}
\item reducing the computer processing time
\item automating consistent decision processes 
\item reducing the volume of output (and thus, the number of decisions)
      required for determining transient stability limits by providing
      only that output which is useful 
\item implementing a high-level user interface with the expert 
      system which would provide an improved representation of the data 
      output for the system user.      
\end{enumerate}

In order to study the feasibility of implementing an expert system 
for transient stability analysis, the research team established a Phase 1 
goal of developing a relatively straightforward prototype system which 
would be able to make a gross determination of stability based on the 
PTI PSS/E output.  If this goal could be successfully implemented -- 
indicating the potential for real benefits in terms of time savings and/or 
efficiency -- a more fully developed system could be considered.
 

\subsection{Knowledge Acquisition}
\label{sec:KnowAck} 

The acquisition of knowledge through interviews with power 
industry planning engineers (i.e. experts) was the primary objective of 
our research team in accomplishing the goal of  implementing the 
prototype expert system for Phase 1.  In order to devise a rule base for 
an expert system it is necessary to develop a knowledge structure 
which models the cognitive processes (decisions) of the human experts.
Therefore, our activities focused on eliciting the knowledge 
of domain experts (i.e. the planning engineers) as it applies to transient 
stability assessment so that it may be codified in the proposed expert 
system.    

Knowledge acquisition has traditionally been a problematic, but 
critically important aspect of developing an expert system.  Expert 
systems, in order to be successful,  must be able to duplicate the 
heuristic knowledge of human experts.  However, this knowledge will 
be of a qualitative as well as a quantitative nature; as in the case of the 
present application, including knowledge of unique operating conditions, 
weather, environmental variables, experience, etc.  Expert systems are 
well suited for representing such knowledge, but the challenge of fully 
eliciting it during the acquisition process has proven problematic.  We 
believe that the methodology employed in our research is ideally suited 
capture this type of heuristic knowledge.  
      
% Methods

In general, two broadly defined types of knowledge acquisition 
approaches were used when interviewing power industry planning 
engineers (experts): (1) an {\it unstructured technique} with either 
individuals or groups of experts; and, 
(2) a {\it structured technique} with selected individuals.

\bigskip
\noindent
{\bf (1) Unstructured Technique} 
\medskip

All of the interviews included an opportunity for the experts (in 
groups or individually) to lead the conversation and simply talk about 
the task (transient stability analysis).  This ``open'' form of interviewing 
helped to diminish barriers or tension between interviewer and the 
expert and aided in enhancing a freely flowing discussion.  This 
approach was useful in gaining a global perspective on the issues 
involved with transient stability analysis.  Moreover, the perspective 
gained originates with the expert and is expected to reflect the real 
concerns and factors involved with the task.

The IST team interviewed the following individuals during 
Phase 1 of the project:
\begin{description}
\item [\quad] {\bf NMPC} \\ 
 	Larry Eng -- Manager, Transmission Planning  \\
	(Larry Eng was interviewed on several occasions and was our
	primary liaison with NMPC)
\item [\quad] {\bf NYPP} \\ 
 	Ricardo Gonzales -- Associate Senior Engineer/Security Analysis \\
 	Jonathan Mayo -- Senior Engineer/Security Analysis \\
 	Robert Waldele -- Supervisor/Security Analysis
\item [\quad] {\bf NYPA} \\ 
 	Michael Henderson - Manager, Operations Planning \\
 	Larry Hochverg -- Associate Planning Engineer \\
 	Paul Roman -- Senior System Planning Engineer \\
 	George Dunn -- Associate System Planning Engineer
\end{description}

For interviews with groups (at NYPP and NYPA), an outcome of this 
technique was a ``brainstorming'' effect.  In each of these interviews the 
engineers appeared eager to engage one another in an attempt to define 
and summarize the procedures, obstacles, and overall ramifications of 
conducting transient stability analysis.

It may be generally stated that this unstructured interviewing 
technique was successful, and the interviewers were able to gain 
important insights into the overall, or global, picture.  We would attribute 
this success to the high level of interest in the topic exhibited by the 
engineers, and the earnestness with which they approached the task 
and cooperated with the efforts of the researchers.

\bigskip
\noindent
{\bf (2) Structured Technique} 
\medskip

Toward the project goal of implementing a prototype expert system 
for interpreting the PTI PSS/E graphical outputs, and in large part as a 
result of meetings with Larry Eng, the research team designed an 
instrument which was administered in one-on-one interviews with 
engineers from NMPC, NYPP, and NYPA.  The primary objective of this 
instrument was to capture and model the cognitive steps taken by 
planning engineers as they attempt to determine stable transfer limits 
for system interfaces.  These steps would provide a knowledge 
structure to be incorporated into the prototype expert system rule base.

The instrument was based upon a set of time domain graphical 
outputs produced by the PTI PSS/E software program (with inputs 
selected by Larry Eng).  The time domain simulation produces hard 
copy output displaying post-fault oscillations for certain variables 
associated with transient stability analysis (i.e. machine angle, voltages, 
accelerating powers, etc.) over a specified time duration.  This type of 
off-line program output is regularly used by planning engineers at 
NMPC and other NYPP utilities (as well as NYPP itself) for the purpose of 
determining stable transfer limits for specified interfaces within the 
overall power grid.  Engineers will typically examine these graphical 
outputs looking for patterns which may indicate instability.  The most 
salient indication of instability would be an oscillation which did not 
show adequate damping over the specified time duration.  This would 
indicate that the unit was not readily returning to its normal state 
following the simulated fault condition (the broadly conceived 
definition of system stability).       

Our research team was interested in seeing whether there were 
regular patterns in the way planning engineers examined these outputs 
and interpreted the results.  Therefore, it was decided that a set of 
outputs would be produced which would simulate a condition for a 
specific interface in which the transfer limits would create a full range 
of variance (and therefore a full range of decision points) from 
definitely stable to definitely unstable (as defined by the interpretation 
of the engineers).  These outputs would be shown to various engineers 
within the NYPP system to see how each engineer examined the outputs 
and what conclusions they derived.  The results would then be 
compared to determine whether, in fact, consistent patterns did emerge.

Five sets of output were produced based on the Central-East 
interface.  The transfer limits used were selected according to the 
experience and expertise of Larry Eng: 2810 MW, 2910 MW, 3075 MW, 
3200 MW, and 3400 MW.  The time duration for the simulation was 15 seconds.

A total of 5 engineers were interviewed using the instrument.  They 
were selected from the 8 engineers who were interviewed as part of the 
unstructured group and individual interviews (see above).  They were: 
Larry Eng (NMPC), Ricardo Gonzales (NYPP), Robert Waldele 
(NYPP), Larry Hochverg (NYPA), and Paul Roman (NYPA).  In each of 
these interviews the engineer was presented with the set of 5 outputs 
and instructed to select the transfer limit where the condition appears 
to be clearly stable.  The engineers were asked to ``think aloud'' as they 
examined the outputs so that the interviewer could record {\it decision 
points} as they occurred.  The decision points were coded according to a 
scheme designed to elicit both the parameters upon which the decision 
was made (i.e. the content) as well as the temporal sequence of the 
decision process.  For example, a code of 1.2.3 on a voltage angle 
oscillation as part of the 3200 MW output would indicate that the 
transfer limit (i.e. 3200 MW) was the first one examined by the 
engineer, the parameter of the output (i.e. voltage angle) was the 
second one examined for the 3200 MW transfer limit, and the oscillation 
marked was the third one marked for that parameter (i.e. the third 
voltage angle marked for the 3200 MW transfer limit).  In addition, the 
engineers were instructed to explain {\it why} they were looking at each 
portion of the output (i.e. what were the decision criteria?).  This 
information was also recorded by the interviewer.

There appeared to be a general consensus with regard to selecting a 
stable transfer limit for the Central-East interface.  Of the 5 engineers 
interviewed, 3 chose 3200 MW, 1 chose 2910 MW, and 1 chose 2810 
MW.  It should be noted that the 2 engineers who did not select 3200 
MW implied that they might have raised the limit if they were able to 
manipulate the input variables.  This was especially true for Paul 
Roman of NYPA who labeled 3075 MW and 3200 MW as ``borderline''.

The actual results were as follows - transfer levels shown are those 
actually examined by the engineers:
\begin{description}
\item [\quad] Larry Eng -- NMPC  \\
	3400 MW:  clearly unstable  \\
	3200 MW:  generally stable  

\item [\quad] Ricardo Gonzales -- NYPP  \\
	3400 MW:  clearly unstable  \\
	3200 MW:  generally stable  

\item [\quad] Robert Waldele -- NYPP  \\
	3400 MW:  clearly unstable  \\
	3075 MW:  appears unstable  \\
	2910 MW:  marginal          \\
	2810 MW:  stable

\item [\quad] Paul Roman -- NYPA \\
	3400 MW:  clearly unstable  \\
	3200 MW:  borderline        \\
	3075 MW:  borderline        \\
	2910 MW:  clearly stable    \\
	2810 MW:  clearly stable 

\item [\quad] Larry Hochverg -- NYPA \\
	3400 MW:  clearly unstable  \\
	3200 MW:  stable 
\end{description}
  

A detailed analysis of the interviews (structured and unstructured) 
revealed generally consistent patterns in the decision processes of the 
engineers.  As a general rule the engineers tended to begin with the 
output of the transfer limit above the ``expected'' limit (the expected 
limit might typically be the actual or existing transfer limit for the 
case).  They would continue looking at incrementally lower limits until 
the outputs revealed characteristics indicating irrefutable stability 
(``bracketing'' the level finally selected).  At each level, the engineers 
invariably started with the variable {\it rotor angles}. Instability indicated 
for this variable discounted the associated stability level as a stable 
transfer limit for the case.  In essence, the engineers seemed to apply a 
rule of thumb where, if rotor angle outputs did not indicate a stable 
condition, the level could be categorically classified as unstable and no 
further investigation was necessary.  At a level where stability was 
indicated by the rotor angle graphs, {\it accelerating powers} and 
{\it voltages} 
were also examined along with a more detailed examination of rotor 
angle criteria (such as units especially relevant to the particular 
interface; e.g. geographically) in order to provide confirmatory evidence 
of stability.

Based on the results of the interviews, decision models were derived 
for both the engineer's interpretation of the PTI PSS/E time domain 
output (see Fig.~\ref{fig:exp1})
and the overall process of conducting transient stability analysis 
(see Fig.~\ref{fig:exp2}).  
The former model was incorporated into 
the design of the prototype expert system.



\subsection{The Expert System Implementation of the Decision Model}
\label{sec:ESImp} 

Here we give an outline of the expert system implementation of the 
decision model for Transient Stability Analysis. 
The Cosmic CLIPS expert system shell~\cite{CLIPS} was used for this purpose.
A survey on available expert systems was performed, and the decision was 
made to use the CLIPS system, for a number of reasons including popularity,
ease of use, and easy availability of both documentation and code.

Based on the interviews held with experts a decision model was deduced
by the IST team. The outline of the model is as follows:
%
\begin{enumerate}       
\item Start decision process from the maximum power output (MW) available
\item Check stability of ROTOR ANGLE for the Current MW under test --
      if not stable, repeat (2) for a lower power output 
\item Check stability of VOLTAGES for the Current MW under test -- 
      if not stable, repeat (2) for a lower power output 
\item Check stability of ACCELERATING POWER for the Current MW under
      test -- if not stable, repeat (2) for a lower power output 
\item Fix the HIGHEST SAFE LOAD LIMIT as the power output that is stable
      for all 3 parameters, minus a $10\%$ safety margin
\end{enumerate}       
%
Before discussing the implementation details it is necessary to understand
how the input data to the expert system was modeled, to serve as test cases. 

\bigskip
\noindent
{\bf Test Cases:} 
\medskip

The waveforms that were modeled for serving as test cases are from 
the Central-East Interface. The specific waveforms considered were from 
the following sections of the PSS/E output (in order):
\newcounter{bean}
\begin{list}{(\roman{bean})}{\usecounter{bean}}
\item INTRFC SYS ROTOR ANGLE
\item INTRC SYS VOLTAGES
\item NE ACCELERATING POWERS
\end{list}
The implementation can be extended to analyze more input parameters. 

The modeling of the waveforms associated with the above parameters
were done as follows:
\begin{itemize}
\item  An aggregate behavior of the waveform for 0 -- 7.5 secs (s1) and
       7.5 -- 15 secs (s2) were determined visually from the graphs.
       (NOTE: This will be replaced by actual data from the
       transient stability analysis program when this is available)
\item  The peak-to-peak values at  7.5 (V1) and 15 secs (V2) were also
       determined from the graph.
\item  The behaviors s1 and s2 can be: 
       \begin{itemize}
       \item positive (+) damping 
       \item negative (-) damping 
       \item zero (0) damping  
       \end{itemize}
\end{itemize}


A typical waveform as an Expert system data (fact) will look like this:
\begin{verbatim}
        (MW  3200 parameter 1 waveform 2 positive 1 negative 2)
              |             |          |          |          |
             /              |          |          |           \
            /               |          |          |            \
           /               /           |           \            \
          /               /            |            \            \
  The MW output       Represents      Which    s1=positive    s2=negative
of the associated  the precendence  waveform        &              &     
    waveform       order of the      in the    V1=1 unit       V2=2 units
                     parameter      parameter
\end{verbatim}

The input data (facts) were analyzed in the order specified
by the parameter precedence field, which is made to coincide with the 
requirements of the decision model. 
For example in this specific implementation the parameter precedence number for
\vspace*{2mm}

	INTRFC SYS ROTOR ANGLE  will be 1

	INTRC SYS VOLTAGES  will be 2, and

	NE ACCELERATING POWERS  will be 3.

\vspace*{2mm}
\noindent
which is what the decision model calls for.
	
The decision on stability is determined as follows (this as an extension of
the decision model, and was generally agreed to be a suitable decision
process for demonstration purposes):
%
\begin{centering}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{ccc}
\vspace*{5mm} \\
\hline %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\qquad\qquad  s1 \qquad & \qquad     s2  \qquad 
& \qquad Decision (`S'table/`U'nstable) \qquad\qquad \\
\hline %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\qquad\qquad  -  \qquad & \qquad  -  \qquad & \qquad          U        
\qquad\qquad \\
\qquad\qquad  -  \qquad & \qquad  +  \qquad & \qquad  S if V2 $<$ THRESHOLD else U \qquad\qquad \\
\qquad\qquad  -  \qquad & \qquad  0  \qquad & \qquad          U        
\qquad\qquad \\
\qquad\qquad  +  \qquad & \qquad  -  \qquad & \qquad 		U       
\qquad\qquad \\
\qquad\qquad  +  \qquad & \qquad  +  \qquad & \qquad 		S        
\qquad\qquad \\
\qquad\qquad  +  \qquad & \qquad  0  \qquad & \qquad  S if V2 $<$ THRESHOLD else U \qquad\qquad \\
\qquad\qquad  0  \qquad & \qquad  -  \qquad & \qquad 		U      
\qquad\qquad \\
\qquad\qquad  0  \qquad & \qquad  +  \qquad & \qquad 		S     
\qquad\qquad \\
\qquad\qquad  0  \qquad & \qquad  0  \qquad & \qquad  S if V2 $<$ THRESHOLD else U \qquad\qquad \\
\hline %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace*{5mm} \\
\end{tabular}
\end{centering}
%
The value of the THRESHOLD can be varied for analysis.       
Reducing THRESHOLD will yield a more conservative decision. 
In cases where the decision has to be made on comparing V2 and THRESHOLD,
we could instead analyze increase the simulation time to obtain a further
15 seconds of data to see if a decision can be made independent of the
THRESHOLD.
If all the waveforms in all the parameters are S then the corresponding
MW value will be considered as stable. From that the Highest Safe Load
Limit is fixed as MW -- 0.1MW ($10\%$ below MW).   

The decision model implemented chose 3200 MW as the stable output 
(for THRESHOLD = 1 unit) and
fixed the Highest Safe Load Limit as 2880 MW. This conforms with
decision taken by the majority of the experts (who were not ``conservative'').
The mean time taken for the decision process was 0.28 seconds 
(based on 10 sample runs that were made on the same set of facts).



\subsection{Future Work}
\label{sec:ESNext} 

Analysis of the output data for 
stability assessment is only one part of the process of setting stability 
limits.  The engineer must also choose appropriate contingencies and 
parameter values in order to find the limiting value for a certain 
component, above which the network becomes unstable.  While there is 
likely to be considerable individual variability from engineer to 
engineer in these aspects, we would suggest that the process is not 
likely to be chaotic and consistent decision patterns may be elicited.  
This problem is a logical next step for further investigation.

Although our efforts so far were primarily directed toward the problem of 
analyzing transient stability off-line, we recognize an industry demand 
for expert system technologies on the operations side as well as the 
planning side.  In fact, it seems safe to say that the volume and 
complexity of data and decisions (decisions with economic 
consequences) faced by operators on-line in today's EMS systems makes 
this area the most needful of AI technology.  We would further 
recommend that in the face of information overload, issues such as the
human-computer interface will be critically important for the successful 
implementation of expert system technologies (e.g. in Dynamic Security 
Assessment).  We suggest that this is an important area for further 
study, and we plan to investigate this problem starting in Phase 2
of this project. 

Currently the input to the expert system has been given as a simple
user determination of the rules as applied to graphical output 
(basically the user mimicking the experts). 
Clearly for an on-line system we want to remove the user from this loop, 
and process the numerical output from the transient stability program
to obtain this information. So far we have been unable to do this, since
we do not have access to the source code for the PTI PSS/E program.
However we do have the EPRI ETMSP program, and we will be utilizing
this program in the future. We can then do some post-processing on the
output from this program that will allow us to create a genuine
automated system for interpreting the transient stability data.



\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{On-Line Transient Stability Analysis}
\label{sec:OLTS} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

EPRI has an ambitious plan to develop an on-line
dynamic security assessment (DSA) program which will basically 
incorporate transient stability analysis into the contingency evaluation
function of the security model of an energy management system.
This is a very large project which is not expected to be completed until
1997. 

Clearly we want to avoid duplication of effort with this large EPRI project.
Our project has a similar, but more modest, goal -- to provide an on-line
transient stability program which will give dynamic stability limits for
operators to use within a standard static contingency evaluation program.
However the overlap of the two projects is actually quite small, since 
their main focus is on different areas:
%
\begin{itemize}
\item
Although EPRI have done some research on the use of parallel processing
\cite{Bose93},
they are planning to implement their DSA program on a workstation platform,
and are not actively pursuing the development of a parallel transient 
stability program, which is the main thrust of our project.
\item
Most of the EPRI work is concentrated on improving contingency selection,
screening and ranking techniques~\cite{DSA,Screen}, 
which is not part of the scope of our project.
\item
EPRI are working on improving the performance of the ETMSP program, but
this is being done by optimizing the current sequential program and using
the latest generation of high-performance workstations. \\
We are also investigating improvements in the sequential performance of
the transient stability program, but we are focussing on improved
algorithms such as DAE solvers.
\item
EPRI are attempting to classify system stability using analytic and 
neural-net based methods, whereas we are using an expert system approach.
\item
The main similarity between the two projects is in methods to provide
on-line stability limits. However this is viewed more as a side-issue 
for the DSA project, whereas it is the central issue for our project.
\end{itemize}
%
The NMPC and EPRI projects are thus mainly complementary
rather than duplicative.


\subsection{The EPRI On-Line Dynamic Security Assessment Project}
\label{sec:DSA} 

NPAC representatives attended the EPRI Power Systems Analysis Package
(PSAPAC) User's Group meeting in Dallas on December 1-3, 
and the EPRI meeting on On-Line Transient Stability Assessment in
Toronto on December 6-7, to obtain information on EPRI projects and
programs related to our work.

OVERVIEW OF PROPOSED EPRI DSA PROGRAM 

A diagrammatic description of the proposed EPRI DSA program is given in 
Fig.~\ref{fig:DSA}.

provide the large decrease in computation time required 
for an on-line transient stability program, 

EPRI are hoping that these improvements will be enough to allow on-line 
capability.

Ontario Hydro, in collaboration with EPRI, have developed a
prototype on-line transient stability program (OLTSP) that
is currently being studied in the control room at Ontario Hydro.
The developers were quick to point out that this was not an easy 
undertaking -- Ontario Hydro's prototype OLTSP took many years and many
millions of dollars to develop.


\subsection{Parallel Processing for Dynamic Security Assessment}
\label{sec:ParallelDSA} 

A simple way of exploiting parallelism in an on-line DSA program is to 
parallelize over contingencies, as described in section~\ref{sec:ParContEval}. 
The EPRI DSA program aims to do a contingency screening and ranking
of around 1000 contingencies, choosing approximately 30 of the worst
contingencies for a full transient stability analysis using ETMSP.
In both the screening and
analysis phases the effects of each contingency can be calculated 
independently. This independent parallelism would give almost perfect 
speed-up on a parallel computer or workstation cluster, 
and is the best approach for this type of problem where a large number
of independent contingencies have to be studied.

The prototype EPRI program will allow this independent parallelism if the
program is run on a workstation network.
The program enables the user to specify computers in the network
which can be used, and different contingencies are sent to these machines
when they are idle. This is an example of using a workstation cluster
as a parallel computer (see section~\ref{sec:ParComp}).

This method cannot be used for other sections of the DSA program, however.
For example, the stability limit for an interface is calculated by taking
the worst contingency for that interface, and using an iterative procedure:
\begin{enumerate}
\item Increase the flow across the interface
\item Do a full transient stability calculation 
\item Classify the result as stable or unstable
\item Calculate sensitivities and use these to determine the next flow value
\item Repeat until system becomes unstable
\end{enumerate}
Speeding up this part of the DSA program would require parallelizing the 
transient stability analysis program. 



\subsection{Feasibility of an On-Line Transient Stability Program}
\label{sec:OnlineTSA} 

The EPRI DSA program aims to supercede contingency evaluation, and 
therefore needs to be run about as frequently. The goal is for the 
program to cycle every 15 minutes.
A program providing on-line stability limits would ideally generate
them anew every time a contingency evaluation was run, however this
is not necessary. Less frequently revised limits, for example every hour,
should be adequate. Some level of conservatism may have to be 
introduced to allow for changes in the state of the system during this
time interval, however these dynamic limits would still be much less 
conservative than those currently used, which are based on worst-case
scenarios.

Our idea for an on-line calculation of stability limits is to use an
expert system to basically mimic the process an engineer would follow
to set these limits.

EXPLAIN WHAT IS CURRENTLY DONE.

GIVE ROUGH IDEA OF HOW MUCH SPEED-UP IS REQUIRED TO HAVE A USEFUL
ON-LINE PROGRAM.

Parallelism could be used in the same way as described in the previous
section. It should be possible to study the required contingencies 
independently in parallel, that is, just studying a different contingency
on each processor, using the sequential transient stability program.
However a sequential, iterative procedure is necessary to set the stability
limit for a particular interface, for which speed-up can only be achieved
by using a parallel version of the transient stability program.
This would also be used if the program was run in study mode by an
operator to find a stability limit for a particular case.

We believe it is feasible to develop an on-line system which will 
provide transient stability limits to control room operators. 
Parallel algorithms to speed up transient stability analysis, 
and an expert system knowledge base to classify the results, 
have already been identified, and development of prototype programs
to implement these functions has begun. 
There are three main obstacles to be overcome to develop a program which
will provide on-line stability limits:
%
\begin{enumerate}
\item {\bf Development of a fully parallel version of ETMSP.}\\
      The main difficulty here is that ETMSP is a large and complex
      program, so fully implementing the parallel algorithm could
      take some time. 
\item {\bf Development of an expert system for finding stability limits.}\\
      An expert system to classify the ETMSP result as stable or unstable
      is only part of what is required to solve this problem. 
      We also need some expert system rules which describe how an engineer
      would use the transient stability program to find the stability limit
      for a particular component or interface in the system. This is a 
      much more difficult task.
\item {\bf Integration of the on-line stability program into the control
      room EMS.}\\
      This is quite a large software engineering problem, which would 
      probably require input from STAG, the developers of the EMS software. 
      The on-line stability program would need to access the network data 
      from the EMS
      state estimator, and return stability limits in a form which is 
      usable by the EMS security model programs. It should also interface
      to study mode operation of the EMS programs.
\end{enumerate}
%
Assuming no major unforseen problems (which are always possible in any
research project of this kind), we would expect that items 1 and 2
should be completed by the end of this year. Another few months would
be required for intensive testing and validation of these programs.
Estimating the time required to integrate our program into the control
room EMS is a more difficult problem, since we currently do not have
enough information on the STAG EMS software. Ideally we would interface
with the new STAG EMS for an IBM workstation network, however this software
is not even completed yet. As a very rough estimate, we would expect the
integration to take about 6 months, with another few months for validation
and testing. We would thus estimate that an on-line program to provide
transient stability limits could be available by the summer of 1995.
Further testing and improvements would probably be required after that date.
This is in accord with our original proposal for the time required for
Phase 3 (the development phase) of this project, which we originally 
estimated to be roughly 18 months after the completion of Phase 2 (which
ends in March 1994), which would be September 1995.


Network stability is not limited to transient stability. Other problems
such as voltage stability and small signals stability need to be considered
also. This is not being attempted in the initial implementation of EPRI's 
DSA program.
Like transient stability analysis, these studies are currently only
done off-line by planning engineers. Increased computing power could
also allow these studies to be integrated into an all-encompassing
on-line security program. 
The possibility of applying parallel processing to these other stability
programs is discussed in section~\ref{sec:SecurProg}.

\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Parallelism and Computational Issues}
\label{sec:Parallel} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Many industries, including the utility industry, understand the 
possible benefits of high performance computing, 
and are interested in using parallel computing in the future. 
However they are often adopting a ``wait-and-see'' approach, in the 
belief that parallel software and hardware are not yet stable or 
standardized enough to warrant incorporation into production use 
programs.  

This was true in the 1980's, but the situation has changed in the 90's.
Parallel computer architectures and computing models are now 
well defined, and the larger computer companies such as IBM, DEC,
Hewlett-Packard, Intel, Sun and SGI are all becoming major vendors of
parallel computers.
Standardized parallel languages have been defined which are portable
between different computers with the same hardware architecture, 
and in some cases can be used on any kind of computer,
from workstations to shared memory vector supercomputers to 
massively parallel computers. 
This standardization and portability means that parallel software
can be developed which will be usable on any current and 
future high performance computer.


\subsection{Parallel Computers}
\label{sec:ParComp} 

shared vs distributed memory ??
SIMD vs MIMD or SPMD ??

Describe CM-5

One of the criticisms of parallel computers is that they are viewed as
having ``exotic'', rapidly changing, and often unreliable architectures.
However large mainstream computer companies such as IBM, DEC and 
Hewlett-Packard have all recently introduced parallel computers which
are basically a collection of standard high-performance workstations
connected by either a standard network interface (Ethernet or FDDI), or a 
proprietary high-speed communications network incorporating a fast switch
to route messages between different processors. 
This new wave of parallel computers are generically referred to as
{\it workstation clusters}.
IBM is offering its Scalable POWERparallel system (SP-1) based on the
IBM RS/6000, 
DEC is offering the Gigaswitch to connect its Alpha-based workstations
into a parallel computer, 
and Hewlett-Packard have a joint venture with Convex to provide a
cluster of HP Apollo workstations.

Workstation clusters are becoming very popular as cost-effective
computational engines, which can be used as either a network of 
independent workstations for running sequential programs, or a connected
parallel computer for running parallel applications.
The systems software supports both methods of use.
We believe that these types of workstation clusters and workstation-based
parallel computers are ideally suited to power systems applications, 
and could support both on-line systems and as a company-wide shared 
computing resource.

Price and price-performance compared with mainframes??

NPAC has recently installed an 8 processor IBM SP-1, which will soon 
be upgraded to 16 processors and have the fast switch and communications
network installed. We also have access to the 64 processor SP-1 at
Cornell University, which is scheduled to be upgraded to a 512 processor
machine.
NPAC will also be installing an 8 processor DEC Alpha clusters with a
Gigaswitch connection in the near future.
We are planning to test the performance and usability of these machines,
especially for power systems applications.
Our parallel sparse matrix solvers and the prototype parallel ETMSP
and unit commitment programs will be benchmarked on these machines,
as well as on the 32 processor Thinking Machines CM-5 at NPAC.


Sun and Silicon Graphics (SGI) are also offering parallel shared memory 
parallel computers.
Shared memory machines do not scale well to large numbers of processors, 
so these are limited to 4, 8 or 16 processors, % ???
although shared memory computers with larger numbers of processors are
available.
Methods exist for implementing distributed memory message-passing 
programs on shared memory machines, so as long as programs are written 
in standard, portable parallel languages (see the next section for more
details), they should be able to run on any parallel computer, including
shared memory machines.
The recent EPRI report on parallel transient stability analysis~\cite{Bose93} 
claims that this application is better suited to shared memory machines,
however we would dispute that claim.



\subsection{Parallel Languages}
\label{sec:ParLang} 

Most power systems programs, including the EPRI ETMSP transient
stability program, are written in Fortran 77. We believe many of these
programs could be parallelized using High Performance Fortran 
(HPF)~\cite{HPF}. 
This is a new industry standard data parallel language for high performance
computers of any architecture (vector, parallel or sequential), 
based on the Fortran 90 standard, which includes Fortran 77 as a subset.
The HPF language specification
was recently formalized by the High Performance Fortran Forum, 
a collection of universities, government laboratories, and commercial 
software and hardware vendors such as DEC and IBM.
Commercial HPF compilers for a wide variety of computers are expected to 
become available within the next year.

The EPRI and PTI transient stability programs are written in sequential 
Fortran 77. In porting this to HPF, one would restructure the code
so that the loop over generators was clean (i.e. no data dependencies
between iterations), and then provide some compiler directives indicating
how the generator data should be distributed over different processors, 
and indicating that each iteration of the loop could be done in parallel.
The network equations would be solved by a call to a parallel sparse
matrix library.

Programs for distributed memory parallel computers can also be written
using portable message passing languages based on C or Fortran.
Many such languages are currently available, including 
Express~\cite{Express} and Zipcode~\cite{Zipcode}, which we have
used in our parallel programs.
An industry standard message passing language, known as the 
Message Passing Interface (MPI), is currently being developed, and should
be finalized this year. 
Programs in both Express and Zipcode will be able to run on any machine
supporting MPI, however for ease of maintenance we are planning 
to port our message passing programs to MPI format
when the standard is complete and compilers become available.


\subsection{Parallel Algorithms and Software}
\label{sec:ParAlg} 

High performance computers use advanced hardware architectures, 
such as vector and parallel processors, to obtain improved performance
over standard sequential computer architectures.
Effective utilization of these special architectures requires special
techniques and software.
Although vectorizing and parallelizing compilers do exist for many high 
performance computers, it is very rare for a standard sequential program to 
obtain anywhere near the peak performance of these machines, and in many
cases the speed-ups obtained over sequential workstations are negligible. 

Obtaining high performance from parallel computers almost always requires
different methods (or algorithms) to those used on sequential machines.
A great deal of research has been done over the past 15 years on developing 
efficient parallel implementations of common algorithms.
For some of these algorithms parallel programs already exist, for some
these could be fairly easily implemented, however for many problems 
implementing efficient parallel programs can be a very difficult task.
Software libraries containing efficient parallel algorithms are now
available for many parallel computers, as are portable public domain 
software packages.

Power systems simulations generally involve sparse linear algebra
algorithms, since power networks are very sparsely connected.
Although dense matrix algorithms can be parallelized very efficiently, 
sparse matrix algorithms are one of the most difficult commonly-used
techniques to parallelize, especially when the matrices are very sparse,
as with power systems problems. Parallel sparse matrix algorithms are
discussed in detail in section~\ref{sec:Sparse}. 

SCALABILITY ??

SCALAPACK~\cite{SCALAPACK} is a scalable parallel version of 
LAPACK (Linear Algebra Package)~\cite{LAPACK},
which is the most comprehensive and widely used public domain linear algebra
software. 
LAPACK is the successor to LINPACK, which is used in the ETMSP program.
Much of the development of SCALAPACK is being done by members of the Center 
for Research in Parallel Computing (CRPC), a group of universities and
national government laboratories of which NPAC is a member. 
We therefore have access to the developers of the parallel code,
if any problem-specific tuning of the code is required to improve
efficiency.
% ???

SCALAPACK and other parallel software libraries contain general sparse
matrix and sparse eigenvalue solvers, which are the main algorithms used
in many power system simulations.
As noted in section~\ref{sec:Sparse}, general-purpose routines for sparse
matrix solvers are not the most efficient techniques for power systems 
applications, so we have developed our own specialized sparse matrix solver. 
General eigenvalue solvers should be good enough to be used for power
systems problems, although specialized algorithms may also perform better.



\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Survey of High Performance Computing and Parallel Processing 
         in the Utility Industry}
\label{sec:Survey} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Previous surveys of utility applications which could benefit from 
high-performance computing, such as the TASC report~\cite{TASC}, 
have concentrated solely on identifying applications which are very
compute-intensive, and therefore quite slow on current computers. 
In many cases, speeding-up these programs would be of substantial benefit 
to utilities.
However the TASC report does not emphasize the fact that the hardware used 
in high performance computers, such as vector and parallel processors, 
generally requires special parallel algorithms and software to be utilized 
effectively, as discussed in section~\ref{sec:ParAlg}. 

We have been concentrating on trying to identify the core algorithms used
in different programs of interest to NMPC. This is not always an easy task,
since documentation is not always readily available, or else may not
provide sufficient detail on the methods used in the program. 
We generally need special permissions, or even to sign non-disclosure
agreements, just to look at the user's manual for certain programs, which
may not even have enough detail to determine the basic algorithms used.
We have studied or requested the appropriate documentation, or attempted
to talk to the code developers, for a number of different power systems
programs developed by EPRI, PTI, STAG, and GE. 

In Phase 1 we have concentrated our study mainly on programs concerned
with security and economic models, with only preliminary studies done
in other areas of utility operations.
We will broaden the survey in Phase 2 of this project and report our 
findings for these other areas at the end of that phase.

One general method to note which could easily be applied to many
NMPC computing applications is the idea of {\it job-level parallelism}.
If a large corporation is using a centralized mainframe computer as a
shared resource, than many different programs or jobs will generally be 
running on this machine at any one time,
using a {\it time-sharing} approach whereby the computer switches 
rapidly from one program to another. If there are 10 programs running
at the same time, then each one will take 10 times longer to complete
than if it had the computer to itself, giving a slower turnaround time.

A cluster of workstations, or a parallel machine such as the IBM SP-1,
contains many processors, each one faster than an IBM 3090 mainframe. 
If many programs are being run at one, a different program can be
assigned to a different processor, 
thus greatly reducing the turnaround time from submitting the job to
obtaining the results.
Any program which will run on an IBM workstation will run on a
workstation cluster or an IBM SP-1 in this way, using available systems
software which finds the least loaded processor and runs the program on
that processor. No parallel programming is required to obtain this
``effective speed-up'', or reduced turnaround time, from this job-level
parallelism.



\subsection{Security Programs}
\label{sec:SecurProg} 

OVERVIEW ??

%
\begin{center}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{lcr} 
\hline %-----------------------------------------------------------------
Application & \qquad\quad Algorithms used \qquad\qquad & Parallel algorithms \\
\hline %-----------------------------------------------------------------
State Estimation         & sparse matrix solver & many available \\
Load Flow                & sparse matrix solver & many available \\
Contingency Evaluation   &    load flow         & embarrassingly parallel \\
Optimal Power Flow       &      ???             & ??? \\
Voltage Stability        & sparse eigenvalue    & some available \\
Small Signals Stability  & sparse eigenvalue    & some available \\
Transient Stability      & sparse matrix solver & many available \\
			 & DAE or ODE solver    & exist but need modification \\
Long-term Stability      & sparse matrix solver & many available \\
			 & DAE or ODE solver    & exist but need modification \\
\hline %-----------------------------------------------------------------
\end{tabular}
\end{center}

As can be seen, most of the programs which involve network simulation
are based on sparse matrix algorithms or differential equation solvers.
A number of portable sparse matrix solvers and sparse eigenvalue solvers
have been developed, and the best general algorithms are incorporated
into the SCALAPACK parallel linear algebra software library described in
section~\ref{sec:ParAlg}. An overview of parallel sparse matrix algorithms
is also presented in that section. 
The parallel sparse matrix techniques developed for transient stability 
programs would be equally useful for a parallel load flow program, and
possibly a parallel program for state estimation.

Parallel differential equation solvers for transient stability are discussed
in detail in section~\ref{sec:DAE}. These techniques could also provide
faster and more accurate simulations for long-term stability studies.



\subsection{Economic Programs}
\label{sec:EconProg} 

OVERVIEW ??

%
\begin{center}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{lcr} 
\hline %-----------------------------------------------------------------
Application & \qquad\quad Algorithms used \qquad\qquad & Parallel algorithms \\
\hline %-----------------------------------------------------------------
Unit Commitment        & dynamic programming  & have been developed at NPAC \\
Economic Dispatch      & Lagrange multipliers & could be developed \\
MAPS                   & linear programming   & could be developed \\
		       &      ???             &        ???         \\
\hline %-----------------------------------------------------------------
\end{tabular}
\end{center}

Unit commitment is discussed in section~\ref{sec:UnitComm}, and 
MAPS is discussed in the next section.


\subsection{ENGINEOUS and MAPS}
\label{sec:ENGINEOUS} 

General Electric have developed a program called ENGINEOUS~\cite{ENGINEOUS},
that is used to help engineers develop products for which the best 
design involves optimizing requirements from many different areas.
ENGINEOUS provides an expert system for such multi-disciplinary analysis 
and design engineering, by using software integration techniques
to interface the programs for the different disciplines, and a combination
of optimization tools and an expert system knowledge base to find an
optimal design.
In December we met with Nelson Simons and Rana Mukerji from the 
General Electric Power Systems Engineering Department, concerning
possible common interests between our project and their project on the 
application of ENGINEOUS to power system planning~\cite{EPSP}.

GE are using ENGINEOUS to combine a pricing program, GE's 
Multi-Area Production System (MAPS), with a stability program, PTI's PSS/E.
Their main concern is that both the application programs (MAPS and PSS/E)
are very slow, so running a single scenario can take 5 hours on an
IBM workstation. % ???
This makes it very difficult for the program to optimize over different
scenarios. To obtain good results it is clearly necessary to speed up both
these programs, and using parallel processing seems the obvious solution. 
Although PTI's PSS/E program is currently used for transient stability 
analysis, it would be a straightforward task to integrate the EPRI ETMSP
program instead.

One of the main algorithms used in the MAPS program is linear programming.
Probably the world's leading experts on linear programming optimization
techniques are the group at Rice University, who are part of the Center
for Research in Parallel Computing (CRPC), in which NPAC has a major role.
We would therefore be able to utilize their expertise to investigate 
the possibility of implementing better and faster linear programming 
algorithms for MAPS. In this field using improves sequential algorithms
can often give large speed-ups.

There are also some common areas of interest with the expert system
and user interface aspects of our project.
Their program also needs a method of determining system stability and
stability limits from the output of the transient stability program.
Currently they are using
rather crude methods to do this. Our work on developing a more sophisticated
expert system for this task could be effectively utilized in their program.
They are also grappling with the problem of how best to interface the program
with the planning engineer, for which the Syracuse IST group has substantial
expertise.
Another problem is that the engineer may want to add or change rules in the
knowledge base. Currently this requires knowledge of the Lisp programming
language used in ENGINEOUS. They are interested in utilizing a simpler
method of entering rules using simple English constructions,
as is done by the CLIPS expert system package being used in our project.

We believe a collaborative effort between GE, NMPC, NPAC and possibly EPRI
to parallelize the EPRI ETMSP and GE MAPS programs and improve the 
application of ENGINEOUS to power systems planning would be extremely 
beneficial to all parties. We are hoping to obtain extra funding from
ESEERCO, the New York State Energy Research and Development Authority, 
and possibly EPRI, to expand our project to include work on speeding-up
the MAPS program and the use of ENGINEOUS for power systems planning.



\subsection{An Advanced Integrated On-Line EMS}
\label{sec:EMS} 

The use of parallel computing could not only enable the use of on-line 
transient stability analysis, but also allow for the incorporation of
other stability programs, and possibly pricing optimization 
programs such as unit commitment and economic dispatch, into a single
integrated on-line EMS. The different programs which make up this 
integrated system could be run on different processors or groups of
processors, to allow the system to work in near-real-time.

For example, one processor could be collecting data from the telemetry
units and updating the database; another could be running the state estimator; 
16 processors could be running a parallel transient stability program to set 
stability limits; while another 16 could be running voltage stability and/or
small signals stability checks; 32 processors could be running a 
contingency evaluation; another could be running a pricing optimization
program incorporating unit commitment and economic dispatch; and still 
others could be running a program to try to find a scenario which would 
minimize the cost of handling a particular contingency should it occur.
Other programs to optimize the security and economy of power generation 
and transmission could also be included, as well as processors to handle
requests from operators to run simulations in study-mode.
Results from all these simulations would be passed to a controlling
workstation running a program to manage the overall EMS. An expert system
would probably be necessary to manage this daunting mass of data, 
and to present a simple user interface to the operator.

A single serial computer would not have anywhere near the computing 
capacity to run all these different applications at once with the speed
necessary for an on-line system.
Currently every element of the Operations EMS, from adding telemetry
data to the database, to running a contingency evaluation, to updating
the graphical user interface displaying results to the operator, is
run on the single processor of the mainframe. 
However there is a large amount of job-level parallelism 
in an advanced, integrated on-line EMS, that could be exploited by
simply assigning the many different tasks to different processors of
a parallel computer.



\newpage
\bibliography{power}
% \bibliographystyle{plain}
\bibliographystyle{ieee}


\end{document}


