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

\reversemarginpar

\input{psfig}


\begin{document}

\section{Transient Stability Problem}

The six steps for computing transient stability problem

\begin{enumerate}
\item Obtain the generator data for the generator model.
\item Obtain the pre-fault network state information, e.g. network impedances.
(The overdetermined network will be reduced before forming the admittance
matrix.)
\item Use loadflow analysis to compute initial data of the generators' terminal
voltages and the mechanical power.
\item Obtain the post-fault network state information, e.g.
fault location, fault clearing time and breaker reclosings time. \label{fault}
\item Obtain the timing of the simulated system, e.g. time of switching and
the length of simulation time.
\item Compute the rotor angle using DAE solver. \label{Solve_DAE}
\end{enumerate}

The last three steps are massively parallel for screening fault scenarios.
All the data should be normalize to a common base unit.  A system base of
100 MVA is generally used.
Simulation is started before switching occurs in the network
to safeguard the simulation from inconsistent data.
In addition to the generator equation, we needed the followings
to simulate a complete system

\begin{itemize}
\item load or network constraints
\item the excitation system
\item mechanical torque
\end{itemize}

Detail model of these variables are nonlinear.  However, for short transient
stability simulation, a linearized model can sometime be used.

A linearize model can be used (with care) to simulate the network's
oscillating behavior.

The discontinuity across the network is caused by switching and fault.
Limits of controllers contributes to the local discontinuity for a particular
device.

Reduced network model is used to study the effect of a fault in the local area. 
Remote network is transformed to equivalent circuit and classical generator
model is used for remote generator.


\section{Classical Generator Model}

Assumptions used for classical model

\begin{itemize}
\item The mechanical power input remains constant during the period of the
transient. ($P_{m} =$ Constant)
\item Damping or asynchronous power is negligible (This can be simulated by
adding a damping term).
\item The synchronous machine can be represented (electrically) by a constant
voltage source behind a transient reactance.
(This implies the magnitude of the internal machine voltages
remain constant for the simulation time.)
\item The mechanical angle of the synchronous machine rotor coincides with the
electrical phase angle of the voltage behind transient reactance.
(The delta angle of the internal machine voltages induced by the transient
reactance is the delta rotor angle.)
\item If a local load is fed at the terminal voltage of the machine, it can be
represented by a constant impedance (or admittance) to neutral.
(No load change except the change cause by the faulty component.)
\end{itemize}

The change in delivered power is reflected by the change in current flow.
This model is designed only for the first swing (ie. accurate up to 1 to 2
sec simulation time).

\section{Classical Generator Model for Multimachine System}

Assumptions applying classical model to multimachine system

\begin{itemize}
\item Mechanical power input is constant.
\item Damping or asynchronous power is negligible (This can be simulated by
adding a damping term to each simulated generator).
\item Constant-voltage-behind-transient-reactance model for the synchronous
machines is valid.
\item The mechanical rotor angle of a machine coincides with the angle of the
voltage behind the transient reactance.
\item Loads are represented by passive impedances. (No non-linear load.)
\end{itemize}

  These assumptions only true for studying first swing period.

\section{Classical Generator equations for Multimachine System}

For $i = 1, 2, \cdots, n$

\begin{math}
\begin{array}{lcl}
\frac{2H_{i}}{\omega_{R}}\frac{d \omega_{i}}{dt} + D_{i}\omega_{i}
& = & P_{mi} - \left[ E_{i}^{2} G_{ii}
+ \sum_{j = 1, j \neq i}^{n} E_{i}E_{j}
\cos(\theta_{ij} - \delta_{i} +  \delta_{j}) \right] \\
\frac{d \delta_{i}}{dt} & = & \omega_{i} - \omega_{R} \\
\end{array}
\end{math}

At $(t = 0^{-})P_{mi0} = P_{ei0}$

\begin{math}
\begin{array}{lcl}
P_{mi0} & = & E_{i}^{2}G_{ii0} + \sum_{j = 1, j \neq i}^{n} E_{i}E_{j}Y_{ij0}
\cos(\theta_{ij0} - \delta_{i0} +  \delta_{j0})
\end{array}
\end{math}

$G_{ij}$ is the real part of the $i$ th, $j$ th entry of the $Y$ matrix.

$E_{i}$ is the voltage of the node $i$. In the simulation, this value will
obtain by a load flow calculation.

$D$ is the damping factor.

