Retrospective Cost Model Refinement for On-Line Estimation of ...
AIAA Guidance, Navigation, and Control (GNC) Conference August 19-22, 2013, Boston, MA
AIAA 2013-5193
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
Retrospective Cost Model Refinement for On-line Estimation of Constant and Time-Varying Flight Parameters Jiapeng Zhong ∗, Jin Yan †, Yousaf Rahman ‡, Changhong Wang §, and Dennis S. Bernstein¶ University of Michigan, 1320 Beal Ave., Ann Arbor, MI 48109
In this paper we apply retrospective cost model refinement (RCMR) to parameter estimation. To gain insight into the accuracy and speed of convergence of the RCMR estimates, we consider aircraft dynamics with uncertain entries in the dynamics and input matrices. We consider scenarios that include multiple uncertain parameters, alternative measurements, and sensor noise. All of the examples are discrete time, as required by RCMR. To account for the structure of the uncertain parameters, we develop a Jacobianbased technique for constructing the appropriate feedback structure. We use this approach to account for the appearance of stability derivatives in the dynamics matrix, and we use this technique to estimate airspeed variations in the vicinity of trimmed flight.
I.
Introduction
System identification is concerned with using input-output data to construct empirical models. In many cases, some components of the system are well-modeled, and the goal is to use input-output data to improve estimates of poorly modeled components. These components may be connected in cascade, parallel, or feedback with the well-modeled components, and they may be static or dynamic. This problem is known as model updating, model correction, model calibration, and model refinement [1–6]. The most common model-refinement problem is parameter estimation, where data are used to improve estimates of parameters in a model whose structure is known. Parameter estimation is related to, but distinct from, state estimation, where the states evolve due to external inputs and their interaction with other states. In contrast, an unknown parameter may either be constant or time-varying in a pre-specified manner. The close relationship between parameter estimation and state estimation is evident from the widespread use of state-estimation techniques for parameter estimation. In particular, the extended Kalman filter can be used with a linearized model to propagate state and parameter estimates [7]. Alternatively, techniques developed for nonlinear state estimation can be applied to parameter estimation [8–11]. In the present paper we revisit the retrospective cost approach to model refinement (RCMR) [12–15]. RCMR can be used to estimate the dynamics of a possibly dynamic subsystem in feedback interconnection ∗ Graduate
Student, Control Science and Engineering Department, Harbin Institute of Technology, Harbin, China Student, Aerospace Engineering Department, University of Michigan, Ann Arbor ‡ Graduate Student, Aerospace Engineering Department, University of Michigan, Ann Arbor § Professor, Control Science and Engineering Department, Harbin Institute of Technology, Harbin, China ¶ Professor, Aerospace Engineering Department, University of Michigan, Ann Arbor † Graduate
1 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
with a main subsystem; the unknown subsystem is assumed to be inaccessible in the sense that its inputs and outputs are not measured. A special case of an inaccessible subsystem occurs when the unknown subsystem is static; in this case, inaccessible subsystem identification is equivalent to parameter estimation.
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
The goal of the present paper is to apply RCMR to parameter estimation. To gain insight into the accuracy and speed of convergence of the RCMR estimates, we consider aircraft dynamics with uncertain entries in the dynamics and input matrices. These scenarios include multiple uncertain parameters, alternative measurements, and sensor noise. All of the examples are discrete time, as required by RCMR. To account for the structure of the uncertain parameters, we develop a Jacobian-based technique for constructing the appropriate feedback structure. We use this approach to account for the appearance of stability derivatives in the dynamics matrix, and we use this technique to estimate variations of the airspeed in the vicinity of trimmed flight.
II.
Problem Statement
Consider the MIMO discrete-time main system x(k + 1) = Ax(k) + Bu(k) + D1 w(k), y(k) = Cx(k), y0 (k) = E1 x(k) + v(k),
(1) (2) (3)
where x(k) ∈ Rn , y(k) ∈ Rly , y0 (k) ∈ Rly0 , u(k) ∈ Rlu , w(k) ∈ Rlw , and k ≥ 0. The main system (1)–(3) is interconnected with the unknown subsystem modeled by u(k) = Gs (q)y(k),
(4)
where q is the forward shift operator. The system (1)–(4) represents the true system. We assume that the excitation signal w(k) is known. v(k) denotes measurement noise. Next, we assume a model of the main system of the form ˆx(k) + B ˆu ˆ 1 w(k), x ˆ(k + 1) = Aˆ ˆ(k) + D yˆ(k) = Cˆ x ˆ(k), ˆ1 x yˆ0 (k) = E ˆ(k),
(5) (6) (7)
where x ˆ(k) ∈ Rnˆ , yˆ(k) ∈ Rlyˆ , yˆ0 (k) ∈ Rly0 , u ˆ(k) ∈ Rluˆ . The model of the main system is interconnected with the subsystem model ˆ s (q)ˆ u ˆ(k) = G y (k).
(8)
ˆ s (q) by minimizing a cost function based on the performance The goal is to estimate a subsystem model G variable 4
z(k) = yˆ0 (k) − y0 (k) ∈ Rlz
(9)
ˆ s (q) by retrospectively reconstructing the signal u We estimate G ˆ(k) that minimizes the performance at the current time step. The reconstruction of u ˆ(k) uses minimal modeling information about the true system ˆ s (q). (1)–(3), namely, a limited number of Markov parameters. We then use u ˆ(k) and yˆ(k) to construct G Figure 1 illustrates the model-refinement architecture.
III.
Retrospective Control Model Refinement
ˆ We begin by defining Markov parameters of the main system model G(q). For i ≥ 1, let 4 ˆ1 Aˆi−1 B. ˆ Hi = E
2 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
(10)
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
Figure 1. Model-refinement architecture.
ˆ1 B ˆ and H2 = E ˆ1 AˆB. ˆ Let r be a positive integer. Then, for all k ≥ r, Therefore, H1 = E x ˆ(k) = Aˆr x ˆ(k − r) +
r X
ˆu Aˆi−1 B ˆ(k − i) +
i=1
r X
ˆ 1 w(k − i), Aˆi−1 D
(11)
i=1
and thus ˆ1 Aˆr x z(k) = E ˆ(k − r) +
r X
¯U ˆ1 Aˆi−1 D ˆ 1 w(k − i) − y0 (k) + H ¯ (k − 1), E
(12)
i=1 4 ¯ = where H
h
H1
···
Hr
i
4 ¯ (k − 1) = ∈ Rlz ×rluˆ , and U
h
u ˆT (k − 1) · · ·
u ˆT (k − r)
iT
.
¯ and the components of U ¯ (k − 1) and partition the resulting Next, we rearrange the columns of H matrix and vector so that ¯U ¯ (k − 1) = H0 U 0 (k − 1) + HU (k − 1), H
(13)
where H0 ∈ Rlz ×(rluˆ −lU ) , H ∈ Rlz ×lU , U 0 (k − 1) ∈ Rrluˆ −lU , and U (k − 1) ∈ RlU . Then, we can rewrite (12) as z(k) = S(k) + HU (k − 1),
(14)
where 4
ˆ1 Aˆr x S(k) = E ˆ(k − r) +
r X
ˆ1 Aˆi−1 D ˆ 1 w(k − i) − y0 (k) + H0 U 0 (k − 1). E
i=1
¯ = For example, H
h
H1
H2
H3 H = 0
i
h
,
H1
H2
i
" ,
0
U (k − 1) =
u ˆ(k − 1) u ˆ(k − 2)
3 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
# ,
(15)
and H = H3 , U (k − 1) = u ˆ(k − 3). Next, we rewrite (14) with a delay of kj time steps, where 0 ≤ k1 ≤ k2 ≤ · · · ≤ ks , in the form z(k − kj ) = Sj (k − kj ) + Hj Uj (k − kj − 1),
(16)
where (15) becomes 4 ˆ1 Aˆr x Sj (k − kj ) = E ˆ(k − kj − r) +
r X
ˆ1 Aˆi−1 D ˆ 1 w(k − kj − i) − y0 (k − kj ) + H0 U 0 (k − kj − 1) E j j
(17)
i=1
and (13) becomes
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
¯U ¯ (k − kj − 1) = Hj0 Uj0 (k − kj − 1) + Hj Uj (k − kj − 1), H
(18)
where Hj0 ∈ Rlz ×(rluˆ −lUj ) , Hj ∈ Rlz ×lUj , Uj0 (k − kj − 1) ∈ Rrluˆ −lUj , and Uj (k − kj − 1) ∈ RlUj . Now, by stacking z(k − k1 ), . . . , z(k − ks ), we define the extended performance h iT 4 Z(k) = z T (k − k1 ) · · · z T (k − ks ) ∈ Rslz . (19) Therefore, 4 ˜ ˜U ˜ (k − 1), Z(k) = S(k) +H (20) h iT 4 ˜ ˜ ∈ Rslz ×lU˜ , and U ˜ (k − 1) ∈ RlU˜ . The vector U ˜ (k − 1) where S(k) = ST (k − k1 ) · · · ST (k − ks ) ∈ Rslz , H is formed by stacking U1 (k − k1 − 1), . . . , Us (k − ks − 1) repetitions of components. For " and removing # u ˆ(k − 1) example, with k1 = 0 and k2 = 1, stacking U1 (k − 1) = and U2 (k − 2) = u ˆ(k − 2) results in u ˆ(k − 2) " # u ˆ(k − 1) ˜ consists of the entries of H1 , . . . , Hs arranged according ˜ U (k − 1) = . The coefficient matrix H u ˆ(k − 2) ˜ (k − 1). Furthermore, we assume that the last entry of U ˜ (k − 1) is a component of to the structure of U u ˆ(k − r).
Next, we define the retrospective performance 4
zˆ(k − kj ) = Sj (k − kj ) + Hj Uj∗ (k − kj − 1),
(21)
where the actual past subsystem outputs Uj (k − kj − 1) in (16) are replaced by the surrogate subsystem outputs Uj∗ (k − kj − 1). The extended retrospective performance for (21), which is defined as 4 ˆ Z(k) =
h
zˆT (k − k1 ) · · ·
zˆT (k − ks )
iT
∈ Rslz ,
(22)
is given by ˜ ˜U ˆ ˜ ∗ (k − 1), Z(k) = S(k) +H
(23)
˜ ∗ (k − 1) ∈ RlU˜ are components of U ∗ (k − k1 − 1), . . . , Us∗ (k − ks − 1) ordered in where the components of U 1 ˜ (k − 1). Subtracting (20) from (23) yields the same way as the components of U ˜U ˜U ˆ ˜ (k − 1) + H ˜ ∗ (k − 1). Z(k) = Z(k) − H
(24)
Finally, we define the retrospective cost function 4
˜ ∗ (k − 1), k) = Zˆ T (k)R(k)Z(k), ˆ J(U
(25)
where R(k) ∈ Rslz ×slz is a positive-definite performance weighting. The goal is to determine refined sub˜ ∗ (k − 1) that would have provided better performance than the subsystem outputs U (k) system outputs U ˜ ∗ (k − 1) are subsequently used to that were applied to the system. The refined subsystem outputs values U update the subsystem estimate. 4 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
A.
Cost Function Optimization with Adaptive Regularization To ensure that (25) has a global minimizer, we consider the regularized cost 4 ¯U ˜ ∗ (k − 1), k) = ˆ ˜ ∗T (k − 1)U ˜ ∗ (k − 1), J( Zˆ T (k)R(k)Z(k) + η(k)U
(26)
where η(k) = η¯z T (k)z(k) and η¯ ≥ 0. Substituting (24) into (26) yields ¯U ˜ ∗ (k − 1), k) = U ˜ ∗ (k − 1)T A(k)U ˜ ∗ (k − 1) + B(k)U ˜ ∗ (k − 1) + C(k), J(
(27)
where 4 ˜T ˜ + η(k)Il , A(k) = H R(k)H ˜ U 4
(28)
˜T
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
˜U ˜ (k − 1)], B(k) = 2H R(k)[Z(k) − H
(29)
4
˜U ˜ T R(k)H ˜U ˜ (k − 1) + U ˜ T (k − 1)H ˜ (k − 1). C(k) = Z T (k)R(k)Z(k) − 2Z T (k)R(k)H
(30)
˜ has full column rank or η(k) > 0, then A(k) is positive definite. In this case, J( ¯U ˜ ∗ (k − 1), k) has If either H the unique global minimizer ˜ ∗ (k − 1) = − 1 A−1 (k)B(k). U 2
B.
(31)
Subsystem Modeling The subsystem output u ˆ(k) is given by the exactly proper time-series model of order nc given by u ˆ(k) =
nc X
Mi (k)ˆ u(k − i) +
nc X
Ni (k)ˆ y (k − i),
(32)
i=0
i=1
where, for all i = 0, 1, . . . , nc , Mi (k) ∈ Rluˆ ×luˆ , Ni (k) ∈ Rluˆ ×lyˆ . The subsystem output (32) can be expressed as u ˆ(k) = θ(k)φ(k), where θ(k) ∈ Rluˆ ×[nc (luˆ +lyˆ )+lyˆ ] and φ(k) ∈ Rnc (luˆ +lyˆ )+lyˆ , 4
θ(k) = [M1 (k) · · · Mnc (k) N0 (k) · · · Nnc (k)] , T 4 T φ(k) = u ˆ (k − 1) · · · u ˆT (k − nc ) yˆT (k) · · · yˆT (k − nc ) .
(33)
u ˆ(k) = θ(k)ˆ y (k) = N0 (k)ˆ y (k).
(35)
(34)
If nc = 0, then
C.
Recursive Least Squares Update ˜ ∗ (k − 1) contains u∗ (k − d), and define the retrospective cost function Let d > 0 such that U 4
JR (θ(k)) =
k X
λk−i ||u∗T (k − d) − φT (k − d − 1)θT (k)||2 + λk (θ(k) − θ(0))P −1 (0)(θ(k) − θ(0))T ,
(36)
i=1
where φ(k − d) is given by (34), k · k is the Euclidean norm, and λ(k) ∈ (0, 1] is the forgetting factor. Minimizing the cumulative cost function yields retrospective cost optimization (RCO) θT (k) = θT (k − 1) + P (k − 1)φ(k − d − 1) · [φT (k − d)P (k − 1)φ(k − d − 1) + λ(k)]−1 · (u∗ (k − d) − φT (k − d − 1)θT (k − 1)).
5 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
(37)
The error covariance is updated by P (k) = λ−1 (k)P (k − 1) − λ−1 (k)P (k − 1)φ(k − d − 1) · [φT (k − d − 1)P (k − 1)φ(k − d) + λ(k)]−1 · φT (k − d − 1)P (k − 1).
(38)
We initialize the error covariance matrix as P (0) = βI, where β > 0.
IV.
Mass-Spring Example: Stiffness Estimation
Example IV.1. We consider the mass-spring-damper structure modeled by
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
m¨ q + cq˙ + κq = w,
(39)
where m, c, and κ are the mass, damping, and stiffness, respectively, and w is the force input. We obtain the state-space representation #" # " # " # " 0 q 0 1 q˙ + w, (40) = 1 κ c q˙ −m −m q¨ m " # h i q + v, (41) y= 0 1 q˙ where q and q˙ are the position and velocity, respectively, of the mass. The parameters are chosen as m = 1, c = 5, and κ = 10. The goal is to estimate κ with the initial estimate of κ ˆ (0) = 0. Using Euler discretizion with sample time Ts yields " # " #" # " # x1 (k + 1) 1 Ts x1 (k) 0 = + Ts w(k), s s x2 (k + 1) − κT 1 − cT x2 (k) m m m " # h i x (k) 1 + v(k), y(k) = 0 1 x2 (k) 4
4
(42)
(43)
4
"
4
0 − Tms
#
ˆ ˆ ˆ κC, ˆ where B ˆ = where, x1 (k) = q(kTs ) and x2 (k) = q(kT ˙ s ). We define Aclosed−loop = A + Bˆ h i 4 and Cˆ = 1 0 . The estimate of κ is given by κ ˆ (k) = θ(k). The initial conditions are assumed to be x1 (0) = 0.1 and x2 (0) = 0.01, which are unknown. Let the sensor noise v be Gaussian h white noise i with 2 ˜ mean µv and variance σ . Furthermore, P (0) = 0.6, λ = 0.999, η = 0, µv = 0, and H = H H , which 1
v
2
ˆ The input is w(k) = 0.8 sin(πk/20). RCMR is turned on are the first and second Markov parameters of G. at k = 20 steps. Figure 2 shows the performance z and the estimated κ ˆ. We next consider the effect of sensor noise. We consider the same estimation problem with white 4 σ2
sensor noise. We define SNR = σy¯20 , where σy2¯0 is the variance of the output signal y¯0 . Figure 3 shows the v effect of sensor noise on RCMR. As sensor noise is increased, the accuracy of the parameter identification is reduced.
V.
Estimation of a Repeated Parameter
We next consider the case where one unknown parameter appears in multiple locations within the model. Consider the continuous-time system x˙ = Ac x + Dc w,
6 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
(44)
0.2
20 RCMR estimate Initial estimate True parameter
0.15 Parameter Estimate
Performance z(k)
15 0.1 0.05 0 −0.05
10
5
0 −0.1 −0.15
0
100
200 Time Step (k)
300
−5
400
0
100
200 Time Step (k)
300
400
20
20
20
RCMR estimate Initial estimate True parameter
RCMR estimate Initial estimate True parameter
5
0
15
Parameter Estimate
10
−5
RCMR estimate Initial estimate True parameter
15
Parameter Estimate
15
Parameter Estimate
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
ˆ Figure 2. Example IV.1, performance z(k) and estimated k.
10
5
10
5
0
0
50
100
150
200 250 Time Step (k)
300
350
400
−5
0
0
50
100
(a) SNR=100
150
200 250 Time Step (k)
300
350
400
−5
(b) SNR=10
0
50
100
150
200 250 Time Step (k)
300
350
(c) SNR=1
Figure 3. Example VI.1, the estimate m ˆ of m for various values of SNR. As SNR decreases, the accuracy of the estimate is degraded.
where a parameter α in Ac and Dc is uncertain. Let α = α ˆ + ∆α, where α ˆ is the initial estimate of α. The system model can be written approximately as x˙ = (Ac +
∂Ac ∂Dc ∆α)x + (Dc + ∆α)w. ∂α α=αˆ ∂α α=αˆ
(45)
Discretizing (45) with sampling time Ts yields ∂Dc ∂Ac ∆αTs )x(k) + (Dc Ts + ∆αTs )w(k) x(k + 1) = (I + Ac Ts + ∂α α=αˆ ∂α α=αˆ ∂Ac ∂Dc = (A + ∆αTs )x(k) + (D1 + ∆αTs )w(k). ∂α α=αˆ ∂α α=αˆ
(46)
We use (46) as the true system model with uncertain parameter ∆α. Since our goal is to estimate ∆α, we rewrite (46) w(k) i h ∂Dc c x(k + 1) = Ax + D1 ∂A (47) ∆αx(k) . ˆ Ts ˆ Ts ∂α |α=α ∂α |α=α ∆αw(k)
7 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
400
Next, let ˆ= B
h "
yˆ(k) =
∂Ac ˆ Ts ∂α |α=α
x ˆ(k) w(k)
∂Dc ˆ Ts ∂α |α=α
i
,
(48)
# ,
(49)
ˆ s (q) = θ(k), ∆α(k) = G
(50)
u ˆ(k) = ∆α(k)ˆ y (k).
(51)
Note that (51) shows that u ˆ(k) is a function of yˆ(k). We define
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
4 Yˆ (k) =
h
yˆT (k − k1 ) · · ·
yˆT (k − ks )
iT
∈ Rsly .
(52)
With these signals, (27)-(31) become ∗ ¯ J(∆α (k − 1), k) = A(k)∆α∗2 (k − 1) + B(k)∆α∗ (k − 1) + C(k),
(53)
4
˜ T R(k)H ˜ + η(k)Il ]Yˆ (k − 1), A(k) = Yˆ (k − 1)T [H ˜ U
(54)
4
˜ T R(k)[Z(k) − H ˜U ˜ (k − 1)]Yˆ (k − 1), B(k) = 2H 4
(55)
˜U ˜U ˜ (k − 1) + U (k − 1)H R(k)H ˜ (k − 1), C(k) = Z (k)R(k)Z(k) − 2Z (k)R(k)H ˜ ∗ (k − 1) = ∆α∗ (k − 1)Yˆ (k − 1) = − 1 A−1 (k)B(k)Yˆ (k − 1). U 2 T
˜T
˜T
T
(56) (57)
˜ ∗ (k − 1) and the corresponding parameter in Y˜ (k − 1) to execute the recursive We use one parameter in U least squares update. In the next section we use the same mass-spring system to illustrate the algorithm.
VI.
Mass-Spring Example: Inertia Estimation
Example the mass-spring-damper structure shown in Example IV.1 where Ac = # VI.1." Consider # h i 0 1 0 , D = , and E = . In this example we assume that the parameter m is 0 1 c c κ 1 c −m −m m unknown. We demonstrate the algorithm by choosing κ = 30 and c = 5, and we assume that m = 0.9. We use an initial estimate is m(0) ˆ = 1, so that ∆m(0) = −0.1. "
From (48), we obtain " ˆ= B
0 3
0 0.5
0 −0.1
# ,
(58)
i , which are the first and " # x ˆ (0) 1 ˆ We choose the ramp input w(k) = 0.1k and the initial state second Markov parameters of G. = x ˆ2 (0) " # " # x1 (0) 0 = . The system refinement algorithm is turned on at k = 100 steps. Figure 4 shows the x2 (0) 0 performance z and the estimate m. ˆ ˜ = We choose P (0) = 1, λ = 1, η = 0, µ = 0, σv2 = 0 and H
h
H1
H2
We next consider the effect of sensor noise. We assume µv = 0. Figure 5 shows the estimation performance for several values of σv2 .
8 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
−4
4
x 10
1.1
2 Parameter Estimate
0 Performance z(k)
RCMR estimate Initial estimate True parameter
1.05
−2 −4 −6
1 0.95 0.9
−8 0.85
−10 −12
0
50
100 150 200 Time Step (k)
250
0.8
300
0
50
100 150 200 Time Step (k)
250
300
1.1
1.1 RCMR estimate Initial estimate True parameter
0.95
0.9
0.85
0.8
1.05
Parameter Estimate
1
1.1 RCMR estimate Initial estimate True parameter
1
0.95
0.9
0.85
0
50
100
150 Time Step (k)
200
250
300
0.8
RCMR estimate Initial estimate True parameter
1.05
Parameter Estimate
1.05
Parameter Estimate
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
Figure 4. Example VI.1, the performance z and the estimate m. ˆ
1
0.95
0.9
0.85
0
(a) σv2 = 1 × 10−10
50
100
150 Time Step (k)
200
250
300
(b) σv2 = 1 × 10−9
0.8
0
50
100
150 Time Step (k)
200
250
(c) σv2 = 1 × 10−8
Figure 5. Example VI.1, the estimate m ˆ of m for various values of σv2 . As σv2 increases, the accuracy of the estimate is degraded.
VII.
Aircraft Examples
We now estimate stability parameters under various scenarios of measurements noise. Consider the linearized longitudinal B-767 model u˙ −0.0168 0.1121 0.0003 −0.5608 u −0.0243 0.0519 0.0015 α −0.0634 −0.0005 α˙ −0.0164 −0.7771 0.9945 = + q˙ −0.0417 −3.6595 −0.9544 q −3.6942 0.0243 0 θ˙ 0 0 1.0000 0 θ 0 0
and measurement
"
δe δT
# ,
(59)
where δe is the elevator deflection and δT is the thrust perturbation. We discretize (59) using a zero-order hold with Ts = 0.1s, which yields 0.9983 0.0110 −0.0022 −0.0560 −0.0022 0.0052 0.0907 0.0002 −0.0018 0.9085 −0.0234 0.0001 A= (60) , D1 = , −0.0037 −0.3336 0.8924 −0.3492 0.0023 0.0001 −0.0002 −0.0172 0.0948 1.0000 −0.0178 0.0001 h iT and w = δe δT is the system input. For all examples in this section, RCMR is turned on at k = 100 h iT steps, and x ˆ(0) = x(0) = 0 0 0 0 . We choose δe and δT to be zero-mean white Gaussian noises with variance 0.001. 9 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
300
Example VII.1. Assume an initial guess of A(3, 3) = 0.8000. Theh true value of iA(3, 3) is 0.8924 and ∆A(3, 3) = 0.0924. Our goal is to identify A(3, 3). We choose E1 = 0 0 1 0 so that y¯0 = q h iT h i ˆ 1 = D1 , E ˆ1 = E1 , B ˆ= 0 0 1 0 and yˆ0 = qˆ. In this example, Aˆ = A, D and Cˆ = 0 0 1 0 . 2 ˆ In this h model, nci = 0 and ∆A(3, 3) = θ(k). We choose P (0) = 200, λ = 0.97, η = 0, µv = 0, σv = 0 and ˜ = H H , which are the first and second Markov parameters of G. ˆ The parameters of the model H 1 2 refinement algorithm are chosen such that z(k) is minimized. Figure 6 shows the performance z and the ˆ 3). estimated A(3, 0.02
1 RCMR estimate Initial estimate True parameter
0.015 Parameter Estimate
Performance z(k)
0.005 0 −0.005
0.9
0.85
0.8 −0.01 −0.015
0
100
200 Time Step (k)
300
0.75
400
0
100
200 Time Step (k)
300
400
ˆ 3). Figure 6. Example VII.1, performance z(k) and estimated A(3,
Example VII.2. We next consider the case where multiple parameters in A are possibly uncertain h iT and K = ∆A has the form shown. We assume full state feedback and y = . Thus u α q θ ˆK ˆ C, ˆ where B ˆ = Cˆ = I4 and K(k) ˆ Aˆclosed−loop = Aˆ + B = θ(k), where K=
−0.05 0 0 0
−0.002 −0.04 −0.04 0
0 0 −0.04 0
0 0 0 −0.04
.
(61)
i H1 H2 . Figure 7 shows the performance z(k) and the parameter estimates θ(k). The parameter error kθ(k) − Kk2 converges to zero. ˜ = We choose P (0) = 100I4 , λ = 0.98, η = 1, µv = 0, σv2 = 0 and H
0.14
h
0.08 0.06
0.12
0.04 0.1 0.02 0.08
θ(k)
Parameter Error
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
0.95 0.01
0.06
0 −0.02 −0.04
0.04 −0.06 0.02 0
−0.08 0
500
1000 1500 Time Step (k)
2000
−0.1
0
500
1000 1500 Time Step (k)
2000
Figure 7. Example VII.2, parameter error kθ(k) − Kk2 and estimated θ(k).
Example VII.3. We next consider the case where parameters of the A matrix that are not being ˆ 2) = −0.3000, whereas the true value is −0.3336. We estimated have modeling errors. We assume A(3, 10 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
assume that this modeling error is unknown and that h i our goal is to identify A(3, 3). We choose P (0) = 100, 2 ˜ = H H , which are the first and second Markov parameters of G. ˆ λ = 1, η = 0, µv = 0, σv = 0 and H 1 2 Figure 8 shows that the modelling errors degrade the accuracy of the estimates of the targeted parameters. 0.015
1
0.005 0 −0.005
0.9
0.85
0.8
0
200
400 600 Time Step (k)
800
0.75
1000
0
200
400 600 Time Step (k)
800
1000
ˆ 3). Figure 8. Example VII.3, performance z and estimated A(3,
Example VII.4. We next consider the case where the unknown parameter is time-varying. In this h iT ˆ 1 = D1 , E ˆ1 = E1 , B ˆ = 0 0 1 0 example, we assume A(3, 3) is time-varying, Aˆ = A, D , and h i Cˆ = 0 0 1 0 . Figure 9 shows the performance when A(3, 3) is varying at a constant rate, and Figure 10 shows the performance when A(3, as a sinusoidal signal. In both cases, P (0) = 10, h 3) is varying i 2 ˜ λ = 0.98, η = 0, µv = 0, σv = 0 and H = H1 H2 . −4
2
x 10
0.9
1.5
RCMR estimate Initial estimate True parameter
0.895 Parameter Estimate
1 Performance z(k)
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
−0.01 −0.015
RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
Performance z(k)
0.01
0.5 0 −0.5
0.89 0.885 0.88
−1 0.875
−1.5 −2
0
500
1000 1500 Time Step (k)
2000
0.87
0
500
1000 1500 Time Step (k)
2000
ˆ 3). A(3, 3) is the ramp signal A(3, 3) = 0.8924 − 10−5 k, Figure 9. Example VII.4, performance z and estimated A(3, 2 ˜ P (0) = 10, λ = 0.98, η = 0, µv = 0, σv = 0 and H = [H1 , H2 ].
Example VII.5. We next consider Example VII.1 with sensor noise present. We assume µv = 0. Figure 11 shows the estimation performance for several values of SNR. As the SNR decreases, the accuracy of the parameter estimatition degrades and convergence time increases. Example VII.6. We next consider Example VII.1 with non-zero-mean measurement noise. We assume SNR=100 and compare the estimate for several values of µv . Figure 12 shows the affect on estimation performance as µv increases. As the bias µv increases, the accuracy of the parameter estimation degrades and convergence time increases.
11 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
0.015
0.005 0 −0.005
0.95
0.9
0.85
−0.01 −0.015
RCMR estimate Initial estimate True parameter
1 Parameter Estimate
Performance z(k)
0.01
0
500
1000 1500 Time Step (k)
0.8
2000
0
500
1000 1500 Time Step (k)
2000
1 RCMR estimate Initial estimate True parameter
0.9
0.85
0.8
0.75
0.9
0.85
0.8
0
500
1000 Time Step (k)
1500
0.75
2000
1 RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
0.95 Parameter Estimate
Parameter Estimate
0.95
1 RCMR estimate Initial estimate True parameter
0.9
0.85
0.8
0
(a) SNR=∞
500
1000 Time Step (k)
1500
0.75
2000
RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
1
0.9
0.85
0.8
0
(b) SNR=100
500
1000 Time Step (k)
1500
0.75
2000
0
(c) SNR=10
2000
4000 6000 Time Step (k)
8000
10000
(d) SNR=1
Figure 11. Example.VII.5. (a) P (0) = 200, λ = 0.97. (b) P (0) = 200, λ = 1. (c) P (0) = 100, λ = 1. (d) P (0) = 15, ˜ = [H1 , H2 ]. λ = 1. In all cases, η = 0, and H
1 RCMR estimate Initial estimate True parameter
0.85
0.8
0.9
0.85
0.8
0
500
1000 Time Step (k)
(a) µv = 0
1500
2000
0.75
1 RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
0.9
0.75
0.95 Parameter Estimate
0.95
1 RCMR estimate Initial estimate True parameter
0.9
0.85
0.8
0
500
1000 Time Step (k)
1500
2000
0.75
RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
1
Parameter Estimate
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
ˆ 3). Ad (3, 3) is the sinusoidal signal A(3, 3) = Figure 10. Example VII.4, performance z and estimated A(3, ˜ = [H1 , H2 ]. 0.8924 + 0.08 sin(πk/400), P (0) = 10, λ = 0.98, η = 0, µv = 0, σv2 = 0 and H
0.9
0.85
0.8
0
(b) µv = −0.001
500
1000 Time Step (k)
1500
(c) µv = −0.01
2000
0.75
0
2000
4000 6000 Time Step (k)
8000
10000
(d) µv = −0.1
Figure 12. Example.VII.6. (a) P (0) = 200, λ = 0.97. (b) P (0) = 200, λ = 1. (c) P (0) = 200, λ = 1. (d) P (0) = 20, ˜ = [H1 , H2 ]. λ = 1. In all cases, η = 0, and H
VIII.
Airspeed Estimation
In this section we use RCMR to estimate airspeed. Example VIII.1. Consider the linearized longitudinal transfer functions [16] for a typical business
12 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
by
Xu0 + XTu0
Zu0 +ZTu0 U0 = Mu0 + MTu 0 0
Xα0
Xq0
Zα0 U0
U0 +Zq0 U0
Mα0 + MTα0 0
Mq0 1
−g
0 0 0
u α q θ
Xδe0
Zδe0 U0 + Mδe0 0
δe,
(62)
where U0 is the aircraft speed at which the equations of motion of the aircraft are linearized. X, Z, M are the related coefficients of the aircraft dynamic mode and many of these coefficients are depend on U0 . Using the relations of these coefficients to U0 , (62) can be written as − 4.9950 8.9782 0 −32.1522 0 u u˙ U0 1 − 1255.3650 − 93.8250 − 445.7224 0 α − 42.1968 α˙ U0 U02 U02 U0 + (63) = δe 634.2975 0.6075 q˙ −7.4416 − U0 0 q −17.6737 U0 θ θ˙ 0 0 0 1 0 ˆ0 = 675 The goal is to estimate U0 . The true value of U0 (0) is 625 ft/s, and we use the initial estimate U ft/s. We choose various y¯0 , and compare the performance of the algorithm. Figure 13 shows the performance of RCMR with y¯0 = α.h In this case, the i input w(k) = 0.0008, 2 ˜ P (0) = 0.0003, λ = 0.9999, η = 0, µ = 0, σv = 0 and H = H1 · · · H12 . Figure 14 shows the performance of RCMR with y¯i0 = q. In this case, the input w(k) = 10−7 k, P (0) = 1, λ = 0.99, η = 0, µ = 0, h ˜ = H H . Figure 15 shows the performance of RCMR with y¯0 = θ. In this case, the σv2 = 0 and H 1 2 h i −7 ˜ = H H . input w(k) = 10 k, P (0) = 0.01, λ = 0.999, η = 0, µ = 0, σv2 = 0, and H 1 2 −5
2
x 10
680 RCMR estimate Initial estimate True parameter
670 Parameter Estimate
0 Performance z(k)
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
jet in cruise given u˙ α˙ q˙ θ˙
−2
−4
660 650 640 630
−6 620 −8
0
1000
2000 3000 4000 Time Step (k)
5000
6000
610
0
1000
2000 3000 4000 Time Step (k)
5000
6000
ˆ0 are shown. Figure 13. Airspeed estimation with y¯0 = α. The performance z and the airspeed estimate U
IX.
Conclusions
This paper showed that RCMR can be used to estimate unknown parameters in a state space model. Both constant and time-varying parameters were considered, under various sensor noise levels and choices of measurements. Future work will compare the accuracy of this technique to nonlinear estimation methods.
Acknowledgment This work was supported in part by the National Aeronautics and Space Administration under Cooperative Agreement NNX12AM54A. 13 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
−7
4
x 10
680
Parameter Estimate
Performance z(k)
0 −2 −4 −6 −8
RCMR estimate Initial estimate True parameter
670
2
660 650 640 630 620
0
1000
2000 3000 4000 Time Step (k)
5000
610
6000
0
1000
2000 3000 4000 Time Step (k)
5000
6000
−5
1.5
x 10
700
Parameter Estimate
0.5 0 −0.5 −1 −1.5
RCMR estimate Initial estimate True parameter
680
1 Performance z(k)
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
ˆ0 are shown. Figure 14. Airspeed estimation with y¯0 = q. The performance z and the airspeed estimate U
660 640 620 600 580
0
1000
2000 3000 4000 Time Step (k)
5000
6000
560
0
1000
2000 3000 4000 Time Step (k)
5000
6000
ˆ0 are shown. Figure 15. Airspeed estimation with y¯0 = θ. The performance z and the airspeed estimate U
References 1 Palanthandalam-Madapusi, H., Renk, E. L., and Bernstein, D. S., “Data-Based Model Refinement for Linear and Hammerstein Systems Using Subspace Identification and Adaptive Disturbance Rejection,” Proc. Conf. Contr. Appl., Toronto, Canada, August 2005, pp. 1630–1635. 2 D’Amato, A. M., Arritt, B. J., Banik, J. A., Ardelean, E. V., and Bernstein, D. S., “Structural Health Determination and Model Refinement for a Deployable Composite Boom,” Proc. AIAA SDM Conf., Palm Springs, CA, April 2009, AIAA2009-2373. 3 Aljanaideh, K., Ali, A. A., Holzel, M., Kukreja, S., and Bernstein, D. S., “Sensor-to-Sensor Identification of Hammerstein Systems,” Proc. Conf. Dec. Contr., Maui, HI, December 2012. 4 Haroon, Y. W. L. M., Adams, D. E., and Ferri, A. A., “A Time and Frequency Domain Approach For Identifying Nonlinear Mechanical System Models in the Absence of an Input Measurement,” Sound and Vibration, Vol. 283, 2005, pp. 1137– 1155. 5 Moheimani, S. O. R., “Model Correction for Sampled-data Models of Structures,” Journal of Guidance, Control, and Dynamics, Vol. 24, 2001, pp. 634–637. 6 Danforth, C. M., Kalnay, E., and Miyoshi, T., “Estimating and Correcting Global Weather Model Error,” Monthly Weather Rev., Vol. 135, 2007, pp. 281–299. 7 Ljung, L., “Asymptotic Behavior of the Extended Kalman Filter as a Parameter Estimator for Linear Systems,” IEEE Trans. Autom. Control, Vol. 24, 1979, pp. 36–50. 8 Ristic, B., Arulampalam, S., and Gordon, N., Beyond the Kalman Filter , Artech House, 2004. 9 Julier, S., Uhlmann, J., and Durrant-Whyte, H., “A New Method for the Nonlinear Transformation of Means and Covariances in Filters and Estimators,” IEEE Trans. Autom. Control, Vol. 45, 2000, pp. 477–482. 10 Psiaki, M., “Backward-Smoothing Extended Kalman Filter,” AIAA J. Guid. Contr. Dyn., Vol. 28, 2005, pp. 885–894. 11 Rao, C., Rawlings, J., and Mayne, D., “Constrained State Estimation for Nonlinear Discrete-Time Systems: Stability and Moving Horizon Approximations,” IEEE Trans. Autom. Control, Vol. 48, 2003, pp. 246–258. 12 Santillo, M. A., D’Amato, A. M., and Bernstein, D. S., “System Identification Using a Retrospective Correction Filter for Adaptive Feedback Model Updating,” Proc. Amer. Contr. Conf., St. Louis, MO, June 2009, pp. 4392–4397.
14 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
13 D’Amato, A. M. and Bernstein, D. S., “Linear Fractional Transformation Identification Using Retrospective Cost Optimization,” Proc. SYSID, Saint-Malo, France, July 2009, pp. 450–455. 14 D’Amato, A. M., Wu, A. R., Mitchell, K. S., Kukreja, S. L., and Bernstein, D. S., “Damage Localization for Structural Health Monitoring Using Retrospective Cost Model Refinement,” Proc. AIAA SDM Conf., Orlando, FL, April 2010, AIAA2010-2628. 15 D’Amato, A. M., Ridley, A. J., and Bernstein, D. S., “Retrospective-Cost-Based Adaptive Model Refinement for the Ionosphere and Thermosphere,” Statistical Analysis and Data Mining, Vol. 4, 2011, pp. 446–458. 16 Roskam, J., Airplane Flight Dynamics and Automatic Flight Control, DARcorp, 2001.
15 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
AIAA 2013-5193
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
Retrospective Cost Model Refinement for On-line Estimation of Constant and Time-Varying Flight Parameters Jiapeng Zhong ∗, Jin Yan †, Yousaf Rahman ‡, Changhong Wang §, and Dennis S. Bernstein¶ University of Michigan, 1320 Beal Ave., Ann Arbor, MI 48109
In this paper we apply retrospective cost model refinement (RCMR) to parameter estimation. To gain insight into the accuracy and speed of convergence of the RCMR estimates, we consider aircraft dynamics with uncertain entries in the dynamics and input matrices. We consider scenarios that include multiple uncertain parameters, alternative measurements, and sensor noise. All of the examples are discrete time, as required by RCMR. To account for the structure of the uncertain parameters, we develop a Jacobianbased technique for constructing the appropriate feedback structure. We use this approach to account for the appearance of stability derivatives in the dynamics matrix, and we use this technique to estimate airspeed variations in the vicinity of trimmed flight.
I.
Introduction
System identification is concerned with using input-output data to construct empirical models. In many cases, some components of the system are well-modeled, and the goal is to use input-output data to improve estimates of poorly modeled components. These components may be connected in cascade, parallel, or feedback with the well-modeled components, and they may be static or dynamic. This problem is known as model updating, model correction, model calibration, and model refinement [1–6]. The most common model-refinement problem is parameter estimation, where data are used to improve estimates of parameters in a model whose structure is known. Parameter estimation is related to, but distinct from, state estimation, where the states evolve due to external inputs and their interaction with other states. In contrast, an unknown parameter may either be constant or time-varying in a pre-specified manner. The close relationship between parameter estimation and state estimation is evident from the widespread use of state-estimation techniques for parameter estimation. In particular, the extended Kalman filter can be used with a linearized model to propagate state and parameter estimates [7]. Alternatively, techniques developed for nonlinear state estimation can be applied to parameter estimation [8–11]. In the present paper we revisit the retrospective cost approach to model refinement (RCMR) [12–15]. RCMR can be used to estimate the dynamics of a possibly dynamic subsystem in feedback interconnection ∗ Graduate
Student, Control Science and Engineering Department, Harbin Institute of Technology, Harbin, China Student, Aerospace Engineering Department, University of Michigan, Ann Arbor ‡ Graduate Student, Aerospace Engineering Department, University of Michigan, Ann Arbor § Professor, Control Science and Engineering Department, Harbin Institute of Technology, Harbin, China ¶ Professor, Aerospace Engineering Department, University of Michigan, Ann Arbor † Graduate
1 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
with a main subsystem; the unknown subsystem is assumed to be inaccessible in the sense that its inputs and outputs are not measured. A special case of an inaccessible subsystem occurs when the unknown subsystem is static; in this case, inaccessible subsystem identification is equivalent to parameter estimation.
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
The goal of the present paper is to apply RCMR to parameter estimation. To gain insight into the accuracy and speed of convergence of the RCMR estimates, we consider aircraft dynamics with uncertain entries in the dynamics and input matrices. These scenarios include multiple uncertain parameters, alternative measurements, and sensor noise. All of the examples are discrete time, as required by RCMR. To account for the structure of the uncertain parameters, we develop a Jacobian-based technique for constructing the appropriate feedback structure. We use this approach to account for the appearance of stability derivatives in the dynamics matrix, and we use this technique to estimate variations of the airspeed in the vicinity of trimmed flight.
II.
Problem Statement
Consider the MIMO discrete-time main system x(k + 1) = Ax(k) + Bu(k) + D1 w(k), y(k) = Cx(k), y0 (k) = E1 x(k) + v(k),
(1) (2) (3)
where x(k) ∈ Rn , y(k) ∈ Rly , y0 (k) ∈ Rly0 , u(k) ∈ Rlu , w(k) ∈ Rlw , and k ≥ 0. The main system (1)–(3) is interconnected with the unknown subsystem modeled by u(k) = Gs (q)y(k),
(4)
where q is the forward shift operator. The system (1)–(4) represents the true system. We assume that the excitation signal w(k) is known. v(k) denotes measurement noise. Next, we assume a model of the main system of the form ˆx(k) + B ˆu ˆ 1 w(k), x ˆ(k + 1) = Aˆ ˆ(k) + D yˆ(k) = Cˆ x ˆ(k), ˆ1 x yˆ0 (k) = E ˆ(k),
(5) (6) (7)
where x ˆ(k) ∈ Rnˆ , yˆ(k) ∈ Rlyˆ , yˆ0 (k) ∈ Rly0 , u ˆ(k) ∈ Rluˆ . The model of the main system is interconnected with the subsystem model ˆ s (q)ˆ u ˆ(k) = G y (k).
(8)
ˆ s (q) by minimizing a cost function based on the performance The goal is to estimate a subsystem model G variable 4
z(k) = yˆ0 (k) − y0 (k) ∈ Rlz
(9)
ˆ s (q) by retrospectively reconstructing the signal u We estimate G ˆ(k) that minimizes the performance at the current time step. The reconstruction of u ˆ(k) uses minimal modeling information about the true system ˆ s (q). (1)–(3), namely, a limited number of Markov parameters. We then use u ˆ(k) and yˆ(k) to construct G Figure 1 illustrates the model-refinement architecture.
III.
Retrospective Control Model Refinement
ˆ We begin by defining Markov parameters of the main system model G(q). For i ≥ 1, let 4 ˆ1 Aˆi−1 B. ˆ Hi = E
2 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
(10)
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
Figure 1. Model-refinement architecture.
ˆ1 B ˆ and H2 = E ˆ1 AˆB. ˆ Let r be a positive integer. Then, for all k ≥ r, Therefore, H1 = E x ˆ(k) = Aˆr x ˆ(k − r) +
r X
ˆu Aˆi−1 B ˆ(k − i) +
i=1
r X
ˆ 1 w(k − i), Aˆi−1 D
(11)
i=1
and thus ˆ1 Aˆr x z(k) = E ˆ(k − r) +
r X
¯U ˆ1 Aˆi−1 D ˆ 1 w(k − i) − y0 (k) + H ¯ (k − 1), E
(12)
i=1 4 ¯ = where H
h
H1
···
Hr
i
4 ¯ (k − 1) = ∈ Rlz ×rluˆ , and U
h
u ˆT (k − 1) · · ·
u ˆT (k − r)
iT
.
¯ and the components of U ¯ (k − 1) and partition the resulting Next, we rearrange the columns of H matrix and vector so that ¯U ¯ (k − 1) = H0 U 0 (k − 1) + HU (k − 1), H
(13)
where H0 ∈ Rlz ×(rluˆ −lU ) , H ∈ Rlz ×lU , U 0 (k − 1) ∈ Rrluˆ −lU , and U (k − 1) ∈ RlU . Then, we can rewrite (12) as z(k) = S(k) + HU (k − 1),
(14)
where 4
ˆ1 Aˆr x S(k) = E ˆ(k − r) +
r X
ˆ1 Aˆi−1 D ˆ 1 w(k − i) − y0 (k) + H0 U 0 (k − 1). E
i=1
¯ = For example, H
h
H1
H2
H3 H = 0
i
h
,
H1
H2
i
" ,
0
U (k − 1) =
u ˆ(k − 1) u ˆ(k − 2)
3 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
# ,
(15)
and H = H3 , U (k − 1) = u ˆ(k − 3). Next, we rewrite (14) with a delay of kj time steps, where 0 ≤ k1 ≤ k2 ≤ · · · ≤ ks , in the form z(k − kj ) = Sj (k − kj ) + Hj Uj (k − kj − 1),
(16)
where (15) becomes 4 ˆ1 Aˆr x Sj (k − kj ) = E ˆ(k − kj − r) +
r X
ˆ1 Aˆi−1 D ˆ 1 w(k − kj − i) − y0 (k − kj ) + H0 U 0 (k − kj − 1) E j j
(17)
i=1
and (13) becomes
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
¯U ¯ (k − kj − 1) = Hj0 Uj0 (k − kj − 1) + Hj Uj (k − kj − 1), H
(18)
where Hj0 ∈ Rlz ×(rluˆ −lUj ) , Hj ∈ Rlz ×lUj , Uj0 (k − kj − 1) ∈ Rrluˆ −lUj , and Uj (k − kj − 1) ∈ RlUj . Now, by stacking z(k − k1 ), . . . , z(k − ks ), we define the extended performance h iT 4 Z(k) = z T (k − k1 ) · · · z T (k − ks ) ∈ Rslz . (19) Therefore, 4 ˜ ˜U ˜ (k − 1), Z(k) = S(k) +H (20) h iT 4 ˜ ˜ ∈ Rslz ×lU˜ , and U ˜ (k − 1) ∈ RlU˜ . The vector U ˜ (k − 1) where S(k) = ST (k − k1 ) · · · ST (k − ks ) ∈ Rslz , H is formed by stacking U1 (k − k1 − 1), . . . , Us (k − ks − 1) repetitions of components. For " and removing # u ˆ(k − 1) example, with k1 = 0 and k2 = 1, stacking U1 (k − 1) = and U2 (k − 2) = u ˆ(k − 2) results in u ˆ(k − 2) " # u ˆ(k − 1) ˜ consists of the entries of H1 , . . . , Hs arranged according ˜ U (k − 1) = . The coefficient matrix H u ˆ(k − 2) ˜ (k − 1). Furthermore, we assume that the last entry of U ˜ (k − 1) is a component of to the structure of U u ˆ(k − r).
Next, we define the retrospective performance 4
zˆ(k − kj ) = Sj (k − kj ) + Hj Uj∗ (k − kj − 1),
(21)
where the actual past subsystem outputs Uj (k − kj − 1) in (16) are replaced by the surrogate subsystem outputs Uj∗ (k − kj − 1). The extended retrospective performance for (21), which is defined as 4 ˆ Z(k) =
h
zˆT (k − k1 ) · · ·
zˆT (k − ks )
iT
∈ Rslz ,
(22)
is given by ˜ ˜U ˆ ˜ ∗ (k − 1), Z(k) = S(k) +H
(23)
˜ ∗ (k − 1) ∈ RlU˜ are components of U ∗ (k − k1 − 1), . . . , Us∗ (k − ks − 1) ordered in where the components of U 1 ˜ (k − 1). Subtracting (20) from (23) yields the same way as the components of U ˜U ˜U ˆ ˜ (k − 1) + H ˜ ∗ (k − 1). Z(k) = Z(k) − H
(24)
Finally, we define the retrospective cost function 4
˜ ∗ (k − 1), k) = Zˆ T (k)R(k)Z(k), ˆ J(U
(25)
where R(k) ∈ Rslz ×slz is a positive-definite performance weighting. The goal is to determine refined sub˜ ∗ (k − 1) that would have provided better performance than the subsystem outputs U (k) system outputs U ˜ ∗ (k − 1) are subsequently used to that were applied to the system. The refined subsystem outputs values U update the subsystem estimate. 4 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
A.
Cost Function Optimization with Adaptive Regularization To ensure that (25) has a global minimizer, we consider the regularized cost 4 ¯U ˜ ∗ (k − 1), k) = ˆ ˜ ∗T (k − 1)U ˜ ∗ (k − 1), J( Zˆ T (k)R(k)Z(k) + η(k)U
(26)
where η(k) = η¯z T (k)z(k) and η¯ ≥ 0. Substituting (24) into (26) yields ¯U ˜ ∗ (k − 1), k) = U ˜ ∗ (k − 1)T A(k)U ˜ ∗ (k − 1) + B(k)U ˜ ∗ (k − 1) + C(k), J(
(27)
where 4 ˜T ˜ + η(k)Il , A(k) = H R(k)H ˜ U 4
(28)
˜T
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
˜U ˜ (k − 1)], B(k) = 2H R(k)[Z(k) − H
(29)
4
˜U ˜ T R(k)H ˜U ˜ (k − 1) + U ˜ T (k − 1)H ˜ (k − 1). C(k) = Z T (k)R(k)Z(k) − 2Z T (k)R(k)H
(30)
˜ has full column rank or η(k) > 0, then A(k) is positive definite. In this case, J( ¯U ˜ ∗ (k − 1), k) has If either H the unique global minimizer ˜ ∗ (k − 1) = − 1 A−1 (k)B(k). U 2
B.
(31)
Subsystem Modeling The subsystem output u ˆ(k) is given by the exactly proper time-series model of order nc given by u ˆ(k) =
nc X
Mi (k)ˆ u(k − i) +
nc X
Ni (k)ˆ y (k − i),
(32)
i=0
i=1
where, for all i = 0, 1, . . . , nc , Mi (k) ∈ Rluˆ ×luˆ , Ni (k) ∈ Rluˆ ×lyˆ . The subsystem output (32) can be expressed as u ˆ(k) = θ(k)φ(k), where θ(k) ∈ Rluˆ ×[nc (luˆ +lyˆ )+lyˆ ] and φ(k) ∈ Rnc (luˆ +lyˆ )+lyˆ , 4
θ(k) = [M1 (k) · · · Mnc (k) N0 (k) · · · Nnc (k)] , T 4 T φ(k) = u ˆ (k − 1) · · · u ˆT (k − nc ) yˆT (k) · · · yˆT (k − nc ) .
(33)
u ˆ(k) = θ(k)ˆ y (k) = N0 (k)ˆ y (k).
(35)
(34)
If nc = 0, then
C.
Recursive Least Squares Update ˜ ∗ (k − 1) contains u∗ (k − d), and define the retrospective cost function Let d > 0 such that U 4
JR (θ(k)) =
k X
λk−i ||u∗T (k − d) − φT (k − d − 1)θT (k)||2 + λk (θ(k) − θ(0))P −1 (0)(θ(k) − θ(0))T ,
(36)
i=1
where φ(k − d) is given by (34), k · k is the Euclidean norm, and λ(k) ∈ (0, 1] is the forgetting factor. Minimizing the cumulative cost function yields retrospective cost optimization (RCO) θT (k) = θT (k − 1) + P (k − 1)φ(k − d − 1) · [φT (k − d)P (k − 1)φ(k − d − 1) + λ(k)]−1 · (u∗ (k − d) − φT (k − d − 1)θT (k − 1)).
5 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
(37)
The error covariance is updated by P (k) = λ−1 (k)P (k − 1) − λ−1 (k)P (k − 1)φ(k − d − 1) · [φT (k − d − 1)P (k − 1)φ(k − d) + λ(k)]−1 · φT (k − d − 1)P (k − 1).
(38)
We initialize the error covariance matrix as P (0) = βI, where β > 0.
IV.
Mass-Spring Example: Stiffness Estimation
Example IV.1. We consider the mass-spring-damper structure modeled by
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
m¨ q + cq˙ + κq = w,
(39)
where m, c, and κ are the mass, damping, and stiffness, respectively, and w is the force input. We obtain the state-space representation #" # " # " # " 0 q 0 1 q˙ + w, (40) = 1 κ c q˙ −m −m q¨ m " # h i q + v, (41) y= 0 1 q˙ where q and q˙ are the position and velocity, respectively, of the mass. The parameters are chosen as m = 1, c = 5, and κ = 10. The goal is to estimate κ with the initial estimate of κ ˆ (0) = 0. Using Euler discretizion with sample time Ts yields " # " #" # " # x1 (k + 1) 1 Ts x1 (k) 0 = + Ts w(k), s s x2 (k + 1) − κT 1 − cT x2 (k) m m m " # h i x (k) 1 + v(k), y(k) = 0 1 x2 (k) 4
4
(42)
(43)
4
"
4
0 − Tms
#
ˆ ˆ ˆ κC, ˆ where B ˆ = where, x1 (k) = q(kTs ) and x2 (k) = q(kT ˙ s ). We define Aclosed−loop = A + Bˆ h i 4 and Cˆ = 1 0 . The estimate of κ is given by κ ˆ (k) = θ(k). The initial conditions are assumed to be x1 (0) = 0.1 and x2 (0) = 0.01, which are unknown. Let the sensor noise v be Gaussian h white noise i with 2 ˜ mean µv and variance σ . Furthermore, P (0) = 0.6, λ = 0.999, η = 0, µv = 0, and H = H H , which 1
v
2
ˆ The input is w(k) = 0.8 sin(πk/20). RCMR is turned on are the first and second Markov parameters of G. at k = 20 steps. Figure 2 shows the performance z and the estimated κ ˆ. We next consider the effect of sensor noise. We consider the same estimation problem with white 4 σ2
sensor noise. We define SNR = σy¯20 , where σy2¯0 is the variance of the output signal y¯0 . Figure 3 shows the v effect of sensor noise on RCMR. As sensor noise is increased, the accuracy of the parameter identification is reduced.
V.
Estimation of a Repeated Parameter
We next consider the case where one unknown parameter appears in multiple locations within the model. Consider the continuous-time system x˙ = Ac x + Dc w,
6 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
(44)
0.2
20 RCMR estimate Initial estimate True parameter
0.15 Parameter Estimate
Performance z(k)
15 0.1 0.05 0 −0.05
10
5
0 −0.1 −0.15
0
100
200 Time Step (k)
300
−5
400
0
100
200 Time Step (k)
300
400
20
20
20
RCMR estimate Initial estimate True parameter
RCMR estimate Initial estimate True parameter
5
0
15
Parameter Estimate
10
−5
RCMR estimate Initial estimate True parameter
15
Parameter Estimate
15
Parameter Estimate
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
ˆ Figure 2. Example IV.1, performance z(k) and estimated k.
10
5
10
5
0
0
50
100
150
200 250 Time Step (k)
300
350
400
−5
0
0
50
100
(a) SNR=100
150
200 250 Time Step (k)
300
350
400
−5
(b) SNR=10
0
50
100
150
200 250 Time Step (k)
300
350
(c) SNR=1
Figure 3. Example VI.1, the estimate m ˆ of m for various values of SNR. As SNR decreases, the accuracy of the estimate is degraded.
where a parameter α in Ac and Dc is uncertain. Let α = α ˆ + ∆α, where α ˆ is the initial estimate of α. The system model can be written approximately as x˙ = (Ac +
∂Ac ∂Dc ∆α)x + (Dc + ∆α)w. ∂α α=αˆ ∂α α=αˆ
(45)
Discretizing (45) with sampling time Ts yields ∂Dc ∂Ac ∆αTs )x(k) + (Dc Ts + ∆αTs )w(k) x(k + 1) = (I + Ac Ts + ∂α α=αˆ ∂α α=αˆ ∂Ac ∂Dc = (A + ∆αTs )x(k) + (D1 + ∆αTs )w(k). ∂α α=αˆ ∂α α=αˆ
(46)
We use (46) as the true system model with uncertain parameter ∆α. Since our goal is to estimate ∆α, we rewrite (46) w(k) i h ∂Dc c x(k + 1) = Ax + D1 ∂A (47) ∆αx(k) . ˆ Ts ˆ Ts ∂α |α=α ∂α |α=α ∆αw(k)
7 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
400
Next, let ˆ= B
h "
yˆ(k) =
∂Ac ˆ Ts ∂α |α=α
x ˆ(k) w(k)
∂Dc ˆ Ts ∂α |α=α
i
,
(48)
# ,
(49)
ˆ s (q) = θ(k), ∆α(k) = G
(50)
u ˆ(k) = ∆α(k)ˆ y (k).
(51)
Note that (51) shows that u ˆ(k) is a function of yˆ(k). We define
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
4 Yˆ (k) =
h
yˆT (k − k1 ) · · ·
yˆT (k − ks )
iT
∈ Rsly .
(52)
With these signals, (27)-(31) become ∗ ¯ J(∆α (k − 1), k) = A(k)∆α∗2 (k − 1) + B(k)∆α∗ (k − 1) + C(k),
(53)
4
˜ T R(k)H ˜ + η(k)Il ]Yˆ (k − 1), A(k) = Yˆ (k − 1)T [H ˜ U
(54)
4
˜ T R(k)[Z(k) − H ˜U ˜ (k − 1)]Yˆ (k − 1), B(k) = 2H 4
(55)
˜U ˜U ˜ (k − 1) + U (k − 1)H R(k)H ˜ (k − 1), C(k) = Z (k)R(k)Z(k) − 2Z (k)R(k)H ˜ ∗ (k − 1) = ∆α∗ (k − 1)Yˆ (k − 1) = − 1 A−1 (k)B(k)Yˆ (k − 1). U 2 T
˜T
˜T
T
(56) (57)
˜ ∗ (k − 1) and the corresponding parameter in Y˜ (k − 1) to execute the recursive We use one parameter in U least squares update. In the next section we use the same mass-spring system to illustrate the algorithm.
VI.
Mass-Spring Example: Inertia Estimation
Example the mass-spring-damper structure shown in Example IV.1 where Ac = # VI.1." Consider # h i 0 1 0 , D = , and E = . In this example we assume that the parameter m is 0 1 c c κ 1 c −m −m m unknown. We demonstrate the algorithm by choosing κ = 30 and c = 5, and we assume that m = 0.9. We use an initial estimate is m(0) ˆ = 1, so that ∆m(0) = −0.1. "
From (48), we obtain " ˆ= B
0 3
0 0.5
0 −0.1
# ,
(58)
i , which are the first and " # x ˆ (0) 1 ˆ We choose the ramp input w(k) = 0.1k and the initial state second Markov parameters of G. = x ˆ2 (0) " # " # x1 (0) 0 = . The system refinement algorithm is turned on at k = 100 steps. Figure 4 shows the x2 (0) 0 performance z and the estimate m. ˆ ˜ = We choose P (0) = 1, λ = 1, η = 0, µ = 0, σv2 = 0 and H
h
H1
H2
We next consider the effect of sensor noise. We assume µv = 0. Figure 5 shows the estimation performance for several values of σv2 .
8 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
−4
4
x 10
1.1
2 Parameter Estimate
0 Performance z(k)
RCMR estimate Initial estimate True parameter
1.05
−2 −4 −6
1 0.95 0.9
−8 0.85
−10 −12
0
50
100 150 200 Time Step (k)
250
0.8
300
0
50
100 150 200 Time Step (k)
250
300
1.1
1.1 RCMR estimate Initial estimate True parameter
0.95
0.9
0.85
0.8
1.05
Parameter Estimate
1
1.1 RCMR estimate Initial estimate True parameter
1
0.95
0.9
0.85
0
50
100
150 Time Step (k)
200
250
300
0.8
RCMR estimate Initial estimate True parameter
1.05
Parameter Estimate
1.05
Parameter Estimate
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
Figure 4. Example VI.1, the performance z and the estimate m. ˆ
1
0.95
0.9
0.85
0
(a) σv2 = 1 × 10−10
50
100
150 Time Step (k)
200
250
300
(b) σv2 = 1 × 10−9
0.8
0
50
100
150 Time Step (k)
200
250
(c) σv2 = 1 × 10−8
Figure 5. Example VI.1, the estimate m ˆ of m for various values of σv2 . As σv2 increases, the accuracy of the estimate is degraded.
VII.
Aircraft Examples
We now estimate stability parameters under various scenarios of measurements noise. Consider the linearized longitudinal B-767 model u˙ −0.0168 0.1121 0.0003 −0.5608 u −0.0243 0.0519 0.0015 α −0.0634 −0.0005 α˙ −0.0164 −0.7771 0.9945 = + q˙ −0.0417 −3.6595 −0.9544 q −3.6942 0.0243 0 θ˙ 0 0 1.0000 0 θ 0 0
and measurement
"
δe δT
# ,
(59)
where δe is the elevator deflection and δT is the thrust perturbation. We discretize (59) using a zero-order hold with Ts = 0.1s, which yields 0.9983 0.0110 −0.0022 −0.0560 −0.0022 0.0052 0.0907 0.0002 −0.0018 0.9085 −0.0234 0.0001 A= (60) , D1 = , −0.0037 −0.3336 0.8924 −0.3492 0.0023 0.0001 −0.0002 −0.0172 0.0948 1.0000 −0.0178 0.0001 h iT and w = δe δT is the system input. For all examples in this section, RCMR is turned on at k = 100 h iT steps, and x ˆ(0) = x(0) = 0 0 0 0 . We choose δe and δT to be zero-mean white Gaussian noises with variance 0.001. 9 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
300
Example VII.1. Assume an initial guess of A(3, 3) = 0.8000. Theh true value of iA(3, 3) is 0.8924 and ∆A(3, 3) = 0.0924. Our goal is to identify A(3, 3). We choose E1 = 0 0 1 0 so that y¯0 = q h iT h i ˆ 1 = D1 , E ˆ1 = E1 , B ˆ= 0 0 1 0 and yˆ0 = qˆ. In this example, Aˆ = A, D and Cˆ = 0 0 1 0 . 2 ˆ In this h model, nci = 0 and ∆A(3, 3) = θ(k). We choose P (0) = 200, λ = 0.97, η = 0, µv = 0, σv = 0 and ˜ = H H , which are the first and second Markov parameters of G. ˆ The parameters of the model H 1 2 refinement algorithm are chosen such that z(k) is minimized. Figure 6 shows the performance z and the ˆ 3). estimated A(3, 0.02
1 RCMR estimate Initial estimate True parameter
0.015 Parameter Estimate
Performance z(k)
0.005 0 −0.005
0.9
0.85
0.8 −0.01 −0.015
0
100
200 Time Step (k)
300
0.75
400
0
100
200 Time Step (k)
300
400
ˆ 3). Figure 6. Example VII.1, performance z(k) and estimated A(3,
Example VII.2. We next consider the case where multiple parameters in A are possibly uncertain h iT and K = ∆A has the form shown. We assume full state feedback and y = . Thus u α q θ ˆK ˆ C, ˆ where B ˆ = Cˆ = I4 and K(k) ˆ Aˆclosed−loop = Aˆ + B = θ(k), where K=
−0.05 0 0 0
−0.002 −0.04 −0.04 0
0 0 −0.04 0
0 0 0 −0.04
.
(61)
i H1 H2 . Figure 7 shows the performance z(k) and the parameter estimates θ(k). The parameter error kθ(k) − Kk2 converges to zero. ˜ = We choose P (0) = 100I4 , λ = 0.98, η = 1, µv = 0, σv2 = 0 and H
0.14
h
0.08 0.06
0.12
0.04 0.1 0.02 0.08
θ(k)
Parameter Error
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
0.95 0.01
0.06
0 −0.02 −0.04
0.04 −0.06 0.02 0
−0.08 0
500
1000 1500 Time Step (k)
2000
−0.1
0
500
1000 1500 Time Step (k)
2000
Figure 7. Example VII.2, parameter error kθ(k) − Kk2 and estimated θ(k).
Example VII.3. We next consider the case where parameters of the A matrix that are not being ˆ 2) = −0.3000, whereas the true value is −0.3336. We estimated have modeling errors. We assume A(3, 10 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
assume that this modeling error is unknown and that h i our goal is to identify A(3, 3). We choose P (0) = 100, 2 ˜ = H H , which are the first and second Markov parameters of G. ˆ λ = 1, η = 0, µv = 0, σv = 0 and H 1 2 Figure 8 shows that the modelling errors degrade the accuracy of the estimates of the targeted parameters. 0.015
1
0.005 0 −0.005
0.9
0.85
0.8
0
200
400 600 Time Step (k)
800
0.75
1000
0
200
400 600 Time Step (k)
800
1000
ˆ 3). Figure 8. Example VII.3, performance z and estimated A(3,
Example VII.4. We next consider the case where the unknown parameter is time-varying. In this h iT ˆ 1 = D1 , E ˆ1 = E1 , B ˆ = 0 0 1 0 example, we assume A(3, 3) is time-varying, Aˆ = A, D , and h i Cˆ = 0 0 1 0 . Figure 9 shows the performance when A(3, 3) is varying at a constant rate, and Figure 10 shows the performance when A(3, as a sinusoidal signal. In both cases, P (0) = 10, h 3) is varying i 2 ˜ λ = 0.98, η = 0, µv = 0, σv = 0 and H = H1 H2 . −4
2
x 10
0.9
1.5
RCMR estimate Initial estimate True parameter
0.895 Parameter Estimate
1 Performance z(k)
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
−0.01 −0.015
RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
Performance z(k)
0.01
0.5 0 −0.5
0.89 0.885 0.88
−1 0.875
−1.5 −2
0
500
1000 1500 Time Step (k)
2000
0.87
0
500
1000 1500 Time Step (k)
2000
ˆ 3). A(3, 3) is the ramp signal A(3, 3) = 0.8924 − 10−5 k, Figure 9. Example VII.4, performance z and estimated A(3, 2 ˜ P (0) = 10, λ = 0.98, η = 0, µv = 0, σv = 0 and H = [H1 , H2 ].
Example VII.5. We next consider Example VII.1 with sensor noise present. We assume µv = 0. Figure 11 shows the estimation performance for several values of SNR. As the SNR decreases, the accuracy of the parameter estimatition degrades and convergence time increases. Example VII.6. We next consider Example VII.1 with non-zero-mean measurement noise. We assume SNR=100 and compare the estimate for several values of µv . Figure 12 shows the affect on estimation performance as µv increases. As the bias µv increases, the accuracy of the parameter estimation degrades and convergence time increases.
11 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
0.015
0.005 0 −0.005
0.95
0.9
0.85
−0.01 −0.015
RCMR estimate Initial estimate True parameter
1 Parameter Estimate
Performance z(k)
0.01
0
500
1000 1500 Time Step (k)
0.8
2000
0
500
1000 1500 Time Step (k)
2000
1 RCMR estimate Initial estimate True parameter
0.9
0.85
0.8
0.75
0.9
0.85
0.8
0
500
1000 Time Step (k)
1500
0.75
2000
1 RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
0.95 Parameter Estimate
Parameter Estimate
0.95
1 RCMR estimate Initial estimate True parameter
0.9
0.85
0.8
0
(a) SNR=∞
500
1000 Time Step (k)
1500
0.75
2000
RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
1
0.9
0.85
0.8
0
(b) SNR=100
500
1000 Time Step (k)
1500
0.75
2000
0
(c) SNR=10
2000
4000 6000 Time Step (k)
8000
10000
(d) SNR=1
Figure 11. Example.VII.5. (a) P (0) = 200, λ = 0.97. (b) P (0) = 200, λ = 1. (c) P (0) = 100, λ = 1. (d) P (0) = 15, ˜ = [H1 , H2 ]. λ = 1. In all cases, η = 0, and H
1 RCMR estimate Initial estimate True parameter
0.85
0.8
0.9
0.85
0.8
0
500
1000 Time Step (k)
(a) µv = 0
1500
2000
0.75
1 RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
0.9
0.75
0.95 Parameter Estimate
0.95
1 RCMR estimate Initial estimate True parameter
0.9
0.85
0.8
0
500
1000 Time Step (k)
1500
2000
0.75
RCMR estimate Initial estimate True parameter
0.95 Parameter Estimate
1
Parameter Estimate
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
ˆ 3). Ad (3, 3) is the sinusoidal signal A(3, 3) = Figure 10. Example VII.4, performance z and estimated A(3, ˜ = [H1 , H2 ]. 0.8924 + 0.08 sin(πk/400), P (0) = 10, λ = 0.98, η = 0, µv = 0, σv2 = 0 and H
0.9
0.85
0.8
0
(b) µv = −0.001
500
1000 Time Step (k)
1500
(c) µv = −0.01
2000
0.75
0
2000
4000 6000 Time Step (k)
8000
10000
(d) µv = −0.1
Figure 12. Example.VII.6. (a) P (0) = 200, λ = 0.97. (b) P (0) = 200, λ = 1. (c) P (0) = 200, λ = 1. (d) P (0) = 20, ˜ = [H1 , H2 ]. λ = 1. In all cases, η = 0, and H
VIII.
Airspeed Estimation
In this section we use RCMR to estimate airspeed. Example VIII.1. Consider the linearized longitudinal transfer functions [16] for a typical business
12 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
by
Xu0 + XTu0
Zu0 +ZTu0 U0 = Mu0 + MTu 0 0
Xα0
Xq0
Zα0 U0
U0 +Zq0 U0
Mα0 + MTα0 0
Mq0 1
−g
0 0 0
u α q θ
Xδe0
Zδe0 U0 + Mδe0 0
δe,
(62)
where U0 is the aircraft speed at which the equations of motion of the aircraft are linearized. X, Z, M are the related coefficients of the aircraft dynamic mode and many of these coefficients are depend on U0 . Using the relations of these coefficients to U0 , (62) can be written as − 4.9950 8.9782 0 −32.1522 0 u u˙ U0 1 − 1255.3650 − 93.8250 − 445.7224 0 α − 42.1968 α˙ U0 U02 U02 U0 + (63) = δe 634.2975 0.6075 q˙ −7.4416 − U0 0 q −17.6737 U0 θ θ˙ 0 0 0 1 0 ˆ0 = 675 The goal is to estimate U0 . The true value of U0 (0) is 625 ft/s, and we use the initial estimate U ft/s. We choose various y¯0 , and compare the performance of the algorithm. Figure 13 shows the performance of RCMR with y¯0 = α.h In this case, the i input w(k) = 0.0008, 2 ˜ P (0) = 0.0003, λ = 0.9999, η = 0, µ = 0, σv = 0 and H = H1 · · · H12 . Figure 14 shows the performance of RCMR with y¯i0 = q. In this case, the input w(k) = 10−7 k, P (0) = 1, λ = 0.99, η = 0, µ = 0, h ˜ = H H . Figure 15 shows the performance of RCMR with y¯0 = θ. In this case, the σv2 = 0 and H 1 2 h i −7 ˜ = H H . input w(k) = 10 k, P (0) = 0.01, λ = 0.999, η = 0, µ = 0, σv2 = 0, and H 1 2 −5
2
x 10
680 RCMR estimate Initial estimate True parameter
670 Parameter Estimate
0 Performance z(k)
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
jet in cruise given u˙ α˙ q˙ θ˙
−2
−4
660 650 640 630
−6 620 −8
0
1000
2000 3000 4000 Time Step (k)
5000
6000
610
0
1000
2000 3000 4000 Time Step (k)
5000
6000
ˆ0 are shown. Figure 13. Airspeed estimation with y¯0 = α. The performance z and the airspeed estimate U
IX.
Conclusions
This paper showed that RCMR can be used to estimate unknown parameters in a state space model. Both constant and time-varying parameters were considered, under various sensor noise levels and choices of measurements. Future work will compare the accuracy of this technique to nonlinear estimation methods.
Acknowledgment This work was supported in part by the National Aeronautics and Space Administration under Cooperative Agreement NNX12AM54A. 13 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
−7
4
x 10
680
Parameter Estimate
Performance z(k)
0 −2 −4 −6 −8
RCMR estimate Initial estimate True parameter
670
2
660 650 640 630 620
0
1000
2000 3000 4000 Time Step (k)
5000
610
6000
0
1000
2000 3000 4000 Time Step (k)
5000
6000
−5
1.5
x 10
700
Parameter Estimate
0.5 0 −0.5 −1 −1.5
RCMR estimate Initial estimate True parameter
680
1 Performance z(k)
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
ˆ0 are shown. Figure 14. Airspeed estimation with y¯0 = q. The performance z and the airspeed estimate U
660 640 620 600 580
0
1000
2000 3000 4000 Time Step (k)
5000
6000
560
0
1000
2000 3000 4000 Time Step (k)
5000
6000
ˆ0 are shown. Figure 15. Airspeed estimation with y¯0 = θ. The performance z and the airspeed estimate U
References 1 Palanthandalam-Madapusi, H., Renk, E. L., and Bernstein, D. S., “Data-Based Model Refinement for Linear and Hammerstein Systems Using Subspace Identification and Adaptive Disturbance Rejection,” Proc. Conf. Contr. Appl., Toronto, Canada, August 2005, pp. 1630–1635. 2 D’Amato, A. M., Arritt, B. J., Banik, J. A., Ardelean, E. V., and Bernstein, D. S., “Structural Health Determination and Model Refinement for a Deployable Composite Boom,” Proc. AIAA SDM Conf., Palm Springs, CA, April 2009, AIAA2009-2373. 3 Aljanaideh, K., Ali, A. A., Holzel, M., Kukreja, S., and Bernstein, D. S., “Sensor-to-Sensor Identification of Hammerstein Systems,” Proc. Conf. Dec. Contr., Maui, HI, December 2012. 4 Haroon, Y. W. L. M., Adams, D. E., and Ferri, A. A., “A Time and Frequency Domain Approach For Identifying Nonlinear Mechanical System Models in the Absence of an Input Measurement,” Sound and Vibration, Vol. 283, 2005, pp. 1137– 1155. 5 Moheimani, S. O. R., “Model Correction for Sampled-data Models of Structures,” Journal of Guidance, Control, and Dynamics, Vol. 24, 2001, pp. 634–637. 6 Danforth, C. M., Kalnay, E., and Miyoshi, T., “Estimating and Correcting Global Weather Model Error,” Monthly Weather Rev., Vol. 135, 2007, pp. 281–299. 7 Ljung, L., “Asymptotic Behavior of the Extended Kalman Filter as a Parameter Estimator for Linear Systems,” IEEE Trans. Autom. Control, Vol. 24, 1979, pp. 36–50. 8 Ristic, B., Arulampalam, S., and Gordon, N., Beyond the Kalman Filter , Artech House, 2004. 9 Julier, S., Uhlmann, J., and Durrant-Whyte, H., “A New Method for the Nonlinear Transformation of Means and Covariances in Filters and Estimators,” IEEE Trans. Autom. Control, Vol. 45, 2000, pp. 477–482. 10 Psiaki, M., “Backward-Smoothing Extended Kalman Filter,” AIAA J. Guid. Contr. Dyn., Vol. 28, 2005, pp. 885–894. 11 Rao, C., Rawlings, J., and Mayne, D., “Constrained State Estimation for Nonlinear Discrete-Time Systems: Stability and Moving Horizon Approximations,” IEEE Trans. Autom. Control, Vol. 48, 2003, pp. 246–258. 12 Santillo, M. A., D’Amato, A. M., and Bernstein, D. S., “System Identification Using a Retrospective Correction Filter for Adaptive Feedback Model Updating,” Proc. Amer. Contr. Conf., St. Louis, MO, June 2009, pp. 4392–4397.
14 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
Downloaded by UNIVERSITY OF MICHIGAN on August 24, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5193
13 D’Amato, A. M. and Bernstein, D. S., “Linear Fractional Transformation Identification Using Retrospective Cost Optimization,” Proc. SYSID, Saint-Malo, France, July 2009, pp. 450–455. 14 D’Amato, A. M., Wu, A. R., Mitchell, K. S., Kukreja, S. L., and Bernstein, D. S., “Damage Localization for Structural Health Monitoring Using Retrospective Cost Model Refinement,” Proc. AIAA SDM Conf., Orlando, FL, April 2010, AIAA2010-2628. 15 D’Amato, A. M., Ridley, A. J., and Bernstein, D. S., “Retrospective-Cost-Based Adaptive Model Refinement for the Ionosphere and Thermosphere,” Statistical Analysis and Data Mining, Vol. 4, 2011, pp. 446–458. 16 Roskam, J., Airplane Flight Dynamics and Automatic Flight Control, DARcorp, 2001.
15 of 15 American Institute of Aeronautics and Astronautics Copyright © 2013 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.
