Robust State Estimation and Parameter Estimation for Linear and

Subscriber access provided by UNIV OF LOUISIANA

Process Systems Engineering

Robust State Estimation and Parameter Estimation for Linear and Nonlinear Direct Feed-through Systems with Correlated Disturbances Devyani Varshney, Mani Bhushan, and Sachin C. Patwardhan Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.9b00284 • Publication Date (Web): 08 May 2019 Downloaded from http://pubs.acs.org on May 9, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Robust State Estimation and Parameter Estimation for Linear and Nonlinear Direct Feed-through Systems with Correlated Disturbances Devyani Varshney, Mani Bhushan,∗ and Sachin C. Patwardhan Department of Chemical Engineering, Indian Institute of Technology Bombay E-mail: *[email protected]

Abstract In this work, we propose a generalized parameterized optimal filter (GPOF) for linear stochastic direct feed-through systems with correlated process and measurement disturbances. We consider a generic case where any set of unknown parameters can simultaneously affect the stochastic process model and measurements. GPOF allows unbiased estimates of states to be obtained in presence of unknown inputs or parameters by appropriately choosing the gain matrix during the state update step. We additionally prove the stability and convergence properties of proposed GPOF for linear time invariant systems with correlated process and measurement noises. We also propose a parameter estimation approach which utilizes the biased innovations generated by GPOF. We then extend the results of GPOF to nonlinear stochastic direct feed-through systems with non-additive correlated process and measurement noises. In particular we introduce extended generalized parameterized optimal filter (EGPOF) ∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

for state and parameter estimation. The proposed method is applied to two simulation based case studies. The results indicate that the proposed method can be a useful tool for state and parameter estimation of nonlinear systems with non-additive correlated noises.

1

Introduction

State estimation is the methodology of optimally estimating all the states of a system by using limited noisy measurements in combination with uncertain model dynamics. This classical problem has received considerable attention in the past. 1–4 Parameters of a dynamic mechanistic model of a process system are known with varying degree of accuracy. Some parameters/unmeasured disturbances are either poorly known or can drift slowly and change significantly over a period of time. From the viewpoint of process monitoring and control, it is important to develop state estimation schemes that are robust to such parameter inaccuracies and drifts in unmeasured disturbances (henceforth referred to as unknown parameters) and this area has also received attention in literature. 5,6 In many model based control and fault diagnosis applications, it is also important to estimate the unknown parameters. It is to be noted that the unknown parameters can affect not only the state dynamics but also affect the state to output map. Such systems are called direct feed-through systems. For example, faults such as measurement biases or drifts in the sensor-transmitter model parameters are frequently observed in process plant operation. 7 Biased or faulty sensors can lead to performance deterioration of process monitoring and control schemes. Thus, it is important to develop state estimation schemes which are robust to unknown parameters for direct feed-through systems. The approaches for state and parameter estimation available in literature can be broadly classified under two categories: simultaneous estimation of the states and unknown parameters, and decoupled estimation of the states and unknown parameters. A widely used approach in the simultaneous estimation category is to assume a dynamic model (typi2


Page 2 of 55

Page 3 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cally a random walk model) for the variation of the unknown parameter, and then perform joint state and parameter estimation by augmenting the unknown parameters as additional states. 8 While conceptually straightforward, these approaches have a major drawback of requiring a model for the parameter variation which if incorrect can lead to even biased state estimates. 9 The decoupled approaches involve estimating the states in a robust manner without requiring knowledge of the values of the unknown parameters. This is typically achieved by choosing an appropriate gain satisfying some constraints for the update step to ensure that the estimated states are unbiased despite the unknown parameters. Most of these robust estimation approaches have been applied for linear systems, such as the work of Kitanidis 6 which considered systems without direct feed-through i.e. the unknown inputs did not directly affect the measurements. Subsequently, Darouach et al. 5 extended the work of Kitanidis 6 by parameterizing the set of solutions of the constraints which result from the unbiasedness requirements. The result was further simplified by Keller et al. 10 by using another parameterizing solution. However, the aim of these approaches 5,6,10 was robust state estimation and parameter estimation was not considered. Also, majority of these approaches focus on unknown disturbances in the state dynamics and parameters directly affecting the state to output map (i.e. direct feed-through term) is not considered. For systems with direct feed-through, Palanthandalam-Madapusi and Bernstein 11 extended the work of Kitanidis 6 by requiring satisfaction of additional constraints while deriving the gain to ensure unbiased estimates. Darouach et al. 12 also presented a parameterized filter for a linear system with direct feed-through where they parameterized the set of solutions of the unbiasedness constraints and also considered different filter forms. The main focus of these works 11,12 was also on state estimation in presence of unknown inputs. However, Palanthandalam-Madapusi and Bernstein 11 have additionally considered the estimation of unknown parameters in a separate step by using only the current innovations. Further, there are other approaches in literature for direct feed-through systems which have extended the idea of decoupled 3



estimation to estimate parameters along with states, often in multi-step process without making assumptions about the dynamic variation of the parameters. 13–15 However, these filters reduce to the decoupled state estimation approaches when state estimates are made invariant to the parameter estimates. 14 Most of the existing literature in decoupled approaches is for linear systems, with uncorrelated process disturbances and measurement noises. Further, only a few works have considered nonlinear systems. These include the works of Palanthandalam-Madapusi et al. 16 and Ganesh et al. 17 who have proposed extensions of Kitanidis Kalman filter (KKF) 6 for nonlinear systems with additive noises for state and parameter estimation. In PalanthandalamMadapusi et al., 16 an unscented KKF was proposed for state estimation. The parameters, which appeared linearly in the process model, were estimated in a separate step using only the current innovations. Recently, Varshney et al. 18 extended the work of Ganesh et al. 17 for nonlinear systems and provided statistical characterization of the innovations in a window which was subsequently used for parameter estimation in a separate step. All these works 16–18 consider unknown inputs only in the process model and do not consider direct feed-through. Teixeira et al. 19 have proposed a gain constrained unscented Kalman filter for nonlinear direct feed-through systems. However, the work does not explicitly discuss the construction of the matrices involved in the constraints for nonlinear systems. Recently, Song 20 has proposed an unscented filter for simultaneous input-state estimation for nonlinear systems with non-additive correlated disturbances where identical unknown input affects the process and measurements. However, in presence of both unknown parameters and inputs in model, they augment unknown parameters with states for input-state-parameter estimation problem which increases the dimensionality of the problem. Further, unknown inputs are assumed to appear linearly in the measurement model. The decoupled state estimation approaches present an attractive option for state estimation since they ensure unbiased estimates without requiring a model for the parameter variation. In this work we focus on development of decoupled approaches for linear and 4


Page 4 of 55

Page 5 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nonlinear dynamic models with uncertain and drifting parameters. In particular, we present a generalization of the work of Darouach et al. 12 for state estimation for systems where inputs could affect both the process model and measurements but which additionally considers correlated disturbances and measurement noises. The stability results of Darouach et al. 12 are also generalized to establish the stability of the generic filter thus developed for a linear time invariant system. We additionally present a moving window parameter estimation strategy based on maximizing the likelihood of the innovations. The parameter estimation approach preserves the decoupled and robust nature of state estimation. We subsequently extend the proposed decoupled state and parameter estimation approaches to nonlinear systems. In particular, we extend the generalized linear filter to nonlinear systems with non-additively occurring process disturbances which could also be correlated with measurement noises. The proposed extension is generic in nature and considers nonlinearly occurring process disturbances which may be correlated with measurement noises, and nonlinearly occurring unknown inputs both in process and measurement model. While the proposed extended generic filter enables state estimation, the parameters are estimated in a decoupled step using the likelihood function of the innovations resulting from the nonlinear process and measurement model. The applicability of our decoupled state and parameter estimation approaches is demonstrated on a linear example and nonlinear case studies. Rest of the paper is organized as follows: Section 2.1 describes the proposed generalized parameterized optimal filter (GPOF) for state estimation of linear systems; Section 2.2 presents the stability and convergence properties of proposed GPOF for linear time invariant systems while Section 2.4 discusses parameter estimation approach for the proposed GPOF. The proposed GPOF is extended to nonlinear systems (EGPOF) for state and parameter estimation in Section 3. Section 4 presents two case studies to demonstrate the applicability of the proposed estimator and Section 5 presents some discussions to give additional insights into the proposed approach. Finally, concluding remarks are presented in Section 6.

5



2

Page 6 of 55

Generalized Parameterized Optimal Filter

In this section we propose a generalized Parameterized Optimal filter (GPOF) that incorporates correlated disturbances and measurement noise and which is applicable for the case when any set of unknown inputs could simultaneously affect the stochastic process model and measurements. In particular, we extend the results of Darouach et al. 12 (referred to as Parameterized Optimal filter (POF)) to systems with correlated disturbances. Further, we demonstrate the stability of the developed filter for linear time invariant systems. In later part of the section, we propose a parameter estimation approach using innovations in a moving window framework. Using the stability property of proposed GPOF, we show that the innovations obtained are multivariate Gaussian and use this result to develop a maximum likelihood function approach for estimating the unknown parameters. The applicability of the proposed filter is demonstrated on a simple linear system adapted from Yong et al. 15 Consider the following linear direct feed-through system: ¯ ¯ X (k)θ X (k) + Γ(k)w(k) X(k + 1) = Φ(k)X(k) +G ¯

(1a)

Y(k) = C(k)X(k) + GY (k)θ Y (k) + v(k)       ¯ ¯  0  Q S  w(k) with   .  ∼ N   ,  ST R 0 v(k)

(1b)

In Eq. (1a), X(k) ∈ Rn are the states of the system, θ X (k) ∈ RpX denotes the vector of unknown model parameters while w(k) ¯ ∈ Rm represents process disturbance or process noise ¯ In Eq. and is a zero mean Gaussian white stochastic sequence with covariance matrix Q. (1b), Y(k) ∈ Rr represents the measurement of linear functions of the states, θ Y (k) ∈ RpY denotes the vector of unknown parameters affecting the measurements while v(k) ∈ Rr represents measurement noise which is assumed to be a Gaussian white stochastic sequence with zero mean and covariance matrix R. Further, the stochastic processes w(k) ¯ and v(k) T are correlated as E w(k)[v(k)] ¯ = S. Without any loss of generality, it is assumed that 6


Page 7 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