\section{2-axis Generator Model}

%\begin{math}
%\begin{array}{lcl}
%T'_{qi} \frac{d E'_{di}}{dt} & = & -E'_{di} + (x_{qi} - x'_{qi}) I_{qi} \\
%T'_{di} \frac{dE'_{qi}}{dt} & = & - E'_{qi} - (x_{di} - x'_{di})I_{di} +
%E_{fdi} \\
%\frac{2H_{i}}{\omega_{i}}\frac{d \omega_{i}}{dt} & = &
%P_{mi} - \left[ E'_{di} I_{di} + E'_{qi}I_{qi} \right] + (x'_{qi} -
%x'_{di})I_{di}I_{qi} - D_{i}(\omega_{i} - \omega_{0}) \\
%\frac{d\delta'_{i}}{dt} & = & \omega_{i} - \omega_{0} \\
%\end{array}
%\end{math}

Assumptions used for 2-axis generator model

\begin{itemize}
\item all faults are three phase fault.
\item neglected any subtransient effects
\end{itemize}

\section{The equation for 2-axis Generator Model in state-space from}

%
%   80% original size
\begin{figure}
\begin{picture}(424,492)
\psfig{figure=2-gen-m.eps,width=424pt,height=492pt}
\end{picture}
\caption{2-axis machine model \label{2-gen-m}}
\end{figure}

\begin{math}
\begin{array}{lcl}
\left[
\begin{array}{c}
\dot{i}_{d} \\
\dot{i}_{F} \\
\dot{i}_{D} \\
\dot{i}_{q} \\
\dot{i}_{Q} \\
\dot{\omega} \\
\dot{\delta} \\
\end{array}
\right] & = &
\left[
\begin{array}{cc}
\\
\\
-L^{-1}(R + \omega N) & 0 \\
\\
\\
\begin{array}{ccccc}
-\frac{L_{d} i_{q}}{3 \tau_{j}} & -\frac{k M_{F}i_{q}}{3 \tau{j}} &
-\frac{k M_{D} i_{q}}{3 \tau_{j}} & \frac{L_{q}i_{d}}{3 \tau_{j}} &
\frac{k M_{Q}i_{d}}{3 \tau_{j}} \\
0 & 0 & 0 & 0 & 0 \\
\end{array} &
\begin{array}{cc}
-\frac{D}{3 \tau_{j}} & 0 \\
1 & 0 \\
\end{array} \\
\end{array} \right]
\left[
\begin{array}{c}
i_{d} \\
i_{F} \\
i_{D} \\
i_{q} \\
i_{Q} \\
\omega \\
\delta \\
\end{array}
\right] +
\left[
\begin{array}{c}
\\
\\
-L^{-1}v \\
\\
\\
\frac{T_{m}}{3 \tau_{j}} \\
-1 \\
\end{array}
\right] \\
\end{array}
\end{math}

\begin{math}
\begin{array}{lcl}
L & = &
\left[
\begin{array}{ccccc}
L_{d} & k M_{F} & k M_{D} & 0 & 0 \\
k M_{F} & L_{F} & M_{R} & 0 & 0 \\
k M_{D} & M_{R} & L_{D} & 0 & 0 \\
0 & 0 & 0 & L_{q} & kM_{Q} \\
0 & 0 & 0 & k M_{Q} & L_{Q} \\
\end{array}
\right] \\
\end{array}
\end{math}

\begin{math}
\begin{array}{lcl}
R + \omega N & = &
\left[
\begin{array}{ccccc}
r & 0 & 0 & \omega L_{q} & \omega k M_{q} \\
0 & r_{F} & 0 & 0 & 0 \\
0 & 0 & r_{D} & 0 & 0 \\
-\omega L_{d} & -\omega k M_{F} & -\omega k M_{D} & r & 0 \\
0 & 0 & 0 & 0 & r_{Q} \\
\end{array}
\right]
\end{array}
\end{math}

This set of equation is the most general form of the two-axis model.
Simplified version of the above model is often used, when the model is used
to simulate a specific problem.

In addition to the above system of equation, we needed the followings
to simulate a complete system

\begin{itemize}
\item load or network constraints
\item the excitation system
\item mechanical torque
\end{itemize}

An alternate state-space model, which is called flux linkages model, with
the state variables, $\lambda_{d}$, $\lambda_{F}$, $\lambda_{D}$,
$\lambda_{q}$ and $\lambda_{Q}$.

The relationship between the current and flux linkage is give by

\begin{math}
\begin{array}{lcl}
\left[
\begin{array}{c}
\lambda_{d} \\
\lambda_{F} \\
\lambda_{D} \\
\lambda_{q} \\
\lambda_{Q} \\
\end{array}
\right] & = & 
\left[
\begin{array}{ccccc}
L_{d} & 0 & k M_{F} & k M_{D} & 0 \\
0 & L_{q} & 0 & 0 & k M_{Q} \\
k M_{F} & 0 & L_{F} & M_{R} & 0 \\
k M_{D} & 0 & M_{R} & L_{D} & 0 \\
0 & k M_{Q} & 0 & 0 & L_{Q} \\
\end{array}
\right]
\left[
\begin{array}{c}
i_{d} \\
i_{F} \\
i_{D} \\
i_{q} \\
i_{Q} \\
\end{array}
\right]
\end{array}
\end{math}

The system of equations is for generic generators.
Consequently, some variables may not appropriately used to describe a
specific type of generators.  Those is generally set to zero.

\section{EPRI's ETMSP}

ETMSP uses explicit RK4 at time of switching with step size zero to handle
the discontinuity.  Then ETMSP can switch over to implicit trapezoidal rule.
ETMSP ignores the local discontinuity (This reduce the accuracy of the
simulation).

The assumptions of generator model used in ETMSP

\begin{itemize}
\item Up to three rotor windings at each axis can be represented
\item d and q axes saturation can be separately represented
\item Subtransient saliency is fully represented
\item An interface is provided for excitation system and speed governing
system models
\item Rate of change of stator flux linkages and the effect of speed
variations are neglected in the stator voltage equations (Complexity of the
model is reduced).
\end{itemize}

%
%   75% original size
\begin{figure}
\begin{picture}(460,486)
\psfig{figure=epri-2-gen-m.eps,width=460pt,height=486pt}
\end{picture}
\caption{EPRI detail synchronous machine model \label{epri-2-gen-m}}
\end{figure}

Schematic of the 2 axis generator model used by ETMSP is shown in
figure~\ref{epri-2-gen-m}.

\subsection{State Equation}

\begin{math}
\begin{array}{lcl}
\frac{d \Psi_{f}}{dt} & = & -\omega_{0} r_{f} i_{f} + \omega_{0} V_{f} \\
\frac{d \Psi_{kd1}}{dt} & = & -\omega_{0} r_{kd1} i_{kd1} \\
\frac{d \Psi_{kd2}}{dt} & = & -\omega_{0} r_{kd2} i_{kd2} \\
\frac{d \Psi_{kq1}}{dt} & = & -\omega_{0} r_{kq1} i_{kq1} \\
\frac{d \Psi_{kq2}}{dt} & = & -\omega_{0} r_{kq2} i_{kq2} \\
\frac{d \Psi_{kq3}}{dt} & = & -\omega_{0} r_{kq3} i_{kq3} \\
\frac{d \omega_{m}}{dt} & = & \frac{\omega_{0}}{2H}[T_{m} - (\Psi_{d}i_{q} -
\Psi_{q}i_{d})] \\
\frac{d \delta}{dt} & = & \omega_{m} \\
\end{array}
\end{math}

\begin{tabbing}
$\omega_{m}$ \= - \=   \kill
$t$ \> - \> time (sec) \\
$\delta$ \> - \> rotor angle (rad) \\
$\omega_{m}$ \> - \> rotor speed (rad/sec) \\
$V_{f}$ \> - \> the field voltage determined by the excitation control system
(p.u.) {$V_{f}$ is constant, if no excitation control} \\
$T_{m}$ \> - \> the mechanical torque determined by the speed governing
system. (p.u.) {$T_{m}$is constant, if no speed governor} \\
\end{tabbing}

$\omega_{m}, \delta, \Psi_{f}, \Psi_{kd1}, \Psi_{kd2}, \Psi_{kq1}, \Psi_{kq2},
\Psi_{kq3}$ are the state variable of the generator model.
The damper winding may not present.

\subsection{Stator Equations}

\begin{math}
\begin{array}{lcl}
V_{d} & = & -r_{a} i_{d} - \Psi_{q} \\
V_{q} & = & -r_{a} i_{q} - \Psi_{d} \\
\end{array}
\end{math}

\begin{tabbing}
$V_{d}$ \= - \=  \kill
$V_{d}$ \> - \> machine terminal voltage in machine d reference frame. \\
$V_{q}$ \> - \> machine terminal voltage in machine q reference frame. \\
$i_{d}$ \> - \> machine terminal current in machine d reference frame. \\
$i_{q}$ \> - \> machine terminal current in machine q reference frame. \\
\end{tabbing}

\subsection{Flux Linkage Equation}

\begin{math}
\begin{array}{lcl}
\left[
\begin{array}{c}
\Psi_{d} \\
\Psi_{f} \\
\Psi_{kd1} \\
\Psi_{kd2} \\
\end{array}
\right] & = &
\left[
\begin{array}{cccc}
X_{d} & X_{ad} & X_{ad} & X_{ad} \\
X_{ad} & X_{f} & X_{fad1} & X_{fad2} \\
X_{ad} & X_{fad1} & X_{ad1} & X_{ad12} \\
X_{ad} & X_{fad2} & X_{ad12} & X_{ad2} \\
\end{array}
\right]
\left[
\begin{array}{c}
-i_{d} \\
i_{f} \\
i_{kd1} \\
i_{kd2} \\
\end{array}
\right] \\
\end{array}
\end{math}

\begin{math}
\begin{array}{lcl}
\left[
\begin{array}{c}
\Psi_{q} \\
\Psi_{kq1} \\
\Psi_{kq2} \\
\Psi_{kq2} \\
\end{array}
\right] & = &
\left[
\begin{array}{cccc}
X_{q} & X_{aq} & X_{aq} & X_{aq} \\
X_{aq} & X_{aq1} & X_{aq} & X_{aq} \\
X_{aq} & X_{aq} & X_{aq2} & X_{aq} \\
X_{aq} & X_{aq} & X_{aq} & X_{aq3} \\
\end{array}
\right]
\left[
\begin{array}{c}
-i_{q} \\
i_{kq1} \\
i_{kq2} \\
i_{kq3} \\
\end{array}
\right] \\
\end{array}
\end{math}

\begin{math}
\begin{array}{lcl}
X_{d} & = & X_{ad} + X_{l} \\
X_{f} & = & X_{ad} + X_{fd} + X_{c1} + X_{c2} \\
X_{fad1} & = & X_{ad} + X_{c1} \\
X_{fad2} & = & X_{ad} + X_{c1} + X_{c2} \\
X_{ad1} & = & X_{ad} + X_{kdl} + X_{c1} \\
X_{ad12} & = & X_{ad} + X_{c1} \\
X_{ad2} & = & X_{ad} + Z_{kd2} + X_{c1} + X_{c2} \\
X_{q} & = & X_{aq} + X_{l} \\
X_{aq1} & = & X_{aq} + X_{kq1} \\
X_{aq2} & = & X_{aq} + X_{kq2} \\
X_{aq3} & = & X_{aq} + X_{kq3} \\
\end{array}
\end{math}

$X_{ad}$ and $X_{aq}$ are assumed to be determined by the generator's
magnetic saturation.

The saturation characteristic is partitioned into three zones. The saturation
characteristic is determined by open circuit.

\subsubsection{A linear air gap region}

The system p.u. gives this region a slope of one.

\subsubsection{A nonlinear knee region}

The knee region is model by exponential function
\begin{math}
\begin{array}{lcl}
V_{a} - V & = & A_{sat}e^{[B_{sat}(V - V_{L})]} \\
\end{array}
\end{math}

The knee reactance $X_{s}$ is model by

\begin{math}
\begin{array}{lcl}
X_{s} & = & X_{u}\frac{V}{V_{a}} \\
& = & \frac{X_{u}}{1 + \frac{A_{sat}e^{[B_{sat}(V - V_{L})]}}{V}} \\
\end{array}
\end{math}

\begin{tabbing}
$X_{u}$ \= - \= \kill
$X_{u}$ \= - \= unsaturated magnetising reactance \\
$V_{L}$ \= - \= knee region starting voltage \\
$V_{a}$ \= - \= air gap voltage corresponding to $V$ \\
\end{tabbing}

$V_{a} - V$ is not zero at $V_{L}$.  Hence, discontinuity is introduced
at $V_{L}$.

\subsubsection{A linear heavily saturated region}

The saturated reactance $X_{s}$ is model by
\begin{math}
\begin{array}{lcl}
X_{s} & = & X_{u}\frac{V}{V_{a}} \\
& = & X_{u}\frac{V}{R_{s}(V - V_{m}) + V_{ma}} \\
\end{array}
\end{math}

\begin{math}
\begin{array}{lcl}
V_{ma} & = & V_{m} + A_{sat}e^{[B_{sat}(V - V_{L})]} \\
\end{array}
\end{math}

\begin{tabbing}
$X_{u}$ \= - \= \kill
$X_{u}$ \= - \= unsaturated magnetising reactance \\
$R_{s}$ \= - \= slope for heavily saturated region \\
$V_{m}$ \= - \= knee region ending voltage \\
\end{tabbing}

\subsection{Load Saturation}

When under load, the saturation is approximated from the open circuit data.

\begin{math}
\begin{array}{lcl}
\Psi_{a} & = & \sqrt{\Psi_{ad}^{2} + \Psi_{aq}^{2}}q \\
\Psi_{ad} & = & \Psi_{d} + X_{l} i_{d} \\
\Psi_{aq} & = & \Psi_{q} + X_{l} i_{q} \\
\end{array}
\end{math}

\subsection{Interfacing with the Network}

In order to interface the generator mode to the rest of the network, we need
to transform the actual voltages to synchronous rotating frame.

\begin{math}
\begin{array}{lcl}
\left[
\begin{array}{c}
V_{D} \\
V_{Q} \\
\end{array}
\right] & = &
\left[
\begin{array}{cc}
\cos \delta & - \sin \delta \\
\sin \delta & \cos \delta \\
\end{array}
\right]
\left[
\begin{array}{c}
V_{d} \\
V_{q} \\
\end{array}
\right] \\
\end{array}
\end{math}

\begin{math}
\begin{array}{lcl}
\left[
\begin{array}{c}
i_{D} \\
i_{Q} \\
\end{array}
\right] & = &
\left[
\begin{array}{cc}
\cos \delta & - \sin \delta \\
\sin \delta & \cos \delta \\
\end{array}
\right]
\left[
\begin{array}{c}
i_{d} \\
i_{q} \\
\end{array}
\right] \\
\end{array}
\end{math}

The transformation, $T$, is orthonormal. Hence, the inverse transformation
$T^{t}$.
$i_{D}$, $i_{Q}$, $V_{D}$ and $V_{Q}$ is determined by load-flow calculation.
In directly, $i_{d}$, $i_{q}$, $V_{d}$ and $V_{q}$ is also determined by
load-flow calculation.

The steady state load angle is given by

\begin{math}
\begin{array}{lcl}
\tan(\beta + \delta) & = & \frac{X_{q} P - r_{a} Q}
{\mid V \mid^{2} + r_{a} P + X_{q} Q} \\
\end{array}
\end{math}

\begin{tabbing}
$\tan \beta$ \= = \= \kill
$\tan \beta$ \> = \> $\frac{V_{D}}{V_{Q}}$ \\
$\mid V \mid^{2}$ \> = \> $V_{D}^{2} + V_{Q}^{2}$ \\
$P$ \> = \> generated active power \\
$Q$ \> = \> generated reactive power \\
\end{tabbing}

\section{EPRI's ETMSP Classical Synchronous Machine Model}

%
%   100% original size
\begin{figure}
\begin{picture}(328,201)
\psfig{figure=epri-c-gen-m.eps,width=328pt,height=201pt}
\end{picture}
\caption{EPRI detail synchronous machine model \label{epri-c-gen-m}}
\end{figure}

Schematic of the classical generator model used by ETMSP is shown in
figure~\ref{epri-c-gen-m}

\begin{math}
\begin{array}{lcl}
\frac{d \omega_{m}}{dt} & = & \frac{\omega_{0}}{2H}(T_{m} - V_{D}i_{D}
- V_{Q}i_{D}) \\
\frac{d \delta}{dt} & = & \omega_{m} \\
\end{array}
\end{math}

\subsection{Interfacing with the Network}

\begin{math}
\begin{array}{lcl}
V_{D} & = & E' \cos \delta + X'_{d} I_{Q} \\
V_{Q} & = & E' \sin \delta + X'_{d} I_{D} \\
\end{array}
\end{math}

$V_{D}$ and $V_{Q}$ is determined by load-flow calculation.
Indirectly, $E'$ is also determined by load-flow calculation and it is held
constant over the cause of simulation.

\end{document}
\endinput


