Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2022.017799

ARTICLE

Pattern-Moving-Based Parameter Identification of Output Error Models with Multi-Threshold Quantized Observations

Xiangquan Li1,2, Zhengguang Xu1,*, Cheng Han1 and Ning Li1

1School of Automation and Electrical Engineering, University of Science and Technology Beijing; Key Laboratory of Knowledge Automation for Industrial Processes, Ministry of Education, Beijing, 100083, China
2School of Information Engineering, Jingdezhen University, Jingdezhen, 333000, China
*Corresponding Author: Zhengguang Xu. Email: xzgustb@163.com
Received: 08 June 2021; Accepted: 03 September 2021

1  Introduction

In the metallurgical, petroleum-chemical and steel industries, there are technologically complicated, highly energy-consuming and polluting large-scale equipments such as electrolytic tank, sintering machine, blast furnace cement rotary kiln and so on. The production process of this kind of equipments presents the following characteristics [1]. 1) The complex system mechanism beyonds the accurate description of mathematical and physical equations; 2) Working conditions and quality parameters are in large quantities, and the system moving mode is full of distributiveness, nonlinearity and parameter perturbations; 3) Some physical and chemical processes are in conformity with statistical law of moving. A feasible method of system modeling and control is the pattern recognition technology for these considered processes [2] and most researchers’ practice is to design the corresponding model and controller according to the different pattern class of the system working condition [3,4]. A novel pattern-moving-based system dynamics description method was proposed in [5]. Its basic idea is to take the pattern class as a moving variable, and it is mapped to a computable space by class centers [6,7], interval numbers [8], and cells [9] due to its lack of arithmetic operation attribute. Furthermore, in view of various metric methods of pattern class, the linear autoregressive model with exogenous input (ARX) or interval ARX (IARX) model was established, and the parameter identification algorithm based on least square [6], minimum-variance-based controller [5], optimal controller [10], state-feedback controller [7] and predictive controller [11] were designed.

Although a series of research achievements have been made in the system modeling, parameter identification and control based on pattern moving, the differences of each pattern sample in one pattern class have never been considered. Since a pattern class is a set of pattern samples with the same or similar characteristics, the hybrid metrics, that is, the combination of implicit metric D¯(⋅) and class center explicit metric D(⋅) is proposed. The way of class center explicit metric indicates the statistical attribute of pattern class, and the way of implicit metric denotes the difference of each pattern sample in one pattern class. By using the hybrid metrics, the outputs of pattern-moving-based system dynamics description are in-line with multi-threshold quantized observations. Therefore, this paper will study the pattern-moving-based parameter identification of a kind of linear system models with multi-threshold quantized observations. To our best knowledge, for the model parameter identification based on accurate measured values, a great number of research results were introduced in [12], such as least square method [13], maximum likelihood and prediction error methods [14], unbiased finite impulse response filter algorithm [15], Kalman linear filtering and prediction algorithm [16,17], Bayesian approach [18], and so forth.

The system with multi-threshold quantized observations is considered as a set-valued system which is different from a conventional one with accurate measurement outputs [19,20]. Although the identification of set-valued systems is not easy task, a series of important research results have been achieved in the past decade. For the linear system identification with threshold quantized observations, a full rank input design method (such as repetitive input design) and an empirical distribution function method were proposed in [21–23]. For the identification of Hammerstein and Wiener systems with binary outputs, the methods of proportional full rank signal, joint identifiability and strong full rank periodic signal inputs were proposed in [24,25], which effectively overcame the problems of system nonlinearity and rough binary output information. The parameter identification of Wiener system was also studied in terms of quantized inputs and binary outputs in [26]. Based on the truncated empirical measurement method of probabilistic statistics, a non-truncated empirical measurement method was proposed for the finite impulse response (FIR) system model with binary outputs, and the progressive effectiveness of the algorithm was demonstrated in [27]. A two-segment design method for a class of FIR systems with binary-valued observations was investigated in [28]. The parameter estimation problem of set-valued system was explained comprehensively from the perspective of combining system principle and practical application in [29]. The idea of parameter identification based on an auxiliary model was systematically expounded in [30]. An auxiliary-model-based least squares recursive algorithm was proposed for a quantized control system with communication constraints in [31].

