A Simple Squared Innovations Based Decision Generator for Detection/Selection Problems in Linear Stochastic Control Systems

I. V. Semoushin^{∗}

Department of Mathematical Cybernetics and Informatics Ulyanovsk State University, Ulyanovsk, Russia

Abstract: We consider the linear time-invariant state-space stochastic control system

x(t_{i+1}) = Φ_{θ}x(t_{i}) + Ψ_{θ}u(t_{i}) + w(t_{i}), x ∈ ℝ^{n}  (i)
y(t_{i}) = H_{θ}x(t_{i}) + v(t_{i}), y ∈ ℝ^{m}  (ii)
x̂_{0}(t_{i+1}^{−}) = Φ_{0}x̂_{0}(t_{i}^{+}) + Ψ_{0}u(t_{i}), x̂_{0} ∈ ℝ^{n}  (iii)
x̂_{0}(t_{i}^{+}) = x̂_{0}(t_{i}^{−}) + K_{0}ν(t_{i}), ν(t_{i}) = y(t_{i}) − H_{0}x̂_{0}(t_{i}^{−})  (iv)
u(t_{i}) = f_{R}[x̂_{0}(t_{i}^{+})] = −G_{0}^{⋆}x̂_{0}(t_{i}^{+}), u ∈ ℝ^{q}  (v)

where i ∈ ℤ, with a plant (i) and a sensor (ii), both parameterized by an uncertainty vector parameter θ, and a feedback controller FC, (iii)–(v), designed to cascade the steady-state Kalman filter (iii)–(iv) with a regulator (v) for a nominal value θ_{0} of θ. In (i)–(v), the initial state x(t_{−s}) is given at some t_{−s} ∈ ℝ with expectation E {∥x(t_{−s})∥^{2}} < ∞; 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_{θ} ≥ 0 and E \{v(t_{i})v(t_{i})^{T} \} = R_{θ} > 0 for all t_{i} ∈ ℝ, with \left[{w(t_{i}) \atop v(t_{i})}\right] orthogonal to x(t_{j}) and u(t_{j}) for all j ≥ i; u(t_{i}) is wide-sense stationary and E {∥u(t_{i})∥^{2}} < ∞ for all t_{i} ∈ ℝ.

Assuming that the true value θ^{†} of θ can change abruptly from θ_{0} to some θ_{1} at an instant t_{c} ∈ (t_{0}, t_{i}), we detect t_{c} using the nominal covariance C_{0} of ν(t_{i}) in (iv) and constructing a decision function of the cumulative sum form:

S_{k} = \sqrt{m/(2k)} ∑_{i=1}^{k} [ν^{T} (t_{i})C_{0}^{−1}ν(t_{i})/m − 1] .

Keywords: Empirical decision procedures; Hypothesis testing; Diagnostics; Stochastic systems and control; Stopping times.

Subject Classification: 62C12; 62F03; 62J20; 93E10; 60G40.

1. INTRODUCTION

