\documentclass{sqa_paper}
%\usepackage{txfonts}
\usepackage{apalike}
\usepackage{times}%Wordlike fonts
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{latexsym}
\usepackage{amsmath}
\usepackage[mathscr]{eucal}
\usepackage{texdraw}
\usepackage{rotating}
\usepackage{epsfig}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\hoffset=5mm
\renewcommand{\thefootnote}{\fnsymbol{footnote}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\switchbox #1{\bsegment
\textref h:C v:C
\vtext (0 0){{$\curvearrowleft$}} \move (0.06 0)
\bulbox{$\bullet$}
\move (-0.02 -0.27)
\bulbox{$t_{\mathrm{c}}$}
\move (0.12 0.01)
\ravec(1 0) \move (-0.1 0.4) \forbox{#1}
\move(0 0)
\rmove (-0.25 +0.125) \bulbox{$\bullet$}
\rmove (0 -0.25) \bulbox{$\bullet$} \rmove(+0.02 +0.21)
\lvec(0.1 0)
    \savepos (0.62 0)(*ex *ey)
\esegment \move (*ex *ey)}

\def\sysbox #1{\bsegment
\move(0 0)
\ravec(+1.25 0) \lbox{$\mathscr{S}_{#1}$} \rlvec(+1.25 0)
\rmove(-2.24 +0.75) \ravec(0 -0.44) \rmove(+0.75 +0.44) \ravec(0 -0.44)
\rmove(-1.55 +0.5) \forbox{$w_{#1}(t_{i})$}
\rmove(+0.85 0) \forbox{$v_{#1}(t_{i})$}
\savepos (3.75 0)(*ex *ey)
\esegment \move (*ex *ey)}

\def\forbox #1{\bsegment
\textref h:C v:C \htext (0.375 0){#1} \savepos (0.75 0)(*ex *ey)
\esegment \move (*ex *ey)}

\def\tbox #1{\bsegment
\lvec (0 +0.25) \lvec (0.75 +0.25) \lvec (0.75 -0.25) \lvec (0 -0.25)
\lvec (0 0)
\textref h:C v:C \htext (0.375 0){#1} \savepos (0.75 0)(*ex *ey)
\esegment \move (*ex *ey)}

\def\lbox #1{\bsegment
\lvec (0 +0.3) \lvec (1.25 +0.3) \lvec (1.25 -0.3) \lvec (0 -0.3)
\lvec (0 0)
\textref h:C v:C \htext (0.625 0){#1} \savepos (1.25 0)(*ex *ey)
\esegment \move (*ex *ey)}

\def\bbox #1{\bsegment
\lvec (0 +0.5) \lvec (-1.25 0.5) \lvec (-1.25 -0.5) \lvec (0 -0.5)
\lvec (0 0)
\textref h:C v:C \htext (-0.625 0){#1} \savepos (-1.25 0)(*ex *ey)
\esegment \move (*ex *ey)}

\def\addbox #1{\bsegment
\move(0.15 0)  \lcir r:0.15
\textref h:C v:C \htext (0.15 0){#1} \savepos (0.3 0)(*ex *ey)
\esegment \move (*ex *ey)}

\def\multbox #1{\bsegment
\lvec (0 +0.15) \lvec (0.3 +0.15) \lvec (0.3 -0.15) \lvec (0 -0.15)
\lvec (0 0)
\textref h:C v:C \htext (0.155 0){#1} \savepos (0.3 0)(*ex *ey)
\esegment \move (*ex *ey)}

\def\bulbox #1{\bsegment
\textref h:C v:C \htext (0 0){#1} \savepos (0 0)(*ex *ey)
\esegment \move (*ex *ey)}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Local definitions added by Semoushin
\newcommand{\Rset}{{\mathbb{R}}}
\newcommand{\Zset}{{\mathbb{Z}}}
\newcommand{\Nset}{{\mathbb{N}}}
\newcommand{\Hyp}{{\mathbb{\bf H}}}
%\DeclareMathOperator{\col}{col}
\newcommand{\col}[1]{\mathop{\rm col}[#1]}
\newcommand{\half}{\frac{\scriptstyle 1}{\scriptstyle 2}}
% Remark environment
\newtheorem{remarka}{Remark}%[section]
% Theorem environment
\newtheorem{theo}{Theorem}[section]
% Corollary environment
\newtheorem{coro}{Corollary}%[section]
\newtheorem{alg}{Algorithm}%[section]
\newtheorem{defn}{Definition}%[section]
\newtheorem{lem}{Lemma}[section]
%\newtheorem{ex}{Example}[section]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%завели новый счетчик для ПРИМЕРОВ и установили формат записи ПРИМЕРОВ %%
\newcounter{exc}[section]                                              %%
\newcommand{\ex}% продолжение команды                                  %%
{\par\medskip\noindent\refstepcounter{exc}{\textbf{Example}%%%%%%%%%%%%%%
% продолжение команды:%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\, \textbf{\thesection.\arabic{exc}.}}\; \ }                                    %%
\renewcommand{\theexc}{\thesection.\arabic{exc}}                       %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%завели новый счетчик для REMARKS и установили формат записи ПРИМЕРОВ %%
\newcounter{rmc}[section]                                              %%
\newcommand{\rem}% продолжение команды                                  %%
{\par\medskip\noindent\refstepcounter{rmc}{\textbf{Remark}%%%%%%%%%%%%%%
% продолжение команды:%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\, \textbf{\thesection.\arabic{rmc}.}}\; \ }                                    %%
\renewcommand{\thermc}{\thesection.\arabic{rmc}}                       %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%новая команда типа \subparagraph. Пишите for example так: \pr{Proof}
\newcommand{\pr}[1]{\paragraph{\it #1.}}
%новая команда quod erat demostrandum
\newcommand{\qed}[1]{\hbox to \textwidth{\hfill  #1}}

\newcommand{\keyw}[1]{{\bf #1}}
\newcommand{\wb}{\widebar}
\newcommand{\wt}{\widetilde}
\newcommand{\wh}{\widehat}
\newcommand{\ds}{\displaystyle}
\newcommand{\E}[1]{\mathop{{\rm \bf E}\left\{#1\right\}}\nolimits}
\newcommand{\rank}{\mathop{\rm rank}\nolimits}
\newcommand{\tr}[1]{\mathop{{\rm tr}\{#1\}}\nolimits}
\newcommand{\ud}{\mathrm{d}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Local definitions added by Semoushin
\newcommand{\Nu}{{\rm N}}
\newcommand{\Eta}{{\rm H}}
\newcommand{\diag}[2]{\mathop{\rm diag}[#1, \; #2]}
%%%%%%%%%%%%%%%%%% Gothic fonts %%%%%%%%%%%%%%%%%%%
% Local definitions added by Semoushin            %
%12pt                                             %
\font\got=eufm10  % Text                          %
\font\gots=eufm7  % Script                        %
\font\gotss=eufm5 % Subscript                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\numberwithin{equation}{section}
%%%%%%%%%%%%%%%%%% Gothic letters %%%%%%%%%%%%%%%%%%%%%%%
% Local definitions added by Semoushin
\newcommand{\g}[1]{\mbox{\got #1}}     % Text       %
\newcommand{\gs}[1]{\mbox{\gots #1}}   % Script     %
\newcommand{\gss}[1]{\mbox{\gotss #1}} % Subscript  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%новая команда для центрирования любого текста. Пишите так: \h{Текст}
\newcommand{\h}[1]{\vspace{2pt}\begin{center}\Large \bf #1\end{center}}
\newcommand{\hs}[1]{\vspace{2pt}\begin{center}\large \bf #1\end{center}}
%новая команда вместо \paragraph. Пишите так: \p{Текст}
\newcommand{\p}[1]{\paragraph{\bf #1}}
%новая команда вместо \subparagraph. Пишите так: \s{Текст}
\newcommand{\s}[1]{\subparagraph{\bf #1}}
%новое окружение для заключения текста в рамку шириной, задаваемой в аргументе      %
%пишите так: \begin{myframe}{6cm} Текст Текст Текст \end{myframe}                   %
\newenvironment{myframe}[1]{\begin{tabular}{|p{#1}|}\hline}{\\\hline\end{tabular}}  %
%
\title{A Simple Squared Innovations Based Decision Generator for
Detection/Selection Problems in Linear Stochastic Control Systems
}

\author{I. V. Semoushin\footnote{Address correspondence to I. V. Semoushin, Department of Mathematical Cybernetics and Informatics, Ulyanovsk State University, 42 L. Tolstoy Street, Ulyanovsk 432970, Russia; Fax: +7-842-241-2340; E-mail: \underline{SemushinIV@ulsu.ru}\,.}} 
\address{
Department of Mathematical Cybernetics and Informatics\\
%% \\42 L Tolstoy Street, 
Ulyanovsk State University, Ulyanovsk, Russia
%%\\
%%e-mail: \underline{SemushinIV@ulsu.ru}, web page:
%%\underline{http://staff.ulsu.ru/semoushin} 
} 

\abstract{
We consider the linear time-invariant state-space stochastic
control system
\begin{subequations}
\renewcommand{\theequation}{\roman{equation}}
\begin{eqnarray}
x(t_{i+1})&=&\Phi_{\theta} x(t_i)+\Psi_{\theta} u(t_i)+ w(t_i), \quad
x\in \Rset^n \label{a1} \\
y(t_i) &=&
H_{\theta} x(t_i) + v(t_i), \quad y\in \Rset^m  \label{a2}\\
\hat x_{0}(t_{i+1}^-)&=&{\Phi_{0}}
\hat x_{0}(t_i^+)+ {\Psi_{0}} u(t_i),\quad \hat x_{0}\in \Rset^n
 \label{a3}\\
\hat x_{0}(t_i^+) &=& \hat x_{0}(t_i^-)+
{K_{0}}\nu (t_i), \quad \nu (t_i)=y(t_i) - {H_{0}}\hat x_{0}(t_i^-)
 \label{a4} \\
u(t_i) &=& f_{\mathrm{R}}[\hat x_{0}(t_i^+)] = - {G_{0}^{\star}}
\hat x_{0}(t_i^+), \quad u\in \Rset^q  \label{a5}
\end{eqnarray}
\end{subequations}
where  $i\in \Zset$, with a plant (\ref{a1}) and a sensor
(\ref{a2}), both parameterized by an uncertainty vector
parameter $\theta$, and a feedback controller FC,
(\ref{a3})--(\ref{a5}), designed to cascade the
steady-state Kalman filter (\ref{a3})--(\ref{a4}) with a regulator 
(\ref{a5}) for a
nominal value $\theta_{0}$ of \ $\theta$. In
(\ref{a1})--(\ref{a5}), the initial state $x(t_{-s})$ is
given at some $t_{-s}\in \Rset$ with expectation $\E{
\|x(t_{-s})\|^2 }<\infty$; $w(t_i)$ and $v(t_i)$ are zero-mean
mutually orthogonal wide-sense stationary orthogonal sequences
with $\E{w(t_i)w(t_i)^T}=Q_{\theta} \ge 0$ and
$\E{v(t_i)v(t_i)^T}=R_{\theta}>0$ for all $t_i\in \Rset$, with
$\left[w(t_i)\atop v(t_{i})\right]$ orthogonal to $x(t_j)$ and
$u(t_j)$ for all $j\ge i$\,; $u(t_i)$  is wide-sense stationary and
$\E{\|u(t_i)\|^2}<\infty$ for all $t_i\in \Rset$.

Assuming  that the true value $\theta^{\dagger}$ of \ $\theta$ can
change abruptly from $\theta_{0}$ to some $\theta_{1}$ at an
instant $t_{\mathrm{c}} \in (t_{0}, t_{i})$, we detect
$t_{\mathrm{c}}$ using the nominal covariance $C_{0}$ of $\nu
(t_i)$ in (\ref{a4}) and constructing a decision function
of the cumulative sum form:
\begin{displaymath}
S_{k}=\sqrt{m/(2k)}\sum_{i=1}^{k}\left[\nu^T(t_i)
C_{0}^{-1}\nu (t_i)/m-1\right]\ .
\end{displaymath}
}
\keywords{Empirical decision procedures; Hypothesis testing; Diagnostics;
Stochastic systems and control; Stopping times.} 
\subjectclass{62C12; 62F03; 62J20; 93E10; 60G40.} 
\begin{document}
\setcounter{equation}{0}

\section{\hspace{-2ex}{.}\hspace{1ex}INTRODUCTION}\label{intro}

\noindent
The theory  of abrupt changes is an intensively investigated topic
in time series analysis and identification~\cite{basseville:1993}. Usually
by abrupt changes, it is meant changes in characteristics that
occur very fast with respect to the sampling period of
measurements, if not instantaneously, at unknown time instants.
Besides strong theoretical motivations, the detection of such
changes is a problem of great practical interest.

The subject of abrupt changes basically grew up at the confluence of 
several disciplines, first of all, mathematical statistics and automatic
control theory. In  mathematical statistics, decision tools for detecting 
changes appeared first in the area of quality control, where the Shewhart 
3-sigma control charts were introduced \cite{shewhart:1931}. As an 
alternative to them, Page set up a principle of cumulative sum control 
charts \cite{page:1954}. In succeeding years the CUSUM has been proven to 
be a very effective tool to detect changes in a mean level of a process. 
The expanded CUSUM mathematical theory and its industrial applications 
were given in \cite{johnson:1962}, \cite{johnson:1963}, 
\cite{johnson:1964} and then in \cite{siegmund:1985}, \cite{kenett:1998} 
and \cite{ghosh:1991}. Wide historical notes and references for seminars, 
survey papers, and books related to change detection both in mathematical 
statistics and automatic control theory, can be found 
in~\cite{basseville:1993}.

In~\cite{basseville:1993}, many applications of decision tools for 
detecting changes have been summarized, among them quality control, 
automatic segmentation of signals, fault detection and monitoring in 
industrial plants and navigation systems, seismic data processing, and 
some others. One of possible applications of change detection tools that 
has recently emerged may be in the Replication-Selection Systems (RSS)
\cite{murgu:2002}. 
Genetic algorithms (GAs) used in the RSS require a fitness function that 
informs the GA of what is the best performer within a population of 
performers. Selection of the best performer is very close to 
discriminating between several contesting hypotheses, which is the core of 
any detection problem. By virtue of this fact, it is quite natural to 
search for a suitable fitness function for RSS amongst the well developed 
detection methods.

One of possible  approaches to detect changes that is frequently
used in the engineering literature, is based upon the intuitive
idea of detecting a nonadditive change by monitoring the
innovation sequence (IS) $\nu(t_{i})$ of an autoregression (AR)
model with the aid of a test for its variance $\sigma^{2}$. The
CUSUM
\begin{equation}\label{Sk}
S_{k}=\frac{1}{\sqrt{k}}\sum_{i=1}^{k}s_{i} \quad \mbox{with} \quad s_{i}=\frac{1}{\sqrt{2}}
\left( \frac{\nu(t_{i})^{2}}{\sigma_{0}^{2}} -1\right)
\end{equation}
is asymptotically,  when $k$ goes to infinity, distributed as the
standard Gaussian law:  $S_{k} \sim \mathcal N(0; 1)$ under the
hypothesis $\Hyp_{0}: \; \sigma^{2}=\sigma_{0}^{2}$.

This algorithm,  which is monitoring $s_{i}$ termed the weighted
squared innovations (WSI), was introduced in the AR case
independently in~\cite{jones:1970}, \cite{borodkin:1976},
and~\cite{segen:1980} for the purpose of automatic segmentation of
electroencephalogram signals. This test assumes that the IS is
white, and thus a test for whiteness should be done
first~\cite{mehra:1971}. It is well known that the IS is white within a
filter if the filter is Kalman, i.~e., optimal one~\cite{maybeck:1979}, 
and whiteness of the IS is only necessary, not sufficient condition
for filter optimality~\cite{boozer:1971}. Autocorrelation properties of
the IS after change may be very complex, depending on the change
occurred. Because of this, monitoring squared innovations has received a 
criticism as having its poor behavior in many circumstances.
The criticism, as well as the ease of implementation of this method,
stems from the fact that there no information is used about the after 
change characteristics.

In the present paper we advocate the WSI method by showing that 
monitoring WSI will, in certain RSS systems, serve as a useful tool
due to its tolerable performance and moderately simple computational cost, 
which promises the overall cost reduction when the size of population is 
large.

The paper formulates  and develops the WSI method as applied to
stochastic state estimation and control systems (Sections~2 and
3). In Section~4, the WSI test is shown to be a useful tool of GAs
for the best feedback selection in view of improved performances
and condition-based maintenance. Section~5 of the paper presents
some experimental results of fault detection for an example
relating to inertial navigation system modeling. Section~5
concludes the main content and Appendix completes the paper.

\section{\hspace{-2ex}{.}\hspace{1ex}PROBLEM FORMULATION}\label{problem}

\noindent
Let the closed-loop control system be parameterized by a vector
$\theta \in {\Rset}^{p}$ and the available data~$z=\left[y\atop
u\right]$ be a vector composed  of two parts: the control input $u
\in {\Rset}^q$ and the  measurement output $y \in {\Rset}^m$. The
system is modeled for $i\in \Zset$ \ by  the equations
\begin{eqnarray}
x(t_{i+1})&=&\Phi_{\theta} x(t_i)+\Psi_{\theta} u(t_i)+ w(t_i), \quad
x\in \Rset^n \label{x} \\
y(t_i) &=&
H_{\theta} x(t_i) + v(t_i), \quad y\in \Rset^m  \label{y}\\
\hat x_{0}(t_{i+1}^-)&=&{\Phi_{0}}
\hat x_{0}(t_i^+)+ {\Psi_{0}} u(t_i),\quad \hat x_{0}\in \Rset^n
\label{hatx-}\\
\hat x_{0}(t_i^+) &=& \hat x_{0}(t_i^-)+
{K_{0}}\nu (t_i), \quad \nu (t_i)=y(t_i) - {H_{0}}\hat x_{0}(t_i^-)
\label{hatx+} \\
u(t_i) &=& f_{\mathrm{R}}[\hat x_{0}(t_i^+)] = - {G_{0}^{\star}}
\hat x_{0}(t_i^+), \quad u\in \Rset^q  \label{u}
\end{eqnarray}
with $\{w(\cdot)\}$, $\{v(\cdot)\}$ being zero mean i.i.d.
sequences of covariances $Q_{\theta}$ and $R_{\theta}$
respectively. The state difference equation  (\ref{x}) is
propagated forward from the initial condition $x(t_{-s})$ of the
finite expectation $\E{ \|x(t_{-s})\|^2 }<\infty$. The initial
state is placed at some $t_{-s}\in \Rset$ where $s>0$ determines
what is termed settling time $T_{s}=t_{0}-t_{-s}$ needed to think
of all processes in (\ref{x})--(\ref{u}) as wide-sense stationary
for $i\ge 0$. As usual, equation (\ref{x}) represents a plant,
equation (\ref{y}) a sensor, and (\ref{hatx-})--(\ref{u}) a
feedback. The feedback is composed of a Kalman-like filter
(\ref{hatx-})--(\ref{hatx+}) cascaded with a regulator (\ref{u}).
The regulator is described by a function $f_{\mathrm{R}}[\cdot]$
of the measurement updated estimate $\hat x_{0}(t_i^+)$ or can be
chosen according to the second equality in (\ref{u}) with a matrix
${G_{0}^{\star}}$.

Matrices of the given system (\ref{x})--(\ref{y}) are assumed to
be known as  $\Phi_{0}$, $\Psi_{0}$, $Q_{0}$, $H_{0}$ and $R_{0}$
for a \emph{nominal mode}, i.~e., for a \emph{nominal value}
$\theta_{0}$ of the uncertainty parameter $\theta$. For this mode,
to guarantee existence of the steady-state filter with Kalman gain
$K_{0}$, we assume that the pair $(\Phi_{0}, Q_{0}^{1/2})$ is
\emph{stabilizable},  the pair $(\Phi_{0}, H_{0})$ is
\emph{observable}, and the pair $(\Phi_{0}, \Psi_{0})$ is
\emph{controllable}. Matrix $G_{0}^{\star}$ can be taken then as a
given one or designed to be LQG optimal for the nominal mode of
operation~\cite{maybeck:1982}, \cite{caines:1988}, \cite{mosca:1995}.

Parameter $\theta$ is  subject to an abrupt change at an unknown
time point $t_{\mathrm{c}}\in (t_{0}, t_{i})$. This can be viewed
as a switch of $\theta$ from $\theta_{0}$ to some other unknown
value $\theta_{1}$; to detect the change point, a \emph{decision
generator}, $\mathscr{DG}$, is needed~(Figure~\ref{fig1}).
Formalizing the problem of $\mathscr{DG}$ synthesis, we consider
\begin{equation}
\left.
\begin{aligned}
x_{j}(t_{i+1})&=\Phi_{j} x_{j}(t_i)+\Psi_{j} u(t_i)+ w_{j}(t_i) \; \\
y_{j}(t_i) &= H_{j} x_{j}(t_i) + v_{j}(t_i) \;
\end{aligned}
\right\} \; \mathscr{S}_{0} \;\; (\mbox{for } j=0) \quad \mbox{or} \quad
\mathscr{S}_{1} \;\; (\mbox{for } j=1) \label{xyj}
\end{equation}
as the generalized  description of the two systems:
$\mathscr{S}_{0}$ the system with $\theta=\theta_{0}$;
$\mathscr{S}_{1}$  the system with $\theta=\theta_{1}$. Filter
(\ref{hatx-})--(\ref{hatx+}), which is denoted by
$\mathscr{F}_{0}$ in~Figure~\ref{fig1}, is designed as the Kalman
filter for $\mathscr{S}_{0}$, i.~e., satisfying equations
\begin{equation}
\left.
\begin{aligned}
K_{0} &= \wt {P_{0}}H_{0}^{T}C_{0}^{-1}\ ,     &
C_{0} &= H_{0}\wt {P_{0}}H_{0}^{T} + R_{0}\;
\\
\wh {P_{0}} &= \wt {P_{0}}-\wt {P_{0}}H_{0}^{T}C_{0}^{-1}H_{0}\wt {P_{0}}\ ,
& \wt {P_{0}} &= \Phi_{0}\wh {P_{0}}\Phi_{0}^{T}+ Q_{0}\;
\end{aligned}
\right\} \;\mbox{filter }  \mathscr{F}_{0} \;\; \mbox{for } \mbox{system }
\mathscr{S}_{0} \label{F0}
\end{equation}
and processing data
\begin{equation}
y(t_{i})= \left \{
\begin{aligned}
y_{0}(t_{i}) \;\; &\mbox{if}& t_{i} < t_{\mathrm{c}}   &&(\mbox{data from } \mathscr{S}_{0})\\
y_{1}(t_i) \;\;   &\mbox{if}& t_{i} \ge t_{\mathrm{c}} &&(\mbox{data from } \mathscr{S}_{1})
\end{aligned}
\right. \label{y(t)}
\end{equation}
The problem is to  detect the change point $t_{\mathrm{c}}$ with a
reasonable delay using a suitable decision rule $d_{0}(t_{k})\in
\{0,\ 1 \}$ according to Figure~\ref{fig1}.

\begin{figure}
\input pict1.txt
\caption{The general  framework for the problem and solution. {\sl
Legend:} $\mathscr{S}_{0}$ is the system with $\theta=\theta_{0}$
while $\mathscr{S}_{1}$ is the system with $\theta=\theta_{1}$;
$\mathscr{F}_{0}$ is the Kalman filter for $\mathscr{S}_{0}$; and
$\mathscr{DG}$ is a decision generator;  $q^{-1}$ denotes a one
sample memory (unit delay); \ $\nu (t_{1}^{k})\triangleq \col
\left[ \nu (t_{1}), \nu (t_{2}), \cdots ,  \nu (t_{k})\right]$~.
\label{fig1}}
\end{figure}

\section{\hspace{-2ex}{.}\hspace{1ex}WSI BASED DECISION GENERATING METHOD}\label{method}
\noindent
Many systems, not only  (\ref{x})--(\ref{u}), have the IS, hence
the WSI of $s_{i}$ type in (\ref{Sk}) can be considered. Because
covariance properties of the WSI after change may become very
complex, depending on the change occurred, we simplify matters by
approximating the  covariance kernel function
$\mathcal{K}_{ss}(j)$ of  $s_{i}$ after change equal to the
exponential function:
\begin{equation}\label{corr}
\mathcal{K}_{ss}(j)=D_{s}r^{|j|}\ , \quad     j\in \Zset \ , \quad r=\exp\{-\tau_{\mathrm{s}}/\tau_{\mathrm{c}}\}\ , \quad
D_{s}\triangleq \E{(s_{i}- \E{s_{i}}   )^2}
\end{equation}
where $\tau_{\mathrm{s}}$  is a sampling time interval and
$\tau_{\mathrm{c}}$ a correlation interval of $s_{i}$. Under this
assumption, we state the WSI based  method for detecting changes
as follows.
\begin{theo}\label{theo1}
Let a dynamical system with the IS $\{\nu (t_{i})\}$, $\nu
(t_{i})\in \Rset^{m}$,  have two possible modes of operation
termed ``normal functioning'' (hypothesis $\Hyp_{0}$) and ``a
fault'' (hypothesis $\Hyp_{1}$). Under both hypotheses, $\{\nu
(t_{i})\}$ is a Gaussian sequence with $\E{\nu (t_{i})}=0$. Under
hypothesis $\Hyp_{0}$, $\{\nu (t_{i})\}$ is an independent
sequence with the known covariance $C_{0}=L_{0}L_{0}^T$  of each
element $\nu (t_{i})$ where $L_{0}$ is the square root of  $C_{0}$
determined, for instance, by Cholesky decomposition. Under
hypothesis  $\Hyp_{1}$, $\{\nu (t_{i})\}$ is a correlated sequence
with an unknown covariance $C_{1}\ne C_{0}$ of each element $\nu
(t_{i})$. Define
\begin{equation}\label{mu}
\mu_{i}\triangleq L_{0}^{-1}\nu (t_{i}) \ , \quad
s_{i}\triangleq   \sqrt{\frac{m}{2}}
\left( \frac{1}{m} %(1/m)
\mu_{i}^{T}\mu_{i}-1 \right) \ , \quad
S_{k}\triangleq \frac{1}{\sqrt{k}}
%1/\sqrt{k}
\sum_{i=1}^{k}s_{i}
\end{equation}
and suppose that under  hypothesis $\Hyp_{1}$, the  covariance
kernel function of this $s_{i}$ has been modeled (maybe rather
approximately) by expression
(\ref{corr}).\\
$\hphantom{...}$ Then  asymptotically, when $k$ goes to infinity,
the following properties hold:
\begin{enumerate}
\item The probability  laws  $\mathcal{L} (\cdot)$ of $S_{k}$ under these hypotheses satisfy \footnote{%$^{)}$
The notation $\leadsto$  corresponds to the weak
convergence~\cite{roussas:1972}.}%$^{)}$
\begin{eqnarray}
\Hyp_{0}: \; & \mathcal{L} (S_{k}) & \leadsto \mathcal{N}(0, 1) \label{H0}\\
\Hyp_{1}: \; & \mathcal{L} (S_{k}) & \leadsto
\mathcal{N}(m_{S_{k}}, D_{S_{k}}) \label{H1}
\end{eqnarray}
where  $\mathcal{N}(0; 1)$ stands for the Gaussian (normal)
distribution with zero mean and unit variance,
$\mathcal{N}(m_{S_{k}}, D_{S_{k}})$ is the normal distribution
whose mean $m_{S_{k}}$ and variance $D_{S_{k}}=\sigma_{S_{k}}^{2}$
are given by the below formulae
\begin{equation}\left.\begin{split}
m_{S_{k}}\triangleq & \E{S_{k}}=m_{s}\sqrt{k}\ , \quad
m_{s}\triangleq \E{s_{i}}=(1/\sqrt{2m})\tr{\Delta}\\
D_{S_{k}}\triangleq & \E{(S_{k}-m_{S_{k}})^2}=D_{s}\frac{1+r}{1-r}\ , \quad
D_{s}=1+\frac{2}{m}\tr{\Delta}+\frac{1}{m}\|\Delta \|^{2} \\
\Delta =&L_{0}^{-1}(C_{1}-C_{0})L_{0}^{-T}\ , \quad
\|\Delta \|^{2}\triangleq  \sum_{i, j =1}^{m} | \Delta_{ij} |^{2}
\end{split}\quad\right \}\label{Delta}\end{equation}
\item
If $S_{k}$  of (\ref{mu}) is used as a decision function at a time
instant $k$ in the decision rule
\begin{equation}\label{rule}
\left | S_{k} \right | \begin{array}{c}\Hyp_{1}\\ \gtrless\\\Hyp_{0}\end{array} h
\end{equation}
with $h$ being  a conveniently chosen threshold (for example,
$h=3$), then probability of false alarm, $P_{F}$, and probability
of detection, $P_{D}$, satisfy the following approximate
expressions (they are exact if laws $\mathcal{L}(S_{k})$  in
(\ref{H0}), (\ref{H1}) are taken to be normal):
\begin{equation} \left. \begin{split}
P_{F}&\simeq 1-\phi(h)\\
P_{D}&\simeq  1-\frac{1}{2}\left[ \phi\left(\frac{m_{S_{k}}+h}{\sigma_{S_{k}}}\right) - \phi\left(\frac{m_{S_{k}}-h}{\sigma_{S_{k}}}\right) \right] > \frac{1}{2}\left[ 1+ \phi\left(\frac{m_{S_{k}}-h}{\sigma_{S_{k}}}\right) \right]
\end{split} \;  \right\}\label{P}\end{equation}
through the standard probability integral
\[
\phi(x)\triangleq \frac{2}{\sqrt{2\pi}}\int_{0}^{x}
\exp{\left(-t^{2}/2\right)} \, \ud t
\]
\item If the  threshold $h$ in (\ref{rule}) is chosen as $h=\phi^{-1}(1-\alpha)$ to guarantee a given level $\alpha$ of $P_{F}$ in (\ref{P}), $P_{F}=\alpha$, then
hypothesis  $\Hyp_{1}$ is detected with probability $P_{D}=1-\beta$
(where $0<\beta<1/2$) not later than after
\begin{equation}\label{t-alarm}
k_{a}^{\star}\simeq \left[ \phi^{-1}(1-\alpha)+\sigma_{S_{k}}%\sqrt{D_{s}\frac{1+d}{1-d}}
\phi^{-1}(1-2\beta)\right]^{2}/{m_{s}^{2}}
\end{equation}
disctere time instants where $x=\phi^{-1}(y)$ is the solution  for
$\phi(x)=y$.
\end{enumerate}
\end{theo}
The proof of Theorem~\ref{theo1} is placed in Appendix~A.

Clearly  sequence $S_{k}$ can be used for any system having such
an innovation $\nu (t_{i})$ to detect the changes marked by matrix
$\Delta$ with $\tr{\Delta}\ne 0$, with the aid of rule
(\ref{rule}). In so doing, it is important to keep in mind the
following merits of the rule:
\begin{enumerate}
\item
Rule  (\ref{rule}) contains the same quadratic forms as does the
log-likelihood function
\[
\ln p \left( \nu (t_{1}^{k}) \bigm| \Hyp_{0} \right)=
- \frac{km}{2}\ln (2\pi) - \frac{k}{2}\ln \lvert C_{0}\rvert
- \frac{1}{2}\sum_{i=1}^{k}\nu ^{T}(t_{i})C_{0}^{-1}\nu (t_{i})
\]
with $\nu (t_{1}^{k})\triangleq \col \left[ \nu (t_{1}),
\nu (t_{2}), \cdots ,  \nu (t_{k})\right]$
because
\begin{equation}\label{Sksum}
S_{k}= \sqrt{\frac{m}{2k}}\sum_{i=1}^{k}\left( \frac{1}{m}
\nu ^{T}(t_{i})C_{0}^{-1}\nu (t_{i})-1 \right)
\end{equation}
\item
The  sufficient statistic $\ell(\cdot)\triangleq \ell\left( \nu
(t_{1}^{k}) \right)$ for detecting changes in the system
parameters subject to the assumptions of Theorem~\ref{theo1},
is~\cite[Sec.~2.6.2]{vantrees:1968}
\[
\ell(\cdot)=
\sum_{i=1}^{k}\nu ^{T}(t_{i})C_{0}^{-1}\nu (t_{i}) -
\nu ^{T}(t_{1}^{k}) \boldsymbol{\mathcal{C}}_{1}^{-1}(k)\nu (t_{1}^{k})\ , \;
\boldsymbol{\mathcal{C}}_{1}(k)\triangleq \E{\nu (t_{1}^{k})\nu ^{T}(t_{1}^{k}) \bigm| \Hyp_{1}}
\]
The  evident theoretical drawback of the method is that the
decision function $S_{k}$, (\ref{mu}) or equivalently
(\ref{Sksum}), is not a sufficient statistic; nevertheless it has
found practical use in many applications since its invention, as
working without any information about the model after the change,
for example in checking filter optimality~\cite{semoushin:1979}.
\item
The  type of parametric changes  detectable by rule (\ref{rule})
is restricted to those for which $\tr{\Delta}\ne 0$.
\item
To  strictly satisfy the condition $\tr{\Delta}= 0$ in most real
situations is hardly probable and because of this the restriction
$\tr{\Delta}\ne 0$ should not be considered as critical from the
practical point of view.
\end{enumerate}

%%\subparagraph{Example 1} 
\ex Let $m$ be equal $3$ and
\begin{equation*}
C_{0}=
\left[
\begin{array}{llr}
1 & 2 & 3 \\
2 & 8 & 2 \\
3 & 2 & 14
\end{array}
\right],   \quad
C_{1} = \left[
\begin{array}{llr}
2.0 & 1.8 & 2.4 \\
1.8 & 9.0 & 1.0 \\
2.4 & 1.0 & 13.0
\end{array}
\right]
\end{equation*}
Then
\begin{align*}
L_{0}=
\left[
\begin{array}{lrr}
1 & 0 & 0 \\
2 & 2 & 0 \\
3 & -2 & 1
\end{array}
\right],     &
& \Delta = \left[
\begin{array}{rrr}
1.0 & -1.1 & -5.8 \\
-1.1 & 1.45 & 6.3 \\
-5.8 & 6.3 & 31.0
\end{array}
\right],     &
& m_{s}\approx 13.7\ ,      &
& D_{s}\approx 395.
\end{align*}
Let $\alpha=0.005$ and $\beta=0.005$. We obtain,
$h=\phi^{-1}(1-\alpha)= 2.8$ and $\phi^{-1}(1-2\beta)=2.6$\ .

If $r=0.5$, then $k_{a}^{\star}\simeq 46$ and if $r=0.95$, then
$k_{a}^{\star}\simeq 562$. Hence, starting from the anticipated
(for instance, from the most probable or most dangerous) faults as
the base, and reasoning from the specified quality $\alpha$ and
$\beta$ of the decision rule (\ref{rule}), it is possible to
numerically predict the lapse of time\footnote{$^{)}$ In this
example, the  delay $k_{a}^{\star}$ should be increased by
$k_{\mathrm{norm}}\gtrapprox 30$ needed to take distributions in
(\ref{H0})--(\ref{H1}) to be approximately normal, from the
engineering point  of view.}$^{)}$ $k_{a}^{\star}$ after which the
fault will be detected with confidence not worse than
$P_{D}=1-\beta$.

For the above example,  the condition
of non-observable faults is readily obtained as the equation
\begin{equation}\label{trDelta}
\tr{\Delta}=27a_{11}-11a_{12}-10a_{13}+{\textstyle
\frac{5}{4}}a_{22}+2a_{23}+a_{33}=0
\end{equation}
in the entries $a_{ij}=a_{ji}$ of $A= C_{1} - C_{0}$~.  Keeping
$C_{1}$ within  bounds of positive definiteness, one can choose
$C_{1}$ to match the condition (\ref{trDelta}), for example
\begin{equation*}
C_{1}=
\left[
\begin{array}{llr}
1 & 2 & 3 \\
2 & 8 & 1 \\
3 & 1 & 16
\end{array}
\right],   \quad
A\triangleq C_{1} - C_{0} = \left[ a_{ij} \right] = \left[
\begin{array}{lrr}
0 & 0 & 0 \\
0 & 0 & -1 \\
0 & -1 & 2
\end{array}
\right]
\end{equation*}
This simple example makes it apparent that only certain of the faults
specially chosen to fit the equation $\tr{\Delta}=0$ to a T  cannot
be detected  by the method. \\

Theorem~\ref{theo1} assumes that changes in $\theta$ do not create
a bias (i.e., non-zero expectation) in the
innovation sequence $\{\nu (t_i)\}:$
\begin{alignat}{4}\label{Hyps}
\Hyp_{0}:\ \; &  \nu (t_{i}) \sim \mathcal{N}(0, C_{0})
&\qquad \mbox{and} &\qquad &
\Hyp_{1}:\ \; &  \nu (t_{i}) \sim \mathcal{N}(0, C_{1})
\end{alignat}
If this is not the case, one should turn from (\ref{Hyps}) to
\begin{alignat}{4}\label{Hypsnew}
\Hyp_{0}:\ \; &  \nu (t_{i}) \sim \mathcal{N}(0, C_{0})
&\qquad \mbox{and} &\qquad &
\Hyp_{1}:\ \; &  \nu (t_{i}) \sim \mathcal{N}(m_{1}, C_{1})
\end{alignat}
with some $m_{1}\ne 0$. This entails the
only modification in Theorem~\ref{theo1}, namely, in $\Delta$ of
formulae (\ref{Delta}):
\begin{alignat}{4}\label{Deltanew}
\Delta =L_{0}^{-1}(C_{1}+m_{1}m_{1}^{T} -C_{0})L_{0}^{-T}
&\qquad \mbox{instead of} &\qquad &
\Delta =L_{0}^{-1}(C_{1}-C_{0})L_{0}^{-T}
\end{alignat}

Thus, for  appreciable quantity of faults, the WSI based method
boils down to a straightforward cumulation of $S_{k}$, (\ref{mu}),
through the sequential algorithm
\begin{equation}\label{Skalg}
S_{k}=S_{k-1}\sqrt{1- 1/k} + s_{k}/\sqrt{k}\ , \quad k\in \Nset\ , \quad
S_{0}=0
\end{equation}
followed by  the rule (\ref{rule}). When testing $\Hyp_{1}$
against $\Hyp_{0}$, let us denote $S_{k}$ in (\ref{Skalg}) by
$S_{k}^{(0)}$ thus emphasizing by superscript $^{(0)}$ that the
main hypothesis is $\Hyp_{0}$. Then we have two schemes with the
stopping rule $t_{\mathrm{a}}$ (the alarm time) for the
alternative $\Hyp_{1}$ detection:
\begin{description}
\item[(A)]  decision is made at any current time $t_{k}$:
\[
d_{0}(t_{k})=\left\{
\begin{aligned}
0 \quad & \mbox{if }\; \left| S_{k}^{(0)} \right| < h \ ;  \quad \Hyp_{0}\ \, \mbox{is\, chosen}\\
1 \quad & \mbox{if }\; \left| S_{k}^{(0)} \right| \ge h \ ;  \quad \Hyp_{1}\ \, \mbox{is\, chosen}
\end{aligned}
\right\} , \quad t_{\mathrm{a}}= \min \left\{ t_{k} : \; d_{0}(t_{k})=1  \right\}
\]
\item[(B)] decision is made at the end of a sample number $l=1,2,\dots, L$
each of size $N$:
\begin{align*}
d_{0}(l)&=\left\{
\begin{aligned}
0 \quad & \mbox{if }\; \forall \ k=1,2,\dots,N : \; \left| S_{N(l-1)+k}^{(0)} \right| < h \ ;  \quad \Hyp_{0}\ \, \mbox{is\, chosen}\\
1 \quad & \mbox{if }\; \exists \ k=1,2,\dots,N : \; \left| S_{N(l-1)+k}^{(0)} \right| \ge h \ ;  \quad \Hyp_{1}\ \, \mbox{is\, chosen}
\end{aligned}
\right. \\
t_{\mathrm{a}}&=\tau_{\mathrm{s}} \left[N(l-1) + \min \left\{ k : \; \left| S_{N(l-1)+k}^{(0)} \right| \ge h   \right\} \right]\stackrel{r}{=}
\tau_{\mathrm{s}} N \min  \left\{ l : \; d_{0}(l)=1   \right\}
\end{align*}
where $\stackrel{r}{=}$ denotes the rounded up equality.
\end{description}

This result is of general significance because many  stochastic
dynamical systems have (or can be  supplied with) a whitening
filter that generates a process of innovations.  It fits naturally
into the systems of (\ref{x})--(\ref{u}) type, in which the
whitening filter exists in the Kalman filter form
(\ref{hatx-})--(\ref{hatx+}) with the IS $\{\nu (t_i)\}$.

\section{\hspace{-2ex}{.}\hspace{1ex}WSI BASED FEEDBACK FILTER TESTING AND SELECTION}
\noindent
From the innovation process theory, two facts are
known~\cite{martin:1976}:\\[6pt]
\begin{tabular}{cp{14.5cm}}
{\bf (i)}  & Case~(\ref{Hyps}) holds if changes happen only in
covariances $Q_{\theta}$ and $R_{\theta}$ . \\
{\bf (ii)} & The only cause for
Case~(\ref{Hypsnew}) are changes in matrices $\Phi_{\theta}$,
$\Psi_{\theta}$ and/or $H_{\theta}$ .\\[6pt]
\end{tabular}

\noindent
If Case~(\ref{Hypsnew}) takes
place, (\ref{Deltanew}) with $m_{1}$ calculated by the general
procedure described, for example in~\cite{martin:1976}, should be used.

\ex%
%\subparagraph{Example 2} 
Consider the simple case when $n=m=1$,
$u(t_{i})\equiv 0$, $Q_{1}=Q_{0}$, $R_{1}=R_{0}$, $H_{1}=H_{0}=1$
and $\Phi_{1}\ne \Phi_{0}$ (both $\lvert\Phi_{j}\rvert < 1, \ j=0,1$).
At $t_{i}\gg t_{\mathrm{c}}$ (i.e., under $\Hyp_{1}$ when the change
transient response effects are vanished) from
(\ref{xyj}) (for $j=0$), (\ref{F0}), (\ref{y(t)}) and
(\ref{hatx-})--(\ref{hatx+}), we have
\begin{equation}
\nu (t_{i})=x_{1}(t_{i})+v_{0}(t_{i})-\hat x_{0}(t_{i})=
e(t_{i}) + v_{0}(t_{i})\ ,
\quad
e(t_{i})\triangleq x_{1}(t_{i})-\hat x_{0}(t_{i})   \label{e}
\end{equation}
Denoting
\begin{equation}
\left.
\begin{aligned}
m_{1}(t_{i}) &\triangleq \E{\nu(t_{i})\bigm|\Hyp_{1}} \ ,     &
m_{e}(t_{i}^{-}) &\triangleq \E{e(t_{i})} \;
\\
m_{x_{1}}(t_{i}) &\triangleq \E{x_{1}(t_{i})} \ ,
&
m_{\hat x_{0}}(t_{i}^{\pm}) &\triangleq \E{\hat x_{0}(t_{i}^{\pm})} \;
\end{aligned}
\right\}  \label{def}
\end{equation}
we obtain from (\ref{e}), (\ref{xyj}) (for $j=1$) and (\ref{hatx-})
\begin{equation}
\left.
\begin{aligned}
m_{1}(t_{i}) &= m_{e}(t_{i}^{-}) \ ,     &
m_{e}(t_{i}^{-}) &= m_{x_{1}}(t_{i}) - m_{\hat x_{0}}(t_{i}^{-})\;
\\
m_{x_{1}}(t_{i}) &= \Phi_{1} m_{x_{1}}(t_{i-1})\ ,
&
m_{\hat x_{0}}(t_{i}^{-}) &= \Phi_{0} m_{\hat x_{0}}(t_{i-1}^{+})\;
\end{aligned}
\right\}  \label{ms}
\end{equation}
Taking expectation of (\ref{hatx+}) yields
\[
m_{\hat x_{0}}(t_{i}^{+}) =(1-K_{0})m_{\hat x_{0}}(t_{i}^{-}) +
K_{0}m_{x_{1}}(t_{i})
\]
Taking here a step backwards and then substituting this into the last
formula of (\ref{ms}) give
\[
m_{\hat x_{0}}(t_{i}^{-}) =\Phi_{0}\left[ (1-K_{0})m_{\hat
x_{0}}(t_{i-1}^{-}) + K_{0}m_{x_{1}}(t_{i-1})  \right]
\]
Then we find consequently using (\ref{def}) and (\ref{ms})
\begin{eqnarray*}
m_{e}(t_{i}^{-}) &=& m_{x_{1}}(t_{i}) -
\Phi_{0}\left[ (1-K_{0}) m_{x_{1}}(t_{i-1}) -
(1-K_{0}) m_{e}(t_{i-1}^{-}) +
K_{0}m_{x_{1}}(t_{i-1})  \right]\\
 &=& m_{x_{1}}(t_{i}) -
\Phi_{0}\left[ m_{x_{1}}(t_{i-1}) - (1-K_{0}) m_{e}(t_{i-1}^{-})  \right]\\
 &=& (\Phi_{1}-\Phi_{0}) m_{x_{1}}(t_{i-1}) +
\Phi_{0}(1-K_{0})  m_{e}(t_{i-1}^{-})\\
m_{1}(t_{i})  &=& \Phi_{0}(1-K_{0})  m_{1}(t_{i-1})+
(\Phi_{1}-\Phi_{0}) m_{x_{1}}(t_{i-1})
\end{eqnarray*}
As $0<\Phi_{0}<1$ and $0<\Phi_{0}(1-K_{0}) <1$, there exist the
steady-state values
\[
m_{1}\triangleq  \lim_{t_{i}\to \infty} \{m_{1}(t_{i})\}=
\lim_{t_{i}\to \infty} \{m_{1}(t_{i-1})\} \ , \qquad
m_{x_{1}} \triangleq  \lim_{t_{i}\to \infty} \{ m_{x_{1}}(t_{i-1}) \}
\]
and finally
\[
m_{1}=\frac{\Phi_{1}-\Phi_{0}}{1- \Phi_{0}(1-K_{0})}m_{x_{1}}
\]
\\
%{\bf N.~B.\;} 
\rem By virtue of $0<\Phi_{0}<1$, $m_{x_{1}}=0$. It is
also worthy of note that although $\tr{\Delta}=0$ is not a
necessary or sufficient condition for filter optimality, the
quantity of faults satisfying $\tr{\Delta}=0$ is relatively thin.\\

As for the case of two contesting hypotheses $\Hyp_{0}$ and
$\Hyp_{1}$ in Section~\ref{method}, we can apply one of two
schemes, ({\bf A}) or ({\bf B}), of decision making. However as
differentiated from that case, here every singly taken $\Hyp_{j} \
(j=0,1,\dots, K)$ contests individually against any other
$\Hyp_{i} \  (i=0,1,\dots, K; i\ne j)$ using its own decision
function $S_{k}^{(j)}$ that equals to $S_{k}$ and is formed
according to (\ref{Skalg}) supplied with the superscript $^{(j)}$.
Again, we have two schemes:
\begin{description}
\item[(A)] decision is made at  any current time $t_{k}$. For any
particular  $\Hyp_{j}, \ j\in \{ 0,1,\dots, K \}$, the particular decision rule is
\[
d_{j}(t_{k})=\left\{
\begin{aligned}
0 \quad & \mbox{if }\; \left| S_{k}^{(j)} \right| < h \ ;  \quad \Hyp_{j}\ \, \mbox{is\, supported}\\
1 \quad & \mbox{if }\; \left| S_{k}^{(j)} \right| \ge h \ ;  \quad \Hyp_{j}\ \, \mbox{is\, rejected}
\end{aligned}
\right.
\]
The over-all alarm time is
\[
t_{\mathrm{a}}= \min \left\{ t_{k} : \;
\left(\exists \kappa: d_{\kappa}(t_{k})=0  \;\:  \&
\mathop{\forall}_{j\ne \kappa} d_{j}(t_{k})=1
\right)
  \right\}
\]
with the over-all winner selection rule:
\[
\kappa= \left\{ j\in \{ 0,1,\dots,K \}: d_{j}(t_{\mathrm{a}})=0
\right\}\ ;  \quad  \Hyp_{\kappa}\ \,  \mbox{is\, chosen}
\]

\item[(B)] decision is made at the end of a sample number $l=1,2,\dots,
L$ each of size $N$.  For any particular
$\Hyp_{j}, \ j\in \{ 0,1,\dots, K \}$, the particular decision rule is
\begin{align}
d_{j}(l)&=\left\{
\begin{aligned}
0 \quad & \mbox{if }\; \forall \ k=1,2,\dots,N : \; \left| S_{N(l-1)+k}^{(j)} \right| < h \ ;  \quad \Hyp_{j}\ \,
\mbox{is\, supported}\\
1 \quad & \mbox{if }\; \exists \ k=1,2,\dots,N : \; \left| S_{N(l-1)+k}^{(j)} \right| \ge h \ ;  \quad \Hyp_{j}\ \,
\mbox{is\, rejected}
\end{aligned}
\right. \label{schemeB}\\
t_{\mathrm{a}}&\stackrel{r}{=}
\tau_{\mathrm{s}} N \min  \left\{ l : \; \left(
\exists \kappa :  d_{\kappa}(l)=0  \;\:  \&
\mathop{\forall}_{j\ne \kappa} d_{j}(l)=1  \right)
\right\}\label{schemeB'}
\end{align}
with the over-all winner selection rule:
\begin{equation}\label{winnerB}
\kappa= \left\{ j\in \{ 0,1,\dots,K \}: d_{j}(l)=0
\right\}\ ;  \quad  \Hyp_{\kappa}\ \,  \mbox{is\, chosen}
\end{equation}
\end{description}
This method is realized in a framework corresponding to
Figure~\ref{fig2}.

\begin{figure}
\input pict2.txt
\caption{The best filter selection. {\sl Legend:} $\mathscr{SP}$
stands  for system population; $\mathscr{FP}$  stands for filter
population, $\mathscr{FS}$  for filter selection, $\mathscr{FE}$
for feedback effector, and  $\mathscr{SF}$  stands for the system
feedback. The stopping rule is given by the alarm time
$t_{\mathrm{a}}=\min \Bigl\{ t_{k}: \; \bigl( \exists \kappa : \
d_{\kappa}(t_{k})=0 \;\: \& \mathop{\forall}\limits_{j\ne\kappa}
d_{j}(t_{k})=1 \bigr) \Bigr\}$\ . At this time, $\mathscr{FE}$
uploads the selected $\mathscr{F}_{\kappa}$ into $\mathscr{SF}$ if
by then it is not in there.
%has not yet been there.
\label{fig2}}
\end{figure}

\section{\hspace{-2ex}{.}\hspace{1ex}SOME COMPUTATIONAL EXPERIMENT RESULTS}\label{simul}
\noindent
To have  confidence in the freedom from failures of airborne
equipment or changes in vehicle parameters is of critical
importance for optimal navigation data processing.

As an  example for the problem under study, we take the damped
Shuler loop driven by the exponentially correlated noise, whose
description was first given  in~\cite{gaines:1971}. In compliance
with~\cite{gaines:1971}, consider equations
\begin{equation}
\left.
\begin{aligned}
x_{j}(t_{i+1})&=\Phi_{j} x_{j}(t_i)+ \Gamma_{j} w_{j}(t_i) \; \\
y_{j}(t_i) &= H_{j} x_{j}(t_i) + v_{j}(t_i) \;
\end{aligned}
\right\} \; \mathscr{S}_{0} \;\; (\mbox{for } j=0) \quad \mbox{or} \quad
\mathscr{S}_{j} \;\; (\mbox{for } j=1,\cdots , K) \label{xyN}
\end{equation}
that  have the form of (\ref{xyj}) in the specific case of control
input $u(t_i)$ being zero, with $K+1$ modes of operation, each
mode being associated with a system $\mathscr{S}_{j}$:\,
$\mathscr{S}_{0}$ is the fault free system and every
$\mathscr{S}_{j} \;\; (\mbox{for } j=1,\cdots , K)$ is a faulty
system. System $\mathscr{S}_{0}$ corresponds to
\begin{align}
\Phi_{0}=
\left[
\begin{array}{lllll}
0.75 & -1.74 & -0.3 & 0 & -0.15 \\
0.09 & 0.91 & -0.0005 & 0 & -0.008 \\
0 & 0 & 0.95 & 0 & 0 \\
0 & 0 & 0 & 0.55 & 0 \\
0 & 0 & 0 & 0 & 0.905
\end{array}
\right],     &
& \Gamma_{0} = \left[
\begin{array}{lll}
0 & 0 & 0 \\
0 & 0 & 0 \\
24.64 & 0 & 0 \\
0 & 0.835 & 0 \\
0 & 0 & 1.83
\end{array}
\right]
\label{sys0'}
\end{align}
\begin{align}
Q_{0} &= \left[
\begin{array}{lll}
1 & 0 & 0  \\
0 & 1 & 0  \\
0 & 0 & 1
\end{array}
\right],     &
H_{0} &= \left[
\begin{array}{lllll}
1 & 0 & 0 & 0 & 1 \\
0 & 1 & 0 & 1 & 0
\end{array}
\right],     &
R_{0} &= \left[
\begin{array}{ll}
1 & 0 \\
0 & 1
\end{array}
\right]
\label{sys0''}
\end{align}
Disengaging ourselves at this stage from a physical interpretation of 
the given
characteristics (\ref{sys0'})--(\ref{sys0''}), we shall restrict our
consideration in this section by introducing the set of faulty
systems through the following formal relations:
\begin{equation*}
\left.
\begin{aligned}
\Phi_{j}&=\Phi_{0}\ ,  \quad  \Gamma_{j}=\Gamma_{0}\ , \quad
Q_{j}=Q_{0}\ , \quad  R_{j}=R_{0} \\
H_{j}&=  \left[
\begin{array}{ccccc}
1-e & 0 & 0 & 0 & 1-f \\
0 & 1-g & 0 & 1-h & 0
\end{array}
\right], \quad  \{ \mathit{efgh} \}\triangleq j_{[2]} \;\;
\end{aligned}
\right\} \;\; \mathscr{S}_{j} \;\;  (j=1,\cdots , 15)  \label{Sj}
\end{equation*}
where $\{ \mathit{efgh} \}$ stands for the binary code $j_{[2]}$ of the
fault number $j$. For example, given $j=2$  or  $j=9$, then $\{
\mathit{efgh} \}=\{ 0010 \}$ or, correspondingly $\{ 1001 \}$, and in these
cases
\[
H_{2}=  \left[
\begin{array}{ccccc}
1 & 0 & 0 & 0 & 1 \\
0 & 0 & 0 & 1 & 0
\end{array}
\right], \quad
H_{9}=  \left[
\begin{array}{ccccc}
0 & 0 & 0 & 0 & 1 \\
0 & 1 & 0 & 0 & 0
\end{array}
\right]
\]

The  change $\mathscr{S}_{0}
\stackrel{t_{\mathrm{c}}}{\rightarrowtail} \mathscr{S}_{j} \
(j=1,\cdots , 15)$ occurs at $t_{\mathrm{c}}=2300$, however filter
$\mathscr{F}_{0}$ ``does not know'' about that and continues to
satisfy equations (\ref{F0}). Experimental data for all $K=15$
faulty system cases are shown in Figures~\ref{Case1}--\ref{Case2}
together with the fault free system case, when $\{ \mathit{efgh} \}=\{ 0000\}$. If the value of threshold $h$ in (\ref{rule}) is set to $3$,
then probability of false alarm, $P_{F}$, in the case of
$\mathscr{S}_{0}$ will be not greater than 0.3\% (due to property
(\ref{H0})), and in other cases probability of right detection,
$P_{D}$, comes to 100\% due to monotonicity of increasing
$|\E{S_{k}}|$. The value of delay in right detection depends on
the type of change and, as it can be seen from
Figures~\ref{Case1}--\ref{Case2}, is in the region from $100$ to
$1500$ sampling time intervals.

The  reported simulation results as well as those obtained for
many different samples from the sequences $w(t_i)$ and $v(t_i)$ in
the model (\ref{x})--(\ref{y}), show practicability of the
approach for inertial navigation systems monitoring and thus cause
us to anticipate that it will find a diversity of other
applications. This assumption has received further support in our
work through the similar results  yielded by other experiments.

\begin{figure}[h]
\centering
\includegraphics[scale=0.4]{0415.eps}
\includegraphics[scale=0.4]{2637.eps}
\caption{Behavior  of $S_{k}$ before and after change at
$t_{\mathrm{c}}=2300$; {\sl left:} $j=0,4,1,5$, {\sl right:}
$j=2,6,3,7$.} \label{Case1}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[scale=0.4]{812913.eps}
\includegraphics[scale=0.4]{10141115.eps}
\caption{Behavior  of $S_{k}$ before and after change at
$t_{\mathrm{c}}=2300$; {\sl left:} $j=8,12,9,13$, {\sl right:}
$j=10,14,11,15$.} \label{Case2}
\end{figure}

\section{\hspace{-2ex}{.}\hspace{1ex}SIMULATION TOOL DEVELOPMENT}\noindent 
The first (above shown) experimental plots were made in MATLAB~6 after
the results were obtained using ASPID -- Adaptive System Parameter 
Identification --  the soft package developed in Visual~C 6.0 (service 
pack~6) to carry on research on 
Inertial Navigation System error budget \cite{semoushin:2002}. However, to perform a wide scale investigation of all the proposed 
detection/selection algorithms that would correspond 
Figure~\ref{fig2}, a special software needs to be developed.

Such a tool, a very convenient one named LSCS (Linear Stochastic 
Control System), has been created by Mike Sunoplya, a fourth year Ulyanovsk State University student. It  allows us to 
observe the behavior of $S_k$ before and after the system failure in 
full accordance with Figure~\ref{fig2}. 
%(a screenshot is shown in Fig~\ref{screenshot}).

\begin{figure}%[h]
\centering
\includegraphics[scale=0.75]{0-2{0,2,8,9}.eps}
\caption{$\mathscr{S}_{0}
\stackrel{397}{\rightarrowtail} \mathscr{S}_{2} \bigm|
\mathscr{FP}=\{ \mathscr{F}_{0}, \mathscr{F}_{2}, \mathscr{F}_{8}, \mathscr{F}_{9} \}$. Sample size $N=15$ in Scheme~B, (\ref{schemeB})--(\ref{winnerB}). \textsl{Comment:} In this case the change 
$\mathscr{S}_{0}
\stackrel{t_{\mathrm{c}}}{\rightarrowtail} \mathscr{S}_{j}$ occurs at $t_{\mathrm{c}}=397$ for $j=2$.  Scheme~B, (\ref{schemeB})--(\ref{winnerB}), is used with each $S_k^{(j)}$, $j\in \{0,2,8,9 \}$, reset to zero every $N$ steps of algorithm (\ref{Skalg}). The zoomed version of this plot can be obtained, -- \textsl{cf.}~Figure~\ref{Case3zoom}.} \label{Case3}
\end{figure}

\begin{figure}%[h]
\centering
\includegraphics[scale=0.75]{0-2{0,2,8,9}zoom.eps}
\caption{Zoom of the previous plot of Figure~\ref{Case3} with $\mathscr{S}_{0}
\stackrel{397}{\rightarrowtail} \mathscr{S}_{2} \bigm|
\mathscr{FP}=\{ \mathscr{F}_{0}, \mathscr{F}_{2}, \mathscr{F}_{8}, \mathscr{F}_{9} \}$. Sample size $N=15$ in Scheme~B, (\ref{schemeB})--(\ref{winnerB}). \textsl{Comment:} It can be seen from here that before the change -- for $t_{i}<397$ -- only $S_{k}^{(0)}$ formed by (\ref{Skalg}) in  $\mathscr{F}_{0}$, remains in the vicinity of zero. After the change from $\mathscr{S}_{0}$ to $\mathscr{S}_{2}$ at $t_{i}=397$, behavior of the curves changes drastically: $S_{k}^{(2)}$ formed by (\ref{Skalg}) in  $\mathscr{F}_{2}$, enters the vicinity of zero while $S_{k}^{(0)}$  leaves it as continue to do so all the other curves, -- in this case $S_{k}^{(8)}$ and $S_{k}^{(9)}$ taken from  $\mathscr{F}_{8}$ and  $\mathscr{F}_{9}$, correspondingly. \textsl{Remark:} On insufficiently reduced scale, one curve can melt into another. In this plot, one can easily discriminate $S_{k}^{(8)}$ curve from $S_{k}^{(9)}$ curve between $t_{i}=413$ and $t_{i}=420$. } \label{Case3zoom}
\end{figure}

\begin{figure}%[h]
\centering
\includegraphics[scale=0.75]{0-9{0,2,8,9}.eps}
\caption{$\mathscr{S}_{0}
\stackrel{397}{\rightarrowtail} \mathscr{S}_{9} \bigm|
\mathscr{FP}=\{ \mathscr{F}_{0}, \mathscr{F}_{2}, \mathscr{F}_{8}, \mathscr{F}_{9} \}$. Sample size $N=15$ in Scheme~B, (\ref{schemeB})--(\ref{winnerB}). \textsl{Comment:} In this case the change 
$\mathscr{S}_{0}
\stackrel{t_{\mathrm{c}}}{\rightarrowtail} \mathscr{S}_{j}$ occurs at $t_{\mathrm{c}}=397$ for $j=9$.  Scheme~B, (\ref{schemeB})--(\ref{winnerB}), is used with each $S_k^{(j)}$, $j\in \{0,2,8,9 \}$, reset to zero every $N$ steps of algorithm (\ref{Skalg}). The zoomed version of any plot can be obtained. For this case, it is shown in~Figure~\ref{Case4zoomscreenshot} in the form of screenshot.} \label{Case4}
\end{figure}

\begin{figure}%[h]
\centering
\includegraphics[scale=0.733]{screen{0-9{0,2,8,9}zoomm}.eps}
\caption{Application window. Zoom of the previous plot of Figure~\ref{Case4} with $\mathscr{S}_{0}
\stackrel{397}{\rightarrowtail} \mathscr{S}_{9} \bigm|
\mathscr{FP}=\{ \mathscr{F}_{0}, \mathscr{F}_{2}, \mathscr{F}_{8}, \mathscr{F}_{9} \}$. } \label{Case4zoomscreenshot}
\end{figure}


\noindent
This LSCS application can work in 1 of~3 available modes:\\
\begin{tabular}{p{3cm} p{11.25cm}}
  &\\
  \textsl{normal mode} & application calculates and draws (in real time mode on the work area) the $S_k^{(0)}$ changes, $S_k^{(0)}$ given by the formula (\ref{Skalg}) introduced in Section~\ref{method}; only one filter $\mathscr{F}_{0}$ is available.\\
  &\\
  \textsl{mixed mode} & this mode  supports four different faulty filters $\mathscr{F}_{j}$, $j\in \{1,2,\dots,K\}$ versus one (nominal) filter $\mathscr{F}_{0}$; four curves of $S_k^{(j)}$ changes can be taken from $\mathscr{F}_{j}$ and observed together like in Figures~\ref{Case1}--\ref{Case2}; if we write
  $\mathscr{S}_{0}
\stackrel{2300}{\rightarrowtail} \mathscr{S}_{j} \bigm| j\in
\{ 2,6,3,7\}$, this means that one of four changes
$\mathscr{S}_{0}
\stackrel{t_{\mathrm{c}}}{\rightarrowtail} \mathscr{S}_{j}$ occurs at $t_{\mathrm{c}}=2300$  for $j=2,6,3 \mbox{ \textsl{or} }7$, -- this is the case of Figure~\ref{Case1},~\textsl{right}.  \\
  &\\
\end{tabular}

\noindent
\begin{tabular}{p{3cm} p{11.25cm}}
  &\\
  \textsl{renewal mode} & application calculates and draws (in real time mode on the work area) the $S_{k}^{(j)}$ changes for all $j \in\{0, \dots, K\}$; we identify this by the following notation: if we write 
  $\mathscr{S}_{0}
\stackrel{397}{\rightarrowtail} \mathscr{S}_{9} \bigm|
\mathscr{FP}=\{ \mathscr{F}_{0}, \mathscr{F}_{2}, \mathscr{F}_{8}, \mathscr{F}_{9} \}$ 
this means that the change 
$\mathscr{S}_{0}
\stackrel{t_{\mathrm{c}}}{\rightarrowtail} \mathscr{S}_{j}$ occurs at $t_{\mathrm{c}}=397$ for $j=9$; in this mode Scheme~B, (\ref{schemeB})--(\ref{winnerB}), is used with each $S_k^{(j)}$ reset to zero every $N$ steps of (\ref{Skalg}).\\
  &\\
\end{tabular}
LSCS was developed using the Visual C++.NET  development system and can work under operation system  Windows~98/2000/XP.\; It provides:\\[-7mm]:
\begin{itemize}
    \item on-line plotting\vspace{-0.25cm}
    \item multi-mode support\vspace{-0.25cm}
    \item manual zoom\vspace{-0.25cm}
    \item friendly work area that can be saved to bmp-file\vspace{-0.25cm}
    \item user controlled system matrices $\Phi_0$, $\Gamma_0$ and $H_0$ input and change at any assigned $t_{\mathrm{c}}$\vspace{-0.25cm}
    \item useful options to create the unique style of curves and legend at the work area\vspace{-0.25cm}
    \item handy interface\vspace{-0.25cm}
    \item and more...\vspace{-0.2cm}
\end{itemize}

With this tool switched onto the renewal mode, many
computational experiments were made in this work
(most intersing of them are presented in Figures~\ref{Case3}--\ref{Case4zoomscreenshot}) and many more
are reserved for further research. 

\section{\hspace{-2ex}{.}\hspace{1ex}CONCLUSIONS}\label{conclusion}
\noindent
The weighted squared  innovation method is formulated and tested
as applied to stochastic state estimation and control systems. The findings of this paper may be of considerable practical value thanks to high detecting capacity of the method at not-too-high computational cost. The understanding of the exact capacity for work of the method awaits further investigation. Further studies will probe practical aspects of change point detecting with respect to changes not only in observation matrix $H$ but in other system matrices, too.

\section*{ACKNOWLEDGEMENTS}
\noindent
The authors would like  to thank Prof. Dr.~A.~Butov of Ulyanovsk
State University, Russia, for his interest to this work and useful
discussions.

%\vspace{1.5cm}
\section*{APPENDIX A: A PROOF OF  THEOREM~\ref{theo1} FROM SECTION~\ref{method}}
\noindent
{\bf Lemma A.1.\;}
%\begin{lem}
$\xi \sim \mathcal{N}(0; 1)$ implies  $\E{\xi^{2k}}=(2k-1)!!$\ .
\pr{Proof}
\begin{eqnarray*}
\E{\xi^{2k}} &=& {\textstyle \frac{1}{\sqrt{2\pi}}}
\int\limits_{-\infty}^{+\infty}x^{2k}
\exp\left\{{\textstyle -\frac{x^{2}}{2}}\right\} \,
\ud x ={\textstyle  \frac{2^k}{\sqrt{\pi}}} \int\limits_{0}^{\infty}
t^{k -\frac{1}{2}}\exp\{-t\} \, \ud t \\
&=& {\textstyle \frac{2^k}{\sqrt{\pi}}}
\Gamma \left(k+{\textstyle \frac{1}{2}} \right)
   = {\textstyle \frac{2^k}{\sqrt{\pi}}} \cdot
\left( k- {\textstyle \frac{1}{2}} \right) \cdot
\Gamma\left(k-{\textstyle \frac{1}{2}}\right) \\
   &=&
\left( 2k-1 \right) \cdot \left( 2k-3 \right)
 \cdot \ldots  \cdot  3 \cdot 1  
\end{eqnarray*}\qed{$\Box$} 
%\end{lem}
\pr{Proof of Theorem~\ref{theo1}} Under  hypothesis $\Hyp_{0}$,
$\mu_{i} \sim \mathcal N(0; I)$ is a white Gaussian sequence where
$I$ is $m\times m$ identity matrix.  If $\mu_{ji}$ denotes the
$j$-th element of $\mu_{i}$, it follows that
\[
\xi\triangleq\sum\limits_{i=1}^{k}\sum\limits_{j=1}^{m} \mu_{ji}^{2}\sim
{\mathcal{F}}_{\chi ^{2}(mk)}(x)
\]
where ${\mathcal{F}}_{\chi ^{2}(mk)}(x)$  stands for $\chi^{2}$
distribution with $mk$ degrees of freedom. When $k\to \infty$,
$\xi$ is distributed normally with center $mk$ and variance $2mk$,
${\mathcal{F}}_{\chi ^{2}(mk)}(x) \leadsto \mathcal{N}(mk; 2mk)$.
From (\ref{mu}) we have
\[
S_{k}= \sqrt{mk/2}\left\{\frac{1}{mk}\xi-1 \right\}
\]
therefore $\mathcal{L}(S_{k}) \leadsto \mathcal{N}(0; 1)$.

Under hypothesis $\Hyp_{1}$ we have
\begin{eqnarray*}
C_{\mu}&\triangleq & \E{\mu_{i}\mu_{i}^{T}} = L_{0}^{-1}C_{1}L_{0}^{-T}=
I+L_{0}^{-1}(C_{1}-C_{0})L_{0}^{-T}=I+\Delta\\
m_{s}&\triangleq &\E{s_{k}}= \sqrt{m/2}\left[ (1/m)
\tr{C_{\mu}}-1 \right]=(2m)^{-1/2} \tr{\Delta}
\end{eqnarray*}
Let  $\left[ C_{\mu} \right]_{jk}$ be the $jk$-th element of
$C_{\mu}$, then $\left[ C_{\mu}
\right]_{jk}=\sigma_{j}\sigma_{k}\rho_{jk}$ where $\sigma_{j}^{2},
\, \sigma_{k}^{2}$ are variances of the $j$-th and  $k$-th
elements of vector $\mu_{i}$, correspondingly, and $\rho_{jk}$ is
the correlation between them. From (\ref{mu}) it follows that
\[
D_{s}= {\textstyle \frac{1}{2m}}\left[
\E{\left( \mu_{i}^{T} \mu_{i} \right)^{2}} -
\left( \E{\mu_{i}^{T} \mu_{i}} \right)^{2}
\right]
\]
Lemma  A.1 and straitforward calculations allow us to write
\[
\E{\mu_{i,k}^{4}}=3\sigma_{k}^{4}\ , \quad
\E{\mu_{ji}^{2}\mu_{ki}^{2}}= \sigma_{j}^{2}\sigma_{k}^{2}
\left( 1+2 \rho_{jk}^{2} \right)
\]
where $j$ and $k$ are indices of the corresponding elements of
$\mu_{i}$. It follows from here that
\begin{eqnarray*}
\E{\left( \mu_{i}^{T} \mu_{i} \right)^{2}} &=&
3\sum \limits_{k=1}^{m}\sigma_{k}^{4}+2
\sum \limits_{\substack{j, k=1\\ j<k}}^{m} \sigma_{j}^{2}\sigma_{k}^{2}
\left( 1+2 \rho_{jk}^{2} \right)\\
\left( \E{\mu_{i}^{T} \mu_{i}} \right)^{2} &=&
\sum \limits_{k=1}^{m}\sigma_{k}^{4}+2\sum \limits_{\substack{j, k=1\\
j<k}}^{m} \sigma_{j}^{2}\sigma_{k}^{2}
  \\
D_{s} &=& {\textstyle \frac{1}{m}}
\sum \limits_{j, k=1}^{m}\left[ C_{\mu} \right]_{jk}^{2}=
{\textstyle \frac{1}{m}}\Vert C_{\mu} \Vert ^{2}
\end{eqnarray*}
Because $\left[ C_{\mu} \right]_{jj}=1+\Delta_{jj}$ and
$\left[ C_{\mu} \right]_{\substack{jk\\j\ne k}}=\Delta_{jk}$ where
$\Delta_{jj}$ and $\Delta_{jk}$ denote the corresponding entries of
$\Delta$, then from the above expression we obtain
\[
D_{s} = 1+(2/m)\tr{\Delta}+(1/m)\Vert \Delta \Vert^{2}
\]
Based on definition $s_{i}$ (\ref{mu}), obtain $m_{S_{k}}$
(\ref{Delta}) and variance
\[
D_{S_{k}}= D_{s} + 2\sum \limits_{j=1}^{k-1}(1-j/k)\mathcal{K}_{ss}(j)
\]
Under approximation (\ref{corr}), we have
\[
\lim_{k\to \infty}D_{S_{k}}=D_{s}(1+r)/(1-r)
\]
Due  to the exponential correlated dependence (\ref{corr}), the
mixing condition in Birkhoff-Khinchin theorem
holds~\cite{yoshihawa:1992}, \cite{ibragimov:1965}, hence convergence
(\ref{H1}) and $P_{D}$ in (\ref{P}) are also true. 

\qed{$\Box$}


\bibliographystyle{apalike}
%\small
\bibliography{mybibase}


\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