This work investigates a M-AM-SGRPIA for a kind of linear output error models with multi-threshold quantized observations which is established from the perspective of pattern moving and hybrid metrics. Compared with the existed research results, the main differences and contributions of this paper are summarized as follows:

1) Different from the previous system identification problem of ARX or IARX models based on pattern moving and single metric [6,7], this paper considers hybrid metrics, the model noise distribution and proves the convergence of the designed M-AM-SGRPIA.

2) Compared with the parameter identification algorithms of set-valued system in [21–29], this paper adopts an auxiliary model and designs a M-AM-SGRPIA, which will reduce the estimation error radio to some extern.

3) An auxiliary-model-based least squares recursive algorithm has been designed for a class of linear systems with communication constraints in [31], which is different from the designed algorithm and physical background in this paper.

The outline of this work is organized as follows. Section 2 presents a pattern-moving-based system dynamics description and problem formulation. Section 3 proposes a M-AM-SGRPIA and its convergence is proved by the contraction mapping principle in Section 4. Section 5 demonstrates the feasibility and effectiveness of the M-AM-SGRPIA by two numerical examples. The conclusion comes in Section 6.

Notation: R stands for the real number field; Z+ denotes the positive integer field; λmax[A] is the maximum eigenvalue of matrix A, and λmin[A] is its minimum eigenvalue; f(t)=o(g(t)) represents g(t)>0, limt→∞f(t)g(t)=0; arg⁡maxi(⋅) is the value of i when (⋅) gets the maximum value; ‖⋅‖ is the Euclidean norm; the superscript T denotes the transposition; F(⋅) and F−1(⋅) represent the probability distribution function and its inverse function respectively; E(⋅) is the mathematical expectation; P(⋅) denotes the probability.

2  Preliminary and Problem Formulation

2.1 Pattern-Moving-Based System Dynamics Description

Consider a class of unknown SIMO non-affine nonlinear discrete-time systems as follows:

{y1(k)=f1(y1(k−1),⋯,y1(k−n1),u(k−1),⋯,u(k−m1)),y2(k)=f2(y2(k−1),⋯,y2(k−n2),u(k−1),⋯,u(k−m2)),⋮yq(k)=fq(yn(k−1),⋯,yn(k−nq),u(k−1),⋯,u(k−mq)), (1)

where q>0; yi(k)∈R and u(k)∈R denote the output of fi(⋅) and the whole system's input respectively; ni∈Z+ and mi∈Z+ are the unknown output and input orders; fi(⋅) denotes an unknown nonlinear discrete-time function; i∈{1,⋯,q}.

Assumption 2.1 The input of system (1) is bounded, i.e., a constant M1 exists and satisfies that |u(k)|≤M1.

A pattern-moving-based system dynamics description [32] corresponding to system (1) is proposed in the following three steps:

1) Feature extraction (T(⋅)) and output dimension reduction. A large number of input and output data are collected offline, and the input data set {u(k)} and q-dimensional output vector set {[y1(k),⋯,yq(k)]} are obtained. Drawing on the data-driven modeling ideas [33,34] and through the principal component analysis (PCA) feature extraction [35] of the output data, the first principal component information is obtained, and then one-dimensional principal component information set {y(k)} will be obtained.

2) Classification (M(⋅)) and metrics (D(⋅),D¯(⋅)). Using pattern classification technology to classify the first principal component information, the number of pattern class (N), the class center value (si), class radius (ri) and threshold value (Ci) of each pattern class (dxi) can be obtained, i=[1,…,N]. At any time instant k, the system output vector [y1(k),⋯,yq(k)] corresponds to a working condition and a pattern class dx(k) at the same time. Since the pattern class does not have the arithmetic operation attribute, that is, pattern class 1+ pattern class 2 ≠ pattern class 3, the pattern class variable needs to be measured. Since the pattern class is a set of pattern samples with the same or similar attributes, the method of combining the class center explicit metric D(⋅) and implicit metric D¯(⋅) is adopted, that is, si=D(dxi) and dxi¯=D¯(dxi). The implicit metric values are unknown, but there is a definite relationship between an implicit metric value and a class center explicit metric value by choosing a reasonable classification method, such as |si−dxi¯|≤ri.

