\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{mathrsfs}
\usepackage{graphicx}


\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2017 (2017), No. 298, pp. 1--15.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2017 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2017/298\hfil Dynamical analysis of OMIB]
{Dynamical analysis of one machine to infinite bus power systems under Gauss
type random excitation}

\author[L. Liu, P. Ju, F. Wu \hfil EJDE-2017/298\hfilneg]
{Lei Liu, Ping Ju, Feng Wu}

\address{Lei Liu \newline
 College of Science,
 Hohai University,
Nanjing, 210098, China}
 \email{liulei\_hust@163.com}

\address{Ping Ju (corresponding author)\newline
 College of Energy and Electrical Engineering,
Hohai University, Nanjing, 210098, China}
\email{pju@hhu.edu.cn}

\address{Feng Wu \newline
 College of Energy and Electrical Engineering,
Hohai University, Nanjing, 210098, China}
\email{wufeng@hhu.edu.cn}

\dedicatory{Communicated by Zhaosheng Feng}

\thanks{Submitted October 8, 2016. Published November 30, 2017.}
\subjclass[2010]{34F05, 60H10, 93E03}
\keywords{Power systems; Random excitation; stochastic fluctuation;
\hfill\break\indent stationary distribution; asymptotic pathwise estimation;
asymptotic moment estimation}

\begin{abstract}
 We discuss the asymptotic behavior of the stochastic one machine
 to infinite bus power systems. Using the exponential martingale inequality
 and the Borel-Cantelli Lemma, we obtain asymptotic moment estimation and
 asymptotic pathwise estimation of the stochastic one machine to infinite
 bus systems. Using the ergodic properties, we give a good explanation of the
 fluctuation phenomena. By means of the property of periodicity, Hormander's
 theorem and a detailed balance method, the existence and probability density
function of the stationary distribution on the cylindrical are illustrated.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks


\section{Introduction}

In the previous decades, dynamics of deterministic power systems has
received a lot of attention, see, e.g.
\cite{ju1996b,ju1996a,pan2013,wu2013,wu2007,wu2008},
 and references therein. By nature, a power system is continually
experiencing stochastic disturbances, such as, switching events,
 load level fluctuations, which may have a significant effect on the
operation of power systems and quality of electricity. In recent years,
many researchers have been carried out to study the dynamics of power systems
under random excitations, see, e.g.,
\cite{Hiskens2006,Hockenberr2004,ju1991,liu2011,Timko1983,wangkw2000,Wu1983}.
The random factors in power systems were classified into three categories
in \cite{zhang2012}, namely, random initial values, random parameters,
and random excitation. With the increase integration of the renewable energy
generation system and electric vehicles, and the features of randomness
into the power system, the dynamics of power system under random excitations
 has received more and more attentions (see
\cite{lu2015,Stankovic,Vlachogiannis}).


Since the Gauss white noise is a well-known mathematical interpretation for
the random excitations, the application of It\^{o} stochastic differential
equations (SDEs) has been taken into account in the research on power systems.
The detailed understanding of SDEs can be founded in
\cite{Hasminskii2011,10}. In this paper, we start our analysis by considering
the stochastic one machine to infinite bus (OMIB) power system
(see Figure \ref{fig1}(top))
perturbed with random excitations. Then the stochastic OMIB system can be
described by the It\^{o} equation as follows:
\begin{equation}\label{0}
\begin{gathered}
 {\rm d}\delta=(\omega-1){\rm d}t \\
 {\rm d}\omega=\big(-\frac{D}{M}(\omega-1)+\frac{1}{M}(P_{M}-P_{e})\big){\rm
 d}t+\sigma{\rm d}B(t),
 \end{gathered}
\end{equation}
where $\delta$ is the rotor angle; $\omega$ is the rotor rotating speed;
$P_{e}=P_{\rm max}\sin \delta$ is the electrical power; $P_{M}$ is the
mechanical power and is assumed to be constant, i.e.
$P_{M}=P_{\rm max}\sin \delta_{s}$; $M$ is the inertia constant in per unit and
taken $M=2569.8288$.

It is necessary to reveal how the noise affects the stochastic OMIB systems
\eqref{0}. The mean stability and mean square stability of the corresponding
linear systems of \eqref{0} were analyzed theoretically by Zhang
et al \cite{zhang2012}. Using the stochastic averaging methods,
\cite{zhu2010} has studied the first-passage problem of dynamical power
system of \eqref{0}. Recently, Keyou Wang \cite{keyou} has applied
the Fokker-Planck Equation to study the evolution of the probability
density function of the system \eqref{0}.

However, the asymptotic properties of stochastic OMIB systems \eqref{0}
has not been fully investigated. In addition, if we make a great number
of records of $\omega$ to investigate the dynamic behavior of a
stochastic OMIB power system \eqref{0}, it can be found that a single
record may fluctuate even if the number of records is large.
Then two questions arise naturally:
(1) Can the fluctuation phenomena be given a explanation?
(2)Is there a stationary distribution to the system? As a result,
to solve these two problems is the main motivation of this paper.
The primary contributions of this paper are as follows:
\begin{itemize}
\item The $p$-th moment estimation and asymptotic pathwise of $\omega$ were
obtained using the exponential martingale inequality and Gronwall's inequality;

\item Utilizing the ergodic property and the comparative method, a good
explanation of the fluctuation phenomena is given;

\item With the help of the HormanderĄ¯s theorem and detailed balance method,
 the existence of the stationary distribution on cylindrical is proved.
\end{itemize}

\section{Notation and problem statement} \label{stable}
Throughout this paper, unless otherwise specified, we let
$(\Omega,\mathscr{F},\{\mathscr{F}_{t}\}_{t\geq 0}, \mathbb{P})$ be
a complete probability space with a filtration
$\{\mathscr{F}_{t}\}_{t\geq 0}$ satisfying the usual conditions
(i.e., it is increasing and right continuous while $\mathscr{F}_0$
contains all $\mathbb{P}$-null sets). Let ${ B}(t)$ be an
one-dimensional Brownian motion defined on the probability space.

In the same way as Mao et al.\ \cite{10} did, we can show the
following result on the existence of a global solution.

\begin{lemma} \label{lem2.1}
For any given initial data
$x_0\in R^2$, there is a unique global solution
$\big(\delta(t,x_0),\omega(t,x_0)\big)$ of system \eqref{0}
on $[0,+\infty)$.
\end{lemma}

The dynamics of the deterministic classical OMIB system has been extensively
studied; see \cite{chen2005,kundur}.
These existing literature show clearly that the point $(\delta_{s}, 1)$
is a stable equilibrium point (SEP).
The trajectory will converge to the SEP if the initial
conditions are selected to lie in the attraction region of the SEP.
Applying a random excitation, the dynamic behavior becomes more
complicated than the deterministic classical OMIB system. The simulations
 show that a single record of $\omega(t)$ may fluctuate even if the number
of records is large(see Figure \ref{fig1}(bottom)).

By means of the stochastic analysis techniques and the ergodic properties,
 some estimation and a good explanation of the fluctuation phenomena will
 be provided in Section III and Section IV, respectively.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth, clip = true, bb = 0 0 170 35]{fig1a} % smib1.eps}
\includegraphics[width=0.7\textwidth]{fig1b} \\ % y.eps}
\end{center}
\caption{(top) One machine infinite bus system.
(bottom) Stochastic trajectory of
 $\omega(t)$ generated by the Heun scheme
for time step $h=2^{-8}$ for \eqref{0} with $\sigma=0.6$ on $[0,5000]$.}
\label{fig1}
\end{figure}

\section{Estimation on the stochastic fluctuation}

