Input Design for Online Fault Diagnosis of Nonlinear Systems with

Subscriber access provided by UNIVERSITY OF THE SUNSHINE COAST

Article

Input Design for Online Fault Diagnosis of Nonlinear Systems with Stochastic Uncertainty Joel Anthony Paulson, Marc Martin-Casas, and Ali Mesbah Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b00602 • Publication Date (Web): 31 Jul 2017 Downloaded from http://pubs.acs.org on August 1, 2017

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 free 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 accessible to all readers and 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.

Industrial & Engineering Chemistry Research 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 40

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

Input Design for Online Fault Diagnosis of Nonlinear Systems with Stochastic Uncertainty Joel A. Paulson, Marc Martin-Casas, and Ali Mesbah∗ Department of Chemical and Biomolecular Engineering, University of California, Berkeley E-mail: [email protected] Phone: +1 (510) 642-7998. Fax: +1 (510) 642-4778

Abstract Fault diagnosis is crucial for ensuring stable and reliable operation of high-performance systems in the presence of abnormal events. System uncertainties often render the discrimination between normal and faulty behavior a challenging task. This paper presents an active fault diagnosis (AFD) method for nonlinear systems with stochastic uncertainty. AFD involves the optimal design of system inputs for discriminating between multiple model hypotheses that correspond to various operational scenarios. The proposed AFD method relies on minimizing the probability of error in hypothesis selection subject to hard input and state chance constraints. Moment-based approximations for a bound on the probability of error in hypothesis selection as well as for chance constraint evaluation are introduced in order to derive a tractable surrogate AFD problem that is amenable to online applications. The performance of the AFD method for offline and online fault diagnosis is demonstrated on a continuous bioreactor with multiple operational scenarios. Simulation results demonstrate the effectiveness of the proposed method for online AFD in a practical setting.

1

ACS Paragon Plus Environment


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

Introduction Fault diagnosis allows for maintaining stable, reliable, and profitable operation of technical systems in the presence of component malfunctions, drifting system parameters, and abnormal events. 1 However, system uncertainties, disturbances, and measurement noise, which are ubiquitous in practical applications, often impede fault diagnosis of nonlinear systems. Fault diagnosis is further compounded by the growing complexity of high-performance systems and the increasingly stringent requirements on their operation, which necessitate retaining the system within admissible operational constraints during fault diagnosis. A variety of methods for fault diagnosis have been developed including residual- and observer-based methods, 2–4 set-based methods, 1 and data-based methods. 3 The vast majority of these methods are passive, that is, the status of the system is deduced based on measurements acquired during the nominal system operation and through comparisons with model predictions or historical data. However, system uncertainties and/or corrective actions of a controller, when the system is under feedback control, typically mask the effects of system faults to an extent that can impair reliable fault diagnosis. To address this challenge, auxiliary input signals can be applied to the system for enhancing the diagnosability of faults. 5 This notion has led to the development of so-called active fault diagnosis (AFD) methods in which an optimal input sequence is designed to diagnose faults over a prespecified time horizon while accounting for system uncertainties. 6 When system uncertainties are set-based (i.e., bounded uncertainty descriptions), AFD typically involves designing input sequences that are “robust” to worst-case realizations of the system uncertainty. These AFD methods often look to minimize the energy of the input sequence while guaranteeing separation of the reachable sets of all model hypotheses that correspond to the different fault scenarios, with the focus being on linear systems. 5,7,8 Zonotopes 9 and constrained zonotopes 10 have been adopted for designing active inputs that achieve guaranteed fault diagnosis for linear systems with set-based uncertainty. The latter method was extended to closed-loop input design for guaranteed fault diagnosis. 11 Set-based 2


Page 2 of 40

Page 3 of 40

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


AFD methods for nonlinear systems have also been developed. 12,13 AFD of systems with stochastic uncertainty has received relatively less attention than set-based AFD. In a pioneering work, Zhang 14 presented a probabilistic framework for AFD. The AFD problem for linear systems subject to additive stochastic process and measurement noise has been addressed by Blackmore and coworkers, 15,16 where the AFD criterion is formulated in terms of minimizing an upper bound on the probability of model selection error. This method was later extended to jump Markov linear systems. 17 The AFD problem for nonlinear systems subject to stochastic uncertainty has been addressed in a limited number of works. 18–21 In particular, Mesbah and coworkers 20 applied the generalized polynomial chaos 22 approach to develop a tractable surrogate AFD problem when the nonlinear system model includes probabilistic uncertainty in the parameters and initial conditions. This approach, however, cannot handle time-varying probabilistic uncertainties such as stochastic system and measurement noise. Results from randomized optimization have also been used to provide theoretical guarantees for fault detection in nonlinear systems. 21 This paper presents a probabilistic framework for the AFD of nonlinear systems subject to white-noise stochastic system and measurement noise as well as probabilistic model uncertainty in parameters and initial conditions that stems from our incomplete knowledge of the system. The probabilistic AFD problem is formulated as a chance-constrained optimization problem aimed at minimizing the probability of model selection error, which is quantified in terms of the Bayes risk. State constraints are enforced as a chance constraint to allow for flexibility and reduced conservatism in dealing with the stochastic system uncertainty. Since the probabilistic AFD problem is computationally intractable, the key contribution of this paper is to derive a tractable surrogate for the input design problem that is amenable to online computations in a closed-loop setting. Closed-loop AFD hinges on estimates of the model probabilities. Another contribution of this paper is the use of Bayes rule for developing an efficient recursive algorithm for online estimation of model probabilities. Closed-loop implementation of the AFD problem enables re-planning of the fault diagnosis experiment

3



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

at each measurement sampling time based on new system observations in order to account for unmodeled system uncertainties and disturbances. A major challenge in solving the probabilistic AFD problem arises from the fact that there exists no closed-form expression for evaluating the probability distribution of the predicted states variables as a function of the input sequence. To this end, we derive a moment-based surrogate for the AFD problem by adopting two different uncertainty propagation methods–linearization and unscented transformation, both of which widely used for designing Kalman filters for nonlinear systems based on the propagation of moments of state variables. 23,24 The knowledge of moments is used to obtain a tractable error bound for the Bayes risk 25 as well as a tractable surrogate for the state chance constraint. The paper is organized as follows. First, a general background on the Bayes risk for model selection is provided. The probabilistic AFD problem considered in this work is then presented, followed by the uncertainty propagation techniques. Subsequently, the proposed surrogate for the probabilistic AFD problem is presented, along with a discussion on the offline and online implementation of the AFD method. Lastly, the performance of the proposed AFD method is demonstrated on a continuous bioreactor with multiple fault scenarios including operational and structural changes to the process. ⊤ ⊤ ⊤ Notation. xa:b = [x⊤ denotes a concatenated vector of the sequence of a , xa+1 , . . . , xb ]

values of x from discrete time a to b. diag(·) represents a block diagonal matrix. P (A) denotes the probability of event A. For a random vector X, p(X) denotes the probability ¯ denotes the expected value of X, and cov(X) or density function (pdf) of X, E[X] or X ΣX denotes the covariance matrix of X. For random vectors X and Y , the cross covariance ¯ between X and Y is denoted by cov(X, Y ) = ΣXY = E[(X − X)(Y − Y¯ )]. X | Y denotes that the random vector X is conditioned on realization Y , with conditional pdf p(X | Y ). xˆ denotes the estimated value of variable x.

4


Page 4 of 40

Page 5 of 40

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