3) Establishing the pattern-moving-based system dynamics equations. The input sequence {u(k)}, implicit metric sequence {dx¯(k)} and class center explicit metric sequence {s(k)} are employed to construct the following dynamics equations of the system:

dx¯(k)=f(dx¯(k−1),…,dx¯(k−n),u(k−1),…,u(k−m)); (2)

s(k)=D(M(dx¯(k)))={s1,dx¯(k)∈[s1−r1,s1+r1];s2,dx¯(k)∈[s2−r2,s2+r2];⋮sN,dx¯(k)∈[sN−rN,sN+rN], (3)

where f(⋅) is a unknown system function; m,n denote the input and output orders of system (2) respectively. It exits a classification-metric deviation e(k) between the dx¯(k) and s(k), and |e(k)|=|s(k)−dx¯(k)|≤r1 while s(k)=ci.

If the contribution rate of the first principal component information obtained by feature extraction T(⋅) is greater than 85%, it can be considered that the output information of source system (1) is not lost or loses very little. For many practical SIMO physical systems, the contribution rate of the first principal component information can reach or even exceed 85%. Taking sintering machine as an example, through data analysis, it is found that the contribution rate of first principal component information is far greater than 85% when the temperatures of four positions near the sintering node under the excitation of ignition input are considered and feature extracted [5].

Remark 2.1 In the process of establishing the pattern-moving-based dynamics Eqs. (2) and (3), the condition of parameter configuration of classification method is to ensure that a certain pattern class corresponds to a specific quality index of the product [36]. Furthermore, a physical SISO nonlinear discrete-time system can also be transformed into a pattern-moving-based SISO system, but it does not need the first step of feature extraction process.

2.2 Problem Formulation

Although there inevitably unmodeled dynamics problems, it is common to employ a linear model to approximate the situation that the system (2) is unknown. Choosing a reasonable classification method such as a modified quantized control classification [36], the following linear output error model with multi-threshold quantized observations can be constructed for system (2)–(3).

{dx¯(k)=a1dx¯(k−1)+…+andx¯(k−n)+b1u(k−d−1)+…+bmu(k−d−m);dx¯∧(k)=dx¯(k)+υ(k);s(k)=D(M(dx¯∧(k)))={c1,c1−r1<dx¯∧(k)≤c1+r1;c2,c2−r2<dx¯∧(k)≤c2+r2;⋮cN,cN−rN<dx¯∧(k)≤cN+rN, (4)

where v(k) is the model noise; s(k)∈{c1,c2,⋯,cN} is the class center explicit metric value; D(M(⋅)) denotes a classification-metric process; Ci=ci+ri=ci+1−ri+1 is satisfied under a reasonable classification.

Assumption 2.2 The input and output orders of the model are known and equal, that is, n=m.

Assumption 2.3 v(k) is an independent and identically distributed (i.i.d.) random signal noise, and its probability distribution function F(s) and inverse function F−1(s) are known and continuously differentiable.

Under the Assumption 2.2, the first expression of model (4) can be written as follows:

dx¯(k)=B(q)A(q)u(k)=φT(k)θ, (5)

where A(q)=1−a1q−…−anqn; B(q)=b1q+b2q2+…+bnqn;φ(k)=[dx¯(k−1),⋯,dx¯(k−n),u(k−1),⋯,u(k−n)]T; the parameter vector θ=[a1,⋯,an,b1,⋯,bn]T is to be identified.

If the input sequence {u(k)} and output sequence {dx¯(k)} can be obtained directly, it is known to all that there exists a lot of methods to identify the model parameters. However, one can only get the observation sequence {s(k)} and {u(k)}. Therefore, this paper will design a parameter identification algorithm based on these two sequences.

Remark 2.2 The model orders and cumulative distribution function of v(k) are supposed to be known in Assumption 2.2 and Assumption 2.3, respectively. These two problems are not the focus of this paper. If the model orders are unknown, in practice, the order selection may mainly depend on the experiential knowledge of the system. From a theoretical point of view, the model order selection could be addressed by means of the parameter estimation technique with quantized observations and drawing on the classical order selection methods of cross-validation, information criteria, the F-test and the statistical tests on the residuals. In view of the unknown distribution function of model noise, it is not easy task to identify that. However, one can refer to the parameterization technique proposed in [37].

3  Design of Parameter Identification Algorithm

3.1 Estimation of Implicit Metric Value dx¯(k)

For the convenience of calculation, s(k) is denoted by ss(k)=[ssk1,ssk2,⋯,sskN]T, and the corresponding indicator function of each element sski is depicted as follows:

sski={1,s(k)≤Ci;0,others, (6)

where N is the number of the pattern class; Ci denotes the threshold value and it satisfies that Ci=ci+ri.

The measurement technology of random repeatability test is employed to design the input u(k). The size of input cycle is set to be t(t>2n) and the number of cycles is l. Therefore, the noise free output dx¯(k) is an unmeasurable signal with cycle equal to t. Totally t⋅l rows of vector data {ssT(k)} and a data matrix Q(t⋅l)⋅N is produced. It is known from Assumption 2.3 that the probability distribution of υ(k) is F(s), and it can be derived from the random repeatability test technology that the occurrence probability of event dx¯^(k)≤Ci at time instant k is pi, e.g., the occurrence probability of event sski=1 is pi=P{dx¯^(k)≤Ci}=P{dx¯(k)+υ(k)≤Ci}=F(Ci−dx¯(k)), and then dx¯(k)=Ci−F−1(pi). Due to the influence of υ(k), the corresponding output s(k) obtained from each test may differ under the identical test conditions. Let ξli=1l∑j=0l−1ssk+t⋅ji. It is known from the probability statistics that {ssi(k+t⋅l)} is the result of Bernoulli tests and it satisfies that E(ξli)=pi). Therefore, it follows that ξli is the estimated value of pi. Further, the estimated value of dx¯(k) can be gained that dx¯~(k)=Ci−F−1(ξli) and the estimation error e(k)=dx¯~(k)−dx¯(k). In terms of the estimation error e(k), there exists two Lemmas and one conclusion as follows.

Lemma 3.1 ([31]) For the model (4) and the output sequence {ssi(k+t⋅l)} generated by repeatability test, the estimation error sequence {e(k):e(k)=dx¯~(k)−dx¯(k)} is a time-varying sequence with mean 0 and bounded variance, and the estimation error e(k) converges to 0 according to probability 1.

Lemma 3.2 ([31]) Under the full order input sequence {u(k+t⋅j)},1≤k≤t,0≤j≤l−1 with cycle equal to t(t>2n), the signals to be identified generated by model (4) have the property of persistent excitation, e.g.,

1k∑i=1kφ(k0−i+1)φT(k0−i+1)≥αI, (7)

where α>0; k0 denotes the starting time; I is an identity matrix.

According to Lemma 3.1 and Lemma 3.2, the estimation error sequence {e(k),Γk} is regarded as a martingale difference sequence defined in a probability space {Ω,Γ,P}, and it satisfies the following conditions:

A1) E[e(k)|Γk−1]=0,a.s.;

A2) E[e2(k)|Γk−1]=σe2(k)≤Mσς2<∞,

where {Γk} is an σ algebra sequence generated by {e(k)} and σ¯e2=Mσς2.

3.2 Design of M-AM-SGRA

To accomplish the identification task, it is necessary to determine the probability estimate ξli and find the corresponding threshold Ci firstly. Considering the presence of υ(k) , it is required to determine the corresponding threshold segmentation point Ci of s(k) at the occurrence time of event sski=1. For the matrix Q(t⋅l)⋅N , the corresponding row k is the demonstrative vector ssT(k), and it is known from its definition that if ssT(k), then sskj=1, j≥i. The ξli and threshold value Ci are determined by the following formula:

ξli={arg⁡maxi{ξli:ξli=1l∑j=0l−1ssk+t⋅ji},0<ξli<1,i=1,2,…,N;Num;ξli=0 or ξlj=1,i,j∈[1,N], (8)

1) If ξli∈(0,1), the maximum value of ξli and its corresponding interval number i can be found. Further, the corresponding threshold Ci is determined.

2) If ξli is equal to 0 or 1, then corresponding interval number i can be determined when the first ξli=1. Further, the corresponding threshold Ci is determined. Furthermore, e.g., let ξli=Num=0.99,which is a truncation data.