The aim of this section is to estimate the stochastic
fluctuation for the rotor speed $\omega$.

\subsection{Moment estimation}

\begin{theorem} \label{thm3.1}
For any $p>0$ and initial data ${x}_0=(\delta_{s},1)$, there exists a constant
$M_{p}>0$ such that the global solution $\omega(t)$ of system \eqref{0} has
the property
$$
\sup_{0\leq t<\infty} E|{\omega}(t,x_0)-1|^{p}\leq M_{p}.
$$
\end{theorem}

\begin{proof}
 For any $p>0$, let $V({\omega}, t)=e^{\varepsilon t}|\omega-1|^{p}
=e^{\varepsilon t}\big((\omega-1)^2\big)^{\frac{p}{2}}$, where
$\varepsilon$ is a sufficient small positive number.
Let $\omega(t) = \omega(t; x_0)$ for
simplicity. Applying It\^{o} formula and Young inequality (see
Mao \cite{10}) to $V(\omega,t)$ yields
\begin{align*}
&\mathscr{L}|\omega-1|^{p} \\
&=e^{\varepsilon t}\frac{p}{2}|\omega-1|^{p-2}\big[\big(-\frac{2D}{M}
 +\frac{2}{p}\varepsilon\big)(\omega-1)^2\\
&\quad+\big(\frac{2P_{\rm max}}{M}(\sin\delta_{s}-\sin \delta)+2\big)
 (\omega-1)+(p-1)\sigma^2\big]\\
&\leq e^{\varepsilon t}\frac{p}{2}|\omega-1|^{p-2}
 \big[\big(-\frac{2D}{M}+\frac{2}{p}\varepsilon\big)
 (\omega-1)^2+\big(\frac{4P_{\rm max}}{M}+2\big)|\omega-1|+(p-1)\sigma^2\big]\\
&=e^{\varepsilon t}\frac{p}{2}\big[\big(-\frac{2D}{M}
 +\frac{2}{p}\varepsilon\big)|\omega-1|^{p}
 +\big(\frac{4P_{\rm max}}{M}+2\big)|\omega-1|^{p-1}
 +(p-1)\sigma^2|\omega-1|^{p-2}\big].
\end{align*}
Note that $(-\frac{2D}{M}+\frac{2}{p}\varepsilon)<0$. Then
there exists a constant $K_{P}$ such that
\[
\sup_{0\leq u \leq \infty}\frac{p}{2}
[(-\frac{2D}{M}+\frac{2}{p}\varepsilon)u^{p}
+\frac{4P_{\rm max}}{M}u^{p-1}+(p-1)\sigma^2u^{p-2}]\leq K_{p}.
\]
 This implies
\[
Ee^{\varepsilon t}|\omega(t)-1|^{p}
\leq E\big((\omega-1)^2(0)\big)^{\frac{p}{2}}
+ \int^{t}_0K e^{\varepsilon s}ds\leq E|\omega(0)-1|^{p}
+ \frac{K_{p}}{\varepsilon}( e^{\varepsilon t}-1),
\]
which means
\[
E|\omega(t)-1|^{p}\leq e^{-\varepsilon t}E|\omega(0)-1|^{p}
+\frac{K_{p}}{\varepsilon}(1-e^{\varepsilon -t})
\leq E|\omega(0)-1|^{p}+ \frac{K_{p}}{\varepsilon}:=M_{p}
\]
We can claim that for any $p>0$ there exists a $M_{p}>0$
such that
$$
\sup_{ 0\leq t<\infty} E|\omega(t)-1|^{p}\leq M_{p}.
$$
The proof is complete.
\end{proof}


\subsection{Asymptotic pathwise estimation}

Theorem \ref{thm3.1} shows that the $p$-th moment of $|\omega(t)-1|$ is boundedness.
Now we process to derive some asymptotic pathwise estimation of
$|\omega(t)-1|$ by applying the exponential martingale inequality and
Borel-Cantelli Lemma (see \cite{10}), which shows the solution of
\eqref{0} how to vary in $R^2$ pathwisely.

\begin{theorem} \label{thm3.2}
For any $p>0$ and initial data ${x}_0=(\delta_{s},1)$, the
global solution $\omega(t,x_0)$ of system \eqref{0} has the
property
\begin{equation}\label{3.14}
\lim_{t\to \infty}\sup\frac{1}{t}\int^{t}_0\big(\omega(s, X_0)-1\big)^2ds
<\Big(\sqrt{\frac{\sigma^2}{D^2}+\frac{2P_{\rm max}}{D}}
+\frac{2P_{\rm max}}{D}\Big)^2.
\end{equation}
\end{theorem}

\begin{proof}
 Let $\omega(t) = \omega(t; x_0)$ for
simplicity. Applying It\^{o}'s formula to $\omega^2$ yields
\begin{equation}\label{3.13}
\begin{aligned}
&\big(\omega(t)-1\big)^2\\
&=\big(\omega(0)-1\big)^2+\int^{t}_0
\Big(-\frac{2D}{M}\big(\omega(s)-1\big)^2 \\
&\quad+\big(\frac{2P_{\rm max}}{M}(\sin
\delta_{s}-\sin \delta(s))+\frac{2D}{M}\big)(\omega(s)-1)+\sigma^2\Big)ds+M(t).
\end{aligned}
\end{equation}
where $M(t)=\int^{t}_02(\omega(s))\sigma dB(s)$ is a continuous
local martingale. The quadratic variation of $M(t)$ is
$\langle M(t)\rangle =\int^{t}_04\sigma^2\big(\omega(t)-1\big)^2 d(s)$.
 For any $\alpha\in (0, \frac{D}{M\sigma^2})$ and every integer $k\geq1$,
using the exponential martingale inequality
(cf. Mao \cite[Theorem 1.7.4]{10}) we have
\[
\mathbb{P}\big( \sup_{0\leq t\leq
T}(M(t)-\frac{\alpha}{2}\langle M(t)\rangle )>\frac{1}{\alpha}\log k\big)
\leq \frac{1}{k^2},\quad k = 1, 2, \ldots,
\]
An application of the well-known
Borel-Cantelli lemma then yields that for almost all
$\varpi\in \Omega$ there is a random integer $k_0=k_0(\varpi)\geq1$ such
that
\begin{equation}\label{3.12}
M(t)\leq \frac{2}{\alpha}\log k + 0.5\alpha \langle M(t)\rangle,
\end{equation}
for $t \in [0, k], k\geq k_0$, almost surely.
Substituting \eqref{3.12} into \eqref{3.13},
we obtain that
\begin{align*}
&\big(\omega(t)-1\big)^2\\
&=\big(\omega(0)-1\big)^2+\int^{t}_0
 \Big(-\frac{2D}{M}\big(\omega(s)-1\big)^2+\big(\frac{2P_{\rm max}}{M}(\sin
\delta_{s}-\sin \delta
(s)) \\
&\quad+\frac{2D}{M}\big)(\omega(s)-1)+\sigma^2\Big)ds
 +\frac{2}{\alpha}\log k+0.5\alpha\int^{t}_04\sigma^2\big(\omega(s)-1\big)^2ds\\
&\leq \omega(0)+\int^{t}_0\Big((-\frac{2D}{M}+2\alpha\sigma^2)
\big(\omega(s)-1\big)^2+\big(\frac{4P_{\rm max}}{M}\sin
\delta_{s}+\frac{2D}{M}\big)|\omega(s)-1|\\
&\quad+\sigma^2\Big)ds+\frac{2}{\alpha}\log k,
\end{align*}
for $t \in [0, k], k\geq k_0$, almost surely.
 Setting $h(\alpha)=\frac{2D}{M}-2\alpha\sigma^2$. For any $\eta\in(0,1)$,