The theory of abrupt changes is an intensively investigated topic in time series analysis and identification (Basseville and Nikiforov, 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 (Shewhart, 1931). As an alternative to them, Page set up a principle of cumulative sum control charts (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 (Johnson and Leone, 1962), (Johnson and Leone, 1963), (Johnson and Leone, 1964) and then in (Siegmund, 1985), (Kenett and Zacks, 1998) and (Ghosh and Sen, 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 (Basseville and Nikiforov, 1993).

∗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} .

In (Basseville and Nikiforov, 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) (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) ν(t_{i}) of an autoregression (AR) model with the aid of a test for its variance σ^{2}. The CUSUM

S_{k} = \frac{1}{\sqrt{k}} ∑_{i=1}^{k} s_{i} with s_{i} = \frac{1}{\sqrt{2}} (\frac{ν(t_{i})^{2}}{σ_{0}^{2}} − 1)  (1.1)

is asymptotically, when k goes to infinity, distributed as the standard Gaussian law: S_{k} ∼ N (0; 1) under the hypothesis H_{0} : σ^{2} = σ_{0}^{2}.

This algorithm, which is monitoring s_{i} termed the weighted squared innovations (WSI), was introduced in the AR case independently in (Jones et al., 1970), (Borodkin and Mottl, 1976), and (Segen and Sanderson, 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 (Mehra and Peschon, 1971). It is well known that the IS is white within a filter if the filter is Kalman, i. e., optimal one (Maybeck, 1979), and whiteness of the IS is only necessary, not sufficient condition for filter optimality (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.

2. PROBLEM FORMULATION

Let the closed-loop control system be parameterized by a vector θ ∈ ℝ^{p} and the available data z = \left[{y \atop u}\right] be a vector composed of two parts: the control input u ∈ ℝ^{q} and the measurement output y ∈ ℝ^{m}. The system is modeled for i ∈ ℤ by the equations

x(t_{i+1}) = Φ_{θ}x(t_{i}) + Ψ_{θ}u(t_{i}) + w(t_{i}), x ∈ ℝ^{n}  (2.1)
y(t_{i}) = H_{θ}x(t_{i}) + v(t_{i}), y ∈ ℝ^{m}  (2.2)
x̂_{0}(t_{i+1}^{−}) = Φ_{0}x̂_{0}(t_{i}^{+}) + Ψ_{0}u(t_{i}), x̂_{0} ∈ ℝ^{n}  (2.3)
x̂_{0}(t_{i}^{+}) = x̂_{0}(t_{i}^{−}) + K_{0}ν(t_{i}), ν(t_{i}) = y(t_{i}) − H_{0}x̂_{0}(t_{i}^{−})  (2.4)
u(t_{i}) = f_{R}[x̂_{0}(t_{i}^{+})] = −G_{0}^{⋆}x̂_{0}(t_{i}^{+}), u ∈ ℝ^{q}  (2.5)

with {w(·)}, {v(·)} being zero mean i.i.d. sequences of covariances Q_{θ} and R_{θ} respectively. The state difference equation (2.1) is propagated forward from the initial condition x(t_{−s}) of the finite expectation E {∥x(t_{−s})∥^{2}} < ∞. The initial state is placed at some t_{−s} ∈ ℝ where s > 0 determines what is termed settling time T_{s} = t_{0} − t_{−s} needed to think of all processes in (2.1)–(2.5) as wide-sense stationary for i ≥ 0. As usual, equation (2.1) represents a plant, equation (2.2) a sensor, and (2.3)–(2.5) a feedback. The feedback is composed of a Kalman-like filter (2.3)–(2.4) cascaded with a regulator (2.5). The regulator is described by a function f_{R}[·] of the measurement updated estimate x̂_{0}(t_{i}^{+}) or can be chosen according to the second equality in (2.5) with a matrix G_{0}^{⋆}.

Matrices of the given system (2.1)–(2.2) are assumed to be known as Φ_{0}, Ψ_{0}, Q_{0}, H_{0} and R_{0} for a nominal mode, i. e., for a nominal value θ_{0} of the uncertainty parameter θ. For this mode, to guarantee existence of the steady-state filter with Kalman gain K_{0}, we assume that the pair (Φ_{0}, Q_{0}^{1/2}) is stabilizable, the pair (Φ_{0}, H_{0}) is observable, and the pair (Φ_{0}, Ψ_{0}) is controllable. Matrix G_{0}^{⋆} can be taken then as a given one or designed to be LQG optimal for the nominal mode of operation (Maybeck, 1982), (Caines, 1988), (Mosca, 1995).

Parameter θ is subject to an abrupt change at an unknown time point t_{c} ∈ (t_{0}, t_{i}). This can be viewed as a switch of θ from θ_{0} to some other unknown value θ_{1}; to detect the change point, a decision generator, 𝒟𝒢, is needed (Figure 1). Formalizing the problem of 𝒟𝒢 synthesis, we consider

\begin{matrix}x_{j}(t_{i+1}) & = Φ_{j}x_{j}(t_{i}) + Ψ_{j}u(t_{i}) + w_{j}(t_{i}) \\ y_{j}(t_{i}) & = H_{j}x_{j}(t_{i}) + v_{j}(t_{i}) \end{matrix} \} 𝒮_{0} (for j = 0) or 𝒮_{1} (for j = 1)  (2.6)

as the generalized description of the two systems: 𝒮_{0} the system with θ = θ_{0}; 𝒮_{1} the system with θ = θ_{1}. Filter (2.3)–(2.4), which is denoted by ℱ_{0} in Figure 1, is designed as the Kalman filter for 𝒮_{0}, i. e., satisfying equations

\begin{matrix}K_{0} & = P̃_{0}H_{0}^{T} C_{0}^{−1} , & C_{0} & = H_{0}P̃_{0}H_{0}^{T} + R_{0} \\ P̂_{0} & = P̃_{0} − P̃_{0}H_{0}^{T} C_{0}^{−1}H_{0}P̃_{0} , & P̃_{0} & = Φ_{0}P̂_{0}Φ_{0}^{T} + Q_{0} \end{matrix} \} filterℱ_{0} for system 𝒮_{0}  (2.7)

and processing data

y(t_{i}) = \{ \begin{matrix}y_{0}(t_{i}) & if & t_{i} < t_{c} & (data from 𝒮_{0}) \\ y_{1}(t_{i}) & if & t_{i} ≥ t_{c} & (data from 𝒮_{1}) \end{matrix}  (2.8)

The problem is to detect the change point t_{c} with a reasonable delay using a suitable decision rule d_{0}(t_{k}) ∈ {0, 1} according to Figure 1.

3. WSI BASED DECISION GENERATING METHOD

Many systems, not only (2.1)–(2.5), have the IS, hence the WSI of s_{i} type in (1.1) 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 K_{ss}(j) of s_{i} after change equal to the exponential function:

K_{ss}(j) = D_{s}r^{|j|} , j ∈ ℤ , r = exp{−τ_{s}/τ_{c}} , D_{s} ≜ E {(s_{i} − E {s_{i}})^{2}}  (3.1)

•  𝒮_{0} w_{0}(t_{i})  v_{0}(t_{i}) 𝒮_{1} w_{1}(t_{i})  v_{1}(t_{i})  ↶ • t_{c}y(t_{i}) • •  • ℱ_{0} f_{R}(·) u(t_{i})  𝒟𝒢d_{0}(t_{k}) ν(t_{1}^{k}) x̂_{0}(t_{i}^{+}) •  q^{−1}  u(t_{i−1})

Figure 1. The general framework for the problem and solution. Legend: 𝒮_{0} is the system with θ = θ_{0} while 𝒮_{1} is the system with θ = θ_{1}; ℱ_{0} is the Kalman filter for 𝒮_{0}; and 𝒟𝒢 is a decision generator; q^{−1} denotes a one sample memory (unit delay); ν(t_{1}^{k}) ≜ col[][ν(t_{1}), ν(t_{2}), ⋯, ν(t_{k})] .

where τ_{s} is a sampling time interval and τ_{c} a correlation interval of s_{i}. Under this assumption, we state the WSI based method for detecting changes as follows.

Theorem 3.1 Let a dynamical system with the IS {ν(t_{i})}, ν(t_{i}) ∈ ℝ^{m}, have two possible modes of operation termed “normal functioning” (hypothesis H_{0}) and “a fault” (hypothesis H_{1}). Under both hypotheses, {ν(t_{i})} is a Gaussian sequence with E {ν(t_{i})} = 0. Under hypothesis H_{0}, {ν(t_{i})} is an independent sequence with the known covariance C_{0} = L_{0}L_{0}^{T} of each element ν(t_{i}) where L_{0} is the square root of C_{0} determined, for instance, by Cholesky decomposition. Under hypothesis H_{1}, {ν(t_{i})} is a correlated sequence with an unknown covariance C_{1} ≠ C_{0} of each element ν(t_{i}). Define

μ_{i} ≜ L_{0}^{−1}ν(t_{i}) , s_{i} ≜ \sqrt{\frac{m}{2}} (\frac{1}{m}μ_{i}^{T} μ_{i} − 1) , S_{k} ≜ \frac{1}{\sqrt{k}} ∑_{i=1}^{k} s_{i}  (3.2)

and suppose that under hypothesis H_{1}, the covariance kernel function of this s_{i} has been modeled (maybe rather approximately) by expression (3.1).

Then asymptotically, when k goes to infinity, the following properties hold:

1. The probability laws L(·) of S_{k} under these hypotheses satisfy^{†}

H_{0} : L(S_{k}) ⇝ N (0, 1)  (3.3)
H_{1} : L(S_{k}) ⇝ N (m_{S_{k}}, D_{S_{k}})  (3.4)

where N (0; 1) stands for the Gaussian (normal) distribution with zero mean and unit variance, N (m_{S_{k}}, D_{S_{k}}) is the normal distribution whose mean m_{S_{k}} and variance D_{S_{k}} = σ_{S_{k}}^{2} are given by the below formulae

†The notation ⇝ corresponds to the weak convergence (Roussas, 1972).

\begin{matrix}m_{S_{k}} ≜ & E \{S_{k}\} = m_{s}\sqrt{k} , m_{s} ≜ E \{s_{i}\} = (1/\sqrt{2m}) tr\{Δ\} \\ D_{S_{k}} ≜ & E \{(S_{k} − m_{S_{k}})^{2}\} = D_{s} \frac{1 + r}{1 − r} , D_{s} = 1 + \frac{2}{m} tr\{Δ\} +\frac{1}{m}∥Δ∥^{2} \\ Δ = & L_{0}^{−1}(C_{1} − C_{0})L_{0}^{−T} , ∥Δ∥^{2} ≜ ∑_{i,j=1}^{m} |Δ_{ij}|^{2} \end{matrix} \}  (3.5)

2. If S_{k} of (3.2) is used as a decision function at a time instant k in the decision rule

|S_{k}| \begin{matrix}H_{1} \\ ≷ \\ H_{0} \end{matrix} h  (3.6)

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 L(S_{k}) in (3.3), (3.4) are taken to be normal):

\begin{matrix}P_{F} & ≃ 1 − ϕ(h) \\ P_{D} & ≃ 1 − \frac{1}{2} [ϕ (\frac{m_{S_{k}} + h}{σ_{S_{k}}}) − ϕ (\frac{m_{S_{k}} − h}{σ_{S_{k}}})] > \frac{1}{2} [1 + ϕ (\frac{m_{S_{k}} − h}{σ_{S_{k}}})] \end{matrix} \}  (3.7)