The {e(k),Γk} is considered as a non-stationary martingale difference sequence with mean 0 in Section 3.1. Therefore, one can establish a new identification model by using the estimated value dx¯~(k) such as:

dx¯~(k)=dx¯(k)+e(k)=φT(k)θ+e(k). (9)

Due to the unknown variables dx¯(k−1),⋯,dx¯(k−n) exist in (9), an auxiliary model is established as follows:

dx¯a(k)=φ~T(k)θ~(k), (10)

where θ~(k) is an estimate of θ, and φ~(k)=[dx¯a(k−1),⋯,dx¯a(k−n),u(k−1),⋯,u(k−n)]T is an estimate of φ(k) at time instant k.

By introducing a convergence index ε(12<ε≤1) and using the auxiliary model (10), a modified auxiliary model stochastic gradient recursive algorithm (M-AM-SGRA) for the parameter vector θ is designed as follows:

θ~(k)=θ~(k−1)+φ~(k)rε(k)π(k),θ~(0)=12np0,p0>0; (11)

π(k)=dx¯~(k)−φ~T(k)θ~(k−1); (12)

r(k)=r(k−1)+‖φ~(k)‖2,r(0)=1; (13)

dx¯a(k)=φ~T(k)θ~(k),dx¯a(−i)=1p0,i=0,1,⋯,n; (14)