we have
\begin{align*}
\eta h(\alpha)\int^{t}_0\omega^2(s)ds
&\leq \omega^2(0)
+\int^{t}_0\Big((1-\eta)h(\alpha)\omega^2(s)\\
&\quad +\big(\frac{4P_{\rm max}}{M}\sin \delta_0+\frac{2D}{M}\big)|\omega(s)|
+\sigma^2\Big)ds+\frac{2}{\alpha}\log k.
\end{align*}
For almost all $\varpi\in\Omega$, if $k \geq k_0$
and $k-1 \leq t \leq k,$ simple computations show that
\begin{align*}
&\frac{1}{t}\int^{t}_0\omega^2(s)ds \leq
\frac{1}{\eta h(\alpha)(k-1)}\big(\omega^2(0)
+(4\frac{P^2_{m}}{(1-\eta)h(\alpha)M^2}
+\sigma^2)k+\frac{2}{\alpha}\log k\big).
\end{align*}
Letting $t\to \infty$ ($k\to \infty$) and $\alpha\to 0$ yields
\[
\frac{1}{t}\int^{t}_0\omega^2(s)ds
\leq \frac{P^2_{\rm max}}{\eta(1-\eta)D^2}+\frac{M\sigma^2}{2D \eta}.
\]
Setting $\eta=1/(1+P_{\rm max}/\sqrt{2P^2_{m}+\sigma^2MD})$ yields
\[
\lim_{t\to \infty}\sup\frac{1}{t}\int^{t}_0\big(\omega(s)-1\big)^2ds
<\Big(\sqrt{\frac{P^2_{\rm max}}{D^2}+\frac{\sigma^2M}{2D}}
+\frac{P_{\rm max}M}{D}\Big)^2.
\]
The proof is complete.
\end{proof}



\begin{theorem} \label{thm3.3}
For initial data ${x}_0=(\delta_{s},1)$, the global solution $\omega(t; x_0)$
of system \eqref{0} has the property
$$
\limsup_{t\to \infty} \frac{|\omega(t,x_0)-1|}{\sqrt{\log t}}
\leq\sqrt{\frac{Me}{D}}\sigma.
$$
\end{theorem}

\begin{proof}
Let $\omega(t)=\omega(t; x_0)$ for simplicity.
Simple computation shows that
\begin{align*}
(\omega-1)(\big(-\frac{D}{M}(\omega-1)+\frac{1}{M}(P_{M}-P_{e})\big))
&\leq\frac{2P_{\max}}{M}|(\omega-1)|-\frac{D}{M}(\omega-1)^2 \\
&\leq-(\frac{D}{M}-\frac{P_{\max}}{M}\epsilon)(\omega-1)^2
+\frac{P_{\max}}{M\epsilon}.
\end{align*}
Set
\[
\lambda(\epsilon)=\frac{D}{M}-\frac{P_{\max}}{M}\epsilon,\quad
b=\frac{P_{\max}}{M\epsilon}.
\]
 Applying It\^{o} formula and Young inequality (see Mao \cite{10}) to
$e^{-2\lambda(\epsilon) t}(\omega-1)^2$ yields
\begin{equation}\label{3.11}
\begin{aligned}
e^{-2\lambda(\epsilon)t}(\omega(t)-1)^2
&=(\omega-1)^2(0)+\int^{t}_0e^{-2\lambda(\epsilon)s}(2b+\sigma^2)ds\\
&\quad +2\sigma\int^{t}_0e^{-2\lambda(\epsilon) s}(\omega-1)(s)dB(s),
\end{aligned}
\end{equation}
where $M(t)=2\sigma\int^{t}_0e^{-2\lambda(\epsilon)
s}(\omega-1)(s)dB(s)$ is a continuous local martingale with the
quadratic variation
$$
\langle M(t)\rangle
=4\sigma^2\int^{t}_0e^{-4\lambda(\epsilon) s}|(\omega-1)(s)|^2{\rm d}s.
$$
For each integer $n>0$, $\varepsilon>0$, $\gamma>0$, $\theta>1$, the
exponential martingale inequality (see
\cite[Theorem 7.4 on page44]{Vlachogiannis}) yields
\[
\mathbb{P}\big( \sup_{0\leq t\leq
t_{n+1}}(M(t)-\frac{\gamma_{n}}{2}\langle M(t)\rangle)>\frac{1}{\gamma_{n}}\log
n^{\theta}\big)
\leq \frac{\log n^{\theta}}{\gamma_{n}}, \quad n = 1, 2, \ldots ,
\]
where $t_{n}=n\varepsilon, \gamma_{n}=\gamma
e^{2n\lambda(\epsilon)\varepsilon}$. Noting that
$\sum^{\infty}_{n=1}\frac{1}{n^{\theta}}<\infty$, by the well-known
Borel-Cantelli lemma, there exists $\Omega_0\subset\Omega$ with
$P(\Omega_0)=1$, such that for any $\omega\in \Omega_0$, there
exists an integer $N(\omega)$ such that for all $n\geq N(\omega)$
and $0 \leq t\leq t_{n+1}$
\begin{equation}\label{3.10}
M(t)\leq \frac{\gamma_{n}}{2}\langle M(t)\rangle+\frac{1}{\gamma_{n}}\log
n^{\theta}.
\end{equation}
Setting $h(t)=e^{-2\lambda(\epsilon) t}(\omega(t)-1)^2$,
Substituting \eqref{3.10} into \eqref{3.11}, with probability one,
gives,
\begin{align*}
&h(t)-h(0) \\
&\leq \frac{2b+\sigma^2}{2|\lambda(\epsilon)|}
 e^{-2\lambda(\epsilon) t}+2\sigma^2\gamma_{n}\int^{t}_0
 e^{-2\lambda(\epsilon) s}h(s)ds+\frac{\log n^{\theta}}{\gamma_{n}}\\
&\leq \frac{2b+\sigma^2}{2|\lambda(\epsilon)|}
 e^{-2\lambda(\epsilon)(n+1)\varepsilon}
 +\frac{\log n^{\theta}}{\gamma}e^{-2\lambda(\epsilon) n\varepsilon}
 +2\sigma^2\gamma e^{2n\lambda(\epsilon)\varepsilon}\int^{t}_0
 e^{-2\lambda(\epsilon) s}h(s)ds.