through the standard probability integral

ϕ(x) ≜ \frac{2}{\sqrt{2π}} ∫_{0}^{x} exp (−t^{2}/2) dt

3. If the threshold h in (3.6) is chosen as h = ϕ^{−1}(1 − α) to guarantee a given level α of P_{F} in (3.7), P_{F} = α, then hypothesis H_{1} is detected with probability P_{D} = 1 − β (where 0 < β < 1/2) not later than after

k_{a}^{⋆} ≃ [ϕ^{−1}(1 − α) + σ_{S_{k}}ϕ^{−1}(1 − 2β)]^{2} /m_{s}^{2}  (3.8)

disctere time instants where x = ϕ^{−1}(y) is the solution for ϕ(x) = y.

The proof of Theorem 3.1 is placed in Appendix A.

Clearly sequence S_{k} can be used for any system having such an innovation ν(t_{i}) to detect the changes marked by matrix Δ with tr{Δ} ≠ 0, with the aid of rule (3.6). In so doing, it is important to keep in mind the following merits of the rule:

1. Rule (3.6) contains the same quadratic forms as does the log-likelihood function

ln p (ν(t_{1}^{k}) | H_{0}) = −\frac{km}{2} ln(2π) − \frac{k}{2} ln|C_{0}| − \frac{1}{2} ∑_{i=1}^{k} ν^{T} (t_{i})C_{0}^{−1}ν(t_{i})

with ν(t_{1}^{k}) ≜ col[][ν(t_{1}), ν(t_{2}), ⋯, ν(t_{k})] because

S_{k} = \sqrt{\frac{m}{2k}} ∑_{i=1}^{k} (\frac{1}{m}ν^{T} (t_{i})C_{0}^{−1}ν(t_{i}) − 1)  (3.9)

2. The sufficient statistic ℓ(·) ≜ ℓ (ν(t_{1}^{k})) for detecting changes in the system parameters subject to the assumptions of Theorem 3.1, is (Van-Trees, 1968, Sec. 2.6.2)

ℓ(·) = ∑_{i=1}^{k} ν^{T} (t_{i})C_{0}^{−1}ν(t_{i}) − ν^{T} (t_{1}^{k})C_{1}^{−1}(k)ν(t_{1}^{k}) , C_{1}(k) ≜ E {ν(t_{1}^{k})ν^{T} (t_{1}^{k}) | H_{1}}

The evident theoretical drawback of the method is that the decision function S_{k}, (3.2) or equivalently (3.9), 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 (Semoushin, 1979).

3. The type of parametric changes detectable by rule (3.6) is restricted to those for which tr{Δ} ≠ 0.

4. To strictly satisfy the condition tr{Δ} = 0 in most real situations is hardly probable and because of this the restriction tr{Δ} ≠ 0 should not be considered as critical from the practical point of view.

Example 3.1. Let m be equal 3 and

C_{0} = [ \begin{matrix}1 & 2 & 3 \\ 2 & 8 & 2 \\ 3 & 2 & 14 \end{matrix} ] , C_{1} = [ \begin{matrix}2.0 & 1.8 & 2.4 \\ 1.8 & 9.0 & 1.0 \\ 2.4 & 1.0 & 13.0 \end{matrix} ]