dx¯~(k)=Ci−F−1(ξli),ξli=arg⁡maxi(ξli:ξli=1l∑j=0l−1ssk+t⋅ji)or Num. (15)

Remark 3.1 A truncation method is adopted to get the value ξli. The determination of truncation rules depends on prior information. The truncated empirical measure method was proposed in [25], and the asymptotic efficiency of the algorithm was proved. This paper chooses a reasonable truncation value of ξli by referring to [25].

4  Main Results

For the designed M-AM-SGRA, the following Lemma is given and its convergence is to be proved.

Lemma 4.1 ([31]) For the given M-AM-SGRA (11)–(15), the following inequalities hold.

∑i=1k‖φ~(i)‖2r(i)≤ln⁡r(k),a.s.; (16)

∑i=1∞‖φ~(i)‖2rη(i)<∞,a.s.η>1 (17)

where a.s. means almost surely.

Theorem 4.1 For the model (9) and the corresponding M-AM-SGRA (11)–(15), under the conditions that A1) and A2) are satisfied, if A(q) is strictly positive real, r(k)→∞, letting

R(k)=∑i=1∞φ~(i)φ~T(i)∈R2n⋅2n (18)

and it is assumed that the following inequality holds

rε(k)≥12‖φ~(k)‖2, (19)

thus, the estimation error vector θˇ(k)(θˇ(k)=(θ~(k)−θ)) satisfies that

‖θ~(k)−θ‖2=o{r2−ε(k)λmin[R(k)]},a.s.. (20)

Proof.

It is concluded from (11) that

θˇ(k)=θˇ(k−1)+φ~(k)rε(k)π(k)=θˇ(k−1)+φ~(k)rε(k)[φ~T(k)θ+e(k)−φ~T(k)θ~(k−1)]

=θˇ(k−1)+φ~(k)rε(k)[−φ~T(k)θˇ(k−1)+e(k)]=θˇ(k−1)+φ~(k)rε(k)[−yˇ(k)+e(k)], (21)

where yˇ(k)=φ~T(k)θˇ(k−1).

Taking the norm of both sides of (21), one has

‖θˇ(k)‖2=‖θˇ(k−1)+φ~(k)rε(k)[−yˇ(k)+e(k)]‖2

=‖θˇ(k−1)‖2+2rε(k)θˇT(k−1)φ~(k)[−yˇ(k)+e(k)]+‖φ~(k)‖2r2ε(k)[−yˇ(k)+e(k)]2

=‖θˇ(k−1)‖2−2rε(k)−‖φ~(k)‖2r2ε(k)yˇ2(k)+2[rε(k)−‖φ~(k)‖2]r2ε(k)yˇ(k)e(k)+‖φ~(k)‖2r2ε(k)e2(k)

≤‖θˇ(k−1)‖2−2rε(k)−‖φ~(k)‖2r2ε(k)yˇ2(k)+‖φ~(k)‖2r2ε(k)e2(k)+2rε(k)yˇ(k)e(k). (22)

Since e(k) is uncorrelated with θˇ(k), yˇ(k), r(k), and φ~(k), and it is Γk−1 measurable, according to (16)-(17) and taking the conditional expectation of (22) with respect to Γk−1, one gets