θ X (k) and θ Y (k) can be represented as, θ X (k) = θ¯X + ∆θ X (k) and θ Y (k) = θ¯Y + ∆θ Y (k) where θ¯X and θ¯Y represent the known nominal values of the parameters θ X (k) and θ Y (k) respectively while ∆θ X (k) and ∆θ Y (k) represent the unknown time varying perturbations in θ X (k) and θ Y (k) respectively. Further, without any loss of generality, it is assumed that some common unknown parameters appear in both state and measurement equations (Eq. (1a) and (1b)) simultaneously i.e. θ xy , θ X ∩ θ Y = 6 {∅}. Then, θ X can be partitioned T where θ x ∈ Rpx denotes the exclusive part of θ X in θ X ∪ θ Y and as [θ x ]T [θ xy ]T T pxy where θ y ∈ Rpy denotes the θ xy ∈ R . Similarly, θ Y can be written as [θ y ]T [θ xy ]T exclusive part of θ Y in θ X ∪ θ Y . Thus, pX = px + pxy and pY = py + pxy . We accordingly ¯ X (k) and GY (k) with respect to θ x , θ xy and θ y as follows: partition the matrices G ¯ X (k) = G

Gx (k) Gxy (k)

, GY (k) =

Gy (k) Gyx (k)

(2)

It is to be noted that Gxy (k) 6= Gyx (k) in Eq. (2). Further, the following assumptions have been made for the system under consideration: Assumption 1: Total number of unknown parameters in Eqs. (1a) and (1b) do not exceed the number of states or number of measurements, i.e. px + py + pxy ≤ min(n, r). Before we propose GPOF, we implement a common trick 2 to transform the system into an equivalent system without correlated disturbances. This is achieved by adding a zero term (D(k) [Y(k) − C(k)X(k) − GY (k)θ Y (k) − v(k)]) to Eq. (1a) as follows: ¯ ¯ X (k)θ X (k) + Γ(k)w(k) X(k + 1) = Φ(k)X(k) +G ¯ + D(k) [Y(k) − C(k)X(k) − GY (k)θ Y (k) − v(k)] = Φ(k)X(k) + D(k)Y(k) + GX (k)θ xy (k) + Gx (k)θ x (k) − D(k)Gy (k)θ y (k) + w(k)

(3)

¯ where, Φ(k) = Φ(k) − D(k)C(k), GX (k) = Gxy (k) − D(k)Gyx (k), w(k) = Γ(k)w(k) ¯ − D(k)v(k).

Now D(k) is chosen such that Eq. (1a) gets converted into an equivalent dynamical system 7



Page 8 of 55

with no correlation between process and measurement disturbances. Therefore with choice of D(k) = Γ(k)S[R]−1 , the disturbances in the transformed process model become uncorrelated with measurement noise i.e. E[w(k)[v(k)]T ] = 0. Thus, the reformulated state dynamic equation is given as:

X(k + 1) = Φ(k)X(k) + D(k)Y(k) + GX (k)θ xy (k) + Gx (k)θ x (k) − D(k)Gy (k)θ y (k) + w(k) (4) Further, in Eq. (4), D(k)Y(k) is treated as a known input and w(k) ∼ N (0, Q(k)) as ¯ − S[R]−1 [S]T [Γ(k)]T . process noise with Q(k) = Γ(k) Q

2.1

State Estimation using GPOF