Then

L_{0} = [ \begin{matrix}1 & 0 & 0 \\ 2 & 2 & 0 \\ 3 & −2 & 1 \end{matrix} ] ,  Δ = [ \begin{matrix}1.0 & −1.1 & −5.8 \\ −1.1 & 1.45 & 6.3 \\ −5.8 & 6.3 & 31.0 \end{matrix} ] ,  m_{s} ≈ 13.7 ,  D_{s} ≈ 395.

Let α = 0.005 and β = 0.005. We obtain, h = ϕ^{−1}(1 − α) = 2.8 and ϕ^{−1}(1 − 2β) = 2.6 .

If r = 0.5, then k_{a}^{⋆} ≃ 46 and if r = 0.95, then k_{a}^{⋆} ≃ 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 α and β of the decision rule (3.6), it is possible to numerically predict the lapse of time^{‡}^{)} k_{a}^{⋆} after which the fault will be detected with confidence not worse than P_{D} = 1 − β.

For the above example, the condition of non-observable faults is readily obtained as the equation

tr\{Δ\} = 27a_{11} − 11a_{12} − 10a_{13} + \frac{5}{4}a_{22} + 2a_{23} + a_{33} = 0  (3.10)

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 (3.10), for example

C_{1} = [ \begin{matrix}1 & 2 & 3 \\ 2 & 8 & 1 \\ 3 & 1 & 16 \end{matrix} ] , A ≜ C_{1} − C_{0} = [a_{ij}] = [ \begin{matrix}0 & 0 & 0 \\ 0 & 0 & −1 \\ 0 & −1 & 2 \end{matrix} ]

This simple example makes it apparent that only certain of the faults specially chosen to fit the equation tr{Δ} = 0 to a T cannot be detected by the method.

Theorem 3.1 assumes that changes in θ do not create a bias (i.e., non-zero expectation) in the innovation sequence {ν(t_{i})} :

H_{0} : ν(t_{i}) ∼ N (0, C_{0})  and H_{1} : ν(t_{i}) ∼ N (0, C_{1})  (3.11)

‡^{)} In this example, the delay k_{a}^{⋆} should be increased by k_{norm} ⪆ 30 needed to take distributions in (3.3)–(3.4) to be approximately normal, from the engineering point of view.

If this is not the case, one should turn from (3.11) to

H_{0} : ν(t_{i}) ∼ N (0, C_{0})  and H_{1} : ν(t_{i}) ∼ N (m_{1}, C_{1})  (3.12)

with some m_{1} ≠ 0. This entails the only modification in Theorem 3.1, namely, in Δ of formulae (3.5):

Δ = L_{0}^{−1}(C_{1} + m_{1}m_{1}^{T} − C_{0})L_{0}^{−T}  instead of Δ = L_{0}^{−1}(C_{1} − C_{0})L_{0}^{−T}  (3.13)

Thus, for appreciable quantity of faults, the WSI based method boils down to a straightforward cumulation of S_{k}, (3.2), through the sequential algorithm

S_{k} = S_{k−1}\sqrt{1 − 1/k} + s_{k}/\sqrt{k} , k ∈ ℕ , S_{0} = 0  (3.14)

followed by the rule (3.6). When testing H_{1} against H_{0}, let us denote S_{k} in (3.14) by S_{k}^{(0)} thus emphasizing by superscript^{(0)} that the main hypothesis is H_{0}. Then we have two schemes with the stopping rule t_{a} (the alarm time) for the alternative H_{1} detection:

(A) decision is made at any current time t_{k}:

d_{0}(t_{k}) = \{ \begin{matrix}0 & if |S_{k}^{(0)}| < h ; H_{0} is chosen \\ 1 & if |S_{k}^{(0)}| ≥ h ; H_{1} is chosen \end{matrix}\} , t_{a} = min \{t_{k} : d_{0}(t_{k}) = 1\}

(B) decision is made at the end of a sample number l = 1, 2, …, L each of size N:

d_{0}(l) = \{ \begin{matrix}0 & if ∀ k = 1, 2, …, N : |S_{N(l−1)+k}^{(0)}| < h ; H_{0} is chosen \\ 1 & if ∃ k = 1, 2, …, N : |S_{N(l−1)+k}^{(0)}| ≥ h ; H_{1} is chosen \end{matrix}

t_{a} = τ_{s} [N(l − 1) + min {k : |S_{N(l−1)+k}^{(0)}| ≥ h}] =^{r} τ_{s}N min {l : d_{0}(l) = 1}

where =^{r} denotes the rounded up equality.

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 (2.1)–(2.5) type, in which the whitening filter exists in the Kalman filter form (2.3)–(2.4) with the IS {ν(t_{i})}.

4. WSI BASED FEEDBACK FILTER TESTING AND SELECTION

From the innovation process theory, two facts are known (Martin and Stubberud, 1976):

(i) Case (3.11) holds if changes happen only in covariances Q_{θ} and R_{θ} .
(ii) The only cause for Case (3.12) are changes in matrices Φ_{θ}, Ψ_{θ} and/or H_{θ} .

If Case (3.12) takes place, (3.13) with m_{1} calculated by the general procedure described, for example in (Martin and Stubberud, 1976), should be used.

Example 4.1. Consider the simple case when n = m = 1, u(t_{i}) ≡ 0, Q_{1} = Q_{0}, R_{1} = R_{0}, H_{1} = H_{0} = 1 and Φ_{1} ≠ Φ_{0} (both |Φ_{j}| < 1, j = 0, 1). At t_{i} ≫ t_{c} (i.e., under H_{1} when the change transient response effects are vanished) from (2.6) (for j = 0), (2.7), (2.8) and (2.3)–(2.4), we have

ν(t_{i}) = x_{1}(t_{i}) + v_{0}(t_{i}) − x̂_{0}(t_{i}) = e(t_{i}) + v_{0}(t_{i}) , e(t_{i}) ≜ x_{1}(t_{i}) − x̂_{0}(t_{i})  (4.1)

Denoting

\begin{matrix}m_{1}(t_{i}) & ≜ E \{ν(t_{i}) | H_{1}\} , & m_{e}(t_{i}^{−}) & ≜ E \{e(t_{i})\} \\ m_{x_{1}}(t_{i}) & ≜ E \{x_{1}(t_{i})\} , & mx̂_{0}(t_{i}^{±}) & ≜ E \{x̂_{0}(t_{i}^{±})\} \end{matrix} \}  (4.2)