Background: Model Selection Consider a finite number of possible model hypotheses {Hi } with known prior probabilities P (Hi ). We look to investigate the decision-theory problem of classifying a system observation y, under a given manipulatable input u, as the prediction of one of the model hypotheses {Hi }. It is well-known that the decision rule that minimizes the probability of error P (error) is the Bayes decision rule, which selectes the hypothesis with the largest a posteriori probability. 26 The problem of model selection (aka Bayesian hypothesis selection) can be expressed as follows select Hi such that i ∈ argmaxj P (Hj | y, u).

(1)

Using Bayes rule P (Hj | y, u) = p(y | Hj , u)P (Hj )/p(y | u), the selection criterion (1) yields a number of decision regions Ri defined by

Ri = {y : p(Hi , y | u) > p(Hl , y | u), ∀l 6= i},

(2)

as shown in Figure 1. If y ∈ Ri , then hypothesis Hi is selected. This implies that P (y ∈ Rj , Hi | u) is the probability of incorrectly selecting the model hypothesis j when the true hypothesis is i. By adding all these probabilities across all possible model pairs, the following expression is derived for the probability of selecting the incorrect hypothesis, known as the Bayes risk

P (error) =

XX i

=

XX i

=

j6=i

j6=i

P (y ∈ Rj , Hi | u) P (y ∈ Rj | Hi , u)P (Hi )

XXZ i

j6=i

Rj

5

p(y | Hi , u)P (Hi )dy.


(3)


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 ACS Paragon Plus Environment

Page 6 of 40

Page 7 of 40

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


Rnx × Rnu × Rnw × Rnθ → Rnx and h : Rnx × Rnv × Rnθ → Rny describe the nonlinear system dynamics and measurement functions, respectively. Due to imperfect knowledge of the system, the initial conditions x0 and parameters θ in (4) are modeled as probabilistic uncertainties with a known joint pdf p(x0 , θ). The system noise wk and measurement noise vk are assumed to be independent and identically distributed sequences with known pdfs that are independent of x0 and θ. The system is subject to hard input constraints uk ∈ U and (possibly nonlinear) state constraints xk ∈ X ⊆ Rnx . Since the system evolves as a stochastic process, the state constraints are enforced as a chance constraint of the form

P (xk ∈ X | H) ≥ 1 − ǫ,

k = 1, 2, . . . ,

(5)

where ǫ is the maximum allowed probability of state constraint violation. Assume that a finite number nf of models, with known dynamics of the form (4), exist to describe potential fault scenarios of the system. Denote the set of models by H = {H0 , H1 , . . . , Hnf }, where Hi represents the system information required for specifying the ith model, i.e., Hi = {f [i] , h[i] , p[i] (x0 , θ), p[i] (wk ), p[i] (vk ), X [i] , ǫ[i] }. The superscript [i] is used to distinguish variables in (4) for each model Hi ∈ H. The nominal system model (i.e., no faulty behavior) is denoted by the superscript i = 0, while i = 1, . . . , nf refer to the different fault models. Note that all models must have the same manipulated inputs u and measured outputs y, suggesting that nu and ny are the same for all Hi ∈ H. For notational convenience, the superscripts [i] are dropped in the remainder of the paper unless otherwise stated.

Active Fault Diagnosis The goal of this work is to present an optimization framework for designing an input sequence that minimizes the probability of model selection error in H, while satisfying the input and 7



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

state constraints. The constraints U and X must be defined such that the designed input sequence ensures safe and least intrusive system operation during the AFD experiment. Due to the stochastic nature of system uncertainty in (4), the optimal input sequence should be designed such that the system evolves through a trajectory of state probability distributions that minimizes the Bayes risk. 15 For a finite (prespecified) planning horizon N , the AFD problem can be stated as follows. Problem 1 (Active fault diagnosis). Determine the input sequence u0:N −1 that minimizes the probability of selecting the wrong model P (error) while satisfying input constraints uk ∈ U and chance constraints (5) over the planning horizon N for all possible models Hi ∈ H. The designed input sequence should be applied to the system to observe its behavior and, accordingly, select the model hypothesis that best (in a probabilistic sense) describes the measurements. The model with the highest probability after N time steps should be selected. However, the AFD experiment can be terminated early if any of the models reaches a user-defined probability threshold before time step N . For early termination, the model probabilities must be evaluated online, as discussed later. AFD can be performed in an open-loop fashion where the optimization in Problem 1 is solved only once, offline, yielding an input sequence that is subsequently applied to the system. In offline AFD, system measurements collected during the implementation of the optimal input sequence are not used for re-planning of the AFD experiment. This can reduce the effectiveness of the designed input sequence due to unknown system uncertainties and disturbances. Alternatively, the AFD problem can be solved online in a closed-loop fashion. In this case, Problem 1 is solved repeatedly at every measurement sampling time and only the first part of the input sequence is applied to the system (often referred to as a receding-horizon implementation).

8


Page 8 of 40

Page 9 of 40

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


Challenges of Active Fault Diagnosis Online AFD requires solving the optimization problem faster than the measurement sampling intervals so that the inputs can be implemented in real time. This requirement can be hard to meet for certain applications due to the following challenges inherent to Problem 1: 1. Stochastic uncertainty propagation: Stochastic uncertainties x0 , θ, and w0:N −1 must be propagated through the nonlinear dynamics of each model Hi ∈ H to evaluate the distribution of the concatenated observation vector p(y0:N | Hi , u0:N −1 ). This is required for computing the Bayes risk P (error) in (3). Traditional sample-based approaches to uncertainty propagation, such as Markov-Chain Monte Carlo, 27 are prohibitively expensive and ill-suited for gradient-based optimization methods. 2. Computation of the Bayes risk P (error): There exist no closed-form expressions for the distribution p(y0:N | Hi , u0:N −1 ) or the decision regions Ri , making evaluation of P (error) particularly challenging. This is further compounded by the multivariate integral in the Bayes risk. 3. Chance constraint evaluation: Chance constraints are generally nonconvex and intractable. 28,29 This is because chance constraint evaluation involves solving a multivariate integral of state pdfs, for which no closed-form expression exists. Although construction of the multivariate distribution p(y0:N | Hi , u0:N −1 ) is inherently difficult, p(y0:N | Hi , u0:N −1 ) can be efficiently approximated in terms of its moments. This is the approach taken in this work to obtain an efficient solution to Problem 1. Some wellknown methods proposed in the literature for moment-based approximation of pdfs include (i) linearization, 23 (ii) the unscented transform, 24 and (iii) generalized polynomial chaos (gPC). 22 In this paper, we focus on the first two methods and interested readers are referred to Refs. 20,30,31 for the application of gPC to input design problems. To address the second challenge stated above, a tractable upper bound is derived for the Bayes risk P (error), which 9



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

is computationally amenable to online optimization. Lastly, a moment-based surrogate is derived for approximating the chance constraint (5) to avoid sample-based integration or construction of the state pdfs.

Uncertainty Propagation A closed-form solution for exact propagation of stochastic system uncertainty (i.e., uncertain parameters and/or system noise) through general nonlinear models does not exist. 32 Sampled-based approaches to uncertainty propagation can be expensive for optimization applications since they require many model simulation runs that must be repeated for any candidate decision variables. In this work, we look to propagate moments of p(y0:N | Hi , u0:N −1 ) to avoid sample-based construction of the distribution. To reduce the computational complexity of Problem 1, p(y0:N | Hi , u0:N −1 ) is approximated as a multivariate Gaussian distribution so that it can be fully represented in terms of its first two moments. Two approaches are presented for obtaining closed-form expressions for the moments of p(y0:N | Hi , u0:N −1 ) as a function of the input sequence. The first approach relies on linearization of the nonlinear model equations, which is a common approach used in, for example, the design of the extended Kalman filter. 23 On the other hand, in the second approach linearization is circumvented by using the unscented transformation (UT) for calculating statistics of a random variable that undergoes a nonlinear transformation. 24 The UT method is particularly useful when the system dynamics are highly nonlinear, so that the linearization method provides a poor approximation of the nonlinear dynamics. In the UT method, a set of samples, called sigma points, are selected around the mean of the random variables. The sigma points are then propagated through the nonlinear model equations to construct new mean and covariance estimates for the random variables. A key advantage of the UT method is that Jacobian calculations, which can be cumbersome for complex functions or even impractical for non-smooth/non-differentiable functions, are no longer required.

