Maximum Likelihood Recursive Least Squares ... - Semantic Scholar
Circuits Syst Signal Process (2014) 33:2971–2986 DOI 10.1007/s00034-014-9783-8 SHORT PAPER
Maximum Likelihood Recursive Least Squares Estimation for Multivariable Systems Junhong Li · Feng Ding · Ping Jiang · Daqi Zhu
Received: 9 April 2013 / Revised: 19 March 2014 / Accepted: 19 March 2014 / Published online: 16 April 2014 © Springer Science+Business Media New York 2014
Abstract This paper discusses parameter estimation problems of the multivariable systems described by input–output difference equations. We decompose a multivariable system to several subsystems according to the number of the outputs. Based on the maximum likelihood principle, a maximum likelihood-based recursive least squares algorithm is derived to estimate the parameters of each subsystem. Finally, two numerical examples are provided to verify the effectiveness of the proposed algorithm. Keywords Recursive identification · Multivariable systems · Parameter estimation · Maximum likelihood 1 Introduction The mathematical models can describe the relationships among the system variables [1,17,28]. Mathematical modeling plays an important role in control [14,23,27] J. Li (B) · P. Jiang School of Electrical Engineering, Nantong University, Nantong 226019, People’s Republic of China e-mail: [email protected] P. Jiang e-mail: [email protected] F. Ding School of Internet of Things Engineering, Jiangnan University, Wuxi 214122, People’s Republic of China e-mail: [email protected] D. Zhu Laboratory of Underwater Vehicles and Intelligent Systems, Shanghai Maritime University, Shanghai 201306, People’s Republic of China e-mail: [email protected]
2972
Circuits Syst Signal Process (2014) 33:2971–2986
and other engineering fields [25,26]. Many practical systems in industrial processes are multivariable systems. Compared with the single-input single-output systems, the multivariate systems have many input and output variables [21], the problem of building their mathematical model is more difficult than the single-input single-output case. Different mathematical models have been used to describe the multivariable systems, such as the state-space model, the transfer function matrix model [19]. Some earlier work for identifying multivariable systems focused on the state-space model [8,20,31]. In this area, Guidorzi studied the identification problem of multivariable systems and derived the connections between the state–space model and the input– output difference equation canonical model [10,11], his work provided an effective solution to identify multivariable systems. Recently, many modeling and estimation methods have been proposed to estimate the parameters of multivariable systems [15,18]. For example, Wang et al. [29] presented the hierarchical gradient-based iterative algorithm for multivariable controlled autoregressive system with autoregressive noises using the hierarchical identification principle and the gradient search; Mercère and Bako [22] studied a subspace-based identification algorithm for multivariable state-space systems; Wills et al. [32] proposed gradient-based search algorithms for multivariable systems; Zheng [35] proposed a bias eliminated least squares for multiple-input single-output systems. The maximum likelihood principle is a useful tool in statistics [38] and system identification [12,34]. In the area of maximum likelihood identification, Diversi et al. [7] proposed a maximum likelihood method for the identification of errors-in-variables models corrupted by white and uncorrelated Gaussian noises; Gibson and Ninness [9] developed a robust maximum likelihood method for linear state-space models; Wills et al. [33] proposed a maximum likelihood estimation algorithm of state-space models. Recently, Li and Ding [13] studied maximum likelihood estimation algorithms for Hammerstein nonlinear systems. Iterative or recursive algorithms are important for modeling a system [2–4,36]. Based on the work in [13], the objective of this paper is to study the recursive parameter estimation of multivariable input–output difference equation model using the maximum likelihood principle. The basic idea is to decompose a multivariable system to several multiple-input single-output subsystems, and to derive the recursive maximum likelihood algorithm for each subsystem. The rest of this paper is organized as follows. Section 2 gives the system description of the multivariable systems. Section 3 derives a maximum likelihood-based recursive least squares algorithm for multivariable systems. Section 4 provides two numerical examples to illustrate the proposed method. Finally, we offer some concluding remarks in Sect. 5.
2 The System Description Consider a multivariable system described by the following input–output representation [10]:
Circuits Syst Signal Process (2014) 33:2971–2986
2973
P(z) y(t) = Q(z)u(t) + R(z)v(t), ⎡
p11 (z) p12 (z) · · · p1m (z)
(1)
⎤
⎢ ⎥ ⎢ p21 (z) p22 (z) · · · p2m (z) ⎥ ⎢ ⎥ P(z) := ⎢ . ⎥ ∈ Rm×m , .. .. ⎢ .. ⎥ . . ⎣ ⎦ pm1 (z) pm2 (z) · · · pmm (z) ⎡ ⎤ q11 (z) q12 (z) · · · q1r (z) ⎢ ⎥ ⎢ q21 (z) q22 (z) · · · q2r (z) ⎥ ⎢ ⎥ ∈ Rm×r , Q(z) := ⎢ . .. .. ⎥ ⎢ .. ⎥ . . ⎦ ⎣ qm1 (z) qm2 (z) · · · qmr (z) ⎡ ⎤ d11 (z) d12 (z) · · · d1m (z) ⎢ ⎥ ⎢ d21 (z) d22 (z) · · · d2m (z) ⎥ ⎢ ⎥ ∈ Rm×m , R(z) := ⎢ . .. .. ⎥ ⎢ .. ⎥ . . ⎦ ⎣ dm1 (z) dm2 (z) · · · dmm (z) where u(t) := [u 1 (t), u 2 (t), . . . , u r (t)]T ∈ Rr and y(t) := [y1 (t), y2 (t), . . . , ym (t)]T ∈ Rm are the system input vector and output vector, respectively, v(t) := [v1 (t), v2 (t), . . . , vm (t)]T ∈ Rm is a Gaussian distributed white noise vector with zero mean and variance (σ12 , σ22 , . . . , σm2 ), z is a unit forward shift operator: z y(t) = y(t + 1), and P(z), Q(z), and R(z) are polynomial matrices whose entries are defined as [10] pi j (z) = δi j z n i j − ai j,n i j z n i j −1 − · · · − ai j,2 z − ai j,1 , i, j = 1, 2, . . . , m, (2) qi j (z) = bi j,n i z n i −1 + bi j,n i −1 z n i −2 · · · + bi j,2 z + bi j,1 , i = 1, 2, . . . , m, j = 1, 2, . . . , r, di j (z) = δi j z n i j + di j,n i j z n i j −1 + · · · + di j,2 z + di j,1 , i, j = 1, 2, . . . , m,
ni j
⎧ i= j ⎪ ⎨ ni , 1, i = j i < j , δi j = = min[n i , n j ], . ⎪ 0, i = j ⎩ min[n i + 1, n j ], i > j
(3) (4)
(5)
The ith subsystem in (1) can be equivalently written as m
j=1
pi j (z)y j (t) =
r
j=1
qi j (z)u j (t) +
m
j=1
di j (z)v j (t), i = 1, 2, . . . , m.
(6)
2974
Circuits Syst Signal Process (2014) 33:2971–2986
That is yi (t + n i ) =
m
[ai j,n i j y j (t + n i j − 1)+ai j,n i j −1 y j (t+n i j − 2) + · · · + ai j,1 y j (t)] j=1 r
+ [bi j,n i u j (t + n i − 1)+bi j,n i −1 u j (t + n i − 2)+· · ·+bi j,1 u j (t)] j=1 m
+ [di j,n i j v j (t + n i j − 1) + di j,n i j −1 v j (t + n i j − 2) + · · · i=1
+ di j,1 vi (t)] + vi (t + n i ). Define the parameter vector of subsystem i as θ iT := [ai1 , ai2 , . . . , aim , bi1 , bi2 , . . . , bir , d i1 , d i2 , . . . , d im ] ∈ R1×(2μi +r n i ) , μi := n i1 + n i2 + · · · + n im , ai j := [ai j,n i j , ai j,n i j −1 , . . . , ai j,1 ] ∈ R1×n i j , bi j := [bi j,n i , bi j,n i −1 , . . . , bi j,1 ] ∈ R1×n i , d i j := [di j,n i j , di j,n i j −1 , . . . , di j,1 ] ∈ R1×n i j , ϕ i (t) = [ ¯yi1 (t), ¯yi2 (t), . . . , ¯yim (t), u¯ i1 (t), u¯ i2 (t), . . . , u¯ ir (t), v¯ i1 (t), v¯ i2 (t), . . . , v¯ im (t)]T ∈ R(2μi +r n i ) , ¯yi j (t) := [y j (t − 1), y j (t − 2), . . . , y j (t − n i j )] ∈ R1×n i j , u¯ i j (t) := [u j (t − 1), u j (t − 2), . . . , u j (t − n i j )] ∈ R1×n i , v¯ i j (t) := [v j (t − 1), v j (t − 2), . . . , v j (t − n i j )] ∈ R1×n i j . Then the ith subsystem can be written as the following identification model: yi (t) = θ iT ϕ i (t) + vi (t).
(7)
The objective of this paper is to present new identification methods for determining the parameter vector θ i from the available input–output data ϕ i (t) and yi (t). 3 The Maximum Likelihood-Based Recursive Least Squares Algorithm Referring to the work in [13], we propose the maximum likelihood criterion function for the ith subsystem in (7) as follows: 1 2 vi (k). 2 t
J (θ i , t) =
(8)
k=1
Let θˆ i (t) denote the maximum likelihood estimate of θ i at time t. The Taylor expansion of J (θ i , t) at point θ i = θˆ i (t − 1) is given by
Circuits Syst Signal Process (2014) 33:2971–2986
2975
∂ J (θˆ i (t − 1), t) J (θ i , t) = J (θˆ i (t − 1), t) + [θ i − θˆ i (t − 1)] ∂θ i 1 ∂ 2 J (θˆ i (t −1), t) + [θ i − θˆ i (t −1)]T [θ i − θˆ i (t −1)]+◦(θ i − θˆ i (t −1)2 ). 2 ∂θ i θ iT Then the maximum likelihood estimate at time t can be achieved by the Newton– Raphson method −1 2 ˆ ∂ J (θˆ i (t − 1), t) ˆθ i (t) = θˆ i (t − 1) − ∂ J (θ i (t − 1), t) . T ∂θ i θ i ∂θ i In order to compute as a recursive form:
∂ 2 J (θˆ i (t−1),t) ∂θ i θ iT
and
∂ J (θˆ i (t−1),t) , ∂θ i
(9)
we rewrite the cost function in (8)
1 2 1 vi (k) = J (θ i , t − 1) + vi2 (t). 2 2 t
J (θ i , t) =
(10)
k=1
Computing the derivative of J (θ i , t) at point θˆ i (t − 1) gives ∂ J (θˆ i (t − 1), t) ∂ J (θˆ i (t − 1), t − 1) ∂vi (t) = + vˆi (t) . ∂θ i ∂θ i ∂θ i θˆi (t−1)
(11)
Let ϕˆ i f (t) := − Then we have
∂vi (t) ∈ R(2μi +r n i )×(2μi +r n i ) . ∂θ i θˆi (t−1)
∂ J (θˆ i (t − 1), t) = −vi (t)ϕˆ i f (t). ∂θ i
(12)
Computing the second-order partial derivative of J (θ i , t) with respect to θ i gives ∂vi2 (t) ∂ 2 J (θˆ i (t − 1), t) ∂ 2 J (θˆ i (t − 1), t − 1) T = + ϕˆ i f (t)ϕˆ i f (t) + vi (t) ∂θ i θ iT ∂θ i ∂θ iT ∂θ i ∂θ iT θˆ i (t−1) ≈
∂ 2 J (θˆ i (t − 1), t − 1) + ϕˆ i f (t)ϕˆ iT f (t). ∂θ i ∂θ iT
Let P i−1 (t) :=
∂ 2 J (θˆ i (t − 1), t) . ∂θ i ∂θ iT
(13)
2976
Circuits Syst Signal Process (2014) 33:2971–2986
Equation (13) can be written as P i−1 (t) = P i−1 (t − 1) + ϕˆ i f (t)ϕˆ iT f (t).
(14)
Applying the matrix inversion formula ( A + BC)−1 = A−1 − A−1 B(I + C A−1 B)−1 C A−1 to (14) gives P i (t) = P i (t − 1) −
P i (t − 1)ϕˆ i f (t)ϕˆ iT f (t) P i (t − 1) 1 + ϕˆ iT f (t) P i (t − 1)ϕˆ i f (t)
.
(15)
Define the gain vector L i (t) : = P i (t)ϕˆ i f (t) =
P i (t − 1)ϕˆ i f (t) 1 + ϕˆ iT f (t) P i (t − 1)ϕˆ i f (t)
∈ R(2μi +r n i ) .
(16)
Thus, Eq. (15) can be rewritten as P i (t) = [I − L i (t)ϕˆ iT f (t)] P i (t − 1).
(17)
In the following, we compute the filtered information vector ϕˆ i f (t) of the ith subsystem. Let aˆ i j (t), bˆ i j (t), and dˆ i j (t) be the estimates of ai j , bi j , and d i j at time t, i.e., aˆ i j (t) := [aˆ i j,n i j (t), aˆ i j,n i j −1 (t), . . . , aˆ i j,1 (t)] ∈ R1×n i j , bˆ i j (t) := [bˆi j,n i (t), bˆi j,n i −1 (t), . . . , bˆi j,1 (t)] ∈ R1×n i , dˆ i j (t) := [dˆi j,n i j (t), dˆi j,n i j −1 (t), . . . , dˆi j,1 (t)] ∈ R1×n i j , dˆii (t, z) := z n i + dˆii,n i (t)z n i −1 + dˆii,n i −1 (t)z n i −2 + · · · + dˆii,1 (t). vi (t) in (7) can be equivalently written as vi (t) = yi (t) − θ iT ϕ i (t) ⎡ ⎤ m r m
= yi (t) − ⎣ ai j ¯yiTj (t) + bi j u¯ iTj (t) + d i j v¯ iTj (t)⎦ . j=1
j=1
j=1
(18)
Circuits Syst Signal Process (2014) 33:2971–2986
2977
The partial derivatives of vi (t) with respect to ai j can be computed by ∂ v¯ iT (t) ∂vi (t) ˆ = − ¯yi j (t) − d ii (t − 1) ∂ ai j θˆi (t−1) ∂ ai j θˆi (t−1) = − ¯yi j (t) − [dˆii,n i (t − 1), . . . , dˆii,2 (t − 1), dˆii,1 (t − 1)] ∂vi (t − 1) ∂vi (t − n i + 1) ∂vi (t − n i ) T × ,..., , ∂ ai j ∂ ai j ∂ ai j θˆ i (t−1)
−1 −n i ∂vi (t) ˆ ˆ = −¯yi j (t)−[dii,n i (t − 1)z +· · · + dii,1 (t − 1)z ] . ∂ ai j θˆi (t−1) Then we have ∂vi (t) 1 ¯yi j (t) =: − ˆ¯yi j, f (t). =− ∂ ai j θˆi (t−1) z −n i dˆii (t − 1, z) Similarly, the partial derivatives of vi (t) with respect to bi j and d i j can be written as follows: ∂vi (t) 1 u¯ i j (t) =: −uˆ¯ i j, f (t), =− −n ˆ ∂ bi j θˆi (t−1) z i dii (t − 1, z) ∂vi (t) 1 =− v¯ˆ i j (t) =: −v¯ˆ i j, f (t). −n ˆ i ∂ d i j θˆ (t−1) z dii (t − 1, z) i
Then the filtered information vectors ˆ¯yi j, f (t), uˆ¯ i j, f (t) and vˆ¯ i j, f (t) can be computed by ˆ¯yi j, f (t) = ¯yi j (t) − dˆii,n i (t − 1) ˆ¯yi j, f (t − 1) − · · · − dˆii,1 (t − 1) ˆ¯yi j, f (t − n i ), u¯ˆ i j, f (t) = u¯ i j (t) − dˆii,n i (t − 1)u¯ˆ i j, f (t − 1) − · · · − dˆii,1 (t − 1)uˆ¯ i j, f (t − n i ), v¯ˆ i j, f (t) = vˆ¯ i j (t) − dˆii,n i (t − 1)vˆ¯ i j, f (t − 1) − · · · − dˆii,1 (t − 1)vˆ¯ i j, f (t − n i ). Thus the filtered information vector can be written as ∂vi (t) ϕˆ i f (t) := − ∂θ i θˆi (t−1) ∂vi (t) ∂vi (t) ∂vi (t) ∂vi (t) ∂vi (t) = − , ,..., , , ,..., ∂ ai1 ∂ ai2 ∂ aim ∂ bi1 ∂ bi2 ∂vi (t) ∂vi (t) ∂vi (t) ∂vi (t) T , , ,..., ∂ bir ∂ d i1 ∂ d i2 ∂ d im = [ ¯ˆyi1, f (t), ¯ˆyi2, f (t), . . . , ¯ˆyim, f (t), u¯ˆ i1, f (t), u¯ˆ i2, f (t), . . . , uˆ¯ ir, f (t), vˆ¯ i1, f (t), vˆ¯ i2, f (t), . . . , vˆ¯ im, f (t)]T .
2978
Circuits Syst Signal Process (2014) 33:2971–2986
Let vˆ¯ i j (t) be the estimates of v¯ i j (t) at time t, and let ϕˆ i (t) denote the information vector obtained by replacing v¯ i j (t) with vˆ¯ i j (t), i.e., ϕˆ i (t) = [ ¯yi1 (t), ¯yi2 (t), . . . , ¯yim (t), u¯ i1 (t), u¯ i2 (t), . . . , u¯ ir (t), vˆ¯ i1 (t), vˆ¯ i2 (t), . . . , v¯ˆ im (t)]T ∈ R(2μi +r n i ) , vˆ¯ i j (t) = [vˆ j (t − 1), vˆ j (t − 2), . . . , vˆ j (t)] ∈ R1×n i j . According to the identification model in (7), we have vi (t) = yi (t) − θ iT ϕ i (t). Replacing ϕ i (t) and θ i with their estimates ϕˆ i (t) and θˆ i (t − 1) gives the estimates of vi (t) T vˆi (t) = yi (t) − θˆ i (t − 1)ϕˆ i (t).
According to (9), (12), (16), and (17), we can summarize the maximum likelihoodbased recursive least squares (ML-RLS) algorithm for the ith subsystem of the multivariable system as θˆ i (t) = θˆ i (t − 1) + L i (t)vˆi (t), P i (t − 1)ϕˆ i f (t) , L i (t) = 1 + ϕˆ iT f (t) P i (t − 1)ϕˆ i f (t)
(19) (20)
P i (t) = [I − L i (t)ϕˆ iT f (t)] P i (t − 1), P i (0) = p0 I,
(21)
vˆi (t) = yi (t) − θˆ i (t − 1)ϕˆ i (t), ϕˆ i (t) = [ ¯yi1 (t), ¯yi2 (t), . . . , ¯yim (t), u¯ i1 (t), u¯ i2 (t),
(22)
T
. . . , u¯ ir (t), v¯ˆ i1 (t), vˆ¯ i2 (t), . . . , vˆ¯ im (t)]T , ϕˆ i f (t) = [ ˆ¯yi1, f (t), ˆ¯yi2, f (t), . . . , ˆ¯yim, f (t), uˆ¯ i1, f (t), uˆ¯ i2, f (t), . . . , uˆ¯ ir, f (t), vˆ¯ i1, f (t), vˆ¯ i2, f (t), . . . , vˆ¯ im, f (t)]T , ¯yi j (t) = [y j (t − 1), y j (t − 2), . . . , y j (t − n i )], u¯ i j (t) = [u j (t − 1), u j (t − 2), . . . , u j (t − n i )], vˆ¯ i j (t) = [vˆ j (t − 1), vˆ j (t − 2), . . . , vˆ j (t − n i )], ˆ¯yi j, f (t) = ¯yi j (t)− dˆii,n i (t −1) ˆ¯yi j, f (t − 1)−· · ·− dˆii,1 (t −1) ˆ¯yi j, f (t −n i ),
(23)
(24) (25) (26) (27) (28)
u¯ˆ i j, f (t) = u¯ i j (t)− dˆii,n i (t −1)u¯ˆ i j, f (t −1)−· · ·− dˆii,1 (t −1)uˆ¯ i j, f (t −n i ), (29) v¯ˆ i j, f (t) = vˆ¯ i j (t)− dˆii,n i (t −1)vˆ¯ i j, f (t − 1) − · · · − dˆii,1 (t −1)vˆ¯ i j, f (t −n i ). (30) The steps involved in the ML-RLS algorithm for subsystem i are listed in the following.
Circuits Syst Signal Process (2014) 33:2971–2986
2979
Fig. 1 The flowchart for computing the ML-RLS estimate θˆ i (t)
1. Let t = 1, set the initial values θˆ i (0) = 1(2μi +r n i ) / p0 and P i (0) = p0 I, ˆ¯yii, f (t) = 0, uˆ¯ i j, f (t) = 0, vˆ¯ ii, f (t) = 0 and v¯ i j (t) = 0 for t ≤ 0, p0 is a large number (i.e., p0 = 106 ). 2. Collect the input-output data u(t) and y(t), form ¯yi j (t) by (25), u¯ i j (t) by (26), vˆ¯ i j (t) by (27), form ϕˆ i (t) by (23). 3. Compute ˆ¯yi j, f (t), uˆ¯ i j, f (t), and vˆ¯ i j, f (t) by (28), (29), and (30), respectively, form ϕˆ i f (t) by (24). 4. Compute L i (t) by (20), compute P i (t) by (21) and vˆi (t) by (22). 5. Update the parameter estimate θˆ i (t) by (19). 6. Increase t by 1 and go to step 2. The flowchart of computing the subsystem parameter estimate θˆ i (t) in the ML-RLS algorithms is shown in Fig. 1. A summary of the variables and their dimensions in the ML-RLS algorithm are shown in Table 1 for convenience. The numbers of multiplications and additions in each recursive computation of the ML-RLS algorithm are shown in Table 2.
2980
Circuits Syst Signal Process (2014) 33:2971–2986
Table 1 The dimensions of the variables in the ML-RLS algorithm Item
Variables
Dimensions
1
Input variable
u(t) ∈ Rr
2
Output variable
y(t) ∈ Rm
3
Parameter vectors
ai j ∈ R1×n i j , bi j ∈ R1×n i , d i j ∈ R1×n i j
4
Parameter vector
θ i ∈ R(2μi +r n i )
5
Information vectors
¯yi j (t) ∈ R1×n i j , u¯ i j (t) ∈ R1×n i , v¯ i j (t) ∈ R1×n i j
6
Information vector
ϕ i (t) ∈ R(2μi +r n i )
5
Filtered information vectors
¯yi j, f (t) ∈ R1×n i j , u¯ i j, f (t) ∈ R1×n i , v¯ i j, f (t) ∈ R1×n i j
6
Filtered information vector
ϕ i f (t) ∈ R(2μi +r n i )
7
Covariance matrix
P i (t) ∈ R(2μi +r n i )×(2μi +r n i )
Table 2 The computational efficiency of the ML-RLS algorithm Variables
Number of multiplications
Number of additions
θˆ i (t) = θˆ i (t − 1) + L i (t)vˆi (t)
2μi + r n i
2μi + r n i
T vˆi (t) = yi (t) − θˆ i (t − 1)ϕˆ i (t)
2μi + r n i
2μi + r n i
L i (t) = ζ i (t)/[1 + ϕ iTf (t)ζ i (t)] 2(2μi + r n i ) (2μi + r n i )2
ζ i (t) := P i (t − 1)ϕ i f (t)
2μi + r n i (2μi + r n i − 1)(2μi + r n i )
P i (t) = P i (t − 1) − L i (t)ζ iT (t) (2μi + r n i )2
(2μi + r n i )2
ˆ¯yi j, f (t) = ¯yi j (t) − · · · − dˆii,1 (t − 1) ˆ¯yi j, f (t − n i )
n i2
n i2
uˆ¯ i j, f (t) = u¯ i j (t) − · · · − dˆii,1 (t − 1)uˆ¯ i j, f (t − n i ) v¯ˆ i j, f (t) = v¯ˆ i j (t) − · · · − dˆii,1 (t − 1)v¯ˆ i j, f (t − n i )
n i2
n i2
n i2
n i2
Sum
2(2μi + r n i )2 + 4(2μi +r n i )+3n i2
2(2μi + r n i )2 + 2(2μi + r n i ) + 3n i2
Total flops
4(2μi + r n i )2 + 6(2μi + r n i ) + 6n i2
In order to show the advantages of the ML-RLS algorithm, we give the RELS algorithm for comparison [19]: T θˆ i (t) = θˆ i (t − 1) + L i (t)[yi (t) − θˆ i (t − 1)ϕˆ i (t)], θˆ (0) = 12μi +r n i / p0 , P i (t − 1)ϕˆ i (t) L i (t) = , 1 + ϕˆ iT (t) P i (t − 1)ϕˆ i (t) P i (t) = [I − L i (t)ϕˆ iT (t)] P i (t − 1), P i (0) = p0 I 2μi +r n i , T vˆi (t) = yi (t) − θˆ (t)ϕˆ i (t),
i
ϕˆ i (t) = [ ¯yi1 (t), ¯yi2 (t), . . . , ¯yim (t), u¯ i1 (t), u¯ i2 (t), . . . , u¯ ir (t), vˆ¯ i1 (t), vˆ¯ i2 (t), . . . , vˆ¯ im (t)]T ,
Circuits Syst Signal Process (2014) 33:2971–2986
2981
¯yi j (t) = [y j (t − 1), y j (t − 2), . . . , y j (t − n i )], u¯ i j (t) = [u j (t − 1), u j (t − 2), . . . , u j (t − n i )], vˆ¯ i j (t) = [vˆ j (t − 1), vˆ j (t − 2), . . . , vˆ j (t − n i )]. The numbers of multiplications and additions in each recursive computation of the RELS algorithm are shown in Table 3. Remark Similarly to the RELS algorithm, the ML-RLS algorithm is a recursive algorithm and can be used for online identification. The two algorithms have almost same computational burden, but the ML-RLS algorithm can achieve higher parameter estimation accuracy with data length increases. 4 Numerical Examples Example 1 Consider the multivariable system in (1), where
z −0.18 0.50 0.80 0.49 z − 0.30 0.66 P(z) = , Q(z) = , R(z) = , 0.01 z + 0.19 0.30 0.50 0.70 z − 0.13 u 1 (t) v (t) y1 (t) , u(t) = , v(t) = 1 . y(t) = y2 (t) u 2 (t) v2 (t)
The parameter vectors of this system are θ 1 = [a11,1 , a12,1 , b11,1 , b12,1 , d11,1 , d12,1 ]T = [0.18,−0.50, 0.80, 0.49,−0.30, 0.66]T , θ 2 = [a21,1 , a22,1 , b21,1 , b22,1 , d21,1 , d22,1 ]T = [−0.01,−0.19, 0.30, 0.50, 0.70,−0.13]T, θ1 . = θ2 In simulation, the inputs {u 1 (t)} and {u 2 (t)} are taken as two uncorrelated stochastic signal sequences with zero mean and unit variance, and the noise {v1 (t)} and {v2 (t)} as two Gaussian white noise sequences with zero mean and variances σ12 = 0.502 and σ22 = 0.302 , respectively. Applying the ML-RLS algorithm and the RELS algorithm to estimate the parameters of this system, the parameter estimates and their errors ˆ δ := (t) − / versus t are shown in Tables 4 and 5, the parameter estimates versus t are shown in Figs. 2 and 3 for the two identification algorithms. From Tables 4, 5 and Figs. 2, 3, we can draw the following conclusions. 1. The parameter estimation error of the two algorithms gradually becomes small as t increases. 2. The parameter estimates of the two algorithms converge to their true values. This verifies the effectiveness of the two algorithms. 3. The RELS algorithm has better performance for small data length t, but for large t, the ML-RLS algorithm can obtain higher parameter estimation accuracy.
2982
Circuits Syst Signal Process (2014) 33:2971–2986
Table 3 The computational efficiency of the recursive extended least squares (RELS) algorithm Variables
Number of multiplications
Number of additions
θˆ i (t) = θˆ i (t − 1) + T L i (t)[yi (t)−θˆ i (t −1)ϕˆ i (t)]
2(2μi + r n i )
2(2μi + r n i )
L i (t) = ζ i (t)/[1 + ϕˆ iT (t)ζ i (t)]
2(2μi + r n i )
2μi + r n i
ζ i (t) := P i (t − 1)ϕˆ i (t)
(2μi + r n i − 1)(2μi + r n i )
P i (t) = P i (t − 1) − L i (t)ζ iT (t)
(2μi + r n i )2 (2μi + r n i )2
vˆi (t) = yi (t) − θˆ i (t)ϕˆ i (t)
2μi + r n i
2μi + r n i
Sum
2(2μi + r n i )2 + 5(2μi + r n i )
2(2μi + r n i )2 + 3(2μi + r n i )
Total flops
4(2μi + r n i )2 + 8(2μi + r n i )
T
(2μi + r n i )2
Table 4 The ML-RLS estimates and errors of Example 1 t
100
200
500
1,000
2,000
3,000
True values
a11,1
0.25088
0.15233
0.18478
0.19745
0.19150
0.18123
0.18000
a12,1
−0.47796
−0.44132
−0.48214
−0.49723
−0.51678
−0.50031
−0.50000
b11,1
0.78764
0.77093
0.79562
0.79351
0.80069
0.80793
0.80000
b12,1
0.55179
0.51400
0.52032
0.51363
0.50714
0.50575
0.49000
d11,1
−0.33228
−0.29851
−0.34230
−0.32844
−0.30864
−0.29515
−0.30000
d12,1
0.27637
0.27588
0.48525
0.59703
0.66636
0.66076
0.66000
a21,1
−0.01866
−0.00335
0.01973
−0.00420
0.00255
−0.00424
−0.01000
a22,1
−0.17889
−0.22734
−0.24365
−0.20447
−0.20723
−0.19312
−0.19000
b21,1
0.32892
0.31850
0.31089
0.30847
0.30652
0.30498
0.30000
b22,1
0.46131
0.48062
0.51007
0.50947
0.49683
0.49780
0.50000
d21,1
0.50289
0.58396
0.62662
0.67795
0.69264
0.70257
0.70000
d22,1
−0.18426
−0.07659
0.00862
−0.07252
−0.09546
−0.12559
−0.13000
δ (%)
28.03052
25.83780
15.55223
6.19942
3.16590
1.30100
Example 2 Consider the multivariable system in (1), where ⎡
⎤ ⎡ ⎤ z − 0.66 0.26 −0.30 0.82 0.49 0.13 0.29 z − 0.49 −0.25 ⎦ , Q(z) = ⎣ 0.35 z − 0.27 −0.25 ⎦ , P(z) = ⎣ −0.29 0.89 z − 0.66 0.97 0.89 0.85 ⎡ ⎤ z − 0.23 0.21 0.23 0.11 z − 0.10 −0.10 ⎦ , R(z) = ⎣ −0.15 0.23 z + 0.13 ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ y1 (t) u 1 (t) v1 (t) y(t) = ⎣ y2 (t) ⎦ , u(t) = ⎣ u 2 (t) ⎦ , v(t) = ⎣ v2 (t) ⎦ . y3 (t) u 3 (t) v3 (t) The parameter vectors of the system are
Circuits Syst Signal Process (2014) 33:2971–2986
2983
Table 5 The RELS estimates and errors of Example 1 t
100
200
500
1,000
2,000
3,000
True values
a11,1
0.26979
0.16553
0.17856
0.19388
0.18532
0.17338
0.18000
a12,1
−0.55086
−0.47729
−0.50521
−0.51738
−0.52678
−0.50455
−0.50000 0.80000
b11,1
0.79906
0.77303
0.79541
0.79138
0.80015
0.80736
b12,1
0.54587
0.51208
0.51914
0.51303
0.50521
0.50466
0.49000
d11,1
−0.29052
−0.28952
−0.33279
−0.32531
−0.30361
−0.28987
−0.30000
d12,1
0.31910
0.31184
0.51463
0.63200
0.68262
0.67119
0.66000
a21,1
−0.03664
−0.01313
0.01460
−0.00683
0.00166
−0.00376
−0.01000
a22,1
−0.13035
−0.19882
−0.23128
−0.19590
−0.20203
−0.19031
−0.19000
b21,1
0.32834
0.31529
0.30729
0.30738
0.30564
0.30461
0.30000
b22,1
0.45357
0.47940
0.50894
0.50841
0.49660
0.49758
0.50000
d21,1
0.51231
0.59512
0.63365
0.68342
0.69598
0.70337
0.70000
d22,1
−0.19767
−0.08989
−0.00121
−0.08534
−0.10425
−0.13067
−0.13000
δ (%)
26.24540
23.05107
13.45366
4.39389
3.12132
1.57697
1
b11,1
Parameter estimates
0.8
d
12,1
0.6
b12,1
0.4
a11,1
0.2 0 −0.2
d
11,1
−0.4
a12,1
−0.6 0
500
1000
1500
2000
2500
3000
t Fig. 2 The ML-RLS estimates of subsystem 1 versus t in Example 1
θ 1 = [a11,1 , a12,1 , a13,1 , b11,1 , b12,1 , b13,1 , d11,1 , d12,1 , d13,1 ]T = [0.66, −0.26, 0.30, 0.82, 0.49, 0.13, −0.23, 0.21, 0.23]T , θ 2 = [a21,1 , a22,1 , a23,1 , b21,1 , b22,1 , b23,1 , d21,1 , d22,1 , d23,1 ]T = [−0.29, 0.49, 0.25, 0.35, −0.27, −0.25, 0.11, −0.10, −0.10]T , θ 3 = [a31,1 , a32,1 , a33,1 , b31,1 , b32,1 , b33,1 , d31,1 , d32,1 , d33,1 ]T = [0.29, −0.89, 0.66, 0.97, 0.89, 0.85, −0.15, 0.23, 0.13]T , ⎡ ⎤ θ1 = ⎣ θ2 ⎦ . θ3
2984
Circuits Syst Signal Process (2014) 33:2971–2986
Parameter estimates
0.8
d21,1
0.6
b
0.4
b
22,1 21,1
0.2
a
21,1
0
d22,1
−0.2
a22,1
−0.4 −0.6 0
500
1000
1500
2000
2500
3000
t Fig. 3 The ML-RLS estimates of subsystem 2 versus t in Example 1
0.25
0.2
δ
0.15
0.1
ML−RLS
RELS
0.05
0 0
500
1000
1500
2000
2500
3000
t Fig. 4 The ML-RLS and RELS estimation errors δ versus t in Example 2
In simulation, the inputs {u 1 (t)}, {u 2 (t)}, and {u 3 (t)} are taken as three uncorrelated stochastic signal sequence with zero mean and unit variance, and the noise {v1 (t)}, {v2 (t)}, and {v3 (t)} as three Gauss white noise sequence with zero mean and variances σ12 = σ22 = σ32 = 0.302 . Applying the ML-RLS algorithm and the RELS algorithm to estimate the parameters of this system, the parameter estimates and their errors are shown in Fig. 4. From Fig. 4, we can see that the parameter estimation error of the two algorithms gradually becomes small as t increases. The ML-RLS algorithm has higher parameter estimation accuracy than the RELS algorithm as t increases, thus can effectively identify multivariable systems.
Circuits Syst Signal Process (2014) 33:2971–2986
2985
5 Conclusions This paper studies the parameter estimation problem of multivariable systems. Based on the maximum likelihood principle and the subsystem decomposition, a recursive maximum likelihood estimation algorithm is derived to identify the parameter vector of the system. The simulation results show the effectiveness of the proposed algorithm. The method in this paper can be extended to study identification problems of other linear scalar of multivariable systems with colored noises [5,16] or nonlinear systems with colored noises [6,24,30,37]. Acknowledgments This work was supported by the National Natural Science Foundation of China (61273024, 61305031, 51307089), the Nantong Science and Technology Project (BK2012060), the Industrialization Project for Colleges and Universities of Jiangsu Province (JHB2012-44).
References 1. F. Ding, K.P. Deng, X.M. Liu, Decomposition based Newton iterative identification method for a Hammerstein nonlinear FIR system with ARMA noise. Circuits Syst. Signal Process. 33(x) (2014). doi:10.1007/s00034-014-9772-y 2. F. Ding, Coupled-least-squares identification for multivariable systems. IET Control Theory Appl. 7(1), 68–79 (2013) 3. F. Ding, Combined state and least squares parameter estimation algorithms for dynamic systems. Appl. Math. Model. 38(1), 403–412 (2014) 4. F. Ding, Hierarchical parameter estimation algorithms for multivariable systems using measurement. Info. Sci. (2014). http://dx.doi.org/10.1016/j.ins.2014.02.103 5. F. Ding, State filtering and parameter identification for state space systems with scarce measurements. Signal Process. (2014). doi:10.1016/j.sigpro.2014.03.031 6. F. Ding, X.P. Liu, G. Liu, Identification methods for Hammerstein nonlinear systems. Digit. Signal Process. 21(2), 215–238 (2011) 7. R. Diversi, R. Guidorzi, U. Soverini, Maximum likelihood identification of noisy input–output models. Automatica 43(3), 464–472 (2007) 8. W. Favoreel, B.D. Moor, P.V. Overschee, Subspace state space system identification for industrial processes. J. Process Control 10(20–3), 149–155 (2000) 9. S. Gibson, B. Ninness, Robust maximum-likelihood estimation of multivariable dynamic systems. Automatica 41(10), 1667–1682 (2005) 10. R. Guidorzi, Canonical structures in the identification of multivariable systems. Automatica 11, 261– 374 (1975) 11. R. Guidorzi, R. Diversi, Minimal representations of MIMO time-varying systems and realization of cyclostationary models. Automatica 39(11), 1903–1914 (2003) 12. J.H. Li, Parameter estimation for Hammerstein CARARMA systems based on the Newton iteration. Appl. Math. Lett. 26(1), 91–96 (2013) 13. J.H. Li, F. Ding, Maximum likelihood stochastic gradient estimation for Hammerstein systems with colored noise based on the key term separation technique. Comput. Math. Appl. 62(11), 4170–4177 (2011) 14. Y. Lin, Y. Shi, R. Burton, Modeling and robust discrete-time sliding mode control design for a fluid power electrohydraulic actuator (EHA) system. IEEE/ASME Trans. Mechatron. 18(1), 1–10 (2013) 15. Y.J. Liu, F. Ding, Convergence properties of the least squares estimation algorithm for multivariable systems. Appl. Math. Model. 37(1–2), 476–483 (2013) 16. Y.J. Liu, F. Ding, Y. Shi, An efficient hierarchical identification method for general dual-rate sampleddata systems. Automatica 50(3), 962–973 (2014) 17. Y.J. Liu, F. Ding, Y. Shi, Least squares estimation for a class of non-uniformly sampled systems based on the hierarchical identification principle. Circuits Syst. Signal Process. 31(6), 1985–2000 (2012) 18. Y.J. Liu, Y.S. Xiao, X.L. Zhao, Multi-innovation stochastic gradient algorithm for multiple-input single-output systems using the auxiliary model. Appl. Math. Computat. 215(4), 1477–1483 (2009)
2986
Circuits Syst Signal Process (2014) 33:2971–2986
19. L. Ljung, System Identification: Theory for the User, 2nd edn. (Prentice-Hall, Englewood Cliffs, NJ, 1999) 20. M. Lovera, T. Gustafsson, M. Verhaegen, Recursive subspace identification of linear and non-linear Wiener state–space models. Automatica 36(11), 1639–1650 (2000) 21. J.X. Ma, F. Ding, Recursive relations of the cost functions for the least squares algorithms for multivariable systems. Circuits Syst. Signal Process. 32(1), 83–101 (2013) 22. G. Mercère, L. Bako, Parameterization and identification of multivariable state-space systems: a canonical approach. Automatica 47(8), 1547–1555 (2011) 23. Y. Shi, H.Z. Fang, Kalman filter-based identification for systems with randomly missing measurements in a network environment. Int. J. Control 83(3), 538–551 (2010) 24. B. Sun, D.Q. Zhu, S.X. Yang, A bio-inspired filtered backstepping cascaded tracking control of 7000 m manned submarine vehicle. IEEE Trans. Ind. Electron. 61(7), 3682–3692 (2014) 25. D.N. Vizireanu, A fast, simple and accurate time-varying frequency estimation method for single-phase electric power systems. Measurement 45(5), 1331–1333 (2012) 26. D.N. Vizireanu, S.V. Halunga, Single sine wave parameters estimation method based on four equally spaced samples. Int. J. Electron. 98(7), 941–948 (2011) 27. J. Vörös, Modeling and identification of systems with backlash. Automatica 46(2), 369–374 (2010) 28. D.Q. Wang, F. Ding, Y.Y. Chu, Data filtering based recursive least squares algorithm for Hammerstein systems using the key-term separation principle. Inf. Sci. 222, 203–212 (2013) 29. D.Q. Wang, R. Ding, X.Z. Dong, Iterative parameter estimation for a class of multivariable systems based on the hierarchical identification principle and the gradient search. Circuits Syst. Signal Process. 31(6), 2167–2177 (2012) 30. D.Q. Wang, F. Ding, X.M. Liu, Least squares algorithm for an input nonlinear system with a dynamic subspace state space model. Nonlinear Dyn. 75(1–2), 49–61 (2014) 31. T. Wigren, Recursive prediction error identification and scaling of non-linear state space models using a restricted black box parameterization. Automatica 42(1), 159–168 (2006) 32. A. Wills, B. Ninness, On gradient-based search for multivariable system estimates. IEEE Trans. Automat. Control 53(1), 298–306 (2008) 33. A. Wills, B. Ninness, S. Gibson, Maximum likelihood estimation of state space models from frequency domain data. IEEE Trans. Automat. Control 54(1), 19–33 (2009) 34. A. Wills, T.B. Schä, L. Ljung, B. Ninness, Identification of Hammerstein-Wiener models. Automatica 45(1), 70–81 (2013) 35. W.X. Zheng, Least-squares identification of a class of multivariable systems with correlated disturbances. J. Frankl. Inst. 336(8), 1309–1324 (1999) 36. J. Zheng, K.W.K. Lui, W.-K. Ma, H.C. So, Two simplified recursive Gauss-Newton algorithms for direct amplitude and phase tracking of a real sinusoid. IEEE Signal Proc. Let. 14(12), 984–987 (2007) 37. D.Q. Zhu, H. Huang, S.X. Yang, Dynamic task assignment and path planning of multi-AUV system based on an improved self-organizing map and velocity synthesis method in 3D underwater workspace. IEEE Trans. Cybern. 43(2), 504–514 (2013) 38. A. Žilinskas, J. Žilinskas, A global optimization method based on the reduced simplicial statistical model. Math. Model. Anal. 16(3), 451–460 (2011)
Maximum Likelihood Recursive Least Squares Estimation for Multivariable Systems Junhong Li · Feng Ding · Ping Jiang · Daqi Zhu
Received: 9 April 2013 / Revised: 19 March 2014 / Accepted: 19 March 2014 / Published online: 16 April 2014 © Springer Science+Business Media New York 2014
Abstract This paper discusses parameter estimation problems of the multivariable systems described by input–output difference equations. We decompose a multivariable system to several subsystems according to the number of the outputs. Based on the maximum likelihood principle, a maximum likelihood-based recursive least squares algorithm is derived to estimate the parameters of each subsystem. Finally, two numerical examples are provided to verify the effectiveness of the proposed algorithm. Keywords Recursive identification · Multivariable systems · Parameter estimation · Maximum likelihood 1 Introduction The mathematical models can describe the relationships among the system variables [1,17,28]. Mathematical modeling plays an important role in control [14,23,27] J. Li (B) · P. Jiang School of Electrical Engineering, Nantong University, Nantong 226019, People’s Republic of China e-mail: [email protected] P. Jiang e-mail: [email protected] F. Ding School of Internet of Things Engineering, Jiangnan University, Wuxi 214122, People’s Republic of China e-mail: [email protected] D. Zhu Laboratory of Underwater Vehicles and Intelligent Systems, Shanghai Maritime University, Shanghai 201306, People’s Republic of China e-mail: [email protected]
2972
Circuits Syst Signal Process (2014) 33:2971–2986
and other engineering fields [25,26]. Many practical systems in industrial processes are multivariable systems. Compared with the single-input single-output systems, the multivariate systems have many input and output variables [21], the problem of building their mathematical model is more difficult than the single-input single-output case. Different mathematical models have been used to describe the multivariable systems, such as the state-space model, the transfer function matrix model [19]. Some earlier work for identifying multivariable systems focused on the state-space model [8,20,31]. In this area, Guidorzi studied the identification problem of multivariable systems and derived the connections between the state–space model and the input– output difference equation canonical model [10,11], his work provided an effective solution to identify multivariable systems. Recently, many modeling and estimation methods have been proposed to estimate the parameters of multivariable systems [15,18]. For example, Wang et al. [29] presented the hierarchical gradient-based iterative algorithm for multivariable controlled autoregressive system with autoregressive noises using the hierarchical identification principle and the gradient search; Mercère and Bako [22] studied a subspace-based identification algorithm for multivariable state-space systems; Wills et al. [32] proposed gradient-based search algorithms for multivariable systems; Zheng [35] proposed a bias eliminated least squares for multiple-input single-output systems. The maximum likelihood principle is a useful tool in statistics [38] and system identification [12,34]. In the area of maximum likelihood identification, Diversi et al. [7] proposed a maximum likelihood method for the identification of errors-in-variables models corrupted by white and uncorrelated Gaussian noises; Gibson and Ninness [9] developed a robust maximum likelihood method for linear state-space models; Wills et al. [33] proposed a maximum likelihood estimation algorithm of state-space models. Recently, Li and Ding [13] studied maximum likelihood estimation algorithms for Hammerstein nonlinear systems. Iterative or recursive algorithms are important for modeling a system [2–4,36]. Based on the work in [13], the objective of this paper is to study the recursive parameter estimation of multivariable input–output difference equation model using the maximum likelihood principle. The basic idea is to decompose a multivariable system to several multiple-input single-output subsystems, and to derive the recursive maximum likelihood algorithm for each subsystem. The rest of this paper is organized as follows. Section 2 gives the system description of the multivariable systems. Section 3 derives a maximum likelihood-based recursive least squares algorithm for multivariable systems. Section 4 provides two numerical examples to illustrate the proposed method. Finally, we offer some concluding remarks in Sect. 5.
2 The System Description Consider a multivariable system described by the following input–output representation [10]:
Circuits Syst Signal Process (2014) 33:2971–2986
2973
P(z) y(t) = Q(z)u(t) + R(z)v(t), ⎡
p11 (z) p12 (z) · · · p1m (z)
(1)
⎤
⎢ ⎥ ⎢ p21 (z) p22 (z) · · · p2m (z) ⎥ ⎢ ⎥ P(z) := ⎢ . ⎥ ∈ Rm×m , .. .. ⎢ .. ⎥ . . ⎣ ⎦ pm1 (z) pm2 (z) · · · pmm (z) ⎡ ⎤ q11 (z) q12 (z) · · · q1r (z) ⎢ ⎥ ⎢ q21 (z) q22 (z) · · · q2r (z) ⎥ ⎢ ⎥ ∈ Rm×r , Q(z) := ⎢ . .. .. ⎥ ⎢ .. ⎥ . . ⎦ ⎣ qm1 (z) qm2 (z) · · · qmr (z) ⎡ ⎤ d11 (z) d12 (z) · · · d1m (z) ⎢ ⎥ ⎢ d21 (z) d22 (z) · · · d2m (z) ⎥ ⎢ ⎥ ∈ Rm×m , R(z) := ⎢ . .. .. ⎥ ⎢ .. ⎥ . . ⎦ ⎣ dm1 (z) dm2 (z) · · · dmm (z) where u(t) := [u 1 (t), u 2 (t), . . . , u r (t)]T ∈ Rr and y(t) := [y1 (t), y2 (t), . . . , ym (t)]T ∈ Rm are the system input vector and output vector, respectively, v(t) := [v1 (t), v2 (t), . . . , vm (t)]T ∈ Rm is a Gaussian distributed white noise vector with zero mean and variance (σ12 , σ22 , . . . , σm2 ), z is a unit forward shift operator: z y(t) = y(t + 1), and P(z), Q(z), and R(z) are polynomial matrices whose entries are defined as [10] pi j (z) = δi j z n i j − ai j,n i j z n i j −1 − · · · − ai j,2 z − ai j,1 , i, j = 1, 2, . . . , m, (2) qi j (z) = bi j,n i z n i −1 + bi j,n i −1 z n i −2 · · · + bi j,2 z + bi j,1 , i = 1, 2, . . . , m, j = 1, 2, . . . , r, di j (z) = δi j z n i j + di j,n i j z n i j −1 + · · · + di j,2 z + di j,1 , i, j = 1, 2, . . . , m,
ni j
⎧ i= j ⎪ ⎨ ni , 1, i = j i < j , δi j = = min[n i , n j ], . ⎪ 0, i = j ⎩ min[n i + 1, n j ], i > j
(3) (4)
(5)
The ith subsystem in (1) can be equivalently written as m
j=1
pi j (z)y j (t) =
r
j=1
qi j (z)u j (t) +
m
j=1
di j (z)v j (t), i = 1, 2, . . . , m.
(6)
2974
Circuits Syst Signal Process (2014) 33:2971–2986
That is yi (t + n i ) =
m
[ai j,n i j y j (t + n i j − 1)+ai j,n i j −1 y j (t+n i j − 2) + · · · + ai j,1 y j (t)] j=1 r
+ [bi j,n i u j (t + n i − 1)+bi j,n i −1 u j (t + n i − 2)+· · ·+bi j,1 u j (t)] j=1 m
+ [di j,n i j v j (t + n i j − 1) + di j,n i j −1 v j (t + n i j − 2) + · · · i=1
+ di j,1 vi (t)] + vi (t + n i ). Define the parameter vector of subsystem i as θ iT := [ai1 , ai2 , . . . , aim , bi1 , bi2 , . . . , bir , d i1 , d i2 , . . . , d im ] ∈ R1×(2μi +r n i ) , μi := n i1 + n i2 + · · · + n im , ai j := [ai j,n i j , ai j,n i j −1 , . . . , ai j,1 ] ∈ R1×n i j , bi j := [bi j,n i , bi j,n i −1 , . . . , bi j,1 ] ∈ R1×n i , d i j := [di j,n i j , di j,n i j −1 , . . . , di j,1 ] ∈ R1×n i j , ϕ i (t) = [ ¯yi1 (t), ¯yi2 (t), . . . , ¯yim (t), u¯ i1 (t), u¯ i2 (t), . . . , u¯ ir (t), v¯ i1 (t), v¯ i2 (t), . . . , v¯ im (t)]T ∈ R(2μi +r n i ) , ¯yi j (t) := [y j (t − 1), y j (t − 2), . . . , y j (t − n i j )] ∈ R1×n i j , u¯ i j (t) := [u j (t − 1), u j (t − 2), . . . , u j (t − n i j )] ∈ R1×n i , v¯ i j (t) := [v j (t − 1), v j (t − 2), . . . , v j (t − n i j )] ∈ R1×n i j . Then the ith subsystem can be written as the following identification model: yi (t) = θ iT ϕ i (t) + vi (t).
(7)
The objective of this paper is to present new identification methods for determining the parameter vector θ i from the available input–output data ϕ i (t) and yi (t). 3 The Maximum Likelihood-Based Recursive Least Squares Algorithm Referring to the work in [13], we propose the maximum likelihood criterion function for the ith subsystem in (7) as follows: 1 2 vi (k). 2 t
J (θ i , t) =
(8)
k=1
Let θˆ i (t) denote the maximum likelihood estimate of θ i at time t. The Taylor expansion of J (θ i , t) at point θ i = θˆ i (t − 1) is given by
Circuits Syst Signal Process (2014) 33:2971–2986
2975
∂ J (θˆ i (t − 1), t) J (θ i , t) = J (θˆ i (t − 1), t) + [θ i − θˆ i (t − 1)] ∂θ i 1 ∂ 2 J (θˆ i (t −1), t) + [θ i − θˆ i (t −1)]T [θ i − θˆ i (t −1)]+◦(θ i − θˆ i (t −1)2 ). 2 ∂θ i θ iT Then the maximum likelihood estimate at time t can be achieved by the Newton– Raphson method −1 2 ˆ ∂ J (θˆ i (t − 1), t) ˆθ i (t) = θˆ i (t − 1) − ∂ J (θ i (t − 1), t) . T ∂θ i θ i ∂θ i In order to compute as a recursive form:
∂ 2 J (θˆ i (t−1),t) ∂θ i θ iT
and
∂ J (θˆ i (t−1),t) , ∂θ i
(9)
we rewrite the cost function in (8)
1 2 1 vi (k) = J (θ i , t − 1) + vi2 (t). 2 2 t
J (θ i , t) =
(10)
k=1
Computing the derivative of J (θ i , t) at point θˆ i (t − 1) gives ∂ J (θˆ i (t − 1), t) ∂ J (θˆ i (t − 1), t − 1) ∂vi (t) = + vˆi (t) . ∂θ i ∂θ i ∂θ i θˆi (t−1)
(11)
Let ϕˆ i f (t) := − Then we have
∂vi (t) ∈ R(2μi +r n i )×(2μi +r n i ) . ∂θ i θˆi (t−1)
∂ J (θˆ i (t − 1), t) = −vi (t)ϕˆ i f (t). ∂θ i
(12)
Computing the second-order partial derivative of J (θ i , t) with respect to θ i gives ∂vi2 (t) ∂ 2 J (θˆ i (t − 1), t) ∂ 2 J (θˆ i (t − 1), t − 1) T = + ϕˆ i f (t)ϕˆ i f (t) + vi (t) ∂θ i θ iT ∂θ i ∂θ iT ∂θ i ∂θ iT θˆ i (t−1) ≈
∂ 2 J (θˆ i (t − 1), t − 1) + ϕˆ i f (t)ϕˆ iT f (t). ∂θ i ∂θ iT
Let P i−1 (t) :=
∂ 2 J (θˆ i (t − 1), t) . ∂θ i ∂θ iT
(13)
2976
Circuits Syst Signal Process (2014) 33:2971–2986
Equation (13) can be written as P i−1 (t) = P i−1 (t − 1) + ϕˆ i f (t)ϕˆ iT f (t).
(14)
Applying the matrix inversion formula ( A + BC)−1 = A−1 − A−1 B(I + C A−1 B)−1 C A−1 to (14) gives P i (t) = P i (t − 1) −
P i (t − 1)ϕˆ i f (t)ϕˆ iT f (t) P i (t − 1) 1 + ϕˆ iT f (t) P i (t − 1)ϕˆ i f (t)
.
(15)
Define the gain vector L i (t) : = P i (t)ϕˆ i f (t) =
P i (t − 1)ϕˆ i f (t) 1 + ϕˆ iT f (t) P i (t − 1)ϕˆ i f (t)
∈ R(2μi +r n i ) .
(16)
Thus, Eq. (15) can be rewritten as P i (t) = [I − L i (t)ϕˆ iT f (t)] P i (t − 1).
(17)
In the following, we compute the filtered information vector ϕˆ i f (t) of the ith subsystem. Let aˆ i j (t), bˆ i j (t), and dˆ i j (t) be the estimates of ai j , bi j , and d i j at time t, i.e., aˆ i j (t) := [aˆ i j,n i j (t), aˆ i j,n i j −1 (t), . . . , aˆ i j,1 (t)] ∈ R1×n i j , bˆ i j (t) := [bˆi j,n i (t), bˆi j,n i −1 (t), . . . , bˆi j,1 (t)] ∈ R1×n i , dˆ i j (t) := [dˆi j,n i j (t), dˆi j,n i j −1 (t), . . . , dˆi j,1 (t)] ∈ R1×n i j , dˆii (t, z) := z n i + dˆii,n i (t)z n i −1 + dˆii,n i −1 (t)z n i −2 + · · · + dˆii,1 (t). vi (t) in (7) can be equivalently written as vi (t) = yi (t) − θ iT ϕ i (t) ⎡ ⎤ m r m
= yi (t) − ⎣ ai j ¯yiTj (t) + bi j u¯ iTj (t) + d i j v¯ iTj (t)⎦ . j=1
j=1
j=1
(18)
Circuits Syst Signal Process (2014) 33:2971–2986
2977
The partial derivatives of vi (t) with respect to ai j can be computed by ∂ v¯ iT (t) ∂vi (t) ˆ = − ¯yi j (t) − d ii (t − 1) ∂ ai j θˆi (t−1) ∂ ai j θˆi (t−1) = − ¯yi j (t) − [dˆii,n i (t − 1), . . . , dˆii,2 (t − 1), dˆii,1 (t − 1)] ∂vi (t − 1) ∂vi (t − n i + 1) ∂vi (t − n i ) T × ,..., , ∂ ai j ∂ ai j ∂ ai j θˆ i (t−1)
−1 −n i ∂vi (t) ˆ ˆ = −¯yi j (t)−[dii,n i (t − 1)z +· · · + dii,1 (t − 1)z ] . ∂ ai j θˆi (t−1) Then we have ∂vi (t) 1 ¯yi j (t) =: − ˆ¯yi j, f (t). =− ∂ ai j θˆi (t−1) z −n i dˆii (t − 1, z) Similarly, the partial derivatives of vi (t) with respect to bi j and d i j can be written as follows: ∂vi (t) 1 u¯ i j (t) =: −uˆ¯ i j, f (t), =− −n ˆ ∂ bi j θˆi (t−1) z i dii (t − 1, z) ∂vi (t) 1 =− v¯ˆ i j (t) =: −v¯ˆ i j, f (t). −n ˆ i ∂ d i j θˆ (t−1) z dii (t − 1, z) i
Then the filtered information vectors ˆ¯yi j, f (t), uˆ¯ i j, f (t) and vˆ¯ i j, f (t) can be computed by ˆ¯yi j, f (t) = ¯yi j (t) − dˆii,n i (t − 1) ˆ¯yi j, f (t − 1) − · · · − dˆii,1 (t − 1) ˆ¯yi j, f (t − n i ), u¯ˆ i j, f (t) = u¯ i j (t) − dˆii,n i (t − 1)u¯ˆ i j, f (t − 1) − · · · − dˆii,1 (t − 1)uˆ¯ i j, f (t − n i ), v¯ˆ i j, f (t) = vˆ¯ i j (t) − dˆii,n i (t − 1)vˆ¯ i j, f (t − 1) − · · · − dˆii,1 (t − 1)vˆ¯ i j, f (t − n i ). Thus the filtered information vector can be written as ∂vi (t) ϕˆ i f (t) := − ∂θ i θˆi (t−1) ∂vi (t) ∂vi (t) ∂vi (t) ∂vi (t) ∂vi (t) = − , ,..., , , ,..., ∂ ai1 ∂ ai2 ∂ aim ∂ bi1 ∂ bi2 ∂vi (t) ∂vi (t) ∂vi (t) ∂vi (t) T , , ,..., ∂ bir ∂ d i1 ∂ d i2 ∂ d im = [ ¯ˆyi1, f (t), ¯ˆyi2, f (t), . . . , ¯ˆyim, f (t), u¯ˆ i1, f (t), u¯ˆ i2, f (t), . . . , uˆ¯ ir, f (t), vˆ¯ i1, f (t), vˆ¯ i2, f (t), . . . , vˆ¯ im, f (t)]T .
2978
Circuits Syst Signal Process (2014) 33:2971–2986
Let vˆ¯ i j (t) be the estimates of v¯ i j (t) at time t, and let ϕˆ i (t) denote the information vector obtained by replacing v¯ i j (t) with vˆ¯ i j (t), i.e., ϕˆ i (t) = [ ¯yi1 (t), ¯yi2 (t), . . . , ¯yim (t), u¯ i1 (t), u¯ i2 (t), . . . , u¯ ir (t), vˆ¯ i1 (t), vˆ¯ i2 (t), . . . , v¯ˆ im (t)]T ∈ R(2μi +r n i ) , vˆ¯ i j (t) = [vˆ j (t − 1), vˆ j (t − 2), . . . , vˆ j (t)] ∈ R1×n i j . According to the identification model in (7), we have vi (t) = yi (t) − θ iT ϕ i (t). Replacing ϕ i (t) and θ i with their estimates ϕˆ i (t) and θˆ i (t − 1) gives the estimates of vi (t) T vˆi (t) = yi (t) − θˆ i (t − 1)ϕˆ i (t).
According to (9), (12), (16), and (17), we can summarize the maximum likelihoodbased recursive least squares (ML-RLS) algorithm for the ith subsystem of the multivariable system as θˆ i (t) = θˆ i (t − 1) + L i (t)vˆi (t), P i (t − 1)ϕˆ i f (t) , L i (t) = 1 + ϕˆ iT f (t) P i (t − 1)ϕˆ i f (t)
(19) (20)
P i (t) = [I − L i (t)ϕˆ iT f (t)] P i (t − 1), P i (0) = p0 I,
(21)
vˆi (t) = yi (t) − θˆ i (t − 1)ϕˆ i (t), ϕˆ i (t) = [ ¯yi1 (t), ¯yi2 (t), . . . , ¯yim (t), u¯ i1 (t), u¯ i2 (t),
(22)
T
. . . , u¯ ir (t), v¯ˆ i1 (t), vˆ¯ i2 (t), . . . , vˆ¯ im (t)]T , ϕˆ i f (t) = [ ˆ¯yi1, f (t), ˆ¯yi2, f (t), . . . , ˆ¯yim, f (t), uˆ¯ i1, f (t), uˆ¯ i2, f (t), . . . , uˆ¯ ir, f (t), vˆ¯ i1, f (t), vˆ¯ i2, f (t), . . . , vˆ¯ im, f (t)]T , ¯yi j (t) = [y j (t − 1), y j (t − 2), . . . , y j (t − n i )], u¯ i j (t) = [u j (t − 1), u j (t − 2), . . . , u j (t − n i )], vˆ¯ i j (t) = [vˆ j (t − 1), vˆ j (t − 2), . . . , vˆ j (t − n i )], ˆ¯yi j, f (t) = ¯yi j (t)− dˆii,n i (t −1) ˆ¯yi j, f (t − 1)−· · ·− dˆii,1 (t −1) ˆ¯yi j, f (t −n i ),
(23)
(24) (25) (26) (27) (28)
u¯ˆ i j, f (t) = u¯ i j (t)− dˆii,n i (t −1)u¯ˆ i j, f (t −1)−· · ·− dˆii,1 (t −1)uˆ¯ i j, f (t −n i ), (29) v¯ˆ i j, f (t) = vˆ¯ i j (t)− dˆii,n i (t −1)vˆ¯ i j, f (t − 1) − · · · − dˆii,1 (t −1)vˆ¯ i j, f (t −n i ). (30) The steps involved in the ML-RLS algorithm for subsystem i are listed in the following.
Circuits Syst Signal Process (2014) 33:2971–2986
2979
Fig. 1 The flowchart for computing the ML-RLS estimate θˆ i (t)
1. Let t = 1, set the initial values θˆ i (0) = 1(2μi +r n i ) / p0 and P i (0) = p0 I, ˆ¯yii, f (t) = 0, uˆ¯ i j, f (t) = 0, vˆ¯ ii, f (t) = 0 and v¯ i j (t) = 0 for t ≤ 0, p0 is a large number (i.e., p0 = 106 ). 2. Collect the input-output data u(t) and y(t), form ¯yi j (t) by (25), u¯ i j (t) by (26), vˆ¯ i j (t) by (27), form ϕˆ i (t) by (23). 3. Compute ˆ¯yi j, f (t), uˆ¯ i j, f (t), and vˆ¯ i j, f (t) by (28), (29), and (30), respectively, form ϕˆ i f (t) by (24). 4. Compute L i (t) by (20), compute P i (t) by (21) and vˆi (t) by (22). 5. Update the parameter estimate θˆ i (t) by (19). 6. Increase t by 1 and go to step 2. The flowchart of computing the subsystem parameter estimate θˆ i (t) in the ML-RLS algorithms is shown in Fig. 1. A summary of the variables and their dimensions in the ML-RLS algorithm are shown in Table 1 for convenience. The numbers of multiplications and additions in each recursive computation of the ML-RLS algorithm are shown in Table 2.
2980
Circuits Syst Signal Process (2014) 33:2971–2986
Table 1 The dimensions of the variables in the ML-RLS algorithm Item
Variables
Dimensions
1
Input variable
u(t) ∈ Rr
2
Output variable
y(t) ∈ Rm
3
Parameter vectors
ai j ∈ R1×n i j , bi j ∈ R1×n i , d i j ∈ R1×n i j
4
Parameter vector
θ i ∈ R(2μi +r n i )
5
Information vectors
¯yi j (t) ∈ R1×n i j , u¯ i j (t) ∈ R1×n i , v¯ i j (t) ∈ R1×n i j
6
Information vector
ϕ i (t) ∈ R(2μi +r n i )
5
Filtered information vectors
¯yi j, f (t) ∈ R1×n i j , u¯ i j, f (t) ∈ R1×n i , v¯ i j, f (t) ∈ R1×n i j
6
Filtered information vector
ϕ i f (t) ∈ R(2μi +r n i )
7
Covariance matrix
P i (t) ∈ R(2μi +r n i )×(2μi +r n i )
Table 2 The computational efficiency of the ML-RLS algorithm Variables
Number of multiplications
Number of additions
θˆ i (t) = θˆ i (t − 1) + L i (t)vˆi (t)
2μi + r n i
2μi + r n i
T vˆi (t) = yi (t) − θˆ i (t − 1)ϕˆ i (t)
2μi + r n i
2μi + r n i
L i (t) = ζ i (t)/[1 + ϕ iTf (t)ζ i (t)] 2(2μi + r n i ) (2μi + r n i )2
ζ i (t) := P i (t − 1)ϕ i f (t)
2μi + r n i (2μi + r n i − 1)(2μi + r n i )
P i (t) = P i (t − 1) − L i (t)ζ iT (t) (2μi + r n i )2
(2μi + r n i )2
ˆ¯yi j, f (t) = ¯yi j (t) − · · · − dˆii,1 (t − 1) ˆ¯yi j, f (t − n i )
n i2
n i2
uˆ¯ i j, f (t) = u¯ i j (t) − · · · − dˆii,1 (t − 1)uˆ¯ i j, f (t − n i ) v¯ˆ i j, f (t) = v¯ˆ i j (t) − · · · − dˆii,1 (t − 1)v¯ˆ i j, f (t − n i )
n i2
n i2
n i2
n i2
Sum
2(2μi + r n i )2 + 4(2μi +r n i )+3n i2
2(2μi + r n i )2 + 2(2μi + r n i ) + 3n i2
Total flops
4(2μi + r n i )2 + 6(2μi + r n i ) + 6n i2
In order to show the advantages of the ML-RLS algorithm, we give the RELS algorithm for comparison [19]: T θˆ i (t) = θˆ i (t − 1) + L i (t)[yi (t) − θˆ i (t − 1)ϕˆ i (t)], θˆ (0) = 12μi +r n i / p0 , P i (t − 1)ϕˆ i (t) L i (t) = , 1 + ϕˆ iT (t) P i (t − 1)ϕˆ i (t) P i (t) = [I − L i (t)ϕˆ iT (t)] P i (t − 1), P i (0) = p0 I 2μi +r n i , T vˆi (t) = yi (t) − θˆ (t)ϕˆ i (t),
i
ϕˆ i (t) = [ ¯yi1 (t), ¯yi2 (t), . . . , ¯yim (t), u¯ i1 (t), u¯ i2 (t), . . . , u¯ ir (t), vˆ¯ i1 (t), vˆ¯ i2 (t), . . . , vˆ¯ im (t)]T ,
Circuits Syst Signal Process (2014) 33:2971–2986
2981
¯yi j (t) = [y j (t − 1), y j (t − 2), . . . , y j (t − n i )], u¯ i j (t) = [u j (t − 1), u j (t − 2), . . . , u j (t − n i )], vˆ¯ i j (t) = [vˆ j (t − 1), vˆ j (t − 2), . . . , vˆ j (t − n i )]. The numbers of multiplications and additions in each recursive computation of the RELS algorithm are shown in Table 3. Remark Similarly to the RELS algorithm, the ML-RLS algorithm is a recursive algorithm and can be used for online identification. The two algorithms have almost same computational burden, but the ML-RLS algorithm can achieve higher parameter estimation accuracy with data length increases. 4 Numerical Examples Example 1 Consider the multivariable system in (1), where
z −0.18 0.50 0.80 0.49 z − 0.30 0.66 P(z) = , Q(z) = , R(z) = , 0.01 z + 0.19 0.30 0.50 0.70 z − 0.13 u 1 (t) v (t) y1 (t) , u(t) = , v(t) = 1 . y(t) = y2 (t) u 2 (t) v2 (t)
The parameter vectors of this system are θ 1 = [a11,1 , a12,1 , b11,1 , b12,1 , d11,1 , d12,1 ]T = [0.18,−0.50, 0.80, 0.49,−0.30, 0.66]T , θ 2 = [a21,1 , a22,1 , b21,1 , b22,1 , d21,1 , d22,1 ]T = [−0.01,−0.19, 0.30, 0.50, 0.70,−0.13]T, θ1 . = θ2 In simulation, the inputs {u 1 (t)} and {u 2 (t)} are taken as two uncorrelated stochastic signal sequences with zero mean and unit variance, and the noise {v1 (t)} and {v2 (t)} as two Gaussian white noise sequences with zero mean and variances σ12 = 0.502 and σ22 = 0.302 , respectively. Applying the ML-RLS algorithm and the RELS algorithm to estimate the parameters of this system, the parameter estimates and their errors ˆ δ := (t) − / versus t are shown in Tables 4 and 5, the parameter estimates versus t are shown in Figs. 2 and 3 for the two identification algorithms. From Tables 4, 5 and Figs. 2, 3, we can draw the following conclusions. 1. The parameter estimation error of the two algorithms gradually becomes small as t increases. 2. The parameter estimates of the two algorithms converge to their true values. This verifies the effectiveness of the two algorithms. 3. The RELS algorithm has better performance for small data length t, but for large t, the ML-RLS algorithm can obtain higher parameter estimation accuracy.
2982
Circuits Syst Signal Process (2014) 33:2971–2986
Table 3 The computational efficiency of the recursive extended least squares (RELS) algorithm Variables
Number of multiplications
Number of additions
θˆ i (t) = θˆ i (t − 1) + T L i (t)[yi (t)−θˆ i (t −1)ϕˆ i (t)]
2(2μi + r n i )
2(2μi + r n i )
L i (t) = ζ i (t)/[1 + ϕˆ iT (t)ζ i (t)]
2(2μi + r n i )
2μi + r n i
ζ i (t) := P i (t − 1)ϕˆ i (t)
(2μi + r n i − 1)(2μi + r n i )
P i (t) = P i (t − 1) − L i (t)ζ iT (t)
(2μi + r n i )2 (2μi + r n i )2
vˆi (t) = yi (t) − θˆ i (t)ϕˆ i (t)
2μi + r n i
2μi + r n i
Sum
2(2μi + r n i )2 + 5(2μi + r n i )
2(2μi + r n i )2 + 3(2μi + r n i )
Total flops
4(2μi + r n i )2 + 8(2μi + r n i )
T
(2μi + r n i )2
Table 4 The ML-RLS estimates and errors of Example 1 t
100
200
500
1,000
2,000
3,000
True values
a11,1
0.25088
0.15233
0.18478
0.19745
0.19150
0.18123
0.18000
a12,1
−0.47796
−0.44132
−0.48214
−0.49723
−0.51678
−0.50031
−0.50000
b11,1
0.78764
0.77093
0.79562
0.79351
0.80069
0.80793
0.80000
b12,1
0.55179
0.51400
0.52032
0.51363
0.50714
0.50575
0.49000
d11,1
−0.33228
−0.29851
−0.34230
−0.32844
−0.30864
−0.29515
−0.30000
d12,1
0.27637
0.27588
0.48525
0.59703
0.66636
0.66076
0.66000
a21,1
−0.01866
−0.00335
0.01973
−0.00420
0.00255
−0.00424
−0.01000
a22,1
−0.17889
−0.22734
−0.24365
−0.20447
−0.20723
−0.19312
−0.19000
b21,1
0.32892
0.31850
0.31089
0.30847
0.30652
0.30498
0.30000
b22,1
0.46131
0.48062
0.51007
0.50947
0.49683
0.49780
0.50000
d21,1
0.50289
0.58396
0.62662
0.67795
0.69264
0.70257
0.70000
d22,1
−0.18426
−0.07659
0.00862
−0.07252
−0.09546
−0.12559
−0.13000
δ (%)
28.03052
25.83780
15.55223
6.19942
3.16590
1.30100
Example 2 Consider the multivariable system in (1), where ⎡
⎤ ⎡ ⎤ z − 0.66 0.26 −0.30 0.82 0.49 0.13 0.29 z − 0.49 −0.25 ⎦ , Q(z) = ⎣ 0.35 z − 0.27 −0.25 ⎦ , P(z) = ⎣ −0.29 0.89 z − 0.66 0.97 0.89 0.85 ⎡ ⎤ z − 0.23 0.21 0.23 0.11 z − 0.10 −0.10 ⎦ , R(z) = ⎣ −0.15 0.23 z + 0.13 ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ y1 (t) u 1 (t) v1 (t) y(t) = ⎣ y2 (t) ⎦ , u(t) = ⎣ u 2 (t) ⎦ , v(t) = ⎣ v2 (t) ⎦ . y3 (t) u 3 (t) v3 (t) The parameter vectors of the system are
Circuits Syst Signal Process (2014) 33:2971–2986
2983
Table 5 The RELS estimates and errors of Example 1 t
100
200
500
1,000
2,000
3,000
True values
a11,1
0.26979
0.16553
0.17856
0.19388
0.18532
0.17338
0.18000
a12,1
−0.55086
−0.47729
−0.50521
−0.51738
−0.52678
−0.50455
−0.50000 0.80000
b11,1
0.79906
0.77303
0.79541
0.79138
0.80015
0.80736
b12,1
0.54587
0.51208
0.51914
0.51303
0.50521
0.50466
0.49000
d11,1
−0.29052
−0.28952
−0.33279
−0.32531
−0.30361
−0.28987
−0.30000
d12,1
0.31910
0.31184
0.51463
0.63200
0.68262
0.67119
0.66000
a21,1
−0.03664
−0.01313
0.01460
−0.00683
0.00166
−0.00376
−0.01000
a22,1
−0.13035
−0.19882
−0.23128
−0.19590
−0.20203
−0.19031
−0.19000
b21,1
0.32834
0.31529
0.30729
0.30738
0.30564
0.30461
0.30000
b22,1
0.45357
0.47940
0.50894
0.50841
0.49660
0.49758
0.50000
d21,1
0.51231
0.59512
0.63365
0.68342
0.69598
0.70337
0.70000
d22,1
−0.19767
−0.08989
−0.00121
−0.08534
−0.10425
−0.13067
−0.13000
δ (%)
26.24540
23.05107
13.45366
4.39389
3.12132
1.57697
1
b11,1
Parameter estimates
0.8
d
12,1
0.6
b12,1
0.4
a11,1
0.2 0 −0.2
d
11,1
−0.4
a12,1
−0.6 0
500
1000
1500
2000
2500
3000
t Fig. 2 The ML-RLS estimates of subsystem 1 versus t in Example 1
θ 1 = [a11,1 , a12,1 , a13,1 , b11,1 , b12,1 , b13,1 , d11,1 , d12,1 , d13,1 ]T = [0.66, −0.26, 0.30, 0.82, 0.49, 0.13, −0.23, 0.21, 0.23]T , θ 2 = [a21,1 , a22,1 , a23,1 , b21,1 , b22,1 , b23,1 , d21,1 , d22,1 , d23,1 ]T = [−0.29, 0.49, 0.25, 0.35, −0.27, −0.25, 0.11, −0.10, −0.10]T , θ 3 = [a31,1 , a32,1 , a33,1 , b31,1 , b32,1 , b33,1 , d31,1 , d32,1 , d33,1 ]T = [0.29, −0.89, 0.66, 0.97, 0.89, 0.85, −0.15, 0.23, 0.13]T , ⎡ ⎤ θ1 = ⎣ θ2 ⎦ . θ3
2984
Circuits Syst Signal Process (2014) 33:2971–2986
Parameter estimates
0.8
d21,1
0.6
b
0.4
b
22,1 21,1
0.2
a
21,1
0
d22,1
−0.2
a22,1
−0.4 −0.6 0
500
1000
1500
2000
2500
3000
t Fig. 3 The ML-RLS estimates of subsystem 2 versus t in Example 1
0.25
0.2
δ
0.15
0.1
ML−RLS
RELS
0.05
0 0
500
1000
1500
2000
2500
3000
t Fig. 4 The ML-RLS and RELS estimation errors δ versus t in Example 2
In simulation, the inputs {u 1 (t)}, {u 2 (t)}, and {u 3 (t)} are taken as three uncorrelated stochastic signal sequence with zero mean and unit variance, and the noise {v1 (t)}, {v2 (t)}, and {v3 (t)} as three Gauss white noise sequence with zero mean and variances σ12 = σ22 = σ32 = 0.302 . Applying the ML-RLS algorithm and the RELS algorithm to estimate the parameters of this system, the parameter estimates and their errors are shown in Fig. 4. From Fig. 4, we can see that the parameter estimation error of the two algorithms gradually becomes small as t increases. The ML-RLS algorithm has higher parameter estimation accuracy than the RELS algorithm as t increases, thus can effectively identify multivariable systems.
Circuits Syst Signal Process (2014) 33:2971–2986
2985
5 Conclusions This paper studies the parameter estimation problem of multivariable systems. Based on the maximum likelihood principle and the subsystem decomposition, a recursive maximum likelihood estimation algorithm is derived to identify the parameter vector of the system. The simulation results show the effectiveness of the proposed algorithm. The method in this paper can be extended to study identification problems of other linear scalar of multivariable systems with colored noises [5,16] or nonlinear systems with colored noises [6,24,30,37]. Acknowledgments This work was supported by the National Natural Science Foundation of China (61273024, 61305031, 51307089), the Nantong Science and Technology Project (BK2012060), the Industrialization Project for Colleges and Universities of Jiangsu Province (JHB2012-44).
References 1. F. Ding, K.P. Deng, X.M. Liu, Decomposition based Newton iterative identification method for a Hammerstein nonlinear FIR system with ARMA noise. Circuits Syst. Signal Process. 33(x) (2014). doi:10.1007/s00034-014-9772-y 2. F. Ding, Coupled-least-squares identification for multivariable systems. IET Control Theory Appl. 7(1), 68–79 (2013) 3. F. Ding, Combined state and least squares parameter estimation algorithms for dynamic systems. Appl. Math. Model. 38(1), 403–412 (2014) 4. F. Ding, Hierarchical parameter estimation algorithms for multivariable systems using measurement. Info. Sci. (2014). http://dx.doi.org/10.1016/j.ins.2014.02.103 5. F. Ding, State filtering and parameter identification for state space systems with scarce measurements. Signal Process. (2014). doi:10.1016/j.sigpro.2014.03.031 6. F. Ding, X.P. Liu, G. Liu, Identification methods for Hammerstein nonlinear systems. Digit. Signal Process. 21(2), 215–238 (2011) 7. R. Diversi, R. Guidorzi, U. Soverini, Maximum likelihood identification of noisy input–output models. Automatica 43(3), 464–472 (2007) 8. W. Favoreel, B.D. Moor, P.V. Overschee, Subspace state space system identification for industrial processes. J. Process Control 10(20–3), 149–155 (2000) 9. S. Gibson, B. Ninness, Robust maximum-likelihood estimation of multivariable dynamic systems. Automatica 41(10), 1667–1682 (2005) 10. R. Guidorzi, Canonical structures in the identification of multivariable systems. Automatica 11, 261– 374 (1975) 11. R. Guidorzi, R. Diversi, Minimal representations of MIMO time-varying systems and realization of cyclostationary models. Automatica 39(11), 1903–1914 (2003) 12. J.H. Li, Parameter estimation for Hammerstein CARARMA systems based on the Newton iteration. Appl. Math. Lett. 26(1), 91–96 (2013) 13. J.H. Li, F. Ding, Maximum likelihood stochastic gradient estimation for Hammerstein systems with colored noise based on the key term separation technique. Comput. Math. Appl. 62(11), 4170–4177 (2011) 14. Y. Lin, Y. Shi, R. Burton, Modeling and robust discrete-time sliding mode control design for a fluid power electrohydraulic actuator (EHA) system. IEEE/ASME Trans. Mechatron. 18(1), 1–10 (2013) 15. Y.J. Liu, F. Ding, Convergence properties of the least squares estimation algorithm for multivariable systems. Appl. Math. Model. 37(1–2), 476–483 (2013) 16. Y.J. Liu, F. Ding, Y. Shi, An efficient hierarchical identification method for general dual-rate sampleddata systems. Automatica 50(3), 962–973 (2014) 17. Y.J. Liu, F. Ding, Y. Shi, Least squares estimation for a class of non-uniformly sampled systems based on the hierarchical identification principle. Circuits Syst. Signal Process. 31(6), 1985–2000 (2012) 18. Y.J. Liu, Y.S. Xiao, X.L. Zhao, Multi-innovation stochastic gradient algorithm for multiple-input single-output systems using the auxiliary model. Appl. Math. Computat. 215(4), 1477–1483 (2009)
2986
Circuits Syst Signal Process (2014) 33:2971–2986
19. L. Ljung, System Identification: Theory for the User, 2nd edn. (Prentice-Hall, Englewood Cliffs, NJ, 1999) 20. M. Lovera, T. Gustafsson, M. Verhaegen, Recursive subspace identification of linear and non-linear Wiener state–space models. Automatica 36(11), 1639–1650 (2000) 21. J.X. Ma, F. Ding, Recursive relations of the cost functions for the least squares algorithms for multivariable systems. Circuits Syst. Signal Process. 32(1), 83–101 (2013) 22. G. Mercère, L. Bako, Parameterization and identification of multivariable state-space systems: a canonical approach. Automatica 47(8), 1547–1555 (2011) 23. Y. Shi, H.Z. Fang, Kalman filter-based identification for systems with randomly missing measurements in a network environment. Int. J. Control 83(3), 538–551 (2010) 24. B. Sun, D.Q. Zhu, S.X. Yang, A bio-inspired filtered backstepping cascaded tracking control of 7000 m manned submarine vehicle. IEEE Trans. Ind. Electron. 61(7), 3682–3692 (2014) 25. D.N. Vizireanu, A fast, simple and accurate time-varying frequency estimation method for single-phase electric power systems. Measurement 45(5), 1331–1333 (2012) 26. D.N. Vizireanu, S.V. Halunga, Single sine wave parameters estimation method based on four equally spaced samples. Int. J. Electron. 98(7), 941–948 (2011) 27. J. Vörös, Modeling and identification of systems with backlash. Automatica 46(2), 369–374 (2010) 28. D.Q. Wang, F. Ding, Y.Y. Chu, Data filtering based recursive least squares algorithm for Hammerstein systems using the key-term separation principle. Inf. Sci. 222, 203–212 (2013) 29. D.Q. Wang, R. Ding, X.Z. Dong, Iterative parameter estimation for a class of multivariable systems based on the hierarchical identification principle and the gradient search. Circuits Syst. Signal Process. 31(6), 2167–2177 (2012) 30. D.Q. Wang, F. Ding, X.M. Liu, Least squares algorithm for an input nonlinear system with a dynamic subspace state space model. Nonlinear Dyn. 75(1–2), 49–61 (2014) 31. T. Wigren, Recursive prediction error identification and scaling of non-linear state space models using a restricted black box parameterization. Automatica 42(1), 159–168 (2006) 32. A. Wills, B. Ninness, On gradient-based search for multivariable system estimates. IEEE Trans. Automat. Control 53(1), 298–306 (2008) 33. A. Wills, B. Ninness, S. Gibson, Maximum likelihood estimation of state space models from frequency domain data. IEEE Trans. Automat. Control 54(1), 19–33 (2009) 34. A. Wills, T.B. Schä, L. Ljung, B. Ninness, Identification of Hammerstein-Wiener models. Automatica 45(1), 70–81 (2013) 35. W.X. Zheng, Least-squares identification of a class of multivariable systems with correlated disturbances. J. Frankl. Inst. 336(8), 1309–1324 (1999) 36. J. Zheng, K.W.K. Lui, W.-K. Ma, H.C. So, Two simplified recursive Gauss-Newton algorithms for direct amplitude and phase tracking of a real sinusoid. IEEE Signal Proc. Let. 14(12), 984–987 (2007) 37. D.Q. Zhu, H. Huang, S.X. Yang, Dynamic task assignment and path planning of multi-AUV system based on an improved self-organizing map and velocity synthesis method in 3D underwater workspace. IEEE Trans. Cybern. 43(2), 504–514 (2013) 38. A. Žilinskas, J. Žilinskas, A global optimization method based on the reduced simplicial statistical model. Math. Model. Anal. 16(3), 451–460 (2011)