E[‖θˇ(k)‖2|Γk−1]≤‖θˇ(k−1)‖2−2rε(k)−‖φ~(k)‖2r2ε(k)yˇ2(k)+‖φ~(k)‖2r2ε(k)Mσς2 (23)

According to Lemma 4.1, the sum of the third term on the right side of (23) from t=1 to t=∞ is finite. Based on the condition (19), it can be concluded from the martingale convergence Theorem (Lemma d.5.3 in [39]) that ‖θˇ(k)‖2 uniformly converges to a finite random variable C and it is noted as

limk→∞‖θˇ(k)‖2=C<∞,a.s.. (24)

It can be also derived from the martingale convergence Theorem that

∑k=1∞2rε(k)−‖φ~(k)‖2r2ε(k)yˇ2(k)<∞,a.s.. (25)

Letting ε=1, it is easy to conclude that (19) is tenable to get and one has

∑k=1∞2rε(k)−‖φ~(k)‖2r2ε(k)yˇ2(k)=[∑k=1∞2r(k)−‖φ~(k)‖2r2(k)yˇ2(k)]ε=1=∑k=1∞r(k)+r(k−1)r2(k)yˇ2(k)

≥∑k=1∞1r(k)yˇ2(k)≥[∑k=1∞1r2−κ(k)yˇ2(k)]12<κ≤1 (26)

One can sequentially get

∑k=1∞[1r(k)yˇ2(k)]≤∑k=1∞[2rε(k)−‖φ~(k)‖2r2ε(k)yˇ2(k)], (27)

[∑k=1∞1r2−k(k)yˇ2(k)]12<κ≤1≤∑k=1∞[2rε(k)−‖φ~(k)‖2r2ε(k)yˇ2(k)]. (28)

Further, using the Kronecker Lemma (Lemma d.5.3 in [38]) for inequalities (27), (28) to yield

limk→∞1r(k)∑i=1kyˇ2(i)=0, (29)

and

limk→∞1r2−ε(k)∑i=1kyˇ2(i)=0. (30)

It is known from the above proof process that the parameter estimation error is continuous and bounded. Using (21), one obtains

θˇ(k)=θˇ(k−i)+∑j=0i−1φ~(k−j)rε(k−j)[−yˇ(k−j)+e(k−j)],i≥1.

Replacing k with k−i in yˇ=φ~T(k)θˇ(k−1) yields

φ~T(k−i)θˇ(k)=yˇ(k−i)+φ~T(k−i)∑j=0iφ~(k−j)rε(k−j)[−yˇ(k−j)+e(k−j)]. (31)

Taking the square on both sides of (31) and adopting the inequality (a+b)2≤2(a2+b2), one gets

[φ~T(k−i)θˇ(k)]2={yˇ(k−i)+φ~T(k−i)∑j=0iφ~(k−j)rε(k−j)[−yˇ(k−j)+e(k−j)]}2

≤2yˇ2(k−i)+2‖φ~(k−i)‖2‖∑j=0iφ~(k−j)rε(k−j)[−yˇ(k−j)+e(k−j)]‖2. (32)

Since e(k−j) is uncorrelated with yˇ(k−j), r(k−j), and φ~(k−j), and it is Γk−1 measurable, according to (16), (17) and taking the conditional expectation on both sides of (32) with respect to Γk−1, one has

E{[φ~T(k−i)θˇ(k)]2|Γk−1}≤2yˇ2(k−i)+2‖φ~(k−i)‖2∑j=0i{‖φ~(k−i)‖2r2ε(k−j)[yˇ2(k−i)+Mσς2]}. (33)

Taking the sum from i=0 to i=k−1 on both sides of (33) and dividing by r2−ε(k) yield

E{∑i=0k−1[φ˜ T(k−i)θˇ(k)]2|Γk−1}1r2−ε(k)≤2r2−ε(k)∑i=0k−1yˇ2(k−i)+2∑i=1k−1‖φ˜(k−i)‖2r2−ε(k)∑j=0i[‖φ˜(k−j)‖2r2ε(k−j)yˇ2(k−j)]+2∑i=1k−1‖φ˜(k−i)‖2r2−ε(k)∑j=0i[‖φ˜(k−j)‖2r2ε(k−j)Mσς2],