\end{align*}
Then Gronwall's inequality implies
\[
h(t)\leq(h(0)+\frac{2b+\sigma^2}{2|\lambda(\epsilon)|}
e^{-2\lambda(\epsilon)(n+1)\varepsilon}\frac{\log
n^{\theta}}{\gamma}e^{-2\lambda(\epsilon)
n\varepsilon})\exp(\frac{\sigma^2}{|\lambda(\epsilon)|}\gamma
e^{2n\lambda(\epsilon) \varepsilon}e^{-2\lambda(\epsilon) t}).
\]
Consequently, for almost all $\omega \in \Omega_0$, if $k \geq N$
and $n\varepsilon \leq t \leq (n+1)\varepsilon$,
\begin{align*}
&(\omega(t)-1)^2 \\
&\leq\Big((\omega(0)-1)^2e^{-2\lambda(\epsilon)
t}+\frac{\log
n^{\theta}}{\gamma}e^{-2\lambda(\epsilon)(
n\varepsilon-t)}+\frac{2b+\sigma^2}{2|\lambda(\epsilon)|}e^{-2\lambda(\epsilon)((n+1)\varepsilon-t)}\Big)\\
&\quad\times\exp(\frac{\sigma^2}{|\lambda(\epsilon)|}\gamma
e^{2n\lambda(\epsilon) \varepsilon}e^{-2\lambda(\epsilon) t})\\
&\leq\Big((\omega(0)-1)^2e^{-2\lambda(\epsilon)
t}+\frac{2b+\sigma^2}{2|\lambda(\epsilon)|}e^{-2\lambda(\epsilon)((n+1)\varepsilon-t)}+\frac{\theta(\log t-\log \varepsilon)}{\gamma}e^{-2\lambda(\epsilon)(
n\varepsilon-t)}\Big)\\
&\quad\times\exp(\frac{\sigma^2}{|\lambda(\epsilon)|}\gamma
e^{2n\lambda(\epsilon) \varepsilon}e^{-2\lambda(\epsilon) t}).
\end{align*}
We therefore have
\begin{align*}
\frac{(\omega(t)-1)^2}{\log t}
&\leq((\omega(0)-1)^2e^{-2\lambda(\epsilon)
t}+\frac{2b+\sigma^2}{2|\lambda(\epsilon)|}
e^{-2\lambda(\epsilon)((n+1)\varepsilon-t)}\\
&\quad+\frac{\theta(\log t-\log \varepsilon)}{\log t\gamma}e^{-2\lambda(\epsilon)(
n\varepsilon-t)})\exp(\frac{\sigma^2}{|\lambda(\epsilon)|}\gamma
e^{-2\lambda(\epsilon)(t-n\varepsilon)}).
\end{align*}
Letting $t\to \infty $ (forcing $k\to \infty$)
and $\epsilon\to 0$ yields
$$
\limsup_{t\to \infty}\frac{(\omega(t)-1)^2}{\log t}
\leq \frac{\theta}{\gamma}\exp(\frac{\sigma^2\gamma M}{D}
e^{-2\frac{D}{M}\varepsilon}).
$$
Letting $\theta\to 1$, $\varepsilon\to 0$, and then
setting $\gamma=\frac{D}{M\sigma^2}$ yields
$$
\limsup_{t\to \infty}\frac{|\omega(t)|}{\sqrt{\log t}}
\leq \sqrt{\frac{Me}{D}}\sigma.
$$
The proof is complete.
\end{proof}

\subsection{Numerical examples}

An stochastic OMIB power system is used in this paper to study power system dynamics
under random excitation. In the OMIB, the transformer reactance
$x_{T_1}=0.138$ pu and
$x_{T_2}=0.122$ pu; the line reactance $x_1=0.234$ pu; the
generator transient reactance $x'_d= 0.295$ pu; the inertia
constant $M=2569.8288$ pu; and the damping coefficient $D=2.0$ pu.
The initial system power $P_0=1.0$ pu, $Q_0=0.2$ pu, voltages
behind the reactance ${E}'=1.41$ pu and rotor angle
$\delta_{s}=34.46$. Per unit system: $S_{B}=220$ MVA,
$U_{B(220)}=209$ kV.

To conform the analytical results above, we use Heun method (see \cite{kloeden})
to simulate the solutions of system \eqref{0} with given initial value and
the parameters. For a certain positive integer $K$,
we have $h=T/K$, $\delta_0=\delta(0)$, $\omega_0=\omega(0)$,
$\Delta B_{i}=B((i+1)h)-B(ih)\sim\sqrt{h}N(0,1)$. The corresponding
discretization equations are
\begin{gather*}
\hat{\delta}=\delta_{i}+h(\omega_{i}-1), \\
\hat{\omega}=\omega_{i}+h\big(-\frac{D}{M}(\omega_{i}-1)
+\frac{P_{\max}}{M}(\sin(\delta_{s})-\sin(\delta_{i}))\big)
+\sigma \Delta B_{i}, \\
\delta_{i+1}=\delta_{i}+0.5h(\omega_{i}+\hat{\omega}), \\
\begin{aligned}
\omega_{i+1}
&=\omega_{i}+0.5h\Big(\big(-\frac{D}{M}(\omega_{i}-1)
 +\frac{P_{\max}}{M}(\sin(\delta_{s})-\sin(\delta_{i}))\big)\\
&\quad +(-\frac{D}{M}(\hat{\omega}-1)
 +\frac{P_{\max}}{M}\big(\sin(\delta_{s})-\sin(\hat{\delta})\big))\Big)
 +\sigma \Delta B_{i}, \quad i=0, \cdot\cdot\cdot,K.
\end{aligned}
\end{gather*}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig2a} % avgy(pf).eps
\includegraphics[width=0.7\textwidth]{fig2b} % y(ybgj).eps
\end{center}
\caption{(top) Stochastic trajectory of
 $\frac{1}{t}\int^{t}_0\omega^2(s)ds$ generated by the Heun scheme
for time step $h=2^{-8}$ for \eqref{0} with $\sigma=0.6$ on $[0,5000]$.
(bottom) Stochastic trajectory of
 $\frac{|\omega(t)|}{\sqrt{\log t}}$ generated by the Heun scheme
for time step $h=2^{-8}$ for \eqref{0} with $\sigma=0.6$ on $[0,5000]$.}
\label{fig2}
\end{figure}

\section{Stochastic fluctuation}

The main aim of this section is to illustrate the fluctuation phenomena.
 At this end, we consider two auxiliary