we obtain from (4.1), (2.6) (for j = 1) and (2.3)

\begin{matrix}m_{1}(t_{i}) & = m_{e}(t_{i}^{−}) , & m_{e}(t_{i}^{−}) & = m_{x_{1}}(t_{i}) − mx̂_{0}(t_{i}^{−}) \\ m_{x_{1}}(t_{i}) & = Φ_{1}m_{x_{1}}(t_{i−1}) , & mx̂_{0}(t_{i}^{−}) & = Φ_{0}mx̂_{0}(t_{i−1}^{+}) \end{matrix} \}  (4.3)

Taking expectation of (2.4) yields

mx̂_{0}(t_{i}^{+}) = (1 − K_{0})mx̂_{0}(t_{i}^{−}) + K_{0}m_{x_{1}}(t_{i})

Taking here a step backwards and then substituting this into the last formula of (4.3) give

mx̂_{0}(t_{i}^{−}) = Φ_{0} [(1 − K_{0})mx̂_{0}(t_{i−1}^{−}) + K_{0}m_{x_{1}}(t_{i−1})]

Then we find consequently using (4.2) and (4.3)

m_{e}(t_{i}^{−}) = m_{x_{1}}(t_{i}) − Φ_{0} [(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})]
= m_{x_{1}}(t_{i}) − Φ_{0} [m_{x_{1}}(t_{i−1}) − (1 − K_{0})m_{e}(t_{i−1}^{−})]
= (Φ_{1} − Φ_{0})m_{x_{1}}(t_{i−1}) + Φ_{0}(1 − K_{0})m_{e}(t_{i−1}^{−})
m_{1}(t_{i}) = Φ_{0}(1 − K_{0})m_{1}(t_{i−1}) + (Φ_{1} − Φ_{0})m_{x_{1}}(t_{i−1})

As 0 < Φ_{0} < 1 and 0 < Φ_{0}(1 − K_{0}) < 1, there exist the steady-state values

m_{1} ≜ lim_{t_{i}→∞}{m_{1}(t_{i})} = lim_{t_{i}→∞}{m_{1}(t_{i−1})} , m_{x_{1}} ≜ lim_{t_{i}→∞}{m_{x_{1}}(t_{i−1})}

and finally

m_{1} = \frac{Φ_{1} − Φ_{0}}{1 − Φ_{0}(1 − K_{0})}m_{x_{1}}

Remark 4.1. By virtue of 0 < Φ_{0} < 1, m_{x_{1}} = 0. It is also worthy of note that although tr{Δ} = 0 is not a necessary or sufficient condition for filter optimality, the quantity of faults satisfying tr{Δ} = 0 is relatively thin.

As for the case of two contesting hypotheses H_{0} and H_{1} in Section 3, we can apply one of two schemes, (A) or (B), of decision making. However as differentiated from that case, here every singly taken H_{j} (j = 0, 1, …, K) contests individually against any other H_{i} (i = 0, 1, …, K; i ≠ j) using its own decision function S_{k}^{(j)} that equals to S_{k} and is formed according to (3.14) supplied with the superscript^{(j)}. Again, we have two schemes:

(A) decision is made at any current time t_{k}. For any particular H_{j}, j ∈ {0, 1, …, K}, the particular decision rule is

d_{j}(t_{k}) = \{ \begin{matrix}0 & if |S_{k}^{(j)}| < h ; H_{j} is supported \\ 1 & if |S_{k}^{(j)}| ≥ h ; H_{j} is rejected \end{matrix}

The over-all alarm time is

t_{a} = min {t_{k} : (∃κ : d_{κ}(t_{k}) = 0 & ∀_{j≠κ} d_{j}(t_{k}) = 1)}

𝒮𝒫  ℱ𝒫  ℱ𝒮 ℱℰ • 𝒮ℱ u(t_{i})  y(t_{i})  (d_{0}, ⋯, d_{K})  κ (ℱ_{0}, ⋯, ℱ_{K}) ℱ_{κ}  κ = {j ∈ {0, …, K} : d_{j}(t_{a}) = 0}

Figure 2. The best filter selection. Legend: 𝒮𝒫 stands for system population; ℱ𝒫 stands for filter population, ℱ𝒮 for filter selection, ℱℰ for feedback effector, and 𝒮ℱ stands for the system feedback. The stopping rule is given by the alarm time t_{a} = min{t_{k} : (∃κ : d_{κ}(t_{k}) = 0 & ∀_{j≠κ} d_{j}(t_{k}) = 1)} . At this time, ℱℰ

uploads the selected ℱ_{κ} into 𝒮ℱ if by then it is not in there.

with the over-all winner selection rule:

κ = {j ∈ {0, 1, …, K} : d_{j}(t_{a}) = 0} ; H_{κ} is chosen

(B) decision is made at the end of a sample number l = 1, 2, …, L each of size N. For any particular H_{j}, j ∈ {0, 1, …, K}, the particular decision rule is

d_{j}(l) = \{ \begin{matrix}0 & if ∀ k = 1, 2, …, N : |S_{N(l−1)+k}^{(j)}| < h ; H_{j} is supported \\ 1 & if ∃ k = 1, 2, …, N : |S_{N(l−1)+k}^{(j)}| ≥ h ; H_{j} is rejected \end{matrix}  (4.4)

t_{a} =^{r} τ_{s}N min {l : (∃κ : d_{κ}(l) = 0 & ∀_{j≠κ} d_{j}(l) = 1)}  (4.5)

with the over-all winner selection rule:

κ = {j ∈ {0, 1, …, K} : d_{j}(l) = 0} ; H_{κ} is chosen  (4.6)

This method is realized in a framework corresponding to Figure 2.

5. SOME COMPUTATIONAL EXPERIMENT RESULTS

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 (Gaines, 1971). In compliance with (Gaines, 1971), consider equations

\begin{matrix}x_{j}(t_{i+1}) & = Φ_{j}x_{j}(t_{i}) + Γ_{j}w_{j}(t_{i}) \\ y_{j}(t_{i}) & = H_{j}x_{j}(t_{i}) + v_{j}(t_{i}) \end{matrix} \} 𝒮_{0} (for j = 0) or 𝒮_{j} (for j = 1, ⋯, K)  (5.1)

that have the form of (2.6) 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 𝒮_{j}: 𝒮_{0} is the fault free system and every 𝒮_{j} (for j = 1, ⋯, K) is a faulty system. System 𝒮_{0} corresponds to

Φ_{0} = [ \begin{matrix}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{matrix} ] , Γ_{0} = [ \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 24.64 & 0 & 0 \\ 0 & 0.835 & 0 \\ 0 & 0 & 1.83 \end{matrix} ] (5.2)

Q_{0} = [ \begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} ] ,  H_{0} = [ \begin{matrix}1 & 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 & 0 \end{matrix} ] ,  R_{0} = [ \begin{matrix}1 & 0 \\ 0 & 1 \end{matrix} ] (5.3)