10


Page 10 of 40

Page 11 of 40

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


Next, we show how the linearization and UT methods can be used for propagating the moments of p(y0:N | Hi , u0:N −1 ).

Augmenting States with Parameters ⊤ ⊤ To simplify presentation, a new “augmented” state vector is defined as zk = x⊤ , k , θk where the system dynamics (4) are augmented with artificial “parameter dynamics” θk+1 = θk

that enforce the time-invariant nature of the unknown parameters. The augmented system dynamics are described by 







xk+1  f (xk , uk , wk , θk )   = θk θk+1 yk = h(xk , vk , θk )

→

zk+1 = fz (zk , uk , wk )

(6a)

→

yk = hz (zk , vk ),

(6b)

where fz : Rnz × Rnu × Rnw → Rnz and hz : Rnz × Rnv → Rny are defined appropriately, with nz = nx + nθ and p(z0 ) = p(x0 , θ) denotes the pdf of the initial augmented state.

Linearization The nonlinear dynamics of the augmented system (6a) can be linearized using a first-order Taylor series

zk+1 ≈ fz (ˆ zk , uk , w) ¯ + ∂zk fz (ˆ zk , uk , w) ¯ (zk − zˆk ) + ∂wk fz (ˆ zk , uk , w) ¯ (wk − w), ¯ | | {z } {z } Ak

(7)

Ek

where zˆk ≈ z¯k = E[zk ] denotes the approximated mean of the augmented states and w¯ = E[wk ]. Taking the expectation of (7) yields the following recursion for the mean

zˆk+1 = fz (ˆ zk , uk , w), ¯

11

zˆ0 = z¯0 .


(8)


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 12 of 40

Define δzk = zk − zˆk , so that (7) can be rewritten as a linear time-varying (LTV) system δzk+1 = Ak δzk + Ek δwk ,

(9)

ˆ z = E[δzk δz ⊤ ] ≈ Σz as the where δz0 has a known pdf and δwk = wk − w. ¯ Define Σ k k k approximate covariance of the augmented states. A recursive expression for this covariance can now be directly derived from (9) ⊤ ˆ z = Ak Σ ˆ z A⊤ Σ k + E k Σw E k , k+1 k

ˆ z0 = Σz0 . Σ

(10)

The covariance of the augmented states over the planning horizon N (including time crosscorrelation) can be approximated by stacking (9) into δz0:N = Aδz0 + Eδw0:N −1 , where 



I       A 0        A1 A0  , A=   . .   .   Q   N A   i=1 N −i  Q  N A i=0 N −i

such that





0 0 ··· ··· 0       E 0 0 0     ..   ...  A1 E 0 E1 .   E=   . . . .  A2 A1 E0 . A2 E1 E2 .      . . .   .. .. .. 0   Q  Q N −1 N −1 A E · · · A E E A E N −i 1 N −1 N −2 N −1 N −i 0 i=2 i=1

ˆ z = AΣz0 A⊤ + Ediag(Σw , . . . , Σw )E⊤ . Σ 0:N

(11)

Note that A and E are functions of the input sequence u0:N −1 (since the linearization is performed locally with respect to an input trajectory) as well as the moments z¯0 and Σz0 of the initial augmented states. Likewise, the measurement function (6b) can be linearized

12


Page 13 of 40

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


around zˆk and v¯ = E[vk ]

z , v¯)(v − v¯), yk ≈ hz (ˆ zk , v¯) + ∂zk hz (ˆ z , v¯)(z − zˆk ) + ∂vk hz (ˆ | {z k } k {z k } k | Ck

(12)

Mk

which leads to the following mean and covariance approximations for the system outputs (13a)

yˆk = hz (ˆ zk , v¯) ˆ y = Ck Σ ˆ z C ⊤ + M k Σv M ⊤ . Σ k k k k

(13b)

Stacking (12) over the planning horizon N leads to δy0:N = Cδz0:N + Mδv0:N with C = diag(C0 , . . . , CN ) and M = diag(M0 , . . . , MN ). The predicted covariance of the outputs over the horizon N is then given by ˆ y = CΣ ˆ z C⊤ + Mdiag(Σv , . . . , Σv )M⊤ . Σ 0:N 0:N

(14)

ˆy Note that Σ is a symmetric matrix that accounts for cross-correlation in the system 0:N outputs at different times, i.e., 



ˆ y ≈ Σy Σ 0:N 0:N

 Σy0 Σy0 y1 · · · Σy0 yN      Σy y Σ · · · Σ y1 y1 yN   10 . = . .. ..  ...  .. . .      Σ y N y 0 Σ y N y1 · · · Σ yN

(15)

The accuracy of the above discussed uncertainty propagation method based on model linearization can be improved by retaining the higher-order terms in (7). However, retaining more terms in the Taylor series expansion can increase the computational burden of the uncertainty propagation substantially and, therefore, should only be done when the increased accuracy is considerable. Including higher-order terms in the Taylor series expansion allows for accurate computation of higher-order moments (such as skewness and kurtosis). Al13



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

though the higher-order moments are not used in the proposed AFD method in this paper, they lead to a more precise approximation of the distribution of the outputs (e.g., when distributions are highly skewed) at the expense of a higher computational cost. Remark 1. The dynamics are linearized only with respect to the states and the system noise (not the inputs) for improved accuracy. In this way, the dynamics are locally linearized around any candidate input sequence that captures variability in the covariance of the predicted outputs through changes in the Jacobian matrices.

Unscented Transformation UT method. The general notion of the UT method is first described for propagation of an n-dimensional random vector x ∈ Rn through a nonlinear function y = g(x). Let x have a known mean x¯ and covariance matrix Σx . Define a set of 2n + 1 sigma points {Xi }2n i=0 X0 = x¯

p Xi = x¯ + (n + λ)Σx i p (n + λ)Σx Xi = x¯ −

i = 1, . . . , n i = n + 1, . . . , 2n,

i

(m) 2n }i=0