stochastic differential equations
\begin{equation}\label{3.5}
\begin{gathered}
 {\rm d}y_1=(-\frac{D}{M}(y_1-1)+\frac{2P_{\max}}{M}\sin \delta_{s}{\rm
 d}t+\sigma{\rm d}B(t),\\
 y_1(0)=\omega(0).
 \end{gathered}
\end{equation}
and
\begin{equation}\label{3.6}
\begin{gathered}
 {\rm d}y_2=(-\frac{D}{M}(y_2-1)-\frac{2P_{\max}}{M}\sin \delta_{s}{\rm
 d}t+\sigma{\rm d}B(t),\\
 y_1(0)=\omega(0).
 \end{gathered}
\end{equation}
It\^{o}'s formula implies
\begin{gather*}
\omega(t)=e^{-\frac{D}{M}t}\big\{\omega(0)
 +\int^{t}_0\big[e^{\frac{D}{M}s}\big(\frac{P_{\max}}{M}(\sin \delta_{s}
 -\sin\delta)\big)+\frac{D}{M}\big]ds
 +\int^{t}_0e^{\frac{D}{M}s}\sigma dB(s)\big\},\\
y_1(t)=e^{-\frac{D}{M}t}\big\{\omega(0)
+\int^{t}_0\big[e^{\frac{D}{M}s}\big(\frac{2P_{\max}}{M}\sin \delta_{s}
+\frac{D}{M}\big]ds
+\int^{t}_0e^{\frac{D}{M}s}\sigma dB(s)\big\},\\
y_2(t)=e^{-\frac{D}{M}t}\big\{\omega(0)
+\int^{t}_0e^{\frac{D}{M}s}\big(\frac{D}{M}-\frac{2P_{\max}}{M}\sin \delta_{s} \big)ds
+\int^{t}_0e^{\frac{D}{M}s}\sigma dB(s)\big\}.
\end{gather*}
We therefore have that
\begin{equation}\label{4.5}
y_2(t)\leq \omega(t)\leq y_1(t).
\end{equation}
Next we have a well known lemma
(see Hasminskii \cite[pp. 106-125]{Hasminskii2011}).
 Let $x(t)$ be a homogeneous Markov process in $R^{n}$ described by
the following stochastic differential equation:
\begin{equation}\label{3.4}
 dx(t)=f(x)dt+g(x)dB(t).
\end{equation}
the drift coefficients and diffusion coefficients of system \eqref{3.4}
are $a(x)=f(x),\sigma(x)=g(x)g^{T}(x)$ respectively.

\begin{lemma}[\cite{Hasminskii2011}] \label{lem4.1}
 We assume that there is a bounded open subset $G \subset R$:

(i) in the domain $G$ and some neighborhood therefore, the
smallest diffusion matrix $g^2(x)$ is bounded away
from zero.

(ii) If $x\in R\setminus G$, the mean time $\tau$ at which a
path issuing from $x$ reaches the set $G$ is finite, and $\sup_{x\in
K\setminus G} E_{x}\tau<+\infty$ for every compact subset $K \in R$.and
throughout this paper we set $\inf\emptyset=\infty$.
Then the Markov process $x(t)$ of system\eqref{3.4} is ergodic
and positive recurrent.
\end{lemma}

\begin{lemma} \label{lem4.2}
 The solution process $(y_1(t))$ of system \eqref{3.5} and $(y_2(t))$
of system \eqref{3.6} are ergodic and positive recurrent.
\end{lemma}

\begin{proof}
Without loss of generality, we only provide
the proof for system \eqref{3.5}, the proof for system \eqref{3.6}
is similar. To validate condition (i) and (ii), it suffices to prove
that there exists some neighborhood $U$ and a nonnegative
$C^2$-function $V$ such that $\sigma^2$ is uniformly elliptical
in $U$ and $\mathscr{L}V\leq -1$ for any $x \in R \setminus U$,
(details refer to \cite{mao2011} p. 400). Applying It\^{o}'s
formula to $V(y_1)=y^2_1$ we have
\[
\mathscr{L}V(y_1)=-\frac{2D}{M}y_1^2
+(\frac{4P_{\max}}{M}\sin\delta_{s}+2\frac{D}{M})y_1+\sigma^2.
\]
Note from $\frac{2D}{M}>0,$ that there is a sufficiently large $N$, such that
$$
\mathscr{L}V(y_1)\leq -1,\quad{\rm \forall} |y_1|>N; \quad
\inf_{|y_1|>N}\lambda_{min}(\sigma^2)=\sigma^2>0.
$$
This immediately implies condition (i) and (ii) in Lemma \ref{lem4.1}.
\end{proof}

Now, we use the recurrence of $(y_1(t))$ and $(y_2(t))$ to
explain such fluctuation phenomena. For fixed $\alpha_1$ and
$\alpha_2$ such that $\alpha_1>\alpha_2$, Let $\alpha_1,
\alpha_2$ denote the higher and lower levels respectively. Set
$\tau^{i}_0=\inf\{t\geq0: y_{i}(t)\geq \alpha_1\},
\tau^{i}_1=\inf\{t\geq \tau^{i}_0: y_{i}(t)\leq \alpha_2\},
i=1,2.$ and $\tau_0=\inf\{t\geq0: \omega(t)\geq \alpha_1\},
\tau_1=\inf\{t\geq \tau_0: \omega(t)\leq \alpha_2\}$.
For $k\geq1, i=1,2,\dots$, we
define the following sequence of stopping times recursively
\begin{gather*}
\tau^{i}_{2k}=\inf\{t\geq\tau^{i}_{2k-1}: y_{i}(t)\leq \alpha_2 \},\quad
 \tau^{i}_{2k+1}=\inf\{t\geq\tau^{i}_{2k}: y_{i}(t)\geq \alpha_1 \},\\
\tau_{2k}=\inf\{t\geq\sigma^{i}_{2k-1}: \omega(t)\leq \alpha_2 \},\quad
 \tau^{i}_{2k+1}=\inf\{t\geq\sigma^{i}_{2k}: \omega(t)\geq \alpha_1 \}.
\end{gather*}
For $i=1,2$, the strong Markov property and Lemma \ref{lem4.2} imply
$$
\tau^{i}_{k}<\infty, \quad k\geq 0 \quad a.s.
$$
For each $k>0$, it therefore follows from \eqref{4.5} that
\begin{gather*}
\omega(\tau^{1}_{2k})\leq y_1(\tau^{1}_{2k})
=\alpha_1, \omega(\tau^{1}_{2k+1})\leq y_1(\tau^{1}_{2k+1})=\alpha_1,\\
\omega(\tau^2_{2k})\geq y_2(\tau^2_{2k})
=\alpha_2,\omega(\tau^2_{2k+1})\geq y_2(\tau^2_{2k+1})=\alpha_2.
\end{gather*}
This implies
\begin{align*}
\tau^{1}_{2k}\leq\tau_{2k}\leq \tau^2_{2k}<\infty, \tau^2_{2k+1}
\leq \tau_{2k+1}\leq \tau^{1}_{2k+1}<\infty, k\geq 0.
\end{align*}
This means that the higher and lower levels of $\omega(t)$ occurs,
 which coincides with the recurrent phenomena
in practice.

\section{Discussion on stationary distribution}

The main aim of this section is to discuss the existence of a
stationary distribution of the system \eqref{0}.
Letting $x=[x_1,x_2]=[\delta, \omega-1]$, the stochastic OMIB
system \eqref{0} can be rewritten as the following vector form
\begin{equation}\label{1}
dx=f(x)dt+g(x)dB(t),
\end{equation}
where
\begin{equation}
 f(x)=\begin{bmatrix}
 f_1(x) \\
 f_2(x)
 \end{bmatrix}
= \begin{bmatrix}
 x_2 \\
 \frac{-D}{M}(x_2)+ \frac{P_{\max}}{M}(\sin{x_1}-\sin\delta_{s})\\
 \end{bmatrix},\quad
 g(x)= \begin{bmatrix}
 0 \\
 \sigma
 \end{bmatrix},
\end{equation}
and the drift coefficients and diffusion coefficients of system \eqref{1} are
\begin{equation}
 a(x)= \begin{bmatrix}
 x_2 \\
 \frac{-D}{M}(x_2)+ \frac{P_{\max}}{M}(\sin{x_1}-\sin\delta_{s})
 \end{bmatrix}, \quad
  \sigma(x)=g(x)g^{T}(x)=\begin{bmatrix}
 0 & 0 \\
 0 & \sigma^2
 \end{bmatrix}.
\end{equation}
Let $P(t,x_0, A)$ denote the probability measure
induced by $x(t,x_0)$; that is,
\begin{align*}
P(t,x_0, A)=\mathbb{P}(x(t,x_0)\in A), \quad A\in \mathscr{B}(E^{n}),
\end{align*}
where $\mathscr{B}(E^{n})$ is the $\sigma$-algebra of all the Borel
sets $A \subset E^{n}$. Now we will show that the solution process $x(t)$
has a transition density function $p(t, x_0, y)$.

\begin{lemma} \label{lem5.1}
The transition probability measure
$P(t,x_0, A)$ of system \eqref{1} has a density $p(t,x_0, y)\in
C^{\infty}((0,\infty)\times R^2\times R^2)$.
\end{lemma}

\begin{proof}
Let us now introduce the notation of Lie bracket. If $a(x)$ and $b(x)$ are
vector fields on $R^{d}$,
then the Lie bracket $[a, b]$ is a vector field given by
$$
[a, b]_{j}(x)\sum_{k=1}^{d}(a_{k}\frac{\partial b_{j}}{\partial x_{k}}(x)-b_{k}
\frac{\partial a_{j}}{\partial x_{k}}(x)), j = 1, 2, \cdot\cdot\cdot,d.
$$
For SMIB system, simple computation shows that
$[f,g](x)=[-\sigma, \frac{\sigma D}{M}]^{T}$. Consequently
 \begin{equation}
 |[f,g](x),g(x)|=
\begin{vmatrix}
 -\sigma & 0 \\
 \frac{\sigma D}{M} & \sigma\\
 \end{vmatrix}
=-\sigma^2\neq0,
\end{equation}
which means that $[f, g], g$ are linearly independent on $R^2$.
Thus for every $(x_1, x_2)\in R^2$, vector $[f, g], g$ span
the space $R^2$. In view of Hormander's Theorem in \cite{Bell}, the
transition probability measure $P(t, x_0, A)$ has a density $p(t,
x_0, y)\in C^{\infty}((0,\infty)\times R^2\times R^2)$.
The proof is complete.
\end{proof}

For the stochastic OMIB system \eqref{1}, the density function
$p(t,x_0, y)$ satisfies the following Fokker-Planck equation (FPE):
\begin{equation}\label{3}
\frac{\partial p(t,x_0,y)}{\partial t}=-\frac{\partial
f_1(y)p(t,x_0, y)}{\partial y_1}-\frac{\partial f_2(y)
p(t,x_0,y)}{\partial y_2}+\frac{\sigma^2}{2}\frac{\partial^2
p(t,x_0,y)}{\partial y^2_2}
\end{equation}
The FPE \eqref{3} is a two-dimensional parabolic partial differential
equation, some specified initial conditions and boundary conditions
are provided by Wang in \cite{keyou}. If the process
$x(t,x_0)$ of \eqref{1} has a stationary distribution, then its
probability density function $p_{s}(y)$ satisfies the following
stationary Fokker-Planck equation
\begin{equation}\label{3.5b}
\frac{\partial f_1(y)p_{s}(y)}{\partial y_1}+\frac{\partial f_2(y)
p_{s}(y)}{\partial y_2}=\frac{\sigma^2}{2}\frac{\partial^2
p_{s}(y)}{\partial y^2_2}.
\end{equation}

Now we try to solve the solve the stationary FPE by the detailed balance method,
 which was first applied by van Kampen in \cite{Kampen}.
This method classify the components of the response vector as either odd or
 even variables, according to the transformation from $x_{j}$ to $\bar{x}_{j}$
upon the time reversal $t \to -t$. The even variables do not change their signs
when time is reversed, whereas the odd variables do. More detailed information,
we refer the reader to \cite{Kampen} and \cite{Graham}.


\begin{lemma} \label{lem5.2}
The stationary Fokker-Planck equation
\eqref{3.5} has a solution as following.
\begin{equation}
p_{s}(y_1,y_2)
=C\exp\Big(-\frac{2D}{M\sigma^2}(\frac{y^2_2}{2}
-\frac{P_{\max}}{M}\sin\delta
y_1-\frac{P_{\max}}{M}\cos y_1)\Big)
\end{equation}
\end{lemma}

\begin{proof}
For the second-order differential equations,
Risken\cite{Risken} showed that the first variable is always even
and the second variable is odd. To obtain conditions of detailed balance
we separate each drift
coefficient $a_{j}(y)$ into reversible and irreversible parts.
$a_{j}(y)=a^{I}_{j}(y) + a^{R}_{j}(y)$, where
$$
a^{I}_{j}(y)=\frac{1}{2}(a_{j}(y)+\Delta_{j} a_{j}(\bar{y})),
a^{R}_{j}(y)=\frac{1}{2}(a_{j}(y)-\Delta_{j} a_{j}(\bar{y})),\quad j=1,2.
$$
where $\bar{y}_{j}=\Delta_{j}y_{j}$, $\Delta_{j}=\pm1$ corresponds
to even $(+)$ or odd $(-)$ variables. Since $p_{s}(y_1,y_2)$ is nonnegative,
it can be expressed in the form
$$
p_{s}(y_1,y_2)=C\exp(-\varphi(y_1, y_2)).
$$
Then according to \cite{Graham} equations to determine
$\varphi$ are
\begin{equation}\label{7}
\begin{gathered}
\sum_{j}2a_{j}^{I}(y)-\sum_{j,k}\frac{\partial
\sigma_{j,k}}{\partial y_{k}}+\sum_{j,k}\sigma_{j,k}\frac{\partial
\varphi} {\partial y_{k}}=0,\\
\sum_{j}\frac{\partial a_{j}^{R}(y)}{\partial
y_{j}}+\sum_{j}a_{j}^{R}(y)\frac{\partial
\varphi} {\partial y_{k}}=0,\\
\sum_{j,k}[\sigma_{jk}(y)-\Delta_{j}\Delta_{k}\sigma_{jk}(\bar{y})]=0.
\end{gathered}
\end{equation}
It is easy to show that the reversible and
irreversible parts of each drift coefficient of system \eqref{1} are
\[
a_1^{I}(y)=0, a_1^{R}(y)=y_2, \quad
a_2^{I}(y)=-\frac{D}{M}y_2, \quad
a_2^{R}(y)=\frac{P_{\max}}{M}(\sin \delta-\sin y_1)
\]
And the corresponding equation \eqref{7} of stochastic OMIB \eqref{1} is given
by
\begin{equation}
\frac{\partial\varphi}{\partial y_2}=\frac{2D}{M\sigma^2}y_2;
\frac{\partial\varphi}{\partial
y_1}=\frac{2DP_{\max}}{M^2\sigma^2}(\sin y_1-\sin \delta).
\end{equation}
Simple computations show that
\begin{equation}
p_{s}(y_1,y_2)
=C\exp\big(-\frac{2D}{M\sigma^2}(\frac{y^2_2}{2}
-\frac{P_{\max}}{M}\sin\delta y_1-\frac{P_{\max}}{M}\cos y_1)\big).
\end{equation}
where $C$ is a constant. The proof is complete.
\end{proof}

Obviously, the solution $p_{s}(y_1,y_2)$ is periodic in $y_1$,
which implies that the integral $\int_{R^2}p_{s}(y_1,y_2)dy_1dy_2$
can not be convergent. As a result, the solution can not be viewed as a
 density function on $R^2$. On the other hand, the periodic might
provide another perspective on this topic. Now we state a lemma
 to indicate the periodic of the solution process.


\begin{lemma} \label{lem5.3}
 Let $x(t,x_0)=(x_1(t,x_0), x_2(t,x_0))$ be the global solution of system
\eqref{1} with initial data $x_0=(x_{01}, x_{02})$. Setting
$x'_0=(x_{01}+2\pi, x_{02})$, for any $t>0$, we can
claim that
$$
x_1(t,x'_0)-x_1(t,x_0)=2\pi, \quad
x_2(t,x'_0)-x_2(t,x_0)=0.
$$
\end{lemma}

\begin{proof}
It follows from the periodicity of $f(x)$ with $x_1$ that the process
$(x_1(t,x_0)+2\pi, x_2(t,x_0))$ is still the solution of system \eqref{1}.
\begin{gather} \label{5.1}
\begin{gathered}
x_1(t,x'_0)=x'_{01}+\int^{t}_0x_2(s,x'_0)ds,   \\
\begin{aligned}
x_2(t,x'_0) & =x'_{02}+\int^{t}_0\Big(\frac{-D}{M}x_2(s,x'_0)+
\frac{P_{\max}}{M}(\sin{x_1(s,x'_0)}-\sin\delta)\Big)ds \\
&\quad+\int^{t}_0\sigma dB(s),
\end{aligned}
\end{gathered}\\
 \label{5.2}
\begin{gathered}
x_1(t,x_0)+2\pi=x_{01}+2\pi+\int^{t}_0x_2(s,x_0)ds, \\
\begin{aligned}
x_2(t,x_0)
&=x_{02}+\int^{t}_0\Big(\frac{-D}{M}x_2(s,x_0)+
\frac{P_{\max}}{M}(\sin(x_1(s,x_0)+2\pi)\\
&\quad -\sin\delta)\Big)ds
 +\int^{t}_0\sigma dB(s).
\end{aligned}
\end{gathered}
\end{gather}
Subtracting \eqref{5.1} from \eqref{5.2} yields
\begin{gather*}
x_1(t,x'_0)-x_1(t,x_0)-2\pi=\int^{t}_0(x_2(s,x'_0)-x_2(s,x_0))ds, \\
\begin{aligned}
x_2(t,x'_0)-x_2(t,x_0)
&=\int^{t}_0\Big(
\frac{-D}{M}(x_2(s,x'_0)-x_2(s,x_0))+\frac{P_{\max}}{M}(\sin(x_1(s,x'_0))\\
&\quad-\sin\big(x_1(s,x_0)+2\pi)\big)\Big)ds.
\end{aligned}
\end{gather*}
Simple computations show that
\begin{align*}
&(x_1(t,x'_0)-x_1(t,x_0)-2\pi)^2+(x_2(t,x'_0)-x_2(t,x_0))^2 \\
&\leq \int^{t}_0\big(\frac{P^2_{\max}}{M^2}+\frac{D}{M}+1\big)
 \Big((x_1(s,x'_0)-x_1(s,x_0)-2\pi)^2+(x_2(s,x'_0) \\
&\quad -x_2(s,x_0))^2\Big)ds.
\end{align*}
The well-known Gronwall inequality yields
$$
x_1(t,x'_0)\equiv x_1(t,x_0)+2\pi, \quad x_2(t,x'_0)\equiv x_2(t,x_0).
$$
Therefore we get the desired assertion.



By Lemma \ref{lem5.3}, the state space can be viewed either in
planar $R^2$ or cylindrical $S^{1}\times R$. The cylindrical
space is a more natural space from the physical perspective. And the
process $x(t,x_0)=(x_1(t,x_0), x_2(t,x_0))$ on planar
$R^2$ can be mapped as a process
$\tilde{x}(t,x_0)=(\tilde{x}_1(t,x_0),
\tilde{x}_2(t,x_0))$ on cylindrical space $S^{1}\times R$, where
\begin{equation}\label{3.3}
\tilde{x}_1(t,x_0)=x_1(t,x_0) \quad\pmod{2\pi}, \quad
\tilde{x}_2(t,x_0)=x_2(t,x_0).
\end{equation}
The corresponding stationary Fokker-Planck equation has the form
\begin{equation}\label{3.9}
\frac{\partial f_1(y)\tilde{p}_{s}(\tilde{y})}{\partial
\tilde{y}_1}+\frac{\partial f_2(\tilde{y})
\tilde{p}_{s}(\tilde{y})}{\partial
\tilde{y}_2}=\frac{\sigma^2}{2}\frac{\partial^2\tilde{
p}_{s}(\tilde{y})}{\partial \tilde{y}^2_2}.
\end{equation}
In this case, we can choose a constant $C_1$ such that
\begin{equation}
C_1\int_{S^{1}\times R}\exp\Big(-\frac{2D}{M\sigma^2}
(\frac{\tilde{y}^2_2}{2}-\frac{P_{\max}}{M}\sin\delta_{s}
\tilde{y}_1-\frac{P_{\max}}{M}\cos \tilde{y}_1)\Big)d\tilde{y}_1d\tilde{y}_2=1.
\end{equation}
That is to say, $\tilde{p}_{s}(\tilde{y})$ can be seen as a density function on
$S^{1}\times R$. Then a interesting question arise naturally:
Is there a stationary distribution to system \eqref{1} on the cylinder?
To illustrate this topic, let us present a lemma which is essential
to the proof.
\end{proof}

\begin{lemma}[\cite{soize}] \label{lem5.4}
 We assume that drift vector $a(x)\in R^{n}$ and diffusion matrix
$\sigma(x)$ of the diffusion process $X$ are continuous on $E^{n}$
and independent of time $t$, and
\begin{itemize}
\item[(i)] The diffusion process $X$ has a transition density
function $p(t,x,y)$.

\item[(ii)] For all $x \in E^{n}$, $j,k \in{1,2,\cdot\cdot\cdot, n}$,
the first order partial derivatives of
$p(t,x,y)$ with respect to $t$; the first order partial derivatives
of $b_{j}(y)p(t,x,y)$ with respect to $y_{j}$, the second order
partial derivatives of $\sigma_{jk}(y(y)p(t,x,y)$ with respect to
$y_{j}$ and $y_{k}$, exist and are all continuously
differentiable.

\item[(iii)] There exists a function $p_{s}(y)$ from $R^{n}$ into
$R$ satisfying the steady state Fokker-Planck equation on $R^{n}$
which can be written for all $x$ in $E^{n}$ as
\begin{equation}\label{6}
-\sum_{i=1}^{n}\frac{\partial a_{i}(y)p_{s}(y)}{\partial
y_{i}}+\sum_{i,j}^{n}\frac{1}{2}\frac{\partial^2b_{ij}(y)
p_{s}(y)}{\partial y_{i}\partial y_{j}}=0.
\end{equation}
with the positivity and normalization condition
$\int_{E^{n}}p_{s}(y)dy=1$.
\end{itemize}
Then the stochastic process $X(t)$ has an invariant measure with
$p_{s}(y)$ as its probability density function.
\end{lemma}


Obviously, the smoothness of $p(t,x_0,y)$ can be guaranteed by
 Lemma \ref{lem5.1}, and probability density function can be obtained by detailed
balance technique in Lemma \ref{lem5.2}. Making use of the property of periodicity,
we can claim that the corresponding probability density function
$\tilde{p}(t, x_0, \tilde{y})$ of process
$\tilde{x}(t,x_0)$ is also sufficiently smooth on cylindrical
space. and the stationary Fokker-Planck equation
of process $\tilde{x}(t,x_0)$ has the following solution on
cylindrical space
\begin{equation}
\tilde{p}_{s}(\tilde{y}_1,\tilde{y}_2)
=C_1\exp\Big(-\frac{2D}{M\sigma^2}(\frac{\tilde{y}^2_2}{2}
-\frac{P_{\max}}{M}\sin\delta_{s}
\tilde{y}_1-\frac{P_{\max}}{M}\cos \tilde{y}_1)\Big).
\end{equation}
where $C_1$ is the normalizing constant. Hence, a simple application of the
 Lemma \ref{lem5.4} implies the main result of this subsection.

\begin{theorem} \label{thm5.1}
For $x_0=(\delta_{s},1)\in S^{1}\times R$, the process
$\tilde{x}(t,x_0)$ on $S^{1}\times R$ has a stationary
distribution.
\end{theorem}

\begin{remark} \label{rmk1} \rm
Denote by $\mu(\cdot)$ the stationary distribution of the process
$\tilde{x}(t,x_0)$ and
$(\tilde{y}_1,\tilde{y}_2)$ be the random variable to which
$\tilde{x}(t,x_0)$ converges in distribution.
The proof of Theorem \ref{thm3.1} implies that the probability density
function of $\tilde{y}_1$ is
\begin{equation}\label{3.2}
p_1(\tilde{y}_1)=C_2\exp\Big(\frac{2DP_{\max}}{M^2\sigma^2}(\sin\delta
\tilde{y}_1+\cos \tilde{y}_1)\Big), \quad
\tilde{y}_1\in[-\pi,\pi),
\end{equation}
where $C_2$ is the normalizing constant. The distribution of
$\tilde{y}_2$ is just the normal distribution with its probability
density function
\begin{align*}
p_2(\tilde{y}_2)=\frac{1}{\sqrt{\frac{M\Pi}{D}}\sigma^2}
\exp\Big(-\frac{2D}{M\sigma^2}(\frac{\tilde{y}^2_2}{2}\big)\Big),
\quad -\infty<\tilde{y}_2<\infty.
\end{align*}
\end{remark}

\subsection*{Conclusion}
In this paper, we have investigated the asymptotic behavior of the stochastic
OMIB systems. Firstly, utilizing stochastic analysis techniques, the
 asymptotic bound properties of $p$th moment and asymptotic pathwise
estimation of the stochastic OMIB systems have been researched. Secondly,
by ergodic property and the strong Markov property, higher and lower
rotor speed levels appear infinitely, which may give a good explanation
of the fluctuation phenomena. Finally, the existence of stationary
distribution on cylinder has been derived by periodicity and some analysis
analysis techniques.

\subsection*{Acknowledgements}
This research was  supported by the National
Science Foundation of China (Grant Nos. 61304070, 61773152, 51137002, 51190102),
by the National Key Basic Research Program of China (973 Program)
(2013CB228204), and by the Fundamental Research Funds for the Central
Universities of China. (Grant No. 2015B19814).


\begin{thebibliography}{99}


\bibitem{Bell} D. R. Bell;
 \emph{The Malliavin Calculus}, Dover Publications, New York, 2006.

\bibitem{chen2005} X. Chen, W. Zhang, W. Zhang;
 Chaotic and subharmonic oscillations of a nonlinear power system,
\emph{IEEE Trans. Circuits Syst. II}, 52(12) (2005), 811-815.

\bibitem{zhu2010} L. Chen, W. Zhu;
 First passage failure of dynamical power systems under
random perturbations, \emph{Science China Technological Sciences,}
53 (9)(2000), 2495-2500.

\bibitem{Graham} R. Graham, H. Hakcn;
 Generalized thermodynamic potential for Markoff systems in detailed balance
and far from thermal equilibrium,
\emph{Zeitschrift f¨šr Physik A Hadrons and nuclei}, 243(3) (1971), 289-302.

\bibitem{ju1996b} P. Ju, E. Handschin, D. Karlsson;
 Nonlinear dynamic load modelling: model and parameter estimation,
\emph{IEEE Trans. on Power Systems}, 11(4) (1996), 1689-1697.

\bibitem{ju1996a} P. Ju, E. Handschin, F. Reyer;
Genetic algorithm aided controller design with application to SVC,
\emph{IEE Proc.-Gen. Transm. Distr}. 143(3) (1996), 258-262.

\bibitem{Hasminskii2011}  R. Hasminskii;
 \emph{Stochastic Stability of Differential Equations},
Springer-Verlag, Berlin, Heidelberg, 2011.

\bibitem{Hiskens2006} I. Hiskens, J. Alseddiqui;
 Sensitivity, approximation, and uncertainty
in power system dynamic simulation, \emph{IEEE Trans. Power Syst.},
21(4) (2006), 1808-1820.

\bibitem{Hockenberr2004} J. Hockenberr, B. Lesieutre;
Evaluation of uncertainty in dynamic
simulations of power system models: The probabilistic collocation
method, \emph{IEEE Trans. Power Syst.}, 19 (3) (2004), 1483-1491.

\bibitem{ju1991} P. Ju, G. Wu, Y. Li;
 Fundamental theorems on probabilistic stability
of power system (in Chinese), \emph{Proc CSEE}, 11(6)(1991), 17-25.

\bibitem{kloeden} P. Kloeden, C. Platen;
\emph{Numerical Solution of Stochastic Differential Equations}.
Berlin etc., Springer-Verlag, 1992.

\bibitem{kundur} P. Kundur;
\emph{Power system and control}, New York, McGraw-Hill, 1994.

\bibitem{liu2011} F. Liu, W. Wei, S. Mei;
On expansion of estimated stability region:
Theory, methodology, and application to power systems,
\emph{Sci China Tech Sci}, 54 (6) (2011), 1394-1406.

\bibitem{lu2015} Z. Lu, W. Zhao, D. Xie, G. Li;
 P-moment stability of power system under small Gauss type random excitation,
\emph{Chaos, Solitons \& Fractals}, 81(2005), 30-37.

\bibitem{10} X. Mao;
 \emph{Stochastic Differential Equations and Applications}, Horwood,
Chichester, 1997.

\bibitem{mao2011} X. Mao;
 Stationary distribution of stochastic population systems,
\emph{Syst. \& Contr. Lett.}, 60 (2011), 398-405.

\bibitem{pan2013} X. Pan, P. Ju, F. Wu, Y. Jin;
 Parameter estimation of drive system in a fixed-speed wind turbine by
 utilizing turbulence excitations,
\emph{IET Generation, Transmission \& Distribution,} 7 (7) (2013), 665-673.

\bibitem{Risken} H. Risken;
 \emph{The Fokkcr-Planck Equation. Methods of Solution and Applications},
Springer, Berlin, 1984.

\bibitem{soize} C. Soize;
\emph{The Fokker-Planck equation for stochastic dynamical systems and Its
explicit Steady Solutions}, World Scientific, London, 1994.

\bibitem{Stankovic} A. Stankovic, B. Lesieutre;
 A probabilistic approach to aggregate induction machine modeling,
\emph{IEEE Trans. Power Syst.,} 11(4)(1996) 1983-1989.

\bibitem{Timko1983} K. Timko, A. Bose, P. M. Anderson;
 Monte Carlo simulation of power system stability,
\emph{IEEE Trans. Power App. Syst.}, PAS-102 (10) (1983) 3453-3459.

\bibitem{Kampen} N. G. Van Kampen;
 Derivation of the phenomenological equations from the master equation: II.
Even and odd variables, \emph{Physical}. 23(6¨C10) (1957), 816-824.

\bibitem{Vlachogiannis} J. Vlachogiannis;
 Probabilistic Constrained Load Flow Considering Integration of Wind
Power Generation and Electric Vehicles,
\emph{IEEE Trans. Power Syst.,} 24(4) (2009), 1808-1817.

\bibitem{wangkw2000} K. Wang, C. Chung, C. Tse, et al.;
 Improved probabilistic method for power system dynamic stability studies,
\emph{IEE P-Gener Transm D}, 147 (1) (2000), 37-43.

\bibitem{keyou} K. Wang, L. Mariesa;
 The Fokker-Planck equation for power system stability probability density
function evolution, \emph{IEEE Trans. Power Syst.}, 28(3)(2013), 2994-3001.

\bibitem{wu2013} F. Wu, P. Ju, X.P. Zhang, H. Huang, C. Qin, J. Fang;
 Modeling, control strategy, and power conditioning for direct-drive wave
 energy conversion to operate with power grid,
\emph{Proceedings of the IEEE}, 101 (4) (2013), 925-941.

\bibitem{Wu1983} F. Wu, Y. Tsai;
 Probabilistic dynamic security assessment of power systems:
Part I, Basic model;
 \emph{IEEE Trans. Circuits Syst.}, 30 (3) (1983), 148-159.

\bibitem{wu2007} F. Wu, X. P. Zhang, K. Godfrey, P. Ju;
 Small signal stability analysis and optimal control of a wind turbine
with doubly fed induction generator,
 \emph{IET Generation, Transmission\& Distribution}, 1(5) (2007), 751-760.

\bibitem{wu2008} F. Wu, X. P. Zhang, P. Ju, M. J. H. Sterling;
Decentralized nonlinear control of wind turbine with doubly fed induction
generator, 	\emph{IEEE Trans. on Power Systems}	, 23 (2) (2008), 613-621.

\bibitem{zhang2012} J. Zhang, P. Ju, Y. Yu, F. Wu;
 Responses and stability of power system under small Gauss type random excitation,
\emph{Science China Technological Sciences}, 55 (7) (2012), 1873-1880.

\end{thebibliography}

\end{document}