Disengaging ourselves at this stage from a physical interpretation of the given characteristics (5.2)–(5.3), we shall restrict our consideration in this section by introducing the set of faulty systems through the following formal relations:

\begin{matrix}Φ_{j} & = Φ_{0} , Γ_{j} = Γ_{0} , Q_{j} = Q_{0} , R_{j} = R_{0} \\ H_{j} & = [ \begin{matrix}1 − e & 0 & 0 & 0 & 1 − f \\ 0 & 1 − g & 0 & 1 − h & 0 \end{matrix} ] , \{efgh\} ≜ j_{[2]} \end{matrix} \} 𝒮_{j} (j = 1, ⋯, 15)

where {efgh} stands for the binary code j_{[2]} of the fault number j. For example, given j = 2 or j = 9, then {efgh} = {0010} or, correspondingly {1001}, and in these cases

H_{2} = [ \begin{matrix}1 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \end{matrix} ] , H_{9} = [ \begin{matrix}0 & 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 & 0 \end{matrix} ]

The change 𝒮_{0} ↣^{t_{c}} 𝒮_{j} (j = 1, ⋯, 15) occurs at t_{c} = 2300, however filter ℱ_{0} “does not know” about that and continues to satisfy equations (2.7). Experimental data for all K = 15 faulty system cases are shown in Figures 3–4 together with the fault free system case, when {efgh} = {0000}. If the value of threshold h in (3.6) is set to 3, then probability of false alarm, P_{F}, in the case of 𝒮_{0} will be not greater than 0.3% (due to property (3.3)), 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 3–4, 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 (2.1)–(2.2), 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.

6. SIMULATION TOOL DEVELOPMENT

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 (Semoushin and Polosenko, 2002). However, to perform a wide scale investigation of all the proposed detection/selection algorithms that would correspond Figure 2, 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 2.

Figure 3. Behavior of S_{k} before and after change at t_{c} = 2300; left: j = 0, 4, 1, 5, right: j = 2, 6, 3, 7.

Figure 4. Behavior of S_{k} before and after change at t_{c} = 2300; left: j = 8, 12, 9, 13, right: j = 10, 14, 11, 15.

This LSCS application can work in 1 of 3 available modes:

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 (3.14) introduced in Section 3; only one filter ℱ_{0} is available.

mixed mode this mode supports four different faulty filters ℱ_{j}, j ∈ {1, 2, …, K} versus one (nominal) filter ℱ_{0}; four curves of S_{k}^{(j)} changes can be taken from ℱ_{j} and

observed together like in Figures 3–4; if we write 𝒮_{0} ↣^{2300} 𝒮_{j} | j ∈ {2, 6, 3, 7},

this means that one of four changes 𝒮_{0} ↣^{t_{c}} 𝒮_{j} occurs at t_{c} = 2300 for j = 2, 6, 3 or 7, – this is the case of Figure 3, right.

Figure 5. 𝒮_{0} ↣^{397} 𝒮_{2} | ℱ𝒫 = {ℱ_{0}, ℱ_{2}, ℱ_{8}, ℱ_{9}}. Sample size N = 15 in Scheme B, (4.4)–(4.6). Comment:

In this case the change 𝒮_{0} ↣^{t_{c}} 𝒮_{j} occurs at t_{c} = 397 for j = 2. Scheme B, (4.4)–(4.6), is used with each S_{k}^{(j)}, j ∈ {0, 2, 8, 9}, reset to zero every N steps of algorithm (3.14). The zoomed version of this plot can be obtained, – cf. Figure 6.