In this section, we develop a generic version of POF for state estimation of direct feed-through systems with correlated disturbances. In particular, it is desired to obtain an unbiased estimate of X(k + 1) using the system of equations given by Eqs. ((4) and (1b)) with ˆ measurements available up to (k + 1)th time instant. Let X(k|k) represent the estimated states at time instant k based on measurements available till time instant k. Further, it is ˆ assumed that the corresponding state estimation error, (k|k) , X(k) − X(k|k) is a zero mean random variable with covariance P(k|k). We now find conditional expectation of X(k + 1) using measurements up to time instant k using Eq. (4) to obtain the prediction equation as follows: ˆ + 1|k) = Φ(k)X(k|k) ˆ X(k + D(k)Y(k) + GX (k)θ¯xy + Gx (k)θ¯x − D(k)Gy (k)θ¯y

(5)

Note that only the nominal values of the parameters are used in prediction equation (Eq.

8


Page 9 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5). The associated prediction error and its covariance can thus be given as, (6a)

ˆ + 1|k) (k + 1|k) , X(k + 1) − X(k

= Φ(k)(k|k) + GX (k)∆θ xy (k) + Gx (k)∆θ x (k) − D(k)Gy (k)∆θ y (k) + w(k) (6b) ˆ where (k|k) , X(k) − X(k|k), P(k + 1|k) , E ((k + 1|k) − E[(k + 1|k)] ) [ ((k + 1|k) − E[(k + 1|k)] ) ]T = Φ(k)P(k|k)Φ(k)T + Q(k).

(6c)

Note that (k + 1|k) is not unbiased i.e. E[(k + 1|k)] 6= 0 even though (k|k) and w(k) are zero mean. This is due to the presence of the non-zero unknown quantity ∆θ x (k), ∆θ xy (k) and ∆θ y (k) on the right hand side of Eq. (6b). Now, we consider the filter of the form, ˆ + 1|k + 1) = X(k ˆ + 1|k) + L(k + 1)e(k + 1) X(k ˆ + 1|k) − GY (k + 1)θ¯Y where, e(k + 1) = Y(k + 1) − C(k + 1)X(k

(7a) (7b)

is the innovation at time instant (k + 1). Using Eqs. ((1b) and (6a)), the innovation at time instant (k + 1) can be written as,

e(k + 1) = C(k + 1)(k + 1|k) + GY (k + 1)∆θ Y (k + 1) + v(k + 1).

(8)

Further, using Eqs. ((7a) and (8)) estimation error at time instant (k + 1) can be expressed as, ˆ + 1|k + 1) (k + 1|k + 1) = X(k + 1) − X(k ˆ + 1|k) − L(k + 1)e(k + 1) = X(k + 1) − X(k = I − L(k + 1)C(k + 1) (k + 1|k) − L(k + 1)GY (k + 1)∆θ Y (k + 1) − L(k + 1)v(k + 1). (9)

9



Page 10 of 55

The corresponding updated covariance is given by: P(k + 1|k + 1) , E ((k + 1|k + 1) − E[(k + 1|k + 1)] ) [ ((k + 1|k + 1) − E[(k + 1|k + 1)] ) ]T T = I − L(k + 1)C(k + 1) P(k + 1|k) I − L(k + 1)C(k + 1) + L(k + 1)R[L(k + 1)]T . (10) We choose the gain matrix L(k +1) to ensure that the state estimates in GPOF are unbiased. For an unbiased estimator we require that E (k + 1|k + 1) = 0. Using Eqs. ((6b) and (9)) and noting that E[w(k)] = 0, E[v(k + 1)] = 0 and E[(k|k)] = 0, it follows: E[(k + 1|k + 1)] = I − L(k + 1)C(k + 1) GX (k)∆θ xy (k) + Gx (k)∆θ x (k) − D(k)Gy (k)∆θ y (k) − L(k + 1)GY (k + 1)∆θ Y (k + 1).

(11)

For E[(k + 1|k + 1)] = 0 to hold regardless of any nonzero ∆θ X (k) and ∆θ Y (k) then requires

(12a)

I − L(k + 1)C(k + 1) GX (k) = 0, I − L(k + 1)C(k + 1) Gx (k) = 0, I − L(k + 1)C(k + 1) D(k)Gy (k) = 0,

(12b)

L(k + 1)GY (k + 1) = 0.

(12d)

(12c)

Eq. (12) together form a system of linear equations which can be written in matrix form as: e e L(k + 1)H(k) = G(k) (13) e where, H(k) = C(k + 1)GX (k) C(k + 1)Gx (k) C(k + 1)D(k)Gy (k) −GY (k + 1) , e and G(k) = GX (k) Gx (k) D(k)Gy (k) [0]n×pY Similar to the filter developed by Darouach and Zasadzinski, 12 we parameterize the solution 10


Page 11 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of Eq. (13) by an arbitrary matrix which can subsequently be specified to obtain a minimumvariance estimator. Further, Eq. (13) has a solution if

e = rank (GX (k)) + rank(Gx (k)) + rank(D(k)Gy (k)) + rank(GY (k + 1)) rank(H)

(14)

Under these conditions, the solution is given by

e L(k + 1) = G(k)Π(k) + K(k + 1)T(k)

(15)

+ e e where Π(k) = (H(k)) ∈ R(pX +pY +py )×r represents any generalized inverse of H(k). Further, + e e T(k) = α(k)(I − H(k) × (H(k)) ) ∈ R(r−(pX +pY +py ))×r , where α(k) is an arbitrary matrix T is non(user specified) which must be chosen such that the matrix [T(k)]T [Π(k)]T

singular which is required to prove stability of GPOF in Section 2.2. Assumption 1 is necessary for this condition to hold. In Eq. (15), K(k + 1) ∈ Rn×(r−(pX +pY +py ) needs to be determined so as to minimize the trace of the error covariance matrix. Substituting Eq. (15) into Eq. (10) gives: ¯ ¯ P(k + 1|k + 1) = (F(k) − K(k + 1)β(k + 1))P(k|k)(F(k) − K(k + 1)β(k + 1))T T ¯ ¯ + Q1 (k) + K(k + 1)Θ(k)[K(k + 1)]T − K(k + 1)[S(k)] − S(k)[K(k + 1)]T T e e where, Q1 (k) = Σ(k)Q(k)[Σ(k)]T + G(k)Π(k)R[Π(k)] [G(k)]T

(16a) (16b)

e ¯ Σ(k) = I − G(k)Π(k)C(k + 1), F(k) = Σ(k)Φ(k), β(k + 1) = T(k)C(k + 1)Φ(k) (16c) T e ¯ + 1) = Σ(k)Q(k)[C(k + 1)]T [T(k)]T − G(k)Π(k)R[T(k)] S(k

Θ(k) = T(k)(C(k + 1)Q(k)[C(k + 1)]T + R)[T(k)]T

(16d) (16e)

The minimum variance estimation problem is reduced to finding the matrix K(k + 1) which minimizes the trace of the covariance matrix Eq. (16a). Therefore, taking derivative of tr(P(k + 1|k + 1)) (Eq. (16a)) with respect to K(k + 1) and equating to zero yields the

11



Page 12 of 55

following result, ¯ ¯ + 1)) β(k + 1)P(k|k)[β(k + 1)]T + Θ(k) −1 K(k + 1) = (F(k)P(k|k)[β(k + 1)]T + S(k

(17)

¯ + 1)[T(k)]T is non-singular where C(k ¯ + 1) = C(k + K(k + 1) will only exist if T(k)C(k 1)(Φ(k)P(k|k)[Φ(k)]T + Q)[C(k + 1)]T + R. Eq. (17) is used to obtain the modified gain which ensures unbiasedness even in presence of parametric plant-model mismatch. The state estimates thus obtained, are robust to the changing unknown parameters, even though the parameter changes are not updated in model equations. Remark 1. Filter given by Eq. (7a) is an extension of the filter described in Section 2.2.1 of Darouach et al. 12 for θ X = θ Y = θ and for uncorrelated process and measurement disturbances, i.e. S = 0. e ¯ + 1) is nonsingular, rank(H(k)) Remark 2. Assume that C(k = pX + pY + py = p˜ and e rank(G(k)) = pX + py , and let pseudo-inverse be given as + T ¯ −1 e e e e ¯ + 1)]−1 Π(k) = (H(k)) = ((H(k)) [C(k + 1)]−1 (H(k))) (H(k))T [C(k

(18)

e ¯ + 1)]−1/2 (H(k)) Further, let the singular value decomposition of [C(k be 

Σlk



  ˜ T e ¯ + 1)]−1/2 (H(k)) ˜ [C(k = U(k)  [V(k)]  0(r−˜p)×˜p

(19)

˜ ˜ where Σlk = diag(σ1 (k), . . . , σp˜(k)) with σ1 (k) ≥ σ2 (k) . . . ≥ σp˜(k) > 0 and U(k) and V(k) T ¯ ˜ are orthogonal matrices. Then for α(k) = 0(r−˜p)×˜p I(r−˜p)×(r−˜p) [U(k)] [C(k + 1)]−1/2 we obtain the filter given by Palanthandalam-Madapusi and Bernstein 11 for θ X = θ Y and S = 0. This can be shown along similar lines as given in Theorem 1 of Darouach et al. 5 which shows connection between their filter and Kitanidis filter 6 for linear systems with unknown parameters affecting only the stochastic process model. 12


Page 13 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


+ e is full row rank, then H(k) e e Remark 3. If matrix H × (H(k)) = Ir×r . In this case, we have

e T(k) = 0 and the gain reduces to L(k + 1) = G(k)Π(k). Here the choice of α(k) does not affect obtained gain. In the next section, we give conditions for convergence and stability properties of GPOF for linear time invariant systems.

2.2

Convergence and Stability of GPOF

In this section we extend the stability and convergence properties given in Darouach et al. 12 to our proposed GPOF for linear time invariant systems which are subsequently utilized in parameter estimation. In this section, we consider the special case where all known matrices are independent of time, and we will use the standard results on the convergence of the Riccati equation ¯ X (k), corresponding to the system that is affected by unknown input. Suppose that Φ(k), G GY (k) and C(k) in Eqs. ((1b) and (4)) are all independent of time instant k, i.e. Φ(k) = Φ, ¯ X (k) = G ¯ X , GY (k) = GY and C(k) = C, then the associated Riccati difference equation G (RDE) from Eq. (16a) can be written as: ¯ − K(k + 1)β)P(k|k)(F ¯ − K(k + 1)β)T + Q1 + K(k + 1)Θ[K(k + 1)]T P(k + 1|k + 1) = (F ¯ T − S[K(k ¯ − K(k + 1)S + 1)]T

(20)

¯s = F ¯ − SΘ ¯ −1 β, Qs = Q1 − SΘ ¯ −1 [S] ¯ T and Ks = K − SΘ ¯ −1 , the algebraic Further, defining F Riccati equation (ARE) can be given as: ¯ s − Ks β)P(F ¯ s − Ks β)T + Ks Θ[Ks ]T + Qs = P (F

13


(21)


Page 14 of 55

e defined in Eq. (13), Assumption 3: We assume that Θ > 0 and for (H) (22)

e = rank(GX ) + rank(Gx ) + rank(DGy ) + rank(GY ) rank(H) We then have the following lemmas: ¯ s , β) is detectable if, and only if, Lemma 1. Under Assumption 3, the pair (F     s s e zI − F −GX −Gx DGy [0]  zI − F −G rank  = rank    C [0] [0] [0] GY C Gy   e G e + rank(Gy ) = n + rank   = n + rank(G) Gy = n + rank(GX ) + rank(Gx ) + rank(DGy ) + rank(GY )

∀z ∈ C, |z| ≥ 1 (23)

e = G where, G X Gx −DGy [0] , Gy = [0] [0] [0] GY Proof. The proof is similar to that of Lemma 4 in Darouach et al. 12 and is hence omitted. ¯ s, Q ¯ s1/2 ) is Lemma 2. Under Assumption 3, the existence of an unreachable mode λ of (F equivalent to the existence of a row vector ω ¯ 6= 0 such that 



1/2 e 0  −Φ + λI G Q ω ¯  = 0. 1/2 λC Gy 0 R

(24)

Proof. The proof is similar to that of Lemma 4 in Darouach et al. 12 and is hence omitted. Then the following theorem can be stated to define stability and convergence properties of the proposed GPOF for linear time invariant systems. Theorem 1. The ARE (Eq. (21)) has a unique stabilizing solution P and subject to P(0), the sequence P(k), generated by the RDE (Eq. (20)), converges exponentially to P if, and 14


Page 15 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


only if,   s e zI − F − G  e + rank(Gy ) rank   = n + rank(G) C Gy

(25)

= n + rank(GX ) + rank(Gx ) + rank(DGy ) + rank(GY ) ∀z ∈ C, |z| ≥ 1

(26)

and   jω 1/2 e 0  −Φ + e I G Q rank  =n+r jω 1/2 Gy 0 R e C

∀ ω ∈ [0, 2π]

(27)

Proof. From the results of Lemma 1 and Lemma 2 and from the known results on the stability and convergence of the standard Kalman filter (Anderson and Moore 3 ), Theorem 2 gives the convergence and the stability of GPOF (Eq. (7a)) for linear time-invariant systems. In the next section we demonstrate the applicability of the proposed GPOF for a linear system adapted from Yong et al. 15

2.3

State Estimation: Illustrative Example

In this example, we consider the state estimation problem for a linear discrete time system adapted from Yong et al. 15 in which the stochastic process model as well as the measurements are affected by the unknown parameters.

15



Page 16 of 55

The various system matrices are:     0.5 2 0 0 0   0 1 0 0 0      0 0.2 1  0 1     0 0 1 0 0     ¯ = 0 Φ ; 0 0.3 0 1  ; C =  0 0 0 1 0         0 0.7 1  0.1 0   0 0 0 0 1 0 0.2 0 0 0.1    1 0 0 0 0     1 0 0 1 0.5 0 0    0 1    −4 −3   ¯ Q = 10 × 0 0.5 1 0 0 ; R = 10 ×  0 0       0 1 0 0 0   0.5 0 0 0 0 0 1 T Γ = I5×5 ; Gy = 0 0 0 0

    1 −0.3         1  0             Gx = 0 ; Gxy =  0  ;         0  0      0 0   1 1  0 0.5  1 1    0 0.3  −5  ; S = 10 ×  1 1  1 0    1 1  0 1 1 1

Gyx

  0   1   =  ; 0     0

 1 1   1 1   1 1    1 1  1 1 (28)

To test the efficacy of the proposed filter, we compare the proposed GPOF with Kalman filter (referred to as KF) which uses nominal values of the unknown parameters. The dynamic variations of the unknown parameters used in this example are as given in Fig. 1. Initial T condition given to the system is X(0) = 0 0 0 0 0 . The nominal values of the e is full row rank and hence unknown parameters are considered to be 0. For this example, H choice of α(k) is immaterial (see Remark 3). The estimator is initialized as,

ˆ X(0|0) = −0.2 −0.1 −0.01 0.01 0.1

T , P(0|0) = 10 × I5×5

(29)

Fig. 2 shows the tracking of the unmeasured state by GPOF. It can be seen that the proposed filter is able to track the true states without any bias whereas Kalman filter results in biased estimation.

16


Page 17 of 55

θx

2

0

-2 0

100

200

300

400

500

600

700

800

900

1000

0

100

200

300

400

500

600

700

800

900

1000

θxy

2

0

-2

Sampling Instant

Figure 1: Profile of unknown parameters used while simulation GPOF

40

KF

True

30 20 10

x(1)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0 -10 -20 -30 -40 0

100

200

300

400

500

600

700

800

900

1000

Sampling Instant

Figure 2: Tracking of unmeasured state, x(1), by GPOF and KF

2.4

Parameter Estimation Using GPOF

In this section we introduce a decoupled approach for estimation of unknown parameters for a system described in Eq. (1). As shown in Section 2.1, the generalized parameterized optimal filter is constructed such that the state estimates are insensitive to unknown variations in the parameters. As a consequence, the innovation sequence contains information about parameter variations which can be exploited to estimate parameters without affecting the state estimates obtained by the GPOF. Taking motivation from Varshney et al., 18 we estimate the parameters in a moving window framework using innovations generated by GPOF. Our proposed parameter estimation approach for direct feed-through systems with 17



Page 18 of 55

correlated disturbances involves the following key concepts: 1. Moving window framework: We extend the moving window based approach given in Varshney et al. 18 to estimate unknown parameters by stacking innovations in a window of length N . 2. Maximum Likelihood estimation of parameters: • Characterization of innovations in a window: We derive density function for characterizing the innovations obtained for direct feed-through systems with correlated disturbances in a window of length N . • Analytical Expression for the estimated parameter : Using the characterized innovations, we construct a Likelihood function for the unknown parameter vector which is then used to estimate the unknown parameters. The moving window approach is discussed in Section 2.4.1 while the maximum likelihood estimation idea is described in Section 2.4.2. Towards this step, the innovations in a window are characterized using Lemma 3 and Assumption 3 in Section 2.4.2. Finally, analytical expression for the parameter estimates based on maximizing the likelihood function are also derived in Section 2.4.2. 2.4.1

Moving Window Approach

In this section we formally introduce the moving window approach in which we stack the innovations in a window of length N . From Eqs. ((6b) and (8)), we can obtain the innovations from GPOF as follows:

e(k + 1) = C(k + 1)GX (k)∆θ xy (k) + C(k + 1)(Gx (k)∆θ x (k) − D(k)Gy (k)∆θ y (k)) + GY (k + 1)∆θ Y (k + 1) + η(k + 1) and η(k + 1) = C(k + 1) Φ(k)(k|k) + w(k) + v(k + 1) 18


(30) (31)

Page 19 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where η(k + 1) is a zero mean random variable with covariance V(k + 1) = C(k + 1)P(k + 1|k)[C(k + 1)]T + R. However, innovations (e(k + 1)) have non-zero mean due to presence of unknown ∆θ X (k) and ∆θ Y (k). We thus propose to use the non-zero mean innovations to estimate the unknown ∆θ X (k) and ∆θ Y (k). The parameters are estimated using the observed innovations in a time window of length N . To proceed further, we introduce the following simplifying assumption: Assumption 2: Variations in the parameters occur at a slow rate and thus parameters are assumed to remain constant over a window [k − N + 2, k + 1] i.e. ∆θ Z (j) = ∆θ¯Z[k−N +2, k+1] ∀j ∈ {k − N + 2, k − N + 3, . . . , k + 1}, T T T T where ∆θ Z (j) , [∆θ x (j)] [∆θ xy (j)] [∆θ y (j)]

(32)

where ∆θ¯Z[k−N +2, k+1] is the constant value over the window [k − N + 2, k + 1]. Under Assumption 2 i.e. ∆θ¯X (k) = ∆θ¯X (k + 1) and ∆θ¯Y (k) = ∆θ¯Y (k + 1), Eq. (30) reduces to:

e(k + 1) = (C(k + 1)GX (k) + Gyx (k + 1))∆θ xy (k) + C(k + 1)Gx (k)∆θ x (k) + (Gy (k + 1) − C(k + 1)D(k)Gy (k))∆θ y (k) + η(k + 1) where η(k + 1) is defined in Eq. (31).

19


(33)


The innovations (Eq. (33)) over a window length N are now stacked together to yield: E (k + 1) =A(k + 1)∆θ¯Z[k−N +2, k+1] + H(k + 1) (34)     e(k − N + 2) η(k − N + 2)      , H(k + 1) =  , where, E (k + 1) =  . . . . . .         e(k + 1) η(k + 1) A(k + 1) = A11 (k + 1) A12 (k + 1) A13 (k + 1) , (35)   C(k − N + 2)Gx (k − N + 1)   ..  A11 (k + 1) =  .     C(k + 1)Gx (k)   C(k − N + 2)GX (k − N + 1) + Gyx (k − N + 2)   ..  A12 (k + 1) =  .     C(k + 1)GX (k) + Gyx (k + 1)   Gy (k − N + 2) − C(k − N + 2)D(k − N + 1)Gy (k − N + 1)   ..  (36) A13 (k + 1) =  .     Gy (k + 1) − C(k + 1)D(k)Gy (k) 2.4.2

Maximum Likelihood framework for Parameter Estimation using GPOF

In this section we show that the innovations obtained from GPOF can be characterized as multivariate Gaussian for linear time invariant systems. This is presented as the following Lemma. Lemma 3. The innovation sequence vector E (k + 1) (Eq. (34)) obtained for linear time invariant system in a window [k − N + 2, k + 1] is multivariate Gaussian. Proof. The detailed proof is given in Section S1 of supporting information document and is similar to the proof in Varshney et al. 18 which showed that innovations obtained from Kitanidis filter are multivariate Gaussian for linear time invariant systems. Here we present 20


Page 20 of 55

Page 21 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


only a broad outline of the proof as follows: Firstly, we show that the estimation error obtained for a linear time invariant system from the proposed filter is Gaussian and zero mean. The Gaussianity of the same is shown by utilizing the stability and convergence properties of GPOF mentioned in Theorem 2. This in turn ensures that H(k + 1) is a Gaussian random variable with zero mean (Eqs 31 and 35). Based on this result, we finally show E (k + 1) (Eq. (34)) to be multivariate Gaussian for linear time invariant systems. However, current work involves linear time varying systems. Thus, we make the following assumption to enable us to characterize innovations obtained from GPOF: Assumption 3: The innovation sequence obtained using GPOF can be approximated as multivariate Gaussian. Using Eq. (34) and Assumption 3, we now formulate a likelihood function for the unknown parameter vector ∆θ¯Z[k−N +2, k+1] as, E (k + 1)|∆θ¯Z[k−N +2, k+1] ). L(∆θ¯Z[k−N +2, k+1] ) = f (E

(37)

∆θ¯Z[k−N +2, k+1] can then be estimated by maximizing log of likelihood function (Eq. (37)). Further, estimate of the unknown parameters can be determined under the following assumption: Assumption 4: A(k + 1) is full column rank at all time instants. Note that Assumption 1 is needed for this assumption to hold. The resulting estimate is 21 : d ∆ θ¯Z[k−N +2, k+1] = [A(k + 1)T [W[k−N +2, k+1] ]−1 A(k + 1)]−1 A(k + 1)T [W[k−N +2, k+1] ]−1E (k + 1). | {z } V(k+1)

(38) W[k−N +2, k+1] in Eq. (38) is covariance matrix of innovations in the window and can be 21



Page 22 of 55

computed as given by Eq. (S2.12) in Section II of supporting information. Further, in Eq. (38), existence of V(k + 1) requires A(k + 1) to be full column rank which is ensured by d Assumption 4. The estimated vector, ∆ θ¯Z[k−N +2, k+1] can be approximated as a Gaussian random vector 21 with: d Mean: E[∆ θ¯Z[k−N +2, k+1] ] = ∆θ¯Z[k−N +2, k+1] , d Covariance: Σ[k−N +2, k+1] , cov(∆ θ¯Z[k−N +2, k+1] ) = [A(k + 1)T [W[k−N +2, k+1] ]−1 A(k + 1)]−1 . (39) From properties of multivariate Gaussian random variables 1 , it follows that the marginal d probability density function of the ith component ∆ θ¯Z[k−N +2, k+1],i is univariate Gaussian given d as ∆ θ¯Z[k−N +2, k+1],i ∼ N (∆θ¯Z[k−N +2, k+1],i , Σ[k−N +2, k+1]i,i ). Thus, a 100(1 − κ)% confidence interval on the true value of ith parameter ∆θ¯Z[k−N +2, k+1],i can be obtained as: d ∆ θ¯Z[k−N +2, k+1],i −zκ/2

q q d Σ[k−N +2, k+1]i,i < ∆θ¯Z[k−N +2, k+1],i < ∆ θ¯Z[k−N +2, k+1],i +zκ/2 Σ[k−N +2, k+1]i,i (40)

where κ is the user specified significance level and zκ/2 is z-score corresponding to right-tailed area of κ/2 of standard normal distribution. The confidence intervals provide likely ranges for the unknown parameters. We next demonstrate the applicability of the proposed parameter estimation approach. 2.4.3

Parameter estimation: Illustrative Example

We implement the parameter estimation technique as discussed in Section 2.4 on the linear system example presented in Section 2.3. We choose window length parameter N = 5 for the same. Fig. 3 shows that the proposed filter is able to track the true parameter variation. It should be noted that the state estimates are still same as given in Fig. 2 since the parameter estimation approach is decoupled from the state estimation approach, and the estimated states are not affected by the estimated parameters. 22


Page 23 of 55

True

GPOF

θx

2 0 -2 0

100

200

300

400

500

600

700

800

900

1000

0

100

200

300

400

500

600

700

800

900

1000

2

θxy

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0 -2

Sampling Instant

Figure 3: Tracking of unknown parameters by GPOF based parameter estimation approach

2.5

Online State and Parameter estimation using GPOF

The procedure of the proposed GPOF for online state and parameter estimation for one time instant (say k th instant) is summarized below: 1. Choose window length (N ) for estimating the unknown parameters. ˆ 2. Given: Initial mean X(k|k), covariance P(k|k) and nominal parameter values θ¯x , θ¯y , θ¯xy . ˆ + 1|k)) using Eq. (5) and the predicted covariance 3. Compute the predicted states (X(k using Eq. (6c). ˆ + 1|k + 1)) 4. On availability of new measurement Y(k + 1), obtain updated states (X(k using Eqs. ((7a) and (15)) and compute corresponding covariance (P(k + 1|k + 1)) using Eq. (10). 5. Compute covariance [W [k−N +2, k+1] ] as mentioned in Section II of supporting information. d 6. Estimate the unknown parameters (∆ θ¯Z[k−N +2, k+1] ) using Eq. (38). 7. Go to Step (2) with k = k + 1. 23



Page 24 of 55

The above online state and parameter estimation procedure was used in the illustrative examples discussed in Sections 2.3 and 2.4.3. Note: For the sake of completeness, a detailed comparative evaluation of GPOF and Kalman filter on the linear system example presented in this section for both state and parameter estimation is presented in Section S7 of supporting information. The procedure and metrics for this evaluation are similar to those presented in Section 4 for nonlinear systems. The comparative results are also qualitatively similar to those presented in Section 4.

3

Extension of GPOF to Nonlinear Systems

GPOF proposed for state and parameter estimation in Section 2.1 is applicable to linear systems. However, most processes are nonlinear in nature. In this section, we further extend the proposed GPOF to nonlinear systems with non-additive correlated disturbances. We additionally extend the parameter estimation methodology to nonlinear systems. Consider a discrete time direct feed-through system consisting of nonlinear process dynamics and nonlinear measurement function which are represented as,

X(k + 1) = F(X(k), u(k), θ X (k), w(k)) ¯ Y(k) = h(X(k), θ Y (k)) + v(k)       ¯ w(k) 0  Q S  with   ∼ N   ,   . T v(k) 0 S R

(41a) (41b)

In Eq. (41a), u(k) ∈ Rm represents the known inputs at time instant k. Rest of the symbol representation is the same as given in Eq. (1). Function F : Rn × Rm × RpX × Rm → Rn in Eq. (41a) represents the nonlinear state dynamics while function h : Rn × RpY → Rr in Eq. (41b) represents the nonlinear observation model. In next section, taking motivation from Varshney et al. 18 on extended version of Kitanidis Kalman filter (KKF) 6 and proposed GPOF in Section 2.1, we propose a generic extended 24


Page 25 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


version of GPOF (referred as EGPOF later) for state estimation of nonlinear, direct feedthrough systems with non-additive correlated disturbances. In particular, it is desired to obtain an unbiased estimate of X(k + 1) using the system of equations given by Eqs. ((41a) and (41b)) with measurements available up to (k + 1)th time instant. Further, as discussed in Section 2.1 earlier, we transform the system represented by Eqs. ((41a) and (41b)) into an equivalent system without correlated disturbances. Towards this end, we employ a common trick (discussed in Section 2.1 earlier) of adding a zero term D(k)[Y(k) − h(X(k), θ Y (k)) − v(k)] to Eq. (41a) as follows:

X(k + 1) = F(X(k), u(k), θ X (k), w(k)) ¯ + D(k)[Y(k) − h(X(k), θ Y (k)) − v(k)]

(42)

T In Eq. (42), matrix D(k) is chosen such that E[w(k)[v(k)] ¯ ] = 0. Choice of matrix D(k) is

discussed further in the following section.

3.1

State Estimation using EGPOF

ˆ Let X(k|k) represent the estimated states at time instant k based on measurements available till time instant k. Further, it is assumed that the corresponding state estimation error, ˆ (k|k) , X(k) − X(k|k) is a zero mean random variable with covariance P(k|k). Using ˆ Taylor series expansion in the neighborhood of X(k|k), θ¯X , θ¯Y and E [w(k)] ¯ = 0 , Eq. (42) can be expressed as follows: ˆ ¯ ˆ ¯ X,U (k)∆θ X (k) + Γ(k)w(k) X(k + 1) ≈ F(X(k|k), u(k), θ¯X , 0) + Φ(k)(X(k) − X(k|k)) +G ¯ ˆ ˆ + D(k)[Y(k) − h(X(k|k), θ¯Y ) − CU (k)(X(k) − X(k|k)) − GY,U (k)∆θ Y (k) − v(k)] (43)

25



¯ where, Φ(k) =

∂F ∂X

, Γ(k) = ˆ (X(k|k),u(k), θ¯X ,0)

∂F ∂w ¯

Page 26 of 55

and CU (k) =

ˆ (X(k|k),u(k), θ¯X ,0)

∂h ∂X

¯ X,U (k) and GY,U are obtained as follows: The system matrices G ¯ X,U (k) = G

=

Gx,U (k) Gxy,U (k)

∂F ∂θ x

∂F ∂θ xy

ˆ (X(k|k),u(k), θ¯X ,0)

. ˆ (X(k|k), θ¯Y )

,

ˆ (X(k|k),u(k), θ¯X ,0)

(44a) GY,U (k) =

Gy,U (k) Gyx,U (k)

=

∂h ∂θ y

ˆ (X(k|k), θ¯Y )

∂h ∂θ xy

.

ˆ (X(k|k), θ¯Y )

(44b)

Note that linearization is carried out at the nominal values of the unknown parameters and hence EGPOF does not require current estimate of θ X (k) or θ Y (k) for state estimation. Now by choosing D(k) = Γ(k)S[R]−1 , the measurement and process noises in this equivalent dynamic representation (Eqs. (43) and (41b)) become uncorrelated. Thus, the reformulated state dynamic equation is given as: ˆ ˆ ˆ X(k + 1) ≈ F(X(k|k), u(k), θ¯X , 0) − D(k)h(X(k|k), θ¯Y ) + Φ(k)(X(k) − X(k|k)) +GX,U (k)∆θ xy (k) + Gx,U (k)∆θ x (k) − D(k)Gy,U (k)∆θ y (k) + w(k) + D(k)Y(k) (45) ¯ where, D(k) = Γ(k)S[R]−1 , Φ(k) = (Φ(k) − D(k)CU (k)), GX,U (k) = (Gxy,U (k) − D(k)Gyx,U (k)), w(k) = (Γ(k)w(k) ¯ − D(k)v(k)) Further, in Eq. (45) D(k)Y(k) is treated as known input and w(k) ∼ N (0, Q(k)) as process ¯ − S[R]−1 [S]T [Γ(k)]T . Similar to Eq. (5), we find predicted noise with Q(k) = Γ(k) Q estimate of X(k + 1) as follows: ˆ + 1|k) = F(X(k|k), ˆ ˆ X(k u(k), θ¯X , 0) − D(k)h(X(k|k), θ¯Y ) + D(k)Y(k)

26


(46)

Page 27 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Using Eq. (6a), the associated prediction error is given by:

(k + 1|k) ≈ Φ(k)(k|k) + GX,U (k)∆θ xy (k) + Gx,U (k)∆θ x (k) − D(k)Gy,U (k)∆θ y (k) + w(k) (47) Corresponding prediction covariance can be computed by Eq. (6c). Now we consider the same filter form given in Eq. (7a) as follows: ˆ + 1|k + 1) = X(k ˆ + 1|k) + L(k + 1)e(k + 1) X(k

(48a)

ˆ + 1|k), θ¯Y ) where, e(k + 1) , Y(k + 1) − h(X(k

(48b)

Using Eqs. ((6a), (41b) and (48b)), innovation at time instant (k + 1) can be defined as:

e(k + 1) ≈ CP (k + 1)(k + 1|k) + GY,P (k + 1)∆θ Y (k + 1) + v(k + 1) ∂h ∂h where, CP (k + 1) = , GY,P (k + 1) = ∂X (X(k+1|k), ∂θ Y (X(k+1|k), ˆ ˆ θ¯Y ) θ¯Y )

(49)

Further, following a similar approach for filter development as given in Section 2.1, the constraints (Eq. (13)) is modified as follows: e e L(k + 1)H(k) = G(k) (50) e where, H(k) = CP (k + 1)GX,U (k) CP (k + 1)Gx,U (k) CP (k + 1)D(k)Gy,U (k) −GY,P (k + 1) , e and G(k) = GX,U (k) Gx,U (k) D(k)Gy,U (k) [0]n×pY Similar to the linear case, Eq. (50) has a solution if the following conditions hold:

e = rank (GX,U (k)) + rank(Gx,U (k)) + rank(D(k)Gy,U (k)) + rank(GY,P (k + 1)) rank(H)

27



Page 28 of 55

Under the above conditions, the associated error covariance matrix and gain L(k + 1) can be easily obtained by Eqs. ((15), (16) and (17)). Note that during gain computation in EGPOF, ˆ the linearized matrices at (k)th instant are computed using X(k|k) and are denoted with ˆ + 1|k) subscript U , whereas linearized matrices at (k + 1)th instant are computed using X(k and are denoted with subscript P . In next section we now extend the parameter estimation approach discussed earlier in Section 2.4 to EGPOF for direct feed-through systems with non-additive correlated disturbances.

3.2

Parameter Estimation using EGPOF

The moving window approach presented in Section 2.4.1 is now extended to EGPOF for parameter estimation. Since the system is nonlinear, using Eqs. ((47) and (49)) and Assumption 2, we can approximate the innovations as follows:

e(k + 1) ≈ (CP (k + 1)GX,U (k) + Gyx,P (k + 1))∆θ xy (k) + CP (k + 1)Gx,U (k)∆θ x (k) + (Gy,P (k + 1) − CP (k + 1)D(k)Gy,U (k))∆θ y (k) + η(k + 1) and η(k + 1) = CP (k + 1) Φ(k)(k|k) + w(k) + v(k + 1)

(52) (53)

where the matrices Gx,U (k), Gy,U (k), Φ(k), GX,U (k), CP (k + 1) and GY,P (k + 1) are obtained using linearization as discussed in Eqs. ((44), (45) and (49)). Further, similar to procedure described in Section 2.4.1, we stack the innovations for window length N . The resulting matrices similar to Eq. (34) for nonlinear systems are given in Section S4 of supporting information. Now corresponding to Section 2.4.2, it is straightforward to estimate the unknown parameter in a moving window framework. The parameter estimation approach remains decoupled and the estimated parameters do not have a bearing on the estimated states.

28


Page 29 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.3

Online State and Parameter estimation using EGPOF

The procedure of the proposed EGPOF for online state and parameter estimation for one time instant (say k th instant) is summarized below: 1. Choose window length (N ) for estimating the unknown parameters. ˆ 2. Given: Initial mean X(k|k), covariance P(k|k) and nominal parameter values θ¯x , θ¯y , θ¯xy . ˆ + 1|k)) using Eq. (46) and corresponding covari3. Compute the predicted states (X(k ance using Eq. (6c) where linearized matrices to compute P(k + 1|k) are obtained as explained in Section 3.1. ˆ + 1|k + 1)) 4. On availability of new measurement Y(k + 1), obtain updated states (X(k using Eqs. ((48a) and (15)) and compute corresponding covariance (P(k + 1|k + 1)) using Eq. (10) (linearized matrices to compute P(k +1|k +1) are obtained as explained in Section 3.1). 5. Compute covariance [W [k−N +2, k+1] ] as mentioned in Section II of supporting information. d 6. Estimate the unknown parameters (∆ θ¯Z[k−N +2, k+1] ) using Eq. (38). 7. Compute the covariance associated with estimated parameters using Eq. (39). 8. Go to Step (2) with k = k + 1.

4

Case Studies

We now present two simulation based case studies corresponding to different scenarios of parameters appearing simultaneously in process model and measurements. The case studies

29



demonstrate the application of the discussed methodology (Section 3.3) for state and parameter estimation of a nonlinear system with correlated process and measurement noises. These are discussed as follows: 1. Case Study I : In this case study we consider the case where different set of unknown parameters affect the stochastic process model and measurements i.e. θ xy = θ X ∩θ Y = {∅}. We consider the Williams-Otto reactor 22 where reaction rate constants are unknown parameters in state dynamics and unknown bias in sensor affects a measurement. Further, we introduce the process disturbances non-additively in the plant dynamics through input flow with correlated process and measurement noises. 2. Case Study II: This case study considers a CSTR with two-phase exothermic reaction. 23 Here flow valve coefficient is the unknown parameter affecting both states and measurements and inlet temperature of coolant is the unknown input affecting only the process dynamics. Thus, this case study considers the case when some common unknown parameters appear in both state and measurement equations i.e. θ xy 6= {∅}. For case study I, we compare EGPOF performance with conventional EKF (labeled as EKF(c)) and simultaneous EKF (labeled as EKF(s)). In EKF(c), nominal values of the unknown parameters are used for state estimation and we do not perform any parameter estimation. This allows us to evaluate the state estimation performance of EGPOF over EKF in presence of unknown parameters. Further, for both the case studies, we compare EGPOF with simultaneous EKF which performs simultaneous state and parameter estimation. We model the dynamic variation of unknown parameters in EKF(s) as random walk model and the noise variance of the random walk model is tuned to ensure good performance of EKF(s). Detailed implementation procedure of EKF(s) for correlated process and measurement noises can be found in Section S4 of supporting information document. Comparison with a well tuned EKF(s) which involves simultaneous state and parameter estimation enables us to benchmark the EGPOF performance. All the computations were performed in Matlab ver30


Page 30 of 55

Page 31 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


sion R2015b running on Windows10 operating system on a computer with 16GB RAM and Intel-core i7 processor. To obtain meaningful comparisons between the filtering approaches, multiple (Nr ) simulation runs were performed with the length of each simulation run being Ns time instants. The simulation runs differed in the realizations of true state initial conditions and process and measurement noises obtained at different time instants. Three performance measures have been used in case studies for comparison: 1. ARMSE (Average Root Mean Square Error): v  u Nr Ns u X X 1 j 2 ˆ t 1 (X(k)j(i) − X(k|k) ARM SE(i) = (i) ) Nr j=1 Ns k=1

(54)

j th ˆ where ARM SE(i), X(k)j(i) and X(k|k) state the ARMSE, true (i) denote for the i

state and estimated state for j th simulation run, respectively. 2. ANEES 4 (Average Normalized Estimation Error Squared) : Nr h iT i h 1 X j j j −1 j j ¯ ˆ ˆ X(k) − X(k|k) P(k|k) X(k) − X(k|k) β(k) = Nr j=1

(55)

j ˆ where X(k)j , X(k|k) and P(k|k)j are true state, estimated state and associated co-

variance matrix for j th run at time instant k. Under the assumption that the state¯ estimation error is N (0, P(k|k)), Nr β(k) is a chi-squared random variable with Nr n degrees of freedom. Thus, the filter is said to be consistent at κ significance level if ¯ Nr β(k) ∈ [ζ1 , ζ2 ] for approximately 100(1 − κ) instances, where ζ1 and ζ2 are thresholds derived from a chi-squared density function with Nr n degrees of freedom at κ significance level.

31



Page 32 of 55

3. Confidence interval for parameter estimation error (Section 2.4.2):

−zκ/2

q q Σ[k−N +2, k+1]i,i < θ(i) < zκ/2 Σ[k−N +2, k+1]i,i

(56)

d where, θ(i) , ∆θ¯Z[k−N +2, k+1]i − ∆ θ¯Z[k−N +2, k+1]i and κ is the user specified significance level. j ˆ ARMSE is also computed for parameters where X(k)j(i) , X(k|k) (i) in Eq. (54) are replaced j j d by ∆θ¯Z[k−N +2, k+1]i and ∆ θ¯Z[k−N +2, k+1]i for ith parameter and j th run respectively.

Remark 4. For all the case studies the matrix α(k) (Eq. (15)) for EGPOF has been chosen T ¯ ˜ ˜ as 0(r−˜p)×˜p I(r−˜p)×(r−˜p) [U(k)] [C(k + 1)]−1/2 , where U(k) is obtained using Eq. (19). Note that this choice of α(k) has been adapted from Darouach et al. 12

4.1

CASE STUDY I: θ X ∩ θ Y = {∅} (Williams-Otto reactor Case Study)

This case study describes the scenario when the unknown parameters affecting the state dynamics and measurements are different. The Williams−Otto reactor is one unit of the Williams−Otto plant model. 22 The reactor is a CSTR with mass holdup W (kg) and temperature TR (K). The reactor is fed with two pure component reactant streams FA and FB . The following three simultaneous reactions involving 6 components occur inside the reactor:

k

k1 = β1 e−6666.7/TR s−1

(57a)

k

k2 = β2 e−8333.3/TR s−1

(57b)

k

k3 = β3 e−11111/TR s−1

(57c)

1 A+B − → C, 2 B+C − → P + E, 3 C +P − → G,

where β1 = β01 × 106 s−1 , β2 = β02 × 108 s−1 and β3 = β03 × 1012 s−1 respectively. 32


Page 33 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The following dynamic mass balances represent the plant behavior: dXA dt dXB W dt dXC W dt dXE W dt dXG W dt dXP W dt W

= FA − (FA + FB )XA − r1

(58a)

= FB − (FA + FB )XB − r1 − r2

(58b)

= −(FA + FB )XC + 2r1 − 2r2 − r3

(58c)

= −(FA + FB )XE + 2r2

(58d)

= −(FA + FB )XG + 1.5r3

(58e)

= −(FA + FB )XP + r2 − 0.5r3

(58f)

where XA , XB , XC , XE , XG , XP are the mass fractions of the respective components and r1 = k1 XA XB W, r2 = k2 XB XC W, r3 = k3 XC XP W are the three reaction rates. Steady flow for F¯A = 1.8270kg/s was considered for the system with W = 2104.7kg. To simulate the dynamics, the manipulated variables, FB and TR , were perturbed by introducing PRBS signals in the frequency range [0, (0.05π/Ts )] (Fig. 4). Further, input FB is corrupted ¯ = 0.0573 (kg/s)2 , which is with white Gaussian noise with zero mean and covariance Q equivalent to standard deviation of 5% steady flow for FB = 4.789kg/s. Further, parameter β03 is considered as a known disturbance and is perturbed as shown in Fig. 4. We have chosen β01 and β02 which are the pre-exponential factors in Arrhenius equations (57), as unknown time varying parameters in state dynamics. Additionally, a time varying bias (yb ) is introduced in the measurement of XB for the current system. Following nominal values of reaction rate frequency factors and bias are considered:

β01 = 1.6599 ,

β02 = 7.2117 ,

yb = 0

(59)

The true values of the parameters are varied as follows: at 300th sampling instant, a negative ramp disturbance is given simultaneously in the frequency factors β01 and β02 . The ramp disturbances are continued till the magnitude of the disturbance variables (β01 and β02 ) 33



Page 34 of 55

reaches 50% of the nominal values. These parameter profiles (Fig. 10, upper and middle panels) are similar to the profiles considered by Yip and Marlin. 24 A time varying bias yb (Fig. 10, lower panel) is added to the measurement of XB .

Figure 4: Input and parameter changes introduced to simulate Williams−Otto reactor for a typical simulation run To incorporate a realistic scenario where process disturbance would be correlated with measurement noise, the flowrate FA is varied as an ARMA (auto-regressive moving average) model as follows: FA = F¯A + δFA (k)

(60a)

where, δFA (k) = 0.999(δFA (k − 1)) + wF (k)

(60b)

and wF (k) ∼ N (0, 4 × 10−6 F¯A ) is the stochastic process noise in the flow rate. Now, δFA (k) is propagated at each time instant using Eq. (60b) (which is considered to be part of the process model) and then follows a zero order hold till the next time instant with a sampling interval of 30 seconds.

34


Page 35 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The measurements of the system are:

y1 (k) = XB (k) + v1 (k),

y2 (k) = XC (k) + v2 (k),

y3 (k) = XE (k) + v3 (k),

y4 (k) = XG (k) + v4 (k),

y5 (k) = XP (k) + v5 (k),

y6 (k) = F¯ + δF (k) + v6 (k)

(61)

where, vi (k) for i = 1, 2, ..., 6 are Gaussian random variables with zero mean and covariance R with R = diag [ (0.00389)2 , (0.00015)2 , (0.0029)2 , (0.0011)2 , (0.0011)2 , 10−4 × F¯A ].

(62)

The correlation between process disturbance and measurement noise is given by E[wF (k)[v6 (k)]T ] = −0.2 × 10−4 F¯A . Further, the initial condition for the plant is assumed to be corrupted with w which is assumed to be Gaussian random variable with zero mean and covariance Qinit (standard deviations are 0.5% of the steady state values of the corresponding states): Qinit = diag [ (0.00044)2 , (0.0019)2 , (0.00008)2 , (0.0014)2 , (0.00053)2 , (0.00054)2 , 0] (63) Thus, the initial condition for plant is given as: X(0) = [0.0874, 0.3896, 0.0153, 0.2907, 0.1075, 0.1095, 1.8270]T + w

(64)

Remark 5. For this case study, θ X ∩ θ Y = {∅} and hence during EGPOF implementation e e Gxy and Gyx have been taken as empty matrices. Therefore, the modified H(k) and G(k) used in EGPOF is as follows: e H(k) = CP (k + 1)Gx,U (k) CP (k + 1)D(k)Gy,U (k) −GY,P (k + 1) , e G(k) = Gx,U (k) D(k)Gy,U (k) [0]n×pY 35



4.1.1

Comparison with Conventional EKF

We first compare EGPOF with EKF(c) to test the efficacy of state estimates obtained by EGPOF in presence of unknown changing parameters. The unknown drifts in the parameters are not updated in the model for both the estimators and instead nominal values of unknown parameters (Eq. (59)) are used. The performances of EGPOF and EKF(c) are compared for Ns = 1900 samples for Nr = 100 simulation runs. The state estimators are initialized with the following initial condition: ˆ X(0|0) = [0.0874, 0.3896, 0.0153, 0.2907, 0.1075, 0.1095, 1.837]T , P(0|0) = 10−2 × I7×7 . (65) Fig. 5 shows the tracking of two filters for unmeasured state. It can be seen that even though EGPOF uses nominal values of the parameters, it is still able to track the true states. However, EKF(c) shows biased estimates due to presence of plant-model mismatch. Further, Table 1 presents the ARMSE values for the unmeasured state estimates. ARMSE obtained from EGPOF is significantly lower than ARMSE obtained from EKF(c) even in the presence of incorrect values of unknown parameters. EGPOF

EKF

True state

0.115 0.11 0.105 0.1 0.095

XA

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 55

0.09 0.085 0.08 0.075 0.07 0.065 0

1

2

3

Time (Sec)

4

5 ×104

Figure 5: Case Study I: Estimates of unmeasured state obtained from EKF(c) and EGPOF for a typical simulation run

36


Page 37 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: ARMSE values for unmeasured states and parameters (Case Study I) Estimator EKF(c) EKF(s) EGPOF 4.1.2

XA β01 β02 yb 0.0099 – – – 0.0011 0.0240 0.0615 0.0016 0.0016 0.0208 0.0493 0.0018

Comparison with Simultaneous EKF

In this section we compare EGPOF with EKF(s). EKF(s) requires assumption of a model associated with dynamic variation of the unknown parameters. In this work, we employ the random walk model in EKF(s) for capturing the parameter drift. The random walk model assumes that the change in parameters at each sampling instant is random in nature as: θ(k + 1) = θ(k) + wθ (k). The noise term (wθ (k)) in the random walk model is assumed to be 0 mean and Gaussian. Appropriate tuning of the covariance matrix (Qp ) of the noise term in the random walk model is critical for successful implementation of simultaneous EKF. 18 After several attempts, we obtained the best tuned matrix (Qp ) for simultaneous EKF approach to be: Qp = 10−2 × diag[ β01 , β02 , 0.04 ]

(66)

The rest of the implementation parameters for the case study are the same as mentioned in Section 4.1.1. Comparison of state estimation performance: Fig. 6 shows the tracking of the unmeasured state by the two filters. It can be seen that even though EGPOF uses nominal values of the parameters, it is still able to track the true state. Also, the state estimates obtained from EGPOF are similar to EKF(s) which uses the estimated values of parameters during state estimation. Further, Table 1 presents the ARMSE values for the unmeasured state. ARMSE of state XA obtained from EGPOF is similar to EKF(s) even though EGPOF uses only nominal value of the unknown parameters whereas EKF(s) uses the most recent estimated values of the parameters during state estimation. 37



EGPOF

0.115

EKF(s)

True state

0.11 0.105 0.1 0.095

XA

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 55

0.09 0.085 0.08 0.075 0.07 0

1

2

3

Time (Sec)

4

5 ×104

Figure 6: Case Study I: Estimates of unmeasured state obtained from EKF(s) and EGPOF for a typical simulation run To check consistency of filters, ANEES (Eq. (55)) is compared for EKF(s) and EGPOF (Fig. 7) at 5% significance level. Square root filtering was implemented to facilitate computation of ANEES values as direct inversion of covariance matrices was not always possible due to presence of significantly small values in the covariance matrices (refer to Remark 6 for details). In Fig. 7, the percentage violation of bounds by EKF(s) is 66% while EGPOF violates the bounds at 71% instants. This may be due to inconsistencies in the actual correlation of state estimation errors and the correlation estimated by the filter. To investigate this further, ANEES was again computed with covariance matrices containing only the diagonal entries of the original matrices. Fig. 8 shows the modified ANEES values. The bound violations by EGPOF have reduced to 21% while EKF(s) violates the bounds at 51% instants. This indicates that estimation error variances computed by EGPOF are more consistent with state estimation errors than EKF(s). Remark 6. The entries in matrix Γ(k) (Eq. (43)) obtained after linearization had small magnitudes. This led to small entries in error covariance matrices which caused numerical precision issues while computing ANEES values. To obtain higher numerical precision, square root filtering was implemented for both EGPOF and EKF(s) where the square roots of the predicted and updated covariance matrices were computed using Householder Trans38


Page 39 of 55

formation. Detailed implementation procedure for the same can be found in Section S5 in supporting information. 16 EKF(s) EGPOF

14 12

¯ β(k)

10 8 6 4 2 0

0

1

2

3

4

5 ×104

Time (Sec)

Figure 7: Case Study I: ANEES for state estimates 10

8

6

¯ β(k)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4 EKF(s) EGPOF

2

0

0

1

2

3

4

5

Time (Sec)

×104

Figure 8: Case Study I: Modified ANEES for state estimates An insight to the working of EGPOF can be obtained by viewing the innovations. Fig. 9 shows innovations obtained from EKF(s) and EGPOF. As the innovation sequences are noisy, they are filtered using MATLAB command f iltf ilt 25 with filter F (z) = 0.1/(1 − 0.9z −1 ) which subjects the innovation sequences to zero-phase reverse and forward digital filtering. It can be seen that the innovations obtained from EGPOF are biased whereas innovations from EKF(s) appear to be zero mean. As discussed in Section 2.4, the biased innovations are then further used to estimate the unknown parameters in the model. 39



EGPOF

Innovation (y1 )

0.06

EKF(s)

0.04 0.02 0 -0.02

0

Innovation (y2 )

1

2

3

4

5 ×104

×10−4

5 0 -5 -10 -15

0

1

2

3

4

5 ×104

Innovation (y3 )

5

×10−3

0 -5 -10

Innovation (y4 )

2

0

1

×10

2

3

4

5 ×104

−3

1 0 -1 -2

0

1

2

3

4

5 ×104

Innovation (y5 )

2

×10

−3

0 -2 -4

0

1

2

3

4

5 ×104

0.02

Innovation (y6 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 55

0.01 0 -0.01

0

1

2

3

Time (Sec)

4

5 ×104

Figure 9: Case Study I: Filtered innovations for a typical simulation run Comparison of parameter estimation performance: Estimates of frequency factors and measurement bias obtained from EGPOF and EKF(s) are shown in Fig. 10 with window length N = 10. It can be seen that the estimates obtained from both filters are able to track the true values of the three parameters. Also, Table 1

40


Page 41 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


shows that the reported ARMSE values of parameters for both the filters are similar. As parameter estimation can be performed only after first N sampling instants in EGPOF, the values of estimated parameters for the first 9 sampling instants are taken to be the corresponding nominal values (Eq. (59)). Thus, the ARMSE computations are based on estimates obtained after 10 sampling instants. Further, Fig. 11 shows the parameter estimation errors with 95% confidence intervals. First 9 sampling instants are removed from the figure as no parameter estimation is performed during that time. The red shaded region is the confidence interval computed by EGPOF and green shaded region is the confidence interval computed by EKF(s). It can be seen that the estimation errors of both the filters lie within the confidence bounds. However, confidence bounds of EKF(s) are larger signifying higher uncertainty in the parameter estimates.

4.2

CASE STUDY II: θ xy = θ X ∩ θ Y 6= {∅} (CSTR Case Study)

The CSTR case study has been used extensively for simultaneous state and parameter estimation. 23 The process involves an exothermic reaction A(l) → B(l) + C(g) . Volume of liquid phase and gaseous phase of the reactor is assumed to be constant. Following dynamic equations represents the plant behavior:

dCA dt dT dt dTc dt dn dt

Fi (CAi − CA ) − rA V Fi rA (−∆H) U A(T − Tc ) = (Ti − T ) + − V ρCp V ρCp Fc U A(T − Tc ) = (Tci − Tc ) + Vj Vj ρj Cpj

(67b)

= rA V − Fvg

(67d)

=

(67a)

(67c)

The pressure in the reactor depends on the number of moles n of vapor. This in turn depends on the rate of the reaction and vent (molar) flow rate Fvg . The vapor is assumed

41



EGPOF EKF(s) True

1.7 1.65 1.6 1.55

β01

1.5 1.45 1.4 1.35 1.3 1.25

0

1

2

3

4

5 ×104

Time (Sec) 7.4

EGPOF EKF(s) True

7.2 7

β02

6.8 6.6 6.4 6.2 6 5.8 5.6

0

1

2

3

4

5 ×104

Time (Sec) 0.045

EGPOF EKF(s) True

0.04 0.035 0.03 0.025

yb

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 55

0.02 0.015 0.01 0.005 0 0

1

2

3

4

5 ×104

Time (Sec)

Figure 10: Case Study I: Estimates of unknown parameters for a typical simulation run to behave ideally: P Vg = nRT . The reaction rate and vent flow rate is given as, rA = CA k0 e−E/RT , Fvg = Kv 42

p

P − Patm


(68)

Page 43 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 11: Estimation errors with 95% confidence bounds for a typical simulation run The measurement equations are:

y1 (k) = T (k) + v1 (k),

y2 (k) = Tc (k) + v2 (k)

y3 (k) = n(k)/Vg + v3 (k) s n(k)RT (k) y4 (k) = Kv − Patm + v4 (k) Vg 43 ACS Paragon Plus Environment

(69a) (69b) (69c)


where, vi (k) for i = 1, 2, 3, 4 are Gaussian random variables with zero mean and covariance R as shown in Table 1 of supporting information. State CA is considered to be unmeasured. To simulate the dynamics, the known input variable Ti in Eq. (67) was perturbed by introducing PRBS signals in the frequency range [0, (0.05π/Ts )] with amplitude of 15◦ R (Fig. 12). The inlet coolant temperature Tci (Fig. 16 (lower panel)) and valve coefficient Kv (Eq. (68)) are considered to be the unknown inputs in the system. Therefore, θ xy = Kv since Kv occurs in the process model (Eq. (67d)) and affects the measurements as well (Eq. (69c)). The nominal values for the unknown parameters are: Kv = 0.0152 lbmol f t lb−0.5 h−1 ,

(70)

Tci = 530 ◦ R

The initial condition for the plant is assumed to be corrupted with w which is assumed to be Gaussian random variable with zero mean and covariance matrix Qinitial . All filter implementation and model parameters are presented in Table S1 of supporting information. We compare EGPOF with EKF(s) approach. The tuned covariance (obtained after signifi550 540

T i(◦ R)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 55

530 520 510

0

5

10

15

Time (Hr)

Figure 12: Input changes introduced to simulate CSTR cant trial and error) associated with the random walk model of unknown parameter used in EKF(s) has been set to: Qp = diag 10−4 × Kv, 10−2 × Tci . Comparison of state estimation performance: Table 2 presents a comparison of ARMSE values of the state variables obtained using EGPOF and EKF(s). From this table it is seen that performance of EGPOF and EKF(s) are similar for CA while EKF(s) has superior performance for state n. This is due to the presence of unknown input (Kv ) in the model (Eq. (67d)) due to which EGPOF gives complete weightage 44


Page 45 of 55

to the noisy measurement of n/Vg (Eq. (69b)) while estimating the state. Fig. 13 shows the tracking of states by estimates obtained from EGPOF and EKF(s). Fig. 14 shows ANEES Table 2: Case Study II: ARMSE values for unmeasured states and parameters Estimator EGPOF EKF(s)

CA n Kv Tci 0.0232 0.6896 0.00086 6.7999 0.0224 0.3174 0.0010 6.2297

EGPOF

EKF(s)

True state

Ca (lbmol/f t3 )

0.45 0.4 0.35 0.3 0.25 0.2 0

5

0

5

10

15

10

15

36 34

n(lbmol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


32 30 28 26

Time (Hr)

Figure 13: Case Study II: Tracking of unmeasured states for a typical simulation run computed for EKF(s) and EGPOF at 5% significance level. Percentage violation of bounds by EGPOF is 11.67% while EKF(s) violates the bounds at 16.33% instants. This shows that computed covariances in EGPOF are more consistent with corresponding estimation errors than EKF(s). The innovations obtained from EGPOF and EKF(s) are subjected to zero-phase forward and reverse digital filtering as in Case Study I (Section 4.1.2) and are displayed in Fig. 15. As expected, innovations obtained from EGPOF are biased which are further used for parameter estimation. Comparison of parameter estimation performance: Fig. 16 shows the parameter and input estimates obtained by EGPOF and EKF(s) with window length N = 5 for EGPOF. It can be seen that both filters are able to track the true 45



14 EKF(s) EGPOF

12

10

¯ β(k)

8

6

4

2

0

0

5

10

15

Time (Hr)

Figure 14: Case Study II: ANEES for state estimates

Innovation (y3 )

EKF(s)

0.5 0 -0.5

Innovation (y2 )

Innovation (y1 )

EGPOF

0

5

10

15

0

5

10

15

0

5

10

15

0

5

10

15

2 0 -2

0.02

0

-0.02

Innovation (y4 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 55

0.4 0.2 0

Time (Hr)

Figure 15: Case Study II: Filtered innovations for a typical simulation run values. Table 2 shows the ARMSE values of the parameters for both filters. It can be seen that both filters report similar values. Also, Fig. 17 provides estimation error with 95%

46


Page 47 of 55

confidence bounds. It can be seen that the estimation errors for both filters lie within the confidence bounds. However, confidence bounds by EKF(s) are larger than EGPOF showing higher uncertainty in parameter estimates. Further, estimates for Tci show relatively high errors during step changes introduced in this parameter due to use of window based approach in EGPOF and random walk model in EKF(s) for parameter estimation. These errors in EGPOF can be reduced by tuning the window length N (not attempted here) used during parameter estimation. EGPOF

0.034

EKF(s)

True

0.032 0.03 0.028

Kv

0.026 0.024 0.022 0.02 0.018 0.016 0.014 0

5

10

15

Time (Hr) EPOF

EKF(s)

True

550 545 540 535

Tci

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


530 525 520 515 510 505

0

5

10

15

Time (Hr)

Figure 16: Case Study II: Estimates of unknown parameters for a typical simulation run

47



Figure 17: Case Study II: Estimation errors with confidence bounds for a typical simulation run

5

Discussions

The two case studies demonstrated superior features of the proposed EGPOF approach over EKF based approaches for state and parameter estimation of nonlinear systems with direct feed-through and non-additive correlated process disturbances and measurement noises. Following are the discussions based on these studies: 1. Advantage for state estimation: • Comparison of EGPOF with EKF(c) shows that the state estimates obtained from EGPOF are unbiased and insensitive to the unknown parameter variations 48


Page 48 of 55

Page 49 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


whereas EKF(c) gives biased estimates when there is parametric plant-model mismatch (Section 4.1.1). • EGPOF allows decoupling of the state and parameter estimation whereas EKF(s) performs simultaneous state and parameter estimation which leads to the use of parameter estimates for update of state estimates. It was observed (results not shown) that unless the covariance matrix for random walk model in EKF(s) was properly tuned, it could lead to incorrect parameter estimates which in turn leads to biased state estimates. This phenomenon has also been reported in literature. 18 For the case studies considered in our work, the covariance could be tuned since the true values of parameters were known. For a practical system, the true values will be unknown and thus tuning the covariance will be a non-trivial task. However, due to decoupling property of EGPOF, effect of incorrect parameter estimates on state estimates is not observed. 2. Advantage for parameter estimation: • Parameter estimation in EGPOF is decoupled from state estimation. This is not the case when a simultaneous EKF is used for joint state and parameter estimation. In EKF(s) the unknown parameters are augmented with the states and estimated. This is turn increases the dimension of the problem which is not the case in EGPOF. • The benefit of EGPOF for parameter estimation can be attributed to the use of window-based approach which does not require any prior knowledge of the past estimates of the parameters. Thus, erroneous estimates of parameters in the past do not affect the current estimate of parameters. • The specification of a model describing the time evolution of the parameters is a non-trivial task for parameter estimation using EKF(s). 18 Often the dynamics governing parameter variation are not known and hence a random walk model is 49



generally assumed in EKF(s). Such a model is not required in case of EGPOF. 3. Effect of modified gain used in EGPOF : EGPOF uses modified gain (Eq. (15)) to obtain state estimates which are robust to the unknown changing parameters. Use of the constraints can however lead to the selection of gains (L(k + 1)) such that the measurements get a significantly higher weightage than the predicted states while computing the estimated states (Eq. (48a)). This in turn can lead to noisy estimated states if the noise variances of the measurements are high. This was observed for state n in Case Study II (Section 4.2). 4. Estimation of Bias: Estimation of bias in a measurement in EGPOF requires the specification of GY matrix (Eq. 30). This in turn requires the user to a priori know the measurements which are corrupted by the bias. For a process with several measurements, this is a non trivial task. Integration of approaches from data reconciliation literature 26 for identifying a biased measurement with EGPOF can be investigated. 5. Nonlinearity and local Gaussian approximation: It should be noted that similar to EKF, EGPOF implicitly assumes the conditional probability densities of the states to be locally Gaussian. Hence, EGPOF is able to perform well even for highly nonlinear systems when local linearization and local Gaussianity assumptions hold. This can be confirmed from the case studies discussed in Section 4. Specifically, the CSTR case study in Section 4.2 has significant nonlinearity due to the presence of term involving the exponential of temperature (which is a state) in the state dynamics. Additionally, the parameters in the two case studies are perturbed significantly. In particular, the true value of flow valve coefficient (Kv ) in Case Study II varies upto −50% of its nominal value. The linearization in EGPOF is performed around the nominal value and not at the true value. In the two case studies, even in presence of nonlinearities in the system as well as the significant perturbations given to the true parameter values, EGPOF is able to accurately track the true states and the unknown parameters. 50


Page 50 of 55

Page 51 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6

Conclusions

In this work, we proposed a generalized parameterized optimal filter (GPOF) for linear direct feed-through systems with correlated disturbances. GPOF ensures that the estimated states are robust in the sense that they are not affected by unknown variations in the parameters. We extended the stability proof given by Darouach et al. 12 to our proposed filter. Using an approach similar to that developed in Varshney et al. 18 we also proposed a decoupled parameter estimation approach which utilizes the innovations generated by GPOF to estimate parameters in a moving window framework. We then extended GPOF to nonlinear systems, and proposed extended generalized parameterized optimal filter (EGPOF) for state and parameter estimation for nonlinear direct feed-through systems with non-additive correlated disturbances. The applicability of GPOF for linear system was demonstrated through an example, while the utility of the proposed EGPOF approach was established by applying it on two case studies. The results demonstrated that EGPOF with nominal values of unknown parameters performs much better than EKF(c) and is similar to a well tuned EKF(s) which performs simultaneous state and parameter estimation.

Acknowledgement We acknowledge financial support for this work from the Ministry of Human Resource Development (MHRD), Government of India.

Supporting Information Available The Supporting Information is available free of charge on the ACS Publications website at DOI: ———— Proof of Lemma 3; Steps for recursive computation of Covariance (W [k−N +2, k+1] ) for parameter estimation; Parameter estimation using EGPOF; Simultaneous EKF for correlated disturbances; Householder Transformation for Square root filtering; Implementation 51



parameters for CSTR study. This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Maybeck, P. S. Stochastic Models, Estimation and Control ; Academic Press, 1979; Vol. 1. (2) Gelb, A. Applied optimal estimation; MIT Press: Cambridge, MA, 1974. (3) Anderson, B.; Moore, J. Optimal Filtering; Prentice-Hall: Englewood Cliffs, NJ, 1979. (4) Crassidis, J. L.; Junkins, J. L. Optimal Estimation of Dynamic Systems, 2nd ed.; Chapman & Hall/CRC, 2011. (5) Darouach, M.; Zasadzinski, M. Unbiased Minimum Variance Estimation for Systems with Unknown Exogenous Inputs. Automatica 1997, 33, 717–719. (6) Kitanidis, P. K. Unbiased-Minimum variance Linear State Estimation. Automatica 1987, 23, 775–778. (7) Venkatasubramanian, V.; Rengaswamy, R.; Yin, K.; Kavuri, S. N. A review of process fault detection and diagnosis part I: Quantitative model-based methods. Computers & Chemical Engineering 2003, 27, 293–311. (8) Simon, D. Optimal State Estimation: Kalman, H∞ and Nonlinear Approaches; WileyInterscience: New York, NY, USA, 2006. (9) Valluru, J.; Lakhmani, P.; Patwardhan, S. C.; Biegler, L. T. Development of moving window state and parameter estimators under maximum likelihood and Bayesian frameworks. Journal of Process Control 2017, 60, 48–67. (10) Keller, J.; Darouach, M.; Caramelle, L. Kalman filter with unknown inputs and robust two-stage filter. International Journal of Systems Science 1998, 29, 41–47. 52


Page 52 of 55

Page 53 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(11) Palanthandalam-Madapusi, H. J.; Bernstein, D. S. Unbiased Minimum-variance Filtering for Input Reconstruction. Proceedings of the 2007 American Control Conference. New York City, USA, 2007. (12) Darouach, M.; Zasadzinski, M.; Boutayeb, M. Extension of minimum variance estimation for systems with unknown inputs. Automatica 2003, 39, 867–876. (13) Gillijns, S.; Moor, B. D. Unbiased minimum-variance input and state estimation for linear discrete-time systems with direct feedthrough. Automatica 2007b, 43, 934–937. (14) Hsieh, C.-S. Extension of unbiased minimum-variance input and state estimation for systems with unknown inputs. Automatica 2009, 45, 2149–2153. (15) Yong, S. Z.; Zhub, M.; Frazzoli, E. A unified filter for simultaneous input and state estimation of linear discrete-time stochastic systems. Automatica 2016, 63, 321–329. (16) Palanthandalam-Madapusi, H. J.; Girard, A.; Bernstein, D. S. Wind-field Reconstruction Using Flight Data. American Control Conference. Seattle, WA, USA, 2008; pp 1863–1868. (17) Ganesh, C.; Ballal, P.; Bhushan, M.; Patwardhan, S. C. 12th International Symposium on Process Systems Engineering and 25th European Symposium on Computer Aided Process Engineering; Computer Aided Chemical Engineering; Elsevier: Copenhagen, Denmark, 2015; Vol. 37; pp 1817 – 1822. (18) Varshney, D.; Bhushan, M.; Patwardhan, S. C. State and Parameter Estimation using Extended Kitanidis Kalman filter. Journal of Process Control, 2019, 76, 98-111. (19) Teixeira, B. O. S.; Chandrasekar, J.; Palanthandalam-Madapusi, H. J.; Torres, L. A. B.; Aguirre, L. A.; Bernstein, D. S. Gain-Constrained Kalman Filtering for Linear and Nonlinear Systems. IEEE Transactions on Signal Processing 2008, 56, 4113–4123.

53



(20) Song, W. Generalized minimum variance unbiased joint input-state estimation and its unscented scheme for dynamic systems with direct feedthrough. Mechanical Systems and Signal Processing 2018, 99, 886–920. (21) Sheldon, S. M. Introduction to Probability and Statistics for Engineers and Scientists, 5th ed.; Academic Press: Boston, 2014. (22) William, T. J.; Otto, R. E. A generalized chemical processing model for the investigation of computer control. Transactions of the American Institute of Electrical Engineers, Part I: Communication and Electronics 1960, 79, 458–473. (23) Vachhani, P.; Rengaswamy, R.; Gangwal, V.; Narasimhan, S. Recursive Estimation in Constrained Nonlinear Dynamical Systems. American Institue of Chemical Engineers Journal 2005, 51, 946–959. (24) Yip, W. S.; Marlin, T. E. Multiple data sets for model updating in real-time operations optimization. Computers & Chemical Engineering 2002, 26, 1345–1362. (25) The MathWorks, Inc., Signal Processing User’s Guide. https://in.mathworks.com/ help/pdf_doc/signal/signal_tb.pdf, 2018; [Online; accessed 14-February-2019]. (26) Narasimhan, S.; Jordache, C. Data Reconciliation & Gross Error Detection: An Intelligent Use of Process Data; Gulf Publishing Co.: Houston, TX, USA, 2000.

54


Page 54 of 55

Page 55 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

55


Robust State Estimation and Parameter Estimation for Linear and

Recommend Documents