where

ψ1(k)=2∑i=1k−1‖φ~(k−i)‖2r2−ε(k)∑j=0i[‖φ~(k−j)‖2r2ε(k−j)Mσς2],

and

ψ2(k)=2∑i=1k−1‖φ~(k−i)‖2r2−ε(k)∑j=0i[‖φ~(k−j)‖2r2ε(k−j)yˇ2(k−j)].

According to Lemma 4.1, the following two inequalities can be obtained.

ψ1(k)=2∑i=1k−1‖φ~(k−i)‖2r2−ε(k)∑j=0i[‖φ~(k−j)‖2r2ε(k−j)Mσς2]

=22r2−ε(k){[r(k−1)−r(0)]‖φ~(k)‖2r2ε(k)Mσς2+∑i=1k−1[[r(k−1)−r(0)]‖φ~(k)‖2r2ε(k)Mσς2]}

≤2r2−ε(k)‖φ~(k)‖2r2ε−1(k)Mσς2+2r2−ε(k)∑i=1k−1[‖φ~(k)‖2r2ε−1(i)Mσς2]

=2Mσς2r2−ε(k)∑i=1k[r2−2ε‖φ~(i)‖2r(i)]≤2Mσς2ln(r(k))rε(k)→0,a.s.,k→∞,12<ε≤1.

ψ2(k)=2∑i=1k−1‖φ~(k−i)‖2r2−ε(k)∑j=0i[‖φ~(k−j)‖2r2ε(k−j)yˇ2(k−j)]

=2‖φ~(k−1)‖2r2−ε(k)[‖φ~(k)‖2r2ε(k)yˇ2(k)+‖φ~(k−1)‖2r2ε(k−1)yˇ2(k−1)]

+2‖φ~(k−2)‖2r2−ε(k)[‖φ~(k)‖2r2ε(k)yˇ2(k)+‖φ~(k−1)‖2r2ε(k−1)yˇ2(k−1)+‖φ~(k−2)‖2r2ε(k−2)yˇ2(k−2)]

+⋯

+2‖φ~(1)‖2r2−ε(k)[‖φ~(k)‖2r2ε(k)yˇ2(k)+‖φ~(k−1)‖2r2ε(k−1)yˇ2(k−1)+…+‖φ~(1)‖2r2ε(1)yˇ2(1)]

=2r2−ε(k)[r(k−1)−r(0)]‖φ~(k)‖2r2ε(k)yˇ2(k)+2r2−ε(k)∑i=1k−1[[r(i)−r(0)]‖φ~(i)‖2r2ε(i)yˇ2(i)]

≤2r2−ε(k)2rε(k)r2ε−1(k)yˇ2(k)+2r2−ε(k)∑i=1k−1[2rε(k)r2ε−1(i)yˇ2(i)]

=4r2−ε(k)∑i=1k[1rε−1(i)yˇ2(i)]≤4r(k)∑i=1kyˇ2(i)→0,a.s.,k→∞,12<ε≤1.

Using the above two inequalities, one has

E{∑i=0k−1[φ~T(k−i)θˇ(k)]2|Γk−1}1r2−ε(k)=E{θˇT(k)R(k)θˇ(k)|Γk−1}r2−ε(k)

≤2r2−ε(k)∑i=0k−1yˇ2(k−i)+ψ1(k)+ψ2(k)→0,a.s. k→∞,12<ε≤1.

Therefore,

E{θˇT(k)R(k)θˇ(k)|Γk−1}r2−ε(k)→0,a.s.,k→∞.

and

‖θ~(k)−θ‖2=o{r2−ε(k)λmin[R(k)]},a.s..

The proof of this Theorem is completed.

5  Numerical Examples

Example 1: Consider a SISO linear output error model with multi-threshold quantized observations which has been established based on pattern-moving and hybrid metrics as follows.

{dx¯∧(k)=0.5dx¯(k−1)+0.3dx¯(k−2)+u(k−1)+0.7u(k−2)+v(k)s(k)=D(M(dx¯∧(k)))={c1,c1−r1≤dx¯∧(k)≤c1+r1;c2,c2−r2<dx¯∧(k)≤c2+r2;⋮c13,c13−r13<dx¯∧(k)≤c13+r13,