Figure 6. Zoom of the previous plot of Figure 5 with 𝒮_{0} ↣^{397} 𝒮_{2} | ℱ𝒫 = {ℱ_{0}, ℱ_{2}, ℱ_{8}, ℱ_{9}}. Sample size N = 15 in Scheme B, (4.4)–(4.6). Comment: It can be seen from here that before the change – for t_{i} < 397 – only S_{k}^{(0)} formed by (3.14) in ℱ_{0}, remains in the vicinity of zero. After the change from 𝒮_{0} to 𝒮_{2} at t_{i} = 397, behavior of the curves changes drastically: S_{k}^{(2)} formed by (3.14) in ℱ_{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 ℱ_{8} and ℱ_{9}, correspondingly. 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.

Figure 7. 𝒮_{0} ↣^{397} 𝒮_{9} | ℱ𝒫 = {ℱ_{0}, ℱ_{2}, ℱ_{8}, ℱ_{9}}. Sample size N = 15 in Scheme B, (4.4)–(4.6). Comment:

In this case the change 𝒮_{0} ↣^{t_{c}} 𝒮_{j} occurs at t_{c} = 397 for j = 9. Scheme B, (4.4)–(4.6), is used with each S_{k}^{(j)}, j ∈ {0, 2, 8, 9}, reset to zero every N steps of algorithm (3.14). The zoomed version of any plot can be obtained. For this case, it is shown in Figure 8 in the form of screenshot.

Figure 8. Application window. Zoom of the previous plot of Figure 7 with 𝒮_{0} ↣^{397} 𝒮_{9} | ℱ𝒫 = {ℱ_{0}, ℱ_{2}, ℱ_{8}, ℱ_{9}}. renewal mode application calculates and draws (in real time mode on the work area) the S_{k}^{(j)} changes for all j ∈ {0, …, K}; we identify this by the following notation: if we

write 𝒮_{0} ↣^{397} 𝒮_{9} | ℱ𝒫 = {ℱ_{0}, ℱ_{2}, ℱ_{8}, ℱ_{9}} this means that the change 𝒮_{0} ↣^{t_{c}} 𝒮_{j}

occurs at t_{c} = 397 for j = 9; in this mode Scheme B, (4.4)–(4.6), is used with each S_{k}^{(j)} reset to zero every N steps of (3.14).

LSCS was developed using the Visual C++.NET development system and can work under operation system Windows 98/2000/XP. It provides:

:

• on-line plotting

• multi-mode support

• manual zoom

• friendly work area that can be saved to bmp-file

• user controlled system matrices Φ_{0}, Γ_{0} and H_{0} input and change at any assigned t_{c}

• useful options to create the unique style of curves and legend at the work area

• handy interface

• and more...

With this tool switched onto the renewal mode, many computational experiments were made in this work (most intersing of them are presented in Figures 5–8) and many more are reserved for further research.

7. CONCLUSIONS

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.

ACKNOWLEDGEMENTS

The authors would like to thank Prof. Dr. A. Butov of Ulyanovsk State University, Russia, for his interest to this work and useful discussions.

APPENDIX A: A PROOF OF THEOREM 3.1 FROM SECTION 3

Lemma A.1. ξ ∼ N (0; 1) implies E {ξ^{2k}} = (2k − 1)!! .

Proof.

E \{ξ^{2k}\} = \frac{1}{\sqrt{2π}} ∫_{−∞}^{+∞}x^{2k} exp \{−\frac{x^{2}}{2} \} dx = \frac{2^{k}}{\sqrt{π}} ∫_{0}^{∞}t^{k− \frac{1}{2}} exp\{−t\} dt

= \frac{2^{k}}{\sqrt{π}} Γ (k + \frac{1}{2}) = \frac{2^{k}}{\sqrt{π}} · (k − \frac{1}{2}) · Γ (k − \frac{1}{2})
= (2k − 1) · (2k − 3) · . . . · 3 · 1
□

Proof of Theorem 3.1. Under hypothesis H_{0}, μ_{i} ∼ N (0; I) is a white Gaussian sequence where I is m × m identity matrix. If μ_{ji} denotes the j-th element of μ_{i}, it follows that

ξ ≜ ∑_{i=1}^{k} ∑_{j=1}^{m} μ_{ji}^{2} ∼ F_{χ^{2}(mk)}(x)

where F_{χ^{2}(mk)}(x) stands for χ^{2} distribution with mk degrees of freedom. When k → ∞, ξ is distributed normally with center mk and variance 2mk, F_{χ^{2}(mk)}(x) ⇝ N (mk; 2mk). From (3.2) we have

S_{k} = \sqrt{mk/2} \{\frac{1}{mk} ξ − 1\}

therefore L(S_{k}) ⇝ N (0; 1).

Under hypothesis H_{1} we have

C_{μ} ≜ E {μ_{i}μ_{i}^{T} } = L_{0}^{−1}C_{1}L_{0}^{−T} = I + L_{0}^{−1}(C_{1} − C_{0})L_{0}^{−T} = I + Δ
m_{s} ≜ E \{s_{k}\} = \sqrt{m/2} [(1/m) tr\{C_{μ}\} −1] = (2m)^{−1/2} tr\{Δ\}

Let [C_{μ}]_{jk} be the jk-th element of C_{μ}, then [C_{μ}]_{jk} = σ_{j}σ_{k}ρ_{jk} where σ_{j}^{2}, σ_{k}^{2} are variances of the j-th and k-th elements of vector μ_{i}, correspondingly, and ρ_{jk} is the correlation between them. From (3.2) it follows that

D_{s} = \frac{1}{2m} [E \{(μ_{i}^{T} μ_{i})^{2}\} − (E \{μ_{i}^{T} μ_{i}\})^{2}]

Lemma A.1 and straitforward calculations allow us to write

E {μ_{i,k}^{4}} = 3σ_{k}^{4} , E {μ_{ji}^{2}μ_{ki}^{2}} = σ_{j}^{2}σ_{k}^{2} (1 + 2ρ_{jk}^{2})

where j and k are indices of the corresponding elements of μ_{i}. It follows from here that

E \{(μ_{i}^{T} μ_{i})^{2}\} = 3 ∑_{k=1}^{m} σ_{k}^{4} + 2 ∑_{\begin{matrix}j,k=1 \\ j<k \end{matrix}}^{m} σ_{j}^{2}σ_{k}^{2} (1 + 2ρ_{jk}^{2})

(E \{μ_{i}^{T} μ_{i}\})^{2} = ∑_{k=1}^{m} σ_{k}^{4} + 2 ∑_{\begin{matrix}j,k=1 \\ j<k \end{matrix}}^{m} σ_{j}^{2}σ_{k}^{2}

D_{s} = \frac{1}{m} ∑_{j,k=1}^{m} [C_{μ}]_{jk}^{2} = \frac{1}{m}∥C_{μ}∥^{2}

Because [C_{μ}]_{jj} = 1 + Δ_{jj} and [C_{μ}]\begin{matrix}jk \\ j≠k \end{matrix} = Δ_{jk} where Δ_{jj} and Δ_{jk} denote the corresponding entries of

Δ, then from the above expression we obtain

D_{s} = 1 + (2/m) tr{Δ} +(1/m)∥Δ∥^{2}

Based on definition s_{i} (3.2), obtain m_{S_{k}} (3.5) and variance

D_{S_{k}} = D_{s} + 2 ∑_{j=1}^{k−1}(1 − j/k)K_{ss}(j)

Under approximation (3.1), we have

lim_{k→∞} D_{S_{k}} = D_{s}(1 + r)/(1 − r)

Due to the exponential correlated dependence (3.1), the mixing condition in Birkhoff-Khinchin theorem holds (Yoshihawa, 1992), (Ibragimov and Linnik, 1965), hence convergence (3.4) and P_{D} in (3.7) are also true.

□

REFERENCES

Basseville, M. and Nikiforov, I. (1993). Detection of Abrupt Changes: Theory and Applications, New

York: Prentice-Hall. Boozer, D. D. (1971). On Innovation Sequence Testing of the Kalman Filter, Teledyne Brown Engineering

1: 251–256. Borodkin, L. I. and Mottl, V. V. (1976). Algorithm for Finding the Jump Times of Random Process

Equation Parameters, Automation and Remote Control 37(6): 23–32. Caines, P. E. (1988). Linear Stochastic Systems, New York: John Wiley & Sons. Gaines, H. T. (1971). Flight Test of a Kalman Filter Doppler Aided Strapdown Inertial Navigator, in Proc.

National Aerospace Electronics Conference, New York. Ghosh, B. K. and Sen, P. K. (1991). Handbook of Sequential Analysis, New York: Marcel Dekker. Ibragimov, I. A. and Linnik, Y. V. (1965). Independent and Stationary Related Values, Moscow: Nauka. Johnson, N. L. and Leone, F. C. (1962). Comulative Sum Control Charts: Mathematical Principles Applied

to Their Construction and Use, I, Industrial Quality Control 18: 15–21. Johnson, N. L. and Leone, F. C. (1963). Comulative Sum Control Charts: Mathematical Principles Applied

to Their Construction and Use, II, Industrial Quality Control 19: 29–36. Johnson, N. L. and Leone, F. C. (1964). Comulative Sum Control Charts: Mathematical Principles Applied

to Their Construction and Use, III, Industrial Quality Control 20: 22–28. Jones, R. H., Crowell, D. H., and Kapuniai, L. E. (1970). Change Detection Model for Serially Correlated

Multivariate Data, Biometrics 26(2): 269–280. Kenett, R. S. and Zacks, S. (1998). Modern Industrial Statistics: Design and Control of Quality and

Reliability, Pacific Grove: Duxbury Press. Martin, W. C. and Stubberud, A. R. (1976). Innovation Process in Identification Problems, in Control

and Dynamic Systems – Advances in Theory and Applications, C. T. Leondes, ed., pp. 121–176, New York: Academic Press.

Maybeck, P. S. (1979). Stochastic Models, Estimation and Control, vol. 1, New York: Academic Press. Maybeck, P. S. (1982). Stochastic Models, Estimation and Control, vol. 3, New York: Academic Press. Mehra, R. K. and Peschon, J. (1971). An Innovations Approach to Fault Detection and Diagnosis in

Dynamic Systems, Automatica 7: 637–640. Mosca, E. (1995). Optimal, Predictive, and Adaptive Control, Englewood Cliffs, New Jersey: Prentice

Hall.

Murgu, A. (2002). Multiclass Service Dynamics Replication Using Evolutionary Selection Games, in Proceedings of the 10^{th} International Symposium on Dynamic Games and Applications (ISDG-2002), St. Petersburg, Russia, 8–11 July 2002, pp. 652–656.

Page, E. S. (1954). Continuous Inspection Schemes, Biometrica 41: 100–114.

Roussas, G. G. (1972). Contiguity of Probability Measures. Some Applications in Statistics, New York: Cambridge University Press.

Segen, J. and Sanderson, A. C. (1980). Detecting Changes in a Time Series, IEEE Transactions on Information Theory IT-26(2): 249–255.

Semoushin, I. V. (1979). Optimality Testing for the Adaptive Kalman Filter by the Use of a Realization of a Scalar Process, Proceedings of National Academy of Sciences U.S.S.R. – Technicheskaya Cybernetika 6: 195–198.

Semoushin, I. V. and Polosenko, V. P. (2002). Iterative Cooperative CAD for Avionics Navigation Systems, in Advances in Concurrent Engineering, Proceedings of the 9^{th} International Conference on Concurrent Engineering, Cranfield, United Kingdom, 27–31 July 2002, R. Jardim-Gonçalves, R. Roy, and A. Steigler-Garção, eds., pp. 753–761, Lisse, The Netherlands: A. A. Balkema.

Shewhart, W. A. (1931). Economic Control of Quality of Manufactures Product, Toronto: Van Nostrand.

Siegmund, D. (1985). Sequential Analysis: Tests and Confidence Intervals, New York: Springer-Verlag.

Van-Trees, H. I. (1968). Detection, Estimation, and Modulation Theory. Part I: Detection, Estimation, and Linear Modulation Theory, New York: John Wiley & Sons.

Yoshihawa, K. (1992). Weakly Dependent Stochastic Sequences and Their Application, vol. 1.4. Tokyo: World Scientific.