- Email: [email protected]
A Novel Approach to the Identification of Exogenous Input of Stochastic Systems Using Pseudomeasurement
Proceedings of the 15th IFAC Symposium on System Identification Saint-Malo, France, July 6-8, 2009
A Novel Approach to the Identification of Exogenous Input of Stochastic Systems Using Pseudomeasurement Akira Ohsumi ∗ Takuro Kimura ∗∗ Michio Kono ∗∗∗ ∗
Graduate School of Engineering, Interdisciplinary Graduate School of Agriculture and Engineering, ∗∗∗ Interdisciplinary Graduate School of Agriculture and Engineering, University of Miyazaki, Kibana, Gakuen, Miyazaki 889-2192, Japan (e-mail: [email protected]; [email protected]) ∗∗
Abstract: In this paper, a novel approach to identify an unknown parameter vector of the exogenous input to a class of stochastic linear systems is proposed. The key of the approach is to introduce an additional information about the unknown parameter vector which is called the pseudomeasurement. Augmenting this pseudomeasurement with the original observation data, the identification of unknown vector as well as the state estimation is performed. The efficacy of the proposed approach is confirmed by simulation studies.
1. INTRODUCTION The idea of pseudomeasurement was first introduced in the state estimation, especially in target tracking in order to rise the estimation accuracy (Tahk and Speyer, 1990; Alouani and Bair, 1991; Ohsumi and Yasuki, 2000). The pseudomeasurement as opposed to the original observation data may be defined as additional information about the system state, and in general the kinematic constraint imposed on the system state is used. This additionally constructed information is used as additional observation data. This idea of pseudomeasurement has been primitively introduced by one of the authors to identify unknown magnitude of a suddenly discharged pollutive load to a river, related to the environmental problem (Ohsumi, Komiyama, Kashiwagi, Watanabe, and Takatsu, 2008). The water quality of the river are evaluated by the major indices BOD (biochemical oxygen demand) and DO (dissolved oxygen) of the river, which are estimated from the observation data made on BOD and DO. The pollutive load which is suddenly discharged by unlawful discharge or unforeseen accident is regarded as an exogenous input with unknown magnitude and further its discharged location is also unknown. Originating from such a research related to the environmental problem, an identification problem has been investigated for identifying the unknown exogenous input and estimating its discharged time for a class of deterministic systems, using the idea of pseudomeasurement (Kimura, Ohsumi, and Kono, 2008). From the system-theoretic point of view, this identification problem is recognized as one of inverse problems. For stochastic systems, the conventional approach will treat the unknown (constant) parameter as a new sate variable and augment it with the original system state
978-3-902661-47-0/09/$20.00 © 2009 IFAC
296
variables. Then, for such augmented system a Kalman filter will be constructed to estimate the augmented system states from the observation data, so the unknown parameter will be inevitably identified. Unfortunately, however, the estimation accuracy for the unknown parameter can not be necessarily attained satisfactorily as shown in simulation studies later. In this paper, as unknown exogenous input to the stochastic system, both stepwise and impulsive inputs and further a quite unknown input are considered. 2. IDENTIFICATION PROBLEM Consider an identification problem for the stochastic system which is described by Itˆo stochastic differential equations: dx(t) = Ax(t) dt + Bξ(t) dt + Cu(t) dt + G dw(t), x(0) = x0 ΣS : (1) dy(t) = Hx(t) dt + R dv(t), where x(t) ∈ Rn , y(t) ∈ Rm (m ≤ n), u(t) ∈ R` ; w(t) ∈ Rd1 and v(t) ∈ Rd2 are mutually independent Wiener processes with covariance matrices Q and I (unit matrix); ξ(t) ∈ Rm is the unknown exogeneous input to be identified; and x0 is the Gaussian random variable. Our problem is to identify the unknown exogenous input ξ(t) as well as the system state x(t) from a couple of input and output data {u(t), y(t)}. All matrices such as A, B and so on are time-invariant. The pair (A, H) is assumed to be observable, and H is of full rank m. As for the exogenous input ξ(t), we can model it such that ξ(t) = b uS (t − t0 ) (2) or ξ(t) = b δ(t − t0 )
(3)
10.3182/20090706-3-FR-2004.0079
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 if we know a priori that it is a constant or impulsive input, where b ∈ Rm is the unknown (constant) magnitude, and uS (·), δ(·) are the Heaviside unit step and the Dirac delta functions, respectively; and t0 (> 0) is some fixed or unknown onset time when the exogenous input is discharged. For simplicity, t0 is assumed to be known until it is noted. In the event of being no a priori information available about the exogenous input, it may be reasonable to assume it by a simple polynomial model of order p [ p ] ∑ ξ(t) = bi φi (t) uS (t − t0 ) (4) i=0
in which {φi (t)} is a set of (scalar) basis functions of time t, while {bi } is a set of unknown constant vectors, bi ∈ Rm . The model (4) can be expressed as ξ(t) = Φ(t) θ uS (t − t0 ),
(5)
where θ = [ bT0 , bT1 , . . . , bTp ]T ∈ R(p+1)m (T denotes transpose),
How can we obtain any information about the unknown parameter vector b or θ? Up to here we have an only information that it is constant. In order to get more information about b (or θ) we make the most of system dynamics. To do this, consider the case of (2) first. Then, the solution is given for t ≥ t0 by t ∫ x(t) = eA(t−t0 ) x(t0 ) + eA(t−τ ) B dτ b t0
eA(t−τ ) Cu(τ ) dτ + {noise}
+
= φT (t) ⊗ Im , (6) (Im : unit matrix of dimension m, ⊗: Kronecker product), and φ(t) = [ φ0 (t), φ1 (t), · · ·, φp (t) ]T . The conventional approach to identify the unknown constant vector can be stated as follows. For the stochastric system ΣS , it is a custom to use the expression for b (or θ) as b(t) (or θ(t)) and to write as ˙ = Gb γb (t), b(0) = b0 b(t) (7)
Consider the case of (2). Augmenting (8) with the system dynamics (1), we have for the augmented state vector z(t) = [ xT (t), bT (t) ]T ∈ Rn+m : dz(t) = A0 z(t) dt + B0 z(t)uS (t − t0 ) dt (9)
(11)
t0
in which the term concerning to the system noise is simply denoted by {noise} because this is not essential to get direct relation between b and x(t). In principle, we could determine b from this if the state {x(·)} is completely known (except {noise}). However, the true state is not available. Keep in mind that only the couple of control input and observation data {u(t), y(t)} are available. So, in order to get the relation between b and {y(t)}, uponmultiplying both sides of (11) by the matrix H dt to obtain
(8)
where γb (t) is the m-vector zero-mean Gaussian random process and wb (t) is the corresponding Wiener process, ∫t viz., wb (t) = 0 γb (τ ) dτ , with covariance matrix Qb . The expression (7) or (8) is introduced to reflect the ambiguity ˙ in the strict constraint b(t) ≡ 0 for all t. The initial value b0 is preassigned and may be set as b0 = 0.
+ C0 u(t) dt + G0 dw0 (t),
3. INTRODUCTION OF PSEUDOMEASUREMENT
∫t
Φ(t) = [ φ0 (t)Im , φ1 (t)Im , · · · , φp (t)Im ]
or in the stochastic differential form, db(t) = Gb dwb (t), b(0) = b0 ,
The identification of b may be performed by the Kalman filter constructed for (9) and (10). Unfortunately, this approach dooms to yield poor and unsatisfactory identification results in b. Unsuccessful results in the conventional approach are caused entirely by the lack of any information about the unknown parameter vector.
Hx(t) dt =
∫t
HeA(t−τ ) B dτ b dt
t0
+ HeA(t−t0 ) x(t0 ) dt +
∫t
HeA(t−τ ) Cu(τ ) dτ dt
t0
+ {noise}, (12) in which the L.H.S. should be equal to dy(t) − Rdv(t). Therefore, we have the relation between b and the observation data, t ∫ dy(t) = HeA(t−τ ) B dτ b dt t0
where w0 (t) = [ wT (t), wbT (t) ]T , and [ ] [ ] [ ] [ ] A0 0 B C G 0 A0 = , B0 = , C0 = , G0 = . 0 0 0 0 0 0 Gb
+
∫t
HeA(t−τ ) Cu(τ ) dτ dt + {noise}
(13)
t0
Then, the observation process is rewritten in terms of z(t) as dy(t) = Hb z(t) dt + R dv(t) (10) with Hb = [ H 0 ].
297
in which the {noise} term reflects now the collection of system noise, observation noise and the random initial state x(t0 ). The term HeA(t−t0 ) x(t0 ) becomes disappearing as time passes as far as the system matrix A is stable. So, this term is regarded as a part of random noise.
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 Subtracting the second term in the R.H.S. from both sides of (13), we have dy(t) −
∫t
By the quite similar way, the pseudomeasurement of the form (16) can be defined for the model (4) (or (5)) to get dyp (t) = Φ0 (t)θ(t) dt + Rp dvp (t) (t ≥ t0 ) (21)
HeA(t−τ ) Cu(τ ) dτ dt
t0
=
∫t
Hδ (t) = HeA(t−t0 ) and dyp (t) is the same as (15).
HeA(t−τ ) B dτ b dt + {noise}.
in which the matrix Φ0 (t) is (14)
∫t
t0
dyp (t) = dy(t) −
∫t
HeA(t−τ ) Cu(τ ) dτ dt.
(15)
(22)
t0
and the noise term is introduced similarly as in (16). In this case the augmented system for z(t) = [ xT (t), θT (t) ]T ∈ Rn+(p+1)m is given by
t0
dz(t) = A0 z(t) dt + Bθ (t)z(t)uS (t − t0 ) dt
Then, (13) can be rewritten, denoting b = b(t), as dyp (t) = Hp (t)b(t) dt + Rp dvp (t) (t ≥ t0 ),
HeA(t−τ ) BΦ(τ ) dτ
Φ0 (t) =
Here, note that the L.H.S. is a known quantity, and hence let it be denoted as dyp (t), i.e.,
(16)
+ C0 u(t) dt + G0 dw0 (t),
(23)
where
where
[
∫t
Bθ (t) =
HeA(t−τ ) B dτ
Hp (t) =
= HA−1 [ eA(t−t0 ) − I ]B in which the second equality holds if A is stable, and the term Rp dvp (t) is introduced with a newly introduced m-vector standard Wiener process vp (t) to describe the {noise} term. At this stage, it should be emphasized that (16) can be regarded just as an observation model for the state b(t) (Ohsumi, 2008). In this sense, the newly defined dyp (t) is called the pseudomeasurement for unknown parameter vector b. Clearly, {dyp (t)} is Yt -measurable (Yt = σ{y(τ ), 0 ≤ τ ≤ t}). In order to adapt such expression as (16) for all t > 0, we assume formally for 0 < t < t0 that dyp (t) = Rp dvp (t)
(0 < t < t0 )
(17)
in which dyp (t) ≡ 0 is set. For the case of impulsive input modeled by (3), noting that ∫t e
]
and other matrices are the same as that in the case of (2).
t0
A(t−τ )
0 BΦ(t) 0 0
b δ(τ − t0 )dτ = e
A(t−τ )
¯ ¯ b ¯¯
0
uS (t − t0 ), (18)
4. IDENTIFICATION OF UNKNOWN INPUT Augmenting the pseudomeasurements (16) and (17) with the original observation process to form the vector y0 (t) = [ y T (t), ypT (t) ]T , we obtain the augmented observation process for the case of (2): dy0 (t) = H0 (t)z(t) dt + R0 dv0 (t) (24) where [ H0 (t) =
H 0 0 χHp (t)
]
in which χ is an indicator taking 0 for 0 < t < t0 and 1 for t ≥ t0 ; v0 (t) = [ v T (t), vpT (t) ]T ; and R0 = block diag {R, Rp }. Then, we have the system formulation: dz(t) = A0 z(t) dt + B0 z(t)uS (t − t0 ) dt + C0 u(t) dt + G0 dw0 (t) dy0 (t) = H0 (t)z(t) dt + R0 dv0 (t).
(25)
τ =t0
Based on (25), the Kalman filter can be constructed to obtain the estimate of b (= b(t)):
we have x(t) = eA(t−t0 ) x(t0 ) + eA(t−t0 ) b
dˆ z (t|t) = Aχ zˆ(t|t) dt + C0 u(t) dt
∫t eA(t−τ ) Cu(τ ) dτ + {noise}
+
(t ≥ t0 ). (19)
P˙ (t|t) = Aχ P (t|t) + P (t|t)ATχ + G0 Q0 GT0
t0
By the similar way as in the case of (2), we have the version dyp (t) = Hδ (t)b(t) dt + Rp dvp (t)
+ P (t|t)H0T (t){R0 R0T }−1 {dy0 (t) − H0 (t)ˆ z (t|t) dt} (26)
(t ≥ t0 ),
(20)
where
298
− P (t|t)H0T (t){R0 R0T }−1 H0 (t)P (t|t), (27) where zˆ(t|t) = E{z(t)|Yt } and P (t|t) = E{[z(t) − zˆ(t|t)][z(t) − zˆ(t|t)]T } are the optimal estimate of z(t)
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 and the estimation error covariance matrix, respectively; Aχ = A0 for 0 < t < t0 and = A1 ≡ A0 + B0 for t ≥ t0 , and Q0 = block diag {Q, Qb }. For the case of (3) the Kalman filter can be also constructed mutatis mutandis. For the case of (4), the system formulation is given as (25) by replacing the matrices B0 and H0 (t) by Bθ (t) and Hθ given below, respectively; [ Hθ (t) =
] H 0 . 0 χΦ0 (t)
The Kalman filter is given by (26) and (27) by replacing H0 (t) in them by Hθ (t). In this case A0 (t) = A0 for 0 < t < t0 and A1 (t) = A0 + Bθ (t) for t ≥ t0 .
5. DETECTION OF ONSET TIME Another important task is to detect the onset time t0 at which an exogeneous input is discharged. To do this, revive t0 to be unknown, and suppose a system which is subjected to no exogeneous input and the same initial condition, dxf (t) = Axf (t) dt + Cu(t) dt + G dw(t), xf (0) = x0 . (28) Note that xf (t) ≡ x(t) until t = t0 because the same initial condition as in (1) is set. After t0 , xf (t) becomes no longer equal to the true x(t). Recalling this fact in mind, let us consider the following innovation process dν(t) = dy(t) − H x ˆf (t|t) dt,
ν(0) = 0
(29)
in which y(t) is the original process in (1) and x ˆf (t|t) is the estimation of xf (t) obtained by the Kalman filter with an appropriate initial estimate x ˆ0 : dˆ xf (t|t) = Aˆ xf (t|t) dt + Cu(t) dt + Pf (t|t)H T {RRT }−1 dν(t)
(30)
P˙f (t|t) = APf (t|t) + Pf (t|t)AT (t) + GQGT T −1
− Pf (t|t)H {RR } T
HPf (t|t).
dx(t) = Ax(t) dt + b {uS (t − t0 ) − uS (t − t1 )} dt + g dw(t), (33) where w(t) is a scalar noise process with variance q = 10: [ ] [ ] 0 1 1 A= , g= ; −1 −2.9 1 and the true value of b is set as b = [3, 2]T . For the observation process (2), H = R = I was set. The user-defined parameters were set as: Gb = I, Qb = 2I for the 2-vector noise process wb (t), Rp = I for the pseudomeasurement process (16) and (17). The initial values for the Kalman filter (26) and (27) were set as zˆ(0|0) = [ 0, 0, 0, 0 ]T and P (0|0) = diag {3, 4, 1, 2}. The simulation studies were performed with discretized version with time partition δt = 0.01 s. The true times t0 and t1 were set as t0 = 20, t1 = 29 s, respectively. To detect t0 and t1 the onset-time detector (32) was used. For the computation of the innovation process (29) the Kalman filter given by (30) and (31) was run with initial values x ˆf (0|0) = [ 0, 0 ]T and Pf (0|0) = diag {3, 4}. The user-defined parameter T0 in the detector (32) was chosen as T0 = 5 δt s. Figure 1 depicts the identification processes for the unknown parameters {bi }. In the figure, solid curves depict the estimation or identification, while the dotted curve the true value. Figure 2 depicts the behavior of the onset-time detector d(t). As seen in the figure, d(t) decreases rapidly after t = 20 s. If we set the threshold to be dthreshold = 0.5, then the estimate of t0 is decided as tˆ0 = 21.55, while for tˆ1 = 31.55 s for the same threshold value. These are sufficiently close to each true values. From Fig. 1 we see that the identification and state estimation is performed well.
(31)
(32)
6.2 Identification of quite-unknown exogenous input
∫t k ν(τ ) k dτ,
Consider the following second-order system in which the input u(t) is assumed to be zero, i.e., u(t) = 0 for all t:
Figure 3 shows the result of identification using the Kalman filter constructed based on the models (9) and (10) in which no idea of pseudomeasurement is used. As suspected, in this case the identification of unknown parameters is quite unsuccessful.
To detect t0 , define a function, d(t) =
6.1 Identification of stepwise exogenous input
t−T0
where k · k is the Euclidean norm and T0 is a short time interval. The onset instant t0 is detected as the time when d(t) exceeds a preassigned threshold. Denote it by tˆ0 and replace t0 in (26) and (27) by it.
The following third-order system is considered in which the input u(t) is assumed to be zero: dx(t) = Ax(t) dt + Bξ(t) dt + G dw(t),
where the variance of the w(t)-process is set as Q = diag {90, 120, 105}; [
6. SIMULATION STUDIES
A=
Owing to the space limit, only the cases of stepwise and quite unknown exogenous input are investigated to grasp the feeling how the use of pseudomeasurement is useful.
299
(34)
] −5 5 10 0 −10 5 , 5 0 −15
[ B=
and onset time t0 was set t0 = 0.
] 1 0 0 1 , 2 0
[ G=
] 60 0 03 0 ; 0 0 4.5
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009
Fig. 1. Identification features of b1 (top) and b2 (bottom).
Fig. 4. Identification of ξ1 (top) and ξ2 (bottom) whose true processes are given by (35). 3 1 1 1 5 1 6 ξ2 (t) = 3 + t − t2 − t3 + t4 + t − t . (35) 8 9 32 120 960 For the observation process, y(t) ∈ R2 , the matrices were set as [ H=
Fig. 2. Behavior of the detector d(t).
] 10 1 , 02 1
[ R=
] 2 0 . 0 1
For identifying the unknown inputs {ξ(t)}, it is assumed that ξi (t) =
p ∑
bik φk (t)
(i = 1, 2).
(36)
k=0
As for the basis functions {φk (t)}, the following simple polynominals of order p = 3 are selected: φk (t) = (t − tj )k for tj ≤ t < tj+1 (k = 0, 1, 2, 3), (37) where {tj } are equi-partitioned times such that t0 = 0 and tj+1 = tj + 1 s (i = 0, 1, 2, · · ·). The polynomial (37) has been modified from the conventionally used form, φk (t) = tk . This modified form is used here to guarantee the condition |φk (t)| < 1 for all t to prevent the divergence as time goes by. The unknown vector θ is given by Fig. 3. Unsuccessful results of identification b1 (top) and b2 (bottom) by the conventional method based on (9) and (10). The true exogenous input ξ(t) ∈ R2 was assumed as follows:
θ = [b10 , b20 , b11 , b21 , b12 , b22 , b13 , b23 ]T . For the model of θ(t), the matrices concerned with the wb (t)-process were assumed as user-defined parameters, Gb = diag {0.4, 0.3, 0.1, 0.2, 0.05, 0.06, 0.003, 0.004}, Qb = diag {8, 6, 2, 4, 1, 1.2, 0.06, 0.08}.
1 1 1 1 1 6 3 t ξ1 (t) = 3 + t − t2 − t3 + t4 + t5 − 2 2 4 24 80 720
Furthermore, as for the pseudomeasurement process (21), the matrix Rp was set as Rp = diag {3, 2}.
300
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 (ii) The pseudomeasurement approach can be used for deterministic systems. For the deterministic systems, observer theory has been developed for the system with unknown input by the name of disturbance observer or unknown input observer (e.g., Medich and Hostetter, 1974). In this theory the disturbance input is modeled by a constant, otherwise the polynomial model (4) with φi (t) = ti . However, the coefficient parameters {bi } are not directly identified, while in the pseudomeasurement approach the unknown input can be recovered surely. Though the adaptive observers can be constructed for identifying unknown input, they have sophisticated structures even in the case of single-input single-output (L¨ uders and Narendra, 1974; Kreisselmeier, 1977). As opposed to adaptive observers, the pseudomeasurement approach is rather simple to construct an observer or Kalman filter to identify the unknown input. Furthermore, impulsive inputs are not identified by the disturbance observer theory. Fig. 5. Identification of ξ1 (top) and ξ2 (bottom) given by (38). The initial value for the Kalman filter (26) and (27) were set as zˆ(0|0) = 0 and P (0|0) = diag {3, 4, 5, 1, 2, 1, 2, 1, 2, 1, 2}. Figure 4 shows the result of identification of the unknown inputs ξ1 (t) and ξ2 (t). The solid curves and dotted ones depict the identification and true values, respectively. Another simulation example is shown in Fig. 5. The system and observation are the same as the above example except the matrix ] [ −50 50 100 0 −100 50 . A= 50 0 −150 The true exogenous inputs are assumed as 1 1 1 1 ξ1 (t) = sin t + cos t, ξ2 (t) = sin t + cos t. (38) 30 40 40 30 These are identified by the Fourier series expansion (35) with p = 16 and k+1 cos t (k : odd) 2 φ0 (t) = 1, φk (t) = (39) sin k t (k : even). 2 As seen in Fig. 5 the identification of the quite-unknown inputs seems to be performed fairly well. 7. CONCLUDING REMARKS The problem to identify unknown exogenous input to the dynamical system is important and attracting much attention to many control engineers. (i) In general, the estimate of unknown input requires higher derivatives of the observation data as inverse system. However, such higher derivatives of observation data is quite unrealistic, particularly for the stochastic systems. The approach proposed in this paper which will be called the pseudomeasurement approach does not need any derivatives of observation data.
301
(iii) As well-known, the fault diagnosis is attracting much interest of process control engineers as one of recent topics. In the process control, the fault signal is regarded as unknown exogenous input and the adaptive observer is used to identify it. The pseudomeasurement approach will have potential ability to identify even such fault signals.
REFERENCES A. T. Alouani and W. D. Bair. Use of a Kinematic Constraint in Tracking Constant Speed, Maneuvering Target. Proc. 30th Conf. on Decision and Control, Brighton, England, 2055–2058, 1991. T. Kimura, A. Ohsumi, and M. Kono. Observer-based Identification of Unknown Exogenous Input via Pseudomeasurement Approach. Proc. SICE Annual Conf. 2008, Chofu, Tokyo, Aug., 2660–2665, 2008. G. Kreisselmeier. Adaptive Observers with Exponential Rate of Convergence. IEEE Trans. on Automatic Control, vol. AC-22, no. 1, 2–8, 1977. G. L¨ uders and K. S. Narendra. A New Canonical Form for an Adaptive Observer. IEEE Trans. on Automatic Control, vol. AC-19, no. 2, 117–119, 1974. J. S. Medich and G. H. Hostetter. Observers for Systems with Unknown and Inaccessible Inputs. Int. J. Control, vol. 19, no. 3, 473–480, 1974. A. Ohsumi. An Outlook on the Use of Pseudomeasurement in System Identification. Proc. 37th SICE Symp. on Control Theory, Kirishima, Sep., 91–96, 2008. A. Ohsumi and S. Yasuki. Tracking of a Maneuvering Target Considering its Kinematic Constraints. Intelligent Automation and Control—Recent Trends in Development and Applications, Vol. 9 (eds: M. Jamshidi, P.Borne, and J. S. Jamshidi), TSI Press, Albuquerque, NM, 407–416, 2000. A. Ohsumi, S. Komiyama, K. Kashiwagi, M. Watanabe, and T. Takatsu. Detection of Polluted Load and Identification of its Discharged Location and Magnitsude for Polluted River. Int. J. of Innovative Computing, Information and Control, vol. 4, no. 1, 63–77, 2008. M. Tahk and J. L. Speyer. Target Tracking Problems Subject to Kinematic Constraints. IEEE Trans. on Automatic Control, vol. 35, no. 3 (1990), 324–326, 1990.
A Novel Approach to the Identification of Exogenous Input of Stochastic Systems Using Pseudomeasurement Akira Ohsumi ∗ Takuro Kimura ∗∗ Michio Kono ∗∗∗ ∗
Graduate School of Engineering, Interdisciplinary Graduate School of Agriculture and Engineering, ∗∗∗ Interdisciplinary Graduate School of Agriculture and Engineering, University of Miyazaki, Kibana, Gakuen, Miyazaki 889-2192, Japan (e-mail: [email protected]; [email protected]) ∗∗
Abstract: In this paper, a novel approach to identify an unknown parameter vector of the exogenous input to a class of stochastic linear systems is proposed. The key of the approach is to introduce an additional information about the unknown parameter vector which is called the pseudomeasurement. Augmenting this pseudomeasurement with the original observation data, the identification of unknown vector as well as the state estimation is performed. The efficacy of the proposed approach is confirmed by simulation studies.
1. INTRODUCTION The idea of pseudomeasurement was first introduced in the state estimation, especially in target tracking in order to rise the estimation accuracy (Tahk and Speyer, 1990; Alouani and Bair, 1991; Ohsumi and Yasuki, 2000). The pseudomeasurement as opposed to the original observation data may be defined as additional information about the system state, and in general the kinematic constraint imposed on the system state is used. This additionally constructed information is used as additional observation data. This idea of pseudomeasurement has been primitively introduced by one of the authors to identify unknown magnitude of a suddenly discharged pollutive load to a river, related to the environmental problem (Ohsumi, Komiyama, Kashiwagi, Watanabe, and Takatsu, 2008). The water quality of the river are evaluated by the major indices BOD (biochemical oxygen demand) and DO (dissolved oxygen) of the river, which are estimated from the observation data made on BOD and DO. The pollutive load which is suddenly discharged by unlawful discharge or unforeseen accident is regarded as an exogenous input with unknown magnitude and further its discharged location is also unknown. Originating from such a research related to the environmental problem, an identification problem has been investigated for identifying the unknown exogenous input and estimating its discharged time for a class of deterministic systems, using the idea of pseudomeasurement (Kimura, Ohsumi, and Kono, 2008). From the system-theoretic point of view, this identification problem is recognized as one of inverse problems. For stochastic systems, the conventional approach will treat the unknown (constant) parameter as a new sate variable and augment it with the original system state
978-3-902661-47-0/09/$20.00 © 2009 IFAC
296
variables. Then, for such augmented system a Kalman filter will be constructed to estimate the augmented system states from the observation data, so the unknown parameter will be inevitably identified. Unfortunately, however, the estimation accuracy for the unknown parameter can not be necessarily attained satisfactorily as shown in simulation studies later. In this paper, as unknown exogenous input to the stochastic system, both stepwise and impulsive inputs and further a quite unknown input are considered. 2. IDENTIFICATION PROBLEM Consider an identification problem for the stochastic system which is described by Itˆo stochastic differential equations: dx(t) = Ax(t) dt + Bξ(t) dt + Cu(t) dt + G dw(t), x(0) = x0 ΣS : (1) dy(t) = Hx(t) dt + R dv(t), where x(t) ∈ Rn , y(t) ∈ Rm (m ≤ n), u(t) ∈ R` ; w(t) ∈ Rd1 and v(t) ∈ Rd2 are mutually independent Wiener processes with covariance matrices Q and I (unit matrix); ξ(t) ∈ Rm is the unknown exogeneous input to be identified; and x0 is the Gaussian random variable. Our problem is to identify the unknown exogenous input ξ(t) as well as the system state x(t) from a couple of input and output data {u(t), y(t)}. All matrices such as A, B and so on are time-invariant. The pair (A, H) is assumed to be observable, and H is of full rank m. As for the exogenous input ξ(t), we can model it such that ξ(t) = b uS (t − t0 ) (2) or ξ(t) = b δ(t − t0 )
(3)
10.3182/20090706-3-FR-2004.0079
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 if we know a priori that it is a constant or impulsive input, where b ∈ Rm is the unknown (constant) magnitude, and uS (·), δ(·) are the Heaviside unit step and the Dirac delta functions, respectively; and t0 (> 0) is some fixed or unknown onset time when the exogenous input is discharged. For simplicity, t0 is assumed to be known until it is noted. In the event of being no a priori information available about the exogenous input, it may be reasonable to assume it by a simple polynomial model of order p [ p ] ∑ ξ(t) = bi φi (t) uS (t − t0 ) (4) i=0
in which {φi (t)} is a set of (scalar) basis functions of time t, while {bi } is a set of unknown constant vectors, bi ∈ Rm . The model (4) can be expressed as ξ(t) = Φ(t) θ uS (t − t0 ),
(5)
where θ = [ bT0 , bT1 , . . . , bTp ]T ∈ R(p+1)m (T denotes transpose),
How can we obtain any information about the unknown parameter vector b or θ? Up to here we have an only information that it is constant. In order to get more information about b (or θ) we make the most of system dynamics. To do this, consider the case of (2) first. Then, the solution is given for t ≥ t0 by t ∫ x(t) = eA(t−t0 ) x(t0 ) + eA(t−τ ) B dτ b t0
eA(t−τ ) Cu(τ ) dτ + {noise}
+
= φT (t) ⊗ Im , (6) (Im : unit matrix of dimension m, ⊗: Kronecker product), and φ(t) = [ φ0 (t), φ1 (t), · · ·, φp (t) ]T . The conventional approach to identify the unknown constant vector can be stated as follows. For the stochastric system ΣS , it is a custom to use the expression for b (or θ) as b(t) (or θ(t)) and to write as ˙ = Gb γb (t), b(0) = b0 b(t) (7)
Consider the case of (2). Augmenting (8) with the system dynamics (1), we have for the augmented state vector z(t) = [ xT (t), bT (t) ]T ∈ Rn+m : dz(t) = A0 z(t) dt + B0 z(t)uS (t − t0 ) dt (9)
(11)
t0
in which the term concerning to the system noise is simply denoted by {noise} because this is not essential to get direct relation between b and x(t). In principle, we could determine b from this if the state {x(·)} is completely known (except {noise}). However, the true state is not available. Keep in mind that only the couple of control input and observation data {u(t), y(t)} are available. So, in order to get the relation between b and {y(t)}, uponmultiplying both sides of (11) by the matrix H dt to obtain
(8)
where γb (t) is the m-vector zero-mean Gaussian random process and wb (t) is the corresponding Wiener process, ∫t viz., wb (t) = 0 γb (τ ) dτ , with covariance matrix Qb . The expression (7) or (8) is introduced to reflect the ambiguity ˙ in the strict constraint b(t) ≡ 0 for all t. The initial value b0 is preassigned and may be set as b0 = 0.
+ C0 u(t) dt + G0 dw0 (t),
3. INTRODUCTION OF PSEUDOMEASUREMENT
∫t
Φ(t) = [ φ0 (t)Im , φ1 (t)Im , · · · , φp (t)Im ]
or in the stochastic differential form, db(t) = Gb dwb (t), b(0) = b0 ,
The identification of b may be performed by the Kalman filter constructed for (9) and (10). Unfortunately, this approach dooms to yield poor and unsatisfactory identification results in b. Unsuccessful results in the conventional approach are caused entirely by the lack of any information about the unknown parameter vector.
Hx(t) dt =
∫t
HeA(t−τ ) B dτ b dt
t0
+ HeA(t−t0 ) x(t0 ) dt +
∫t
HeA(t−τ ) Cu(τ ) dτ dt
t0
+ {noise}, (12) in which the L.H.S. should be equal to dy(t) − Rdv(t). Therefore, we have the relation between b and the observation data, t ∫ dy(t) = HeA(t−τ ) B dτ b dt t0
where w0 (t) = [ wT (t), wbT (t) ]T , and [ ] [ ] [ ] [ ] A0 0 B C G 0 A0 = , B0 = , C0 = , G0 = . 0 0 0 0 0 0 Gb
+
∫t
HeA(t−τ ) Cu(τ ) dτ dt + {noise}
(13)
t0
Then, the observation process is rewritten in terms of z(t) as dy(t) = Hb z(t) dt + R dv(t) (10) with Hb = [ H 0 ].
297
in which the {noise} term reflects now the collection of system noise, observation noise and the random initial state x(t0 ). The term HeA(t−t0 ) x(t0 ) becomes disappearing as time passes as far as the system matrix A is stable. So, this term is regarded as a part of random noise.
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 Subtracting the second term in the R.H.S. from both sides of (13), we have dy(t) −
∫t
By the quite similar way, the pseudomeasurement of the form (16) can be defined for the model (4) (or (5)) to get dyp (t) = Φ0 (t)θ(t) dt + Rp dvp (t) (t ≥ t0 ) (21)
HeA(t−τ ) Cu(τ ) dτ dt
t0
=
∫t
Hδ (t) = HeA(t−t0 ) and dyp (t) is the same as (15).
HeA(t−τ ) B dτ b dt + {noise}.
in which the matrix Φ0 (t) is (14)
∫t
t0
dyp (t) = dy(t) −
∫t
HeA(t−τ ) Cu(τ ) dτ dt.
(15)
(22)
t0
and the noise term is introduced similarly as in (16). In this case the augmented system for z(t) = [ xT (t), θT (t) ]T ∈ Rn+(p+1)m is given by
t0
dz(t) = A0 z(t) dt + Bθ (t)z(t)uS (t − t0 ) dt
Then, (13) can be rewritten, denoting b = b(t), as dyp (t) = Hp (t)b(t) dt + Rp dvp (t) (t ≥ t0 ),
HeA(t−τ ) BΦ(τ ) dτ
Φ0 (t) =
Here, note that the L.H.S. is a known quantity, and hence let it be denoted as dyp (t), i.e.,
(16)
+ C0 u(t) dt + G0 dw0 (t),
(23)
where
where
[
∫t
Bθ (t) =
HeA(t−τ ) B dτ
Hp (t) =
= HA−1 [ eA(t−t0 ) − I ]B in which the second equality holds if A is stable, and the term Rp dvp (t) is introduced with a newly introduced m-vector standard Wiener process vp (t) to describe the {noise} term. At this stage, it should be emphasized that (16) can be regarded just as an observation model for the state b(t) (Ohsumi, 2008). In this sense, the newly defined dyp (t) is called the pseudomeasurement for unknown parameter vector b. Clearly, {dyp (t)} is Yt -measurable (Yt = σ{y(τ ), 0 ≤ τ ≤ t}). In order to adapt such expression as (16) for all t > 0, we assume formally for 0 < t < t0 that dyp (t) = Rp dvp (t)
(0 < t < t0 )
(17)
in which dyp (t) ≡ 0 is set. For the case of impulsive input modeled by (3), noting that ∫t e
]
and other matrices are the same as that in the case of (2).
t0
A(t−τ )
0 BΦ(t) 0 0
b δ(τ − t0 )dτ = e
A(t−τ )
¯ ¯ b ¯¯
0
uS (t − t0 ), (18)
4. IDENTIFICATION OF UNKNOWN INPUT Augmenting the pseudomeasurements (16) and (17) with the original observation process to form the vector y0 (t) = [ y T (t), ypT (t) ]T , we obtain the augmented observation process for the case of (2): dy0 (t) = H0 (t)z(t) dt + R0 dv0 (t) (24) where [ H0 (t) =
H 0 0 χHp (t)
]
in which χ is an indicator taking 0 for 0 < t < t0 and 1 for t ≥ t0 ; v0 (t) = [ v T (t), vpT (t) ]T ; and R0 = block diag {R, Rp }. Then, we have the system formulation: dz(t) = A0 z(t) dt + B0 z(t)uS (t − t0 ) dt + C0 u(t) dt + G0 dw0 (t) dy0 (t) = H0 (t)z(t) dt + R0 dv0 (t).
(25)
τ =t0
Based on (25), the Kalman filter can be constructed to obtain the estimate of b (= b(t)):
we have x(t) = eA(t−t0 ) x(t0 ) + eA(t−t0 ) b
dˆ z (t|t) = Aχ zˆ(t|t) dt + C0 u(t) dt
∫t eA(t−τ ) Cu(τ ) dτ + {noise}
+
(t ≥ t0 ). (19)
P˙ (t|t) = Aχ P (t|t) + P (t|t)ATχ + G0 Q0 GT0
t0
By the similar way as in the case of (2), we have the version dyp (t) = Hδ (t)b(t) dt + Rp dvp (t)
+ P (t|t)H0T (t){R0 R0T }−1 {dy0 (t) − H0 (t)ˆ z (t|t) dt} (26)
(t ≥ t0 ),
(20)
where
298
− P (t|t)H0T (t){R0 R0T }−1 H0 (t)P (t|t), (27) where zˆ(t|t) = E{z(t)|Yt } and P (t|t) = E{[z(t) − zˆ(t|t)][z(t) − zˆ(t|t)]T } are the optimal estimate of z(t)
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 and the estimation error covariance matrix, respectively; Aχ = A0 for 0 < t < t0 and = A1 ≡ A0 + B0 for t ≥ t0 , and Q0 = block diag {Q, Qb }. For the case of (3) the Kalman filter can be also constructed mutatis mutandis. For the case of (4), the system formulation is given as (25) by replacing the matrices B0 and H0 (t) by Bθ (t) and Hθ given below, respectively; [ Hθ (t) =
] H 0 . 0 χΦ0 (t)
The Kalman filter is given by (26) and (27) by replacing H0 (t) in them by Hθ (t). In this case A0 (t) = A0 for 0 < t < t0 and A1 (t) = A0 + Bθ (t) for t ≥ t0 .
5. DETECTION OF ONSET TIME Another important task is to detect the onset time t0 at which an exogeneous input is discharged. To do this, revive t0 to be unknown, and suppose a system which is subjected to no exogeneous input and the same initial condition, dxf (t) = Axf (t) dt + Cu(t) dt + G dw(t), xf (0) = x0 . (28) Note that xf (t) ≡ x(t) until t = t0 because the same initial condition as in (1) is set. After t0 , xf (t) becomes no longer equal to the true x(t). Recalling this fact in mind, let us consider the following innovation process dν(t) = dy(t) − H x ˆf (t|t) dt,
ν(0) = 0
(29)
in which y(t) is the original process in (1) and x ˆf (t|t) is the estimation of xf (t) obtained by the Kalman filter with an appropriate initial estimate x ˆ0 : dˆ xf (t|t) = Aˆ xf (t|t) dt + Cu(t) dt + Pf (t|t)H T {RRT }−1 dν(t)
(30)
P˙f (t|t) = APf (t|t) + Pf (t|t)AT (t) + GQGT T −1
− Pf (t|t)H {RR } T
HPf (t|t).
dx(t) = Ax(t) dt + b {uS (t − t0 ) − uS (t − t1 )} dt + g dw(t), (33) where w(t) is a scalar noise process with variance q = 10: [ ] [ ] 0 1 1 A= , g= ; −1 −2.9 1 and the true value of b is set as b = [3, 2]T . For the observation process (2), H = R = I was set. The user-defined parameters were set as: Gb = I, Qb = 2I for the 2-vector noise process wb (t), Rp = I for the pseudomeasurement process (16) and (17). The initial values for the Kalman filter (26) and (27) were set as zˆ(0|0) = [ 0, 0, 0, 0 ]T and P (0|0) = diag {3, 4, 1, 2}. The simulation studies were performed with discretized version with time partition δt = 0.01 s. The true times t0 and t1 were set as t0 = 20, t1 = 29 s, respectively. To detect t0 and t1 the onset-time detector (32) was used. For the computation of the innovation process (29) the Kalman filter given by (30) and (31) was run with initial values x ˆf (0|0) = [ 0, 0 ]T and Pf (0|0) = diag {3, 4}. The user-defined parameter T0 in the detector (32) was chosen as T0 = 5 δt s. Figure 1 depicts the identification processes for the unknown parameters {bi }. In the figure, solid curves depict the estimation or identification, while the dotted curve the true value. Figure 2 depicts the behavior of the onset-time detector d(t). As seen in the figure, d(t) decreases rapidly after t = 20 s. If we set the threshold to be dthreshold = 0.5, then the estimate of t0 is decided as tˆ0 = 21.55, while for tˆ1 = 31.55 s for the same threshold value. These are sufficiently close to each true values. From Fig. 1 we see that the identification and state estimation is performed well.
(31)
(32)
6.2 Identification of quite-unknown exogenous input
∫t k ν(τ ) k dτ,
Consider the following second-order system in which the input u(t) is assumed to be zero, i.e., u(t) = 0 for all t:
Figure 3 shows the result of identification using the Kalman filter constructed based on the models (9) and (10) in which no idea of pseudomeasurement is used. As suspected, in this case the identification of unknown parameters is quite unsuccessful.
To detect t0 , define a function, d(t) =
6.1 Identification of stepwise exogenous input
t−T0
where k · k is the Euclidean norm and T0 is a short time interval. The onset instant t0 is detected as the time when d(t) exceeds a preassigned threshold. Denote it by tˆ0 and replace t0 in (26) and (27) by it.
The following third-order system is considered in which the input u(t) is assumed to be zero: dx(t) = Ax(t) dt + Bξ(t) dt + G dw(t),
where the variance of the w(t)-process is set as Q = diag {90, 120, 105}; [
6. SIMULATION STUDIES
A=
Owing to the space limit, only the cases of stepwise and quite unknown exogenous input are investigated to grasp the feeling how the use of pseudomeasurement is useful.
299
(34)
] −5 5 10 0 −10 5 , 5 0 −15
[ B=
and onset time t0 was set t0 = 0.
] 1 0 0 1 , 2 0
[ G=
] 60 0 03 0 ; 0 0 4.5
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009
Fig. 1. Identification features of b1 (top) and b2 (bottom).
Fig. 4. Identification of ξ1 (top) and ξ2 (bottom) whose true processes are given by (35). 3 1 1 1 5 1 6 ξ2 (t) = 3 + t − t2 − t3 + t4 + t − t . (35) 8 9 32 120 960 For the observation process, y(t) ∈ R2 , the matrices were set as [ H=
Fig. 2. Behavior of the detector d(t).
] 10 1 , 02 1
[ R=
] 2 0 . 0 1
For identifying the unknown inputs {ξ(t)}, it is assumed that ξi (t) =
p ∑
bik φk (t)
(i = 1, 2).
(36)
k=0
As for the basis functions {φk (t)}, the following simple polynominals of order p = 3 are selected: φk (t) = (t − tj )k for tj ≤ t < tj+1 (k = 0, 1, 2, 3), (37) where {tj } are equi-partitioned times such that t0 = 0 and tj+1 = tj + 1 s (i = 0, 1, 2, · · ·). The polynomial (37) has been modified from the conventionally used form, φk (t) = tk . This modified form is used here to guarantee the condition |φk (t)| < 1 for all t to prevent the divergence as time goes by. The unknown vector θ is given by Fig. 3. Unsuccessful results of identification b1 (top) and b2 (bottom) by the conventional method based on (9) and (10). The true exogenous input ξ(t) ∈ R2 was assumed as follows:
θ = [b10 , b20 , b11 , b21 , b12 , b22 , b13 , b23 ]T . For the model of θ(t), the matrices concerned with the wb (t)-process were assumed as user-defined parameters, Gb = diag {0.4, 0.3, 0.1, 0.2, 0.05, 0.06, 0.003, 0.004}, Qb = diag {8, 6, 2, 4, 1, 1.2, 0.06, 0.08}.
1 1 1 1 1 6 3 t ξ1 (t) = 3 + t − t2 − t3 + t4 + t5 − 2 2 4 24 80 720
Furthermore, as for the pseudomeasurement process (21), the matrix Rp was set as Rp = diag {3, 2}.
300
15th IFAC SYSID (SYSID 2009) Saint-Malo, France, July 6-8, 2009 (ii) The pseudomeasurement approach can be used for deterministic systems. For the deterministic systems, observer theory has been developed for the system with unknown input by the name of disturbance observer or unknown input observer (e.g., Medich and Hostetter, 1974). In this theory the disturbance input is modeled by a constant, otherwise the polynomial model (4) with φi (t) = ti . However, the coefficient parameters {bi } are not directly identified, while in the pseudomeasurement approach the unknown input can be recovered surely. Though the adaptive observers can be constructed for identifying unknown input, they have sophisticated structures even in the case of single-input single-output (L¨ uders and Narendra, 1974; Kreisselmeier, 1977). As opposed to adaptive observers, the pseudomeasurement approach is rather simple to construct an observer or Kalman filter to identify the unknown input. Furthermore, impulsive inputs are not identified by the disturbance observer theory. Fig. 5. Identification of ξ1 (top) and ξ2 (bottom) given by (38). The initial value for the Kalman filter (26) and (27) were set as zˆ(0|0) = 0 and P (0|0) = diag {3, 4, 5, 1, 2, 1, 2, 1, 2, 1, 2}. Figure 4 shows the result of identification of the unknown inputs ξ1 (t) and ξ2 (t). The solid curves and dotted ones depict the identification and true values, respectively. Another simulation example is shown in Fig. 5. The system and observation are the same as the above example except the matrix ] [ −50 50 100 0 −100 50 . A= 50 0 −150 The true exogenous inputs are assumed as 1 1 1 1 ξ1 (t) = sin t + cos t, ξ2 (t) = sin t + cos t. (38) 30 40 40 30 These are identified by the Fourier series expansion (35) with p = 16 and k+1 cos t (k : odd) 2 φ0 (t) = 1, φk (t) = (39) sin k t (k : even). 2 As seen in Fig. 5 the identification of the quite-unknown inputs seems to be performed fairly well. 7. CONCLUDING REMARKS The problem to identify unknown exogenous input to the dynamical system is important and attracting much attention to many control engineers. (i) In general, the estimate of unknown input requires higher derivatives of the observation data as inverse system. However, such higher derivatives of observation data is quite unrealistic, particularly for the stochastic systems. The approach proposed in this paper which will be called the pseudomeasurement approach does not need any derivatives of observation data.
301
(iii) As well-known, the fault diagnosis is attracting much interest of process control engineers as one of recent topics. In the process control, the fault signal is regarded as unknown exogenous input and the adaptive observer is used to identify it. The pseudomeasurement approach will have potential ability to identify even such fault signals.
REFERENCES A. T. Alouani and W. D. Bair. Use of a Kinematic Constraint in Tracking Constant Speed, Maneuvering Target. Proc. 30th Conf. on Decision and Control, Brighton, England, 2055–2058, 1991. T. Kimura, A. Ohsumi, and M. Kono. Observer-based Identification of Unknown Exogenous Input via Pseudomeasurement Approach. Proc. SICE Annual Conf. 2008, Chofu, Tokyo, Aug., 2660–2665, 2008. G. Kreisselmeier. Adaptive Observers with Exponential Rate of Convergence. IEEE Trans. on Automatic Control, vol. AC-22, no. 1, 2–8, 1977. G. L¨ uders and K. S. Narendra. A New Canonical Form for an Adaptive Observer. IEEE Trans. on Automatic Control, vol. AC-19, no. 2, 117–119, 1974. J. S. Medich and G. H. Hostetter. Observers for Systems with Unknown and Inaccessible Inputs. Int. J. Control, vol. 19, no. 3, 473–480, 1974. A. Ohsumi. An Outlook on the Use of Pseudomeasurement in System Identification. Proc. 37th SICE Symp. on Control Theory, Kirishima, Sep., 91–96, 2008. A. Ohsumi and S. Yasuki. Tracking of a Maneuvering Target Considering its Kinematic Constraints. Intelligent Automation and Control—Recent Trends in Development and Applications, Vol. 9 (eds: M. Jamshidi, P.Borne, and J. S. Jamshidi), TSI Press, Albuquerque, NM, 407–416, 2000. A. Ohsumi, S. Komiyama, K. Kashiwagi, M. Watanabe, and T. Takatsu. Detection of Polluted Load and Identification of its Discharged Location and Magnitsude for Polluted River. Int. J. of Innovative Computing, Information and Control, vol. 4, no. 1, 63–77, 2008. M. Tahk and J. L. Speyer. Target Tracking Problems Subject to Kinematic Constraints. IEEE Trans. on Automatic Control, vol. 35, no. 3 (1990), 324–326, 1990.