with the mean weights {Wi (m)

W0

(m)

Wi

(c)

and covariance weights {Wi }2n i=0 given by

λ λ (c) , W0 = + (1 − α2 + β) (n + λ) (n + λ) 1 (c) , i = 1, . . . , 2n. = Wi = 2(n + λ)

=

In the above expressions, (A)i is the ith column of the matrix A; λ = α2 (n+κ)−n is a scaling parameter; α specifies the spread of the sigma points around the mean; κ is another scaling parameter; and β is used to account for prior knowledge of the pdf of x. Some common values for these parameters are α = 1, κ = 2, and β = 0. 24 Note that β = 2 is the optimal choice for β when x has a Gaussian distribution. The UT method relies on propagating the 14


Page 14 of 40

Page 15 of 40

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


sigma points through the nonlinear function y = g(x), i.e.,

Yi = g(Xi ),

i = 0, . . . , 2n,

and approximating the mean and covariance of y based on the weighted sample mean and covariance of Yi y¯ ≈ Σy ≈

2n X

(m)

Wi

i=0

2n X i=0

Yi

(c)

Wi (Yi − y¯)(Yi − y¯)⊤ .

The UT method requires 2n + 1 function evaluations, which scales linearly with respect to n (dimension of the uncertainty). As a result, this method may be substantially cheaper than general sample-based methods (such as Monte Carlo) that typically require orders of magnitude more sample points for accurate computation of y¯ and Σy . 24 Note that the sigma points are chosen heuristically in the UT method, so that its accuracy must be tested for the chosen points. On the other hand, the convergence rate of Monte Carlo methods is √ O(1/ N ), where N denotes the number of samples. The latter expression implies that quadrupling the number samples reduces the error in half regardless of the state dimensions n. 33 Prediction scheme using UT. The UT method can be adopted for approximating the mean and covariance of the predictions of (4) over the planning horizon N . Define the concatenated vector of augmented states, system noise, and measurement noise as ⊤ xak = zk⊤ , wk⊤ , vk⊤ . Over the planning horizon k = 0, 1, . . . , N − 1, we perform the following steps: (i) generate ⊤ ˆ xa = the sigma points Xka = [(Xkz )⊤ , (Xkw )⊤ , (Xkv )⊤ ]⊤ using xâk = zˆk⊤ , w¯ ⊤ , v¯⊤ and Σ k 15



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

(m)

ˆ z , Σw , Σv ) and compute the weights W diag(Σ i k

(c)

, Wi

Page 16 of 40

as described above, (ii) propagate

z the latter sigma points through the dynamics Xk+1 = fz (Xkz , Xkw ), and (iii) compute the

mean and covariance estimates of the augmented states from

zˆk+1 = ˆz = Σ k+1

2n X

i=0 2n X i=0

(m)

Wi

(16a)

z Xi,k+1

(c)

z z Wi (Xi,k+1 − zˆk+1 )(Xi,k+1 − zˆk+1 )⊤ ,

(16b)

where n = nz + nw + nv . Using the sigma points of the augmented states determined above, we can now approximate the mean of the output and the various blocks in the output covariance matrix (15) over the planning horizon. To this end, define the output sigma points as Yk = hz (Xkz , Xkv ) for k = 0, . . . , N . Then, apply the UT method for approximating the mean and covariances as

yˆk =

2n X

(m)

Wi

i=0

ˆy y = Σ k l

2n X i=0

Yi,k ,

(17a)

∀k = 0, 1, . . . , N,

(c)

Wi (Yi,k − yˆk )(Yi,l − yˆl )⊤ ,

∀k, l = 0, 1, . . . , N.

(17b)

Note that only N (N + 1)/2 blocks of the output covariance (15) should be calculated due to its symmetry.

Framework for Active Fault Diagnosis This section presents a computationally tractable surrogate for Problem 1 by using the uncertainty propagation methods discussed in the previous section. The proposed AFD problem relies on a bound on the Bayes risk P (error), which is used as the fault diagnosis criterion for posing the optimization in Problem 1. In the AFD problem, the state chance constraint is approximated using a moment-based surrogate. This section also presents an 16


Page 17 of 40

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


algorithm for receding-horizon implementation of the AFD problem, which enables online re-planning of the fault diagnosis experiment.

Tractable Criterion for Hypothesis Selection Bound on the Bayes risk. There exists no closed-form expression for the Bayes risk due to the multivariate nature of the integral in (3). This poses a key challenge to direct use of the Bayes risk in the AFD Problem 1. A tractable upper bound on the Bayes risk, known as the Bhattacharrya bound , 26 has been derived for the case of multiple hypotheses. 16 For a general system observation y and manipulated input u, the result is summarized as follows. Theorem 1. For hypothesis selection between a finite number of hypotheses, the Bayes risk is upper bounded by XX

1

1

P (Hi ) 2 P (Hj ) 2 Bij (u),

(18)

Z q Bij (u) = p(y | Hi , u)p(y | Hj , u)dy

(19)

P (error) ≤

i

j>i

where

is the Bhattacharrya coefficient between likelihoods of the ith and j th hypotheses. Proof. The key notion is to express (3) in terms of the Bayes risk for each of the pairs of hypotheses, and then bound the Bayes risk between each pair using the established Bhattacharrya bound. 16 The Bhattacharrya coefficient as AFD criterion. The bound (18) alleviates the need for computing the decision regions Ri in (3). However, calculating the Bhattacharrya coefficient Bij hinges on (possibly expensive) high-dimensional numerical integration. In the case of Gaussian system observations, Bij can be evaluated analytically to circumvent the 17



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 18 of 40

numerical integration in (19). This leads to the following approximation [i]

p(y0:N | Hi , u0:N −1 ) ≈ N (¯ y0:N , Σy[i]0:N ),

∀Hi ∈ H,

[i]

[i]

where the mean y¯0:N and covariance Σy0:N can be obtained using the uncertainty propagation methods described in the previous section. Under this approximation, the Bhattacharrya coefficient (19) simplifies to [i]

[i]

[j]

[j]

(20)

Bij (u0:N −1 ) = e−d(¯y0:N ,Σy0:N ,¯y0:N ,Σy0:N ) , where 1 1 d(x, X, y, Y ) = (x − y)⊤ [X + Y ]−1 (x − y) + ln 4 2

det

X+Y 2

p det(X) det(Y )

!

.

(21)

Using the error bound (18) on the Bayes risk P (error), a tractable objective function for the AFD Problem 1 can now be defined as

J(u0:N −1 ) =

nf nf X X

1

1

[i]

[i]

[j]

[j]

P (Hi ) 2 P (Hj ) 2 e−d(¯y0:N ,Σy0:N ,¯y0:N ,Σy0:N ) .

(22)

i=0 j=i+1

Moment-based Chance Constraint Approximation A key feature of the AFD Problem 1 is to enforce the system inputs and states to satisfy their constraints during the fault diagnosis experiment. Input constraints typically arise from practical considerations such as hard bounds on the manipulated variables (i.e., umin ≤ uk ≤ umax ) or the maximum allowable system perturbation for fault diagnosis, which is commonly defined in terms of the “energy” of the input sequence. 34,35 On the other hand, state constraints are generally intended to ensure that the system remains in a “safe” operating region (or in a region for which the model hypotheses have been validated) as well as to ensure that certain system performance requirements are met during the fault diagnosis 18


Page 19 of 40

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


experiment. In the system (4), the state constraints can be enforced only in probability, i.e., the chance constraint (5), because of the stochastic system uncertainty. Without loss of generality, assume that the feasible region X , {x : A⊤ x ≤ b} is a polytope defined by matrix A ∈ Rnx ×nc and vector b ∈ Rnc . Since evaluating the probability distribution p(xk | H, u0:k−1 ) can be computationally prohibitive, we look to derive a momentbased approximation for the chance constraint (5) that is not overly conservative. To this end, p(xk | H, u0:k−1 ) is approximated to have a Gaussian distribution, so that it can be represented in terms of its first two moments. By defining xk = Tk z0:N , where Tk is a matrix of all zeros except for Inx in the position of the k th step states, the latter moments can be directly computed as E[xk | H] = x¯k = Tk z¯0:N and cov(xk | H) = Σxk = Tk Σz0:N Tk⊤ . For any normally distributed random variable x ∼ N (¯ x, Σx ) of dimension n, enforcing x ∈ X with at least probability 1 − ǫ can be expressed exactly in terms of an integral condition over X . However, the integral constraint is difficult to enforce directly since there is no closed-form solution. A simple alternative relaxation can be used instead, which involves 2 finding the largest ellipsoid of radius r, denoted by Er , x : x⊤ Σ−1 , inside of the x x ≤ r

polytope X . If it can be ensured that x¯ ⊕ Er ⊂ X , then P (x ∈ X ) > P (x ∈ x¯ ⊕ Er ) must hold by the axioms of probability. The expression P (x ∈ x¯ ⊕ Er ) can be stated in terms of the chi-squared distribution P (x ∈ x¯ ⊕ Er ) = P ((x − x¯)⊤ Σ−1 ¯) ≤ r2 ) = Fχ2 (r2 ; n), x (x − x

where Fχ2 (·; n) is the cdf of the chi-squared distribution with n degrees of freedom and the radius r is selected so that Fχ2 (r2 ; n) = 1 − ǫ holds. This implies that P (x ∈ X ) can be replaced with the stronger constraint x¯ ⊕ Er ⊂ X , which is equivalent to requiring that th the ellipsoid lies within a collection of half-spaces Hj = {x : a⊤ j x ≤ bj } where aj is the j

column of A and bj is the j th element of b. The constraint x¯ ⊕ Er ⊂ X is equivalent to 36

19




Page 20 of 40

Page 21 of 40

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


constraints that it can handle. The AFD problem (24) is a nonlinear programming (NLP) problem that can be solved efficiently to local optimality using various nonlinear optimization methods such as interior point methods 37 and sequential quadratic programming. 38 The optimal input sequence, denoted by u⋆0:N −1 , will enable discrimination between the model hypotheses in a probabilistic sense since the system uncertainty is of stochastic nature. Remark 2. There are two main optimization approaches for solving (24): (i) a condensed approach in which only the inputs are treated as decision variables, so that the dynamic equations become dense nonlinear equality constraints, and (ii) a large-scale sparse approach in which the inputs and states are treated as decision variables, so that the dynamic equations form a simpler structure. Which approach suits the optimization problem (24) best in light of the complexity of the AFD problem remains to be investigated. ⋆ ⋆ } that correspond to the sys= {y0⋆ , , . . . , yN In offline AFD, once the measurements y0:N

tem evolution under the optimal input sequence u⋆0:N −1 have been collected, the probability of each model hypothesis can be evaluated using Bayes rule 16 p(y ⋆ | Hi , u⋆0:N −1 )P (Hi ) ⋆ , u⋆0:N −1 ) = Pnf 0:N ⋆ P (Hi | y0:N . ⋆ i=0 p(y0:N | Hi , u0:N −1 )P (Hi )

(25)

⋆ The main challenge in (25) arises from computing the likelihood function p(y0:N | Hi , u⋆0:N −1 ).

The likelihood function can be approximated using sample-based methods such as particle filtering. 39 Alternatively, approximating the likelihood function as a Gaussian distribution, which is also done for deriving a bound on the Bayes risk in (18), greatly simplifies the computation of (25). In the next section, an algorithm is presented for online implementation of the AFD method, which relies on solving the optimization problem (24) in tandem with Bayes rule (25) in real time as the fault diagnosis experiment is conducted.

21



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 22 of 40

Online Active Fault Diagnosis The unmodeled system uncertainties and disturbances in the AFD problem (24) can impede effective fault diagnosis when the optimal input sequence u⋆0:N −1 is designed completely offline. On the other hand, as the input sequence is applied to the true (but unknown) system, system measurements can be used to inform AFD in real time via updating the estimates of the augmented states and the probabilities of each model Hi ∈ H. The measurement feedback to the AFD problem (24) is likely to give some degree of robustness to unmodeled system uncertainties and disturbances when designing the input sequence. The receding-horizon implementation of the AFD problem hinges on recursively updating the model probabilities online based on Bayes rule (25) and, subsequently, redesigning the input sequence in (24) using the most recent system information.

Recursive Calculation of Model Probabilities Online Given the system measurements y0:k and inputs u0:k−1 up until sampling time k stacked into a single information vector I k = {y0:k , u0:k−1 }, the probability of each model Hi ∈ H at every k can be evaluated recursively using Bayes rule p(yk | y0:k−1 , u0:k−1 , Hi )P (Hi | I k−1 ) , P (Hi | I ) = Pnf k−1 ) i=0 p(yk | y0:k−1 , u0:k−1 , Hi )P (Hi | I k

P (Hi | I −1 ) = P (Hi ).

(26)

The likelihood function p(yk | y0:k−1 , u0:k−1 , Hi ) in fact represents the so-called evidence (aka innovation), which can be obtained from state estimation for the given model Hi . As is wellknown from the Bayesian estimation literature, 32 no closed-form solution exists for state estimation of nonlinear systems. In this work, the evidence is approximated as a multivariate Gaussian distribution p(yk | y0:k−1 , u0:k−1 , H) ≈ N (¯ yk|k−1 , Σyk|k−1 ), where yk|k−1 = yk | y0:k−1 , u0:k−1 , H is a conditional random variable representing the evidence with mean y¯k|k−1 and covariance Σyk|k−1 . The extended Kalman filter (EKF) 23 and the unscented Kalman filter (UKF) 24 are straightforward 22


Page 23 of 40

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


extensions of the linearization and UT uncertainty propagation methods that can be used for estimating the mean and covariance of the evidence. These filters both rely on a model ˆz prediction and a measurement update step. The prediction step calculates zˆk|k−1 , Σ , k|k−1 ˆy yˆk|k−1 , and Σ using either (8), (10), (13a), and (13b) in the EKF or (16a), (16b), (17a), k|k−1 and (17b) in the UKF. In the update step, the latter model predictions are updated using the system observations. To this end, the optimal filter gain is defined by ˆz ˆ −1 Kk = Σ Σ yk|k−1 . k|k−1 yk|k−1

(27)

The estimates of the mean and covariance of the posterior distribution of the augmented states can then be updated as

zˆk|k = zˆk|k−1 + Kk (yk − yˆk|k−1 )

(28a)

ˆz = Σ ˆz ˆy Σ − Kk Σ Kk⊤ . k|k k|k−1 k|k−1

(28b)

The prediction and update steps in the estimator must be repeated at every time instant k. This procedure yields the estimates of y¯k|k−1 and Σyk|k−1 , which are required for approximating the evidence in (26). The recursive estimation procedure is summarized in Algorithm 1. Note that the convergence of the estimated model probabilities to the true probabilities in the Bayesian estimation framework remains a theoretically open problem for general nonlinear systems.

Online Design of Input Sequence In online AFD, the optimization problem (24) should be re-solved based on newly obtained system information I k at every measurement sampling time k = 0, . . . , N − 1. This amounts to minimizing P (error | I k ) subject to the system constraints. Hence, for online AFD, the [i]

[i]

priors P (Hi ) and initial conditions z¯0 , Σz0 in (24) must be replaced with the updated priors

23



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 24 of 40

Algorithm 1 Recursive estimation of states, parameters, and model probabilities. Require: Specify model H and priors P (H) for all H ∈ H. 1: Offline: ˆz 2: For all H ∈ H, calculate the initial estimates zˆ0|−1 and Σ from p(z0 ) = p(x0 , θ). 0|−1 3: Online: 4: for k = 0, 1, 2, . . . do 5: Get measurements yk from the true system. 6: for all models H ∈ H do ˆz ˆy 7: Calculate yˆk|k−1 , Σ from zˆk|k−1 , Σ with linearization or UT method. k|k−1 k|k−1 k ˆy 8: Compute P (H | I ) from (26) with p(yk | y0:k−1 , u0:k−1 , H) ≈ N (ˆ yk|k−1 , Σ ). k|k−1 ˆ z using (28). 9: Determine posterior state estimate zˆk|k and covariance Σ k|k 10: end for 11: Apply input uk to the true system. 12: for all models H ∈ H do ˆz ˆ z and uk using linearization or UT method. 13: Calculate zˆk+1|k , Σ from zˆk|k , Σ k|k k+1|k 14: end for 15: end for [i] ˆ z[i] , respectively. These quantities are readily P (Hi | I k ) and updated state estimates zˆk|k , Σ k|k

available from the recursive evaluation of Bayes rule (26) detailed in Step 9 of Algorithm 1. Hence, the online AFD problem that is solved at each sampling time k = 0, . . . , N − 1 takes the form

min

uk:N −1|k

nf nf X X

i=0 j=i+1

1

1

P (Hi | I k ) 2 P (Hj | I k ) 2 e

[i]

[i]

[j]

[j]

−d(¯ yk:N |k ,Σyk:N |k ,¯ yk:N |k ,Σyk:N |k )

(29)

[i] ˆ z[i] , subject to: mean and covariance dynamics initialized with zˆk|k and Σ k|k q [i] [i] ut|k ∈ U, a⊤ ¯t+1|k + r a⊤ j x j Σxt+1|k aj ≤ bj ,

∀i = 0, . . . , nf , ∀j = 1, . . . , nc , ∀t = k, . . . , N − 1.

The receding-horizon implementation of the AFD problem (29) follows the same recursive estimation procedure detailed in Algorithm 1 with the input uk = u⋆k|k in Step 11 being the first element of the optimal input sequence u⋆k:N −1|k .

24


Page 25 of 40

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


Simulation Case Study System Description The performance of the proposed active fault diagnosis framework is demonstrated on a biochemical reaction occurring in a continuous well-stirred bioreactor. Assuming constant volume, the continuous-time process dynamics are described by 40,41 (30a)

dX = (−DX + µX) dt + σX dwX (t) 1 dS = D(Sf − S) − µX dt + σS dwS (t) YX/S

(30b)

dP = (−DP + (αµ + β)X) dt + σP dwP (t),

(30c)

where X, S, and P are, respectively, the concentration of biomass, substrate, and product in the bioreactor; D is the dilution rate (equal to the inlet feed divided by the reactor volume); Sf is the concentration of substrate in the inlet feed; YX/S is the yield of biomass per substrate consumed; α and β are the yield parameters for the production of P ; and wX (t), wS (t), and wP (t) are independent, zero-mean unit-variance Wiener processes scaled by standard deviations σX , σS , and σP , respectively. The quantity µ represents the (substrate-dependent) growth rate of biomass that is given by

µ=

µmax S , KM + S

(31)

where µmax denotes the maximum growth rate; and KM is an affinity constant. The dilution rate D is the only manipulated input of the process. The nominal parameter values are listed in Table 1. The initial conditions for the states are assumed to be X(0) = 7.038 g/L, S(0) = 2.404 g/L, and P (0) = 24.87 g/L.

25



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 26 of 40

Table 1: Nominal parameters of the continuous bioreactor. 40,41 Variable Nominal Value Units YX/S 0.4 g/g α 2.2 g/g β 0.2 hr−1 µmax 0.48 hr−1 KM 1.2 g/L Sf 20 g/L σX 0.25 g/L σS 0.25 g/L σP 0.25 g/L

Normal Operation and Fault Scenarios Normal operation is described by (30) with parameters listed in Table 1. The nominal model is the same as that used for normal operation, however, the parameter µmax is treated as an uncertainty since it cannot be perfectly estimated from data. The prior distribution µmax ∼ N (¯ µmax , Σµmax ) is assumed to be normal with mean µ ¯max = 0.6 hr−1 and variance Σµmax = 0.05 hr−1 . Note that the mean of µmax is different from that of its true value. Under continuous operation of the bioeractor, faults must be swiftly identified in order to prevent product loss, which can lead to negative economic consequences. To demonstrate the efficacy of the proposed AFD framework, two different fault scenarios relative to normal operation are considered. First, a change in the model structure is considered in which inhibition of the biochemical reaction takes place due to the presence of excess substrate (commonly referred to as substrate inhibition). In this fault scenario, the specific growth rate µ is described by

µ= KM

µmax S +S 1+

S KI

,

(32)

instead of (31), where KI is the rate constant that specifies the strength of inhibition. The parameters for the substrate inhibition fault model are the same as those used for the nominal model with KI = 1 g/L. The second fault scenario is an operation fault in which the substrate 26


Page 27 of 40

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


concentration in the inlet feed drops by 20%, i.e., Sf ← 0.8Sf , resulting in Sf = 16 g/L. Three model hypotheses based on the model structure (30) are used to represent the normal operation and the two fault scenarios. The continuous-time models are discretized to obtain discrete-time models of the form (4) with states x = [X, S, P ]⊤ and input u = D. Measurements of S and P are taken every 0.2 hr. Both measurements are corrupted with zero-mean Gaussian white noise with standard deviations of 0.03 g/L and 0.2 g/L, respectively. The initial states are not known exactly and are represented by the prior x0 ∼ N (¯ x0 , Σx0 ) with x¯0 = [6.5, 2.75, 26]⊤ g/L and Σx0 = diag(0.5, 0.5, 0.5) g2 /L2 . Throughout the case study, the model probabilities are estimated online using the Bayesian estimation method described in Algorithm 1.

Performance Evaluation of Offline AFD The dynamics of the three operation scenarios are fairly similar, suggesting that fault diagnosis might be challenging in the presence of uncertainty in the maximum growth rate as well as disturbances and measurement noise. The difficulty in fault diagnosis under a nominal input sequence 41 (i.e., constant dilution rate at 0.15 hr−1 ) is demonstrated in Figure 3. The figure shows the evolution of the probability of each model, P (Hi | I k ), for 10 Monte Carlo runs when the normal operation model, substrate inhibition model, and substrate offset model are active. In all three cases, the substrate inhibition model initially has the highest probability of being correct, which is clearly inaccurate. In fact, under normal operation, the estimated probability of the nominal model being true is less than 50% over the entire operation, as shown in the top plot in Figure 3. To effectively discriminate between the three model hypotheses, we look to determine the optimal dilution rate trajectory, denoted by D⋆ , that minimizes the bound on the Bayes risk P (error) by solving the AFD problem (24). The optimization problem was solved using the Matlab R2016 (Version 9.1) constrained optimization routine fmincon on a 2015 MacBook Pro with 8 GB of RAM and a 3.1 GHz Intel i7 processor. Figure 4 depicts ten Monte 27



1 Normal operation Substr. inhibition Substr. offset

0.9 0.8

Model probability

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.5

1

1.5

2

Time (hr) 1 0.9

Model probability

0.8 Normal operation Substr. inhibition Substr. offset

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.5

1

1.5

2

Time (hr) 1 Normal operation Substr. inhibition Substr. offset

0.9 0.8 0.7

Model probability

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 28 of 40

0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.5

1

1.5

2

Time (hr)

Figure 3: Ten Monte Carlo runs of the evolution of probability of each model being active using nominal input sequence under normal operation (top), substrate inhibition (middle), and substrate concentration offset in the feed (bottom).

28


Page 29 of 40

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


Carlo runs of the evolution of the model probabilities P (Hi | I k ) with D⋆ applied to the process over a span of two hours. In all three cases, the active model is correctly identified reasonably fast (less than 1 hour). Note that subsequent measurements further improve upon the confidence in the true model since the probability of the active model gradually increases. This illustrates the significance of the input design for fault diagnosis in contrast to the nominal input sequence. Note that the efficiency and accuracy of gradients required for solving the optimization problem has a strong influence on the performance of the optimizer. For simplicity of implementation, this case study relied on the finite difference method for evaluating the gradients in the optimization routine fmincon. It is well-known, however, that finite differences can be slow and inaccurate as the size of the problem grows. Automatic differentiation (AD) is an emerging alternative that exploits the fact that any differentiable function can be decomposed into elementary operations (e.g., addition, multiplication, exponential functions) on which the chain rule can be applied. 42 This enables computation of derivatives up to machine precision at a cost that is on the same order of magnitude as evaluating the original function when sparsity is exploited.

Comparing Offline to Online AFD The above results revealed the efficacy of offline AFD, where the optimal input sequence was designed only once, offline. As such, when a deviation from normal operation is detected, for example using process monitoring methods, 43,44 the optimal sequence D⋆ should be applied to the process in order to diagnose the current mode of operation using measurements collected during the experiment. When the tractable AFD problem (24) can be solved within the measurement sampling intervals, it will often be advantageous to re-design the input sequence online based on the most recent system observations. Due to regular adaptation of the model probabilities and state/parameter estimates, online fault diagnosis can partly account for the effects of unmodeled system uncertainties and disturbances, which can be detrimental to the 29



1 0.9

Model probability

0.8

Normal operation Substr. inhibition Substr. offset

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.5

1

Time (hr)

1.5

2

1 0.9

Normal operation Substr. inhibition Substr.offset

Model probability

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.5

1

1.5

2

Time (hr) 1 0.9 Normal operation Substr. inhibition Substr. offset

0.8 0.7

Model probability

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 30 of 40

0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.5

1

1.5

2

Time (hr)

Figure 4: Ten Monte Carlo runs of the evolution of probability of each model being active using optimal input sequence D⋆ under normal operation (top), substrate inhibition (middle), and substrate concentration offset in the feed (bottom).

30


Page 31 of 40

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


effectiveness of input sequences designed offline. For the case that the second fault scenario is active, Figure 5 shows the input sequences designed using the online AFD method and the offline AFD method under the same uncertainty realizations. A shrinking-horizon implementation is utilized for online AFD in which the final time of the planning horizon is fixed at all iterations of the Algorithm 1. The corresponding model probability evolutions are shown in Figure 6. Figure 5 suggests that the online optimized dilution rate sequence exhibits significant fluctuations between the upper and lower bounds of the manipulated input, whereas the offline computed input sequence exhibits a smoother behavior. The online and offline input sequences are significantly different due to the fact that the model probabilities are changing during the fault diagnosis experiment. For example, the substrate inhibition model can be ruled out after the first few sampling times as its probability goes to zero (see Figure 6). This allows the online AFD method to direct its efforts at discrimination between the normal operation and substrate feed offset models, in contrast to considering all three models as in offline AFD. Figure 6 suggests that the correct model is identified faster in online AFD, while leading to higher probability of the true model than the case of offline AFD. In addition, the pulselike dilution rate sequence computed online is likely to be less intrusive to the process since less substrate is fed to the bioreactor (see Figure 5). It should be noted that the average CPU time required to solve the optimization problem at each measurement sampling time was approximately 5 seconds, which is significantly lower than the 12-min sampling interval between the measurements. 41

Online AFD in Practice We now demonstrate how the proposed online AFD method can be effectively applied in practice. Consider three distinct operational stages of a process: (i) normal process operation during which a performance loss can be detected using standard performance monitoring methods, (ii) online AFD to diagnose the current process operation mode, and (iii) recovery of 31



0.8 Online Offline

0.7

Dilution rate (hr-1 )

0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.5

1

1.5

2

Time (hr)

Figure 5: Comparison of the optimal input sequences for online AFD (blue) and offline AFD (red) computed under the same uncertainty realizations when the second fault scenario is active. 1 0.9 0.8

Model probability

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 32 of 40

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.5

1

1.5

2

Time (hr)

Figure 6: Evolution of the probabilities for the normal operation (black), substrate inhibition (red), and substrate concentration offset in the feed (blue) models under the optimal input sequences of online AFD (solid line) and offline AFD (dashed line) shown in Figure 5.

32


Page 33 of 40

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 process to desired conditions. For the bioreactor case study at hand, the three operational stages are illustrated in Figure 7, which shows the product concentration profile, normal and fault model probabilities evaluated online, and the dilution rate sequence. The first three hours of operation correspond to normal operation wherein the dilution rate is fixed and the product concentration remains relatively constant. Towards the end of the third hour of operation, the product concentration starts to decrease (Figure 7a), which indicates that a process change/fault might have occurred. Note that the probability of the normal operation model decreases during this period, but changes relatively slowly (Figure 7c). The proposed online AFD method is then executed during the next two hours of operation (from 3-5 hours). During this period, the dilution rate is computed by solving the AFD problem (24) to generate informative measurements that can aid identifying the model that corresponds to the correct mode of operation with high confidence (Figure 7d). The latter model can then be used to recover the operation. In this case study, the dilution rate is set to a value that is known to yield a steady-state product concentration at the setpoint of 27 g/L. In general, advanced control strategies, such as model predictive control, can be used in this stage to enable high-performance system operation. Twenty Monte Carlo runs were performed to evaluate the performance of online AFD under various uncertainty realizations (not shown here). In all of the Monte Carlo runs, the correct model was identified in only few sampling times, suggesting the robustness of the online AFD method to system uncertainties. In addition, the product concentration eventually returned to the desired setpoint in all runs. To evaluate the effectiveness of the proposed moment-based surrogate for chance constraint approximation, a chance constraint on the product concentration was included in the AFD problem. Even though up to 15% violation of the state constraint P > 23 g/L was allowed in the chance constraint, only a maximum of 10% constraint violation was observed in the Monte Carlo runs, suggesting a non-conservative chance constraint approximation.

33



P (g/l)

30

(a)

25

(b)

Model active

20

Constraint Set point

0

1

(d)

Model probability

(c)

2

3

3

4

5

6

7

8

4

5

6

7

8

Substrate offset

2 Normal operation

1 0

D (hr-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

Page 34 of 40

1

2

3

1 Normal operation Substr. inhibition Substrate offset

0.5 0

0

1

2

3

0

1

2

3

4

5

6

7

8

4

5

6

7

8

0.5

0

Time (hr)

Figure 7: One run of online AFD during system operation. Different segments of time correspond to different operational stages: 0-3 hours is normal operation (fault occurs at start of second hour), 3-5 hours is the fault diagnosis period during which online AFD is performed, and 5-8 hours is the performance recovery stage. The plots show (a) the product concentration, (b) the active model in different stages of operation, (c) the online model probabilities, and (d) the dilution rate profile.

Conclusions A computationally tractable method is presented for active fault diagnosis of nonlinear systems subject to stochastic model uncertainty, disturbances, and measurement noise. The proposed method aims to design input sequences that can optimally discriminate between uncertain models of normal and faulty operation in a probabilistic sense, while ensuring safe and least intrusive system operation during the fault diagnosis experiment through enforcement of input and state chance constraints. Tractable approximations for chance constraints and the Bayes risk for hypothesis selection are introduced to enable online solution of the active fault diagnosis problem. The efficacy of the proposed method for offline and online active fault diagnosis is demonstrated on a continuous bioreactor. The simulation results indicate that the designed input sequences achieve correct fault diagnosis in all the considered cases despite the presence of stochastic system uncertainty. The online implementation of the active fault diagnosis method, which relies on repeated design of the input sequences 34


Page 35 of 40

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


based on system observations, shows superior performance compared to that of the offline implementation. The simulation results also demonstrate that the proposed active fault diagnosis method can be seamlessly integrated into practical settings.

References (1) Nikoukhah, R. Guaranteed active failure detection and isolation for linear dynamical systems. Automatica 1998, 34, 1345–1358. (2) Patton, R. J.; Chen, J. Observer-based fault detection and isolation: Robustness and applications. Control Engineering Practice 1997, 5, 671–682. (3) Chiang, L. H.; Russell, E. L.; Braatz, R. D. Fault Detection and Diagnosis in Industrial Systems; Springer, 2000. (4) Ding, S. Model-based Fault Diagnosis Techniques: Design Schemes, Algorithms, and Tools; Springer, 2008. (5) Campbell, S. L.; Nikoukhah, R. Auxiliary Signal Design for Failure Detection; Princeton University Press, 2015. (6) Blanke, M.; Kinnaert, M.; Lunze, J.; Staroswiecki, M.; Schröder, J. Diagnosis and Fault-tolerant Control ; Springer, Berlin, 2006. (7) Sampath, M.; Lafortune, S.; Teneketzis, D. Active diagnosis of discrete-event systems. IEEE Transactions on Automatic Control 1998, 43, 908–929. (8) Ashari, A. E.; Nikoukhah, R.; Campbell, S. L. Auxiliary signal design for robust active fault detection of linear discrete-time systems. Automatica 2011, 47, 1887–1895. (9) Scott, J. K.; Findeisen, R.; Braatz, R. D.; Raimondo, D. M. Input design for guaranteed fault diagnosis using zonotopes. Automatica 2014, 50, 1580–1589. 35



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

(10) Scott, J. K.; Raimondo, D. M.; Marseglia, G. R.; Braatz, R. D. Constrained zonotopes: A new tool for set-based estimation and fault detection. Automatica 2016, 69, 126–136. (11) Raimondo, D. M.; Marseglia, G. R.; Braatz, R. D.; Scott, J. K. Closed-loop input design for guaranteed fault diagnosis using set-valued observers. Automatica 2016, 74, 107–117. (12) Andjelkovic, I.; Sweetingham, K.; Campbell, S. L. Active fault detection in nonlinear systems using auxiliary signals. Proceedings of the American Control Conference. 2008; pp 2142–2147. (13) Paulson, J. A.; Raimondo, D. M.; Findeisen, R.; Braatz, R. D.; Streif, S. Guaranteed active fault diagnosis for uncertain nonlinear systems. Proceedings of the European Control Conference. 2014; pp 926–931. (14) Zhang, X. J. Auxiliary Signal Design in Fault Detection and Diagnosis; Springer, 1989. (15) Blackmore, L.; Williams, B. Finite horizon control design for optimal model discrimination. Proceedings of the 45th IEEE Conference on Decision and Control and European Control Conference. 2005; pp 3795–3802. (16) Blackmore, L.; Williams, B. Finite horizon control design for optimal discrimination between several models. Proceedings of the 45th IEEE Conference on Decision and Control. 2006; pp 1147–1152. (17) Blackmore, L.; Rajamanoharan, S.; Williams, B. C. Active estimation for jump Markov linear systems. IEEE Transactions on Automatic Control 2008, 53, 2223–2236. (18) Šimandl, M.; Punčochář, I. Active fault detection and control: Unified formulation and optimal design. Automatica 2009, 45, 2052–2059. (19) Mirzaee, A.; Salahshoor, K. Fault diagnosis and accommodation of nonlinear systems

36


Page 36 of 40

Page 37 of 40

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


based on multiple-model adaptive unscented Kalman filter and switched MPC and H-infinity loop-shaping controller. Journal of Process Control 2012, 22, 626–634. (20) Mesbah, A.; Streif, S.; Findeisen, R.; Braatz, R. D. Active fault diagnosis for nonlinear systems with probabilistic uncertainties. Proceedings of the IFAC World Congress Conference. 2014; pp 7079–7084. (21) Mohajerin Esfahani, P.; Lygeros, J. A tractable fault detection and isolation approach for nonlinear systems with probabilistic performance. IEEE Transactions on Automatic Control 2016, 61, 633–647. (22) Xiu, D.; Karniadakis, G. E. The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM Journal of Scientific Computation 2002, 24, 619–644. (23) Ljung, L. Asymptotic behavior of the extended Kalman filter as a parameter estimator for linear systems. IEEE Transactions on Automatic Control 1979, 24, 36–50. (24) Wan, E. A.; van der Merwe, R. The unscented Kalman filter for nonlinear estimation. Proceedings of the IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium. 2000; pp 153–158. (25) Kailath, T. The divergence and Bhattacharyya distance measures in signal selection. IEEE Transactions on Communication Technology 1967, 15, 52–60. (26) Cover, T.; Hart, P. Nearest neighbor pattern classification. IEEE Transactions on Information Theory 1967, 13, 21–27. (27) Cowles, M. K.; Carlin, B. P. Markov chain Monte Carlo convergence diagnostics: A comparative review. Journal of the American Statistical Association 1996, 91, 883–904. (28) Calafiore, G. C.; El Ghaoui, L. On distributionally robust chance-constrained linear programs. Journal of Optimization Theory and Applications 2006, 130, 1–22.

37



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

(29) Nemirovski, A.; Shapiro, A. Convex approximations of chance constrained programs. SIAM Journal on Optimization 2006, 17, 969–996. (30) Streif, S.; Petzke, F.; Mesbah, A.; Findeisen, R.; Braatz, R. D. Optimal experimental design for probabilistic model discrimination using polynomial chaos. Proceedings of the IFAC World Congress Conference. 2014; pp 4103–4109. (31) Mesbah, A.; Streif, S. A probabilistic approach to robust optimal experiment design with chance constraints. Proceedings of the 9th International Symposium on Advanced Control of Chemical Processes (ADCHEM). 2015; pp 100–105. (32) Chen, Z. Bayesian filtering: From Kalman filters to particle filters, and beyond. Statistics 2003, 182, 1–69. (33) Caflisch, R. E. Monte carlo and quasi-monte carlo methods. Acta Numerica 1998, 7, 1–49. (34) Goodwin, G. C.; Payne, R. L. Dynamic system identification: Experiment design and data analysis. Mathematics in Science and Engineering 1977, 136 . (35) Pronzato, L. Optimal experimental design and some related control problems. Automatica 2008, 44, 303–325. (36) van Hessem, D. H.; Scherer, C. W.; Bosgra, O. H. LMI-based closed-loop economic optimization of stochastic process operation under state and input constraints. Proceedings of the 40th IEEE Conference on Decision and Control. 2001; pp 4228–4233. (37) Mehrotra, S. On the implementation of a primal-dual interior point method. SIAM Journal on Optimization 1992, 2, 575–601. (38) Nocedal, J.; Wright, S. J. Sequential Quadratic Programming; Springer, 2006.

38


Page 38 of 40

Page 39 of 40

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


(39) Gordon, N. J.; Salmond, D. J.; Smith, A. F. M. Novel approach to nonlinear/nonGaussian Bayesian state estimation. IEE Proceedings on Radar and Signal Processing. 1993; pp 107–113. (40) Agrawal, P.; Koshy, G.; Ramseier, M. An algorithm for operating a fed-batch fermentor at optimum specific-growth rate. Biotechnology and Bioengineering 1989, 33, 115–125. (41) Henson, M. A.; Seborg, D. E. Nonlinear control strategies for continuous fermenters. Chemical Engineering Science 1992, 47, 821–835. (42) Andersson, J.; Åkesson, J.; Diehl, M. Recent Advances in Algorithmic Differentiation; Springer, 2012; pp 297–307. (43) MacGregor, J.; Cinar, A. Monitoring, fault diagnosis, fault-tolerant control and optimization: Data driven methods. Computers & Chemical Engineering 2012, 47, 111– 120. (44) Zagrobelny, M.; Ji, L.; Rawlings, J. B. Quis custodiet ipsos custodes? Annual Reviews in Control 2013, 37, 260–270.

39




Page 40 of 40

Input Design for Online Fault Diagnosis of Nonlinear Systems with

Recommend Documents