# Semiadaptive Fault Diagnosis via Variational Bayesian Mixture Factor Analysis with Application to Wastewater Treatment

Authors: Hongjun Xiao; Yiqi Liu; Daoping Huang
Journal: Journal of Control Science and Engineering
 (2016)

---

## Abstract

Mainly due to the hostile environment in wastewater plants (WWTPs), the reliability of sensors with respect to important qualities is often poor. In this work, we present the design of a semiadaptive fault diagnosis method based on the variational Bayesian mixture factor analysis (VBMFA) to support process monitoring. The proposed method is capable of capturing strong nonlinearity and the significant dynamic feature of WWTPs that seriously limit the application of conventional multivariate statistical methods for fault diagnosis implementation. The performance of proposed method is validated through a simulation study of a wastewater plant. Results have demonstrated that the proposed strategy can significantly improve the ability of fault diagnosis under fault-free scenario, accurately detect the abrupt change and drift fault, and even localize the root cause of corresponding fault properly.


---

## Body

## 1. Introduction

Due to stringent environmental regulations, wastewater treatment plants (WWTPs) are always challenged to meet new effluent quality requirements [1]. In such context, the use of online sensors for control and optimization in the wastewater treatment plants is indispensable. Measurements from instrumentation shall be available 24 hours a day and 7 days a week. However, even reliable instrumentation can fail during operation, which can have serious consequences if the instrumentation is used in the closed form. The huge amount of correlated data adds more complexity to this problem, thus making manual expert-based quality evaluation infeasible [2]. The development of methods that allow us to extract useful information for automatic detection of faults and identifying their root cause is urgently needed.This vision has motivated inspiring efforts toward using multivariate statistical process monitoring (MSPM) methods, such as principal component analysis (PCA) [3, 4] and partial least-squares (PLS) [5, 6], for fault diagnosis in the wastewater treatment process, even though they had been used for many years in the chemical process industries [3, 7, 8]. The goal of such methods is to project a high dimensional measurement space onto a space with significantly fewer dimensions. In PCA, Hotelling statistics, T2, and the sum of squared residuals, SPE, or Q statistic are formulated as tools to identify the deviations of normal process. The T2 statistic is a measure of the variation in the PCA model, whereas the Q statistic is a measure of the amount of variation not captured by the PCA model [9]. However, these methods adhere to an assumption that process noises for each variable are normally distributed. The presence of variations such as shifting operational conditions, feed flow fluctuations, control strategy changes, and non-Gaussian disturbances in the WWTPs often deviates from such an assumption. Independent component analysis (ICA) provides an alternative to make data more comprehensible and useful for non-Gaussian fault diagnosis [10–12]. Other complementary MSPM approaches include Fisher discriminant analysis (FDA) [13], canonical variable analysis (CVA) [14], and support vector machines (SVM) [15]. Nevertheless, WWTPs pose an additional challenge to the development and application of these methods. The problem is highly uncertain environment. To alleviate such problem, a probabilistic PCA based method was proposed with the assumption that both of process variables and noises follow Gaussian distribution [16, 17]. The issue of PPCA is that it assumes the same noise level for all measurement variables with potentially underfitting the fault diagnosis model. Factor analysis provides a potential way to address such a problem using different noise variance for each variable [18]. Unfortunately, parameter identification of FA using the least-squares method or maximum likelihood for optimization potentially results in an overfitting problem. Therefore, variational Bayesian learning was introduced to identify corresponding parameters for FA or PCA and further used to process monitoring [19]. These methods work reasonably well with synthetic data sets but they have difficulties in handling wastewater plant dynamics (e.g., seasonal or diurnal changes) and nonlinearity (e.g., temperature-dependent kinetics). One of common and plausible ways is to extend the conventional methods to a mixture form with the capability of approaching global nonlinearity by some local linearity models. However, there are several drawbacks of the mixture model based methods [20]: (1) they cannot determine the effective factor number for each factor analyzer automatically through the modeling process; (2) the importance of each retained factor in the model structure cannot be easily differentiated; (3) singularity and overfitting problems maybe caused when the number of modeling data samples is limited. Variational Bayesian mixture factor analysis offers a good alternative to overcome these problems. Additionally, the VBMFA method is not only able to determine the dimensionality automatically, but also capable of avoiding overfitting problem [21].It is common sense that batch processes are in fact a special type of nonlinear processes. Therefore, in this paper, unlike other mixture models for batch multimode process monitoring in general [22–25], the VBMFA-based process monitoring method is proposed to address the significant nonlinear problems in the WWTPs. This is more practical because systems are always nonclassifiable clearly for most cases. By doing so, the corresponding T2 and SPE statistics are formulated as tools for fault diagnosis and identification. T2 and SPE statistics from local factor analyzers are assimilated by weights, respectively. In order to account for changing process behavior, such weights can change over time, thereby offering more reliable way for fault diagnosis and identification. Therefore, the implementation of adaptive mechanism for the proposed method is not by their mean and variance but rather depending on their weight update. In other words, the resulting online VBMFA model is obtained by weighting adaptively several factor analyzers that are already trained offline. Therefore, the resulting method is termed as a semiadaptive fault diagnosis method. In summary, the contributions of this paper is given as follows: (1) typical variables are generated to account for the fluctuations of diurnal wastewater changes in workdays and weekends, serving as data pretreatment and an initial step to break up dynamics of the wastewater process; (2) the VBMFA is introduced for fault diagnosis; (3) an adaptive probabilistic weighted strategy is adapted to mix different local factor analyzers; (4) an adaptive control limit is proposed to monitor a WWTP with strong nonlinearity and dynamic feature.The rest of this paper is organized as follows. In Section2, an overview of issues in the WWTPs is given, together with a description of the simulation systems involved in this work. Section 3 will deal with the description of the FA and VBMFA. Section 4 provides a detailed demonstration for the development of the VBMFA-based fault diagnosis strategy, which is followed by the case study of the corresponding process monitoring scheme in the next section. Finally, some conclusions are made.


## 2. A Simulated Platform for Process Monitoring

### 2.1. WWTP

A wastewater treatment plant is designed to degrade wastewater and return clean water to the river. One of the most popular processes available for dealing with wastewater is the activated sludge. However, the complexity of the physical, chemical, and biological phenomena associated with treatment units always frustrates the reliability of WWTPs dramatically. Moreover, the reactors exhibit strong nonlinear dynamics and have coupling effects among the variables. More care should also be taken for fault diagnosis model construction; more details are shown as follows:(1)
Modeling Challenges for Process Monitor. The overall WWTP model consists of the hydraulic model and the biochemical model. The hydraulic model is related to reactor behavior, flow rates, and recirculation, whereas the biochemical model is mainly associated with the activated sludge and describes microbial growth, death, and nutrient consumption. Such huge number of biological processes discourage the process model construction.(2)
Uncertainties. The main goal of a WWTP is to maximize the reduction of wastewater while meeting the environmental requirement from government. However, uncertainties such as influent flow rate, water quality, and contamination loads add more complexity for process monitor model building.
(a)
Uncertainty about the Initial Conditions. The large variety and complexity of all present biological species in WWTPs make it identify the role of each one in a process. Moreover, even if it is possible to determine a total initial biomass concentration, this does not generally correspond to the concentration of the effectively active species.(b)
Uncertainty in Process Kinetics. Again, the large variety complexity of biological species presented in WWTPs makes it impossible to determine all kinetic parameters for all biological reactions. Moreover, bacteria can be sensitive to environmental changes and even if the kinetics could be globally described, their parameters vary with time as a function of the time-varying distribution of the biological population.(c)
Uncertainty in the Input Contractions. In contrast with food and pharmaceutical biological processes [26, 27], where feed concentrations are usually well known, the inputs of a WWTP have to be considered as disturbances. Input concentrations depend on several random factors. For instance, in municipal WWTPs, the organic load varies during the day according to the level of human activity. Also, they themselves are strongly affected by weather conditions and seasonal change. Such fluctuations often result in the degradation of the performance or even plant failure.


### 2.2. BSM1 (Benchmark Simulation Model 1) Platform

The simulation platform used is the BSM1 developed by the second IWA Respirometry Task Group, offering an unbiased benchmarking system for comparing various control strategies (http://www.benchmarkwwtp.org/) [28]. In this study, a relatively simple plant layout, which is most commonly used for removal of organic matter, nitrification, and denitrification, is selected in the simulation benchmark (see Figure 1). The plant consists of five compartment biological tanks (5999 m3) and a secondary settler (6000 m3), designing to treat an average flow of 20,000 m3/day with an average biodegradable chemical oxygen demand (COD) concentration of 300 mg/L. The last three compartments of the bioreactor are aerated, whereas the others are not. The nonreactive secondary settler is modeled with a series of 10 layers with an area of 1500 m2 and a depth of 4 m. There are two internal recycles: the nitrate internal recycle from the last tank to the first tank and the activated sludge recycle from the underflow of the secondary settler to the front end of the plant. The feed from the final biological reactor is connected to the sixth layer from the bottom of the secondary settler. The scenario with dry weather condition is tested in this case study. The simulated BSM1 comprises 14 days of influent data recorded at 15-minute intervals from a real plant. In order to make sure that the methodology is able to translate to a real WWTP easily in the future, the BSM1 model and proposed method were run in separate computer and connected by two I/O cards.Figure 1
The schematic of wastewater plant for BSM1.


## 2.1. WWTP

A wastewater treatment plant is designed to degrade wastewater and return clean water to the river. One of the most popular processes available for dealing with wastewater is the activated sludge. However, the complexity of the physical, chemical, and biological phenomena associated with treatment units always frustrates the reliability of WWTPs dramatically. Moreover, the reactors exhibit strong nonlinear dynamics and have coupling effects among the variables. More care should also be taken for fault diagnosis model construction; more details are shown as follows:(1)
Modeling Challenges for Process Monitor. The overall WWTP model consists of the hydraulic model and the biochemical model. The hydraulic model is related to reactor behavior, flow rates, and recirculation, whereas the biochemical model is mainly associated with the activated sludge and describes microbial growth, death, and nutrient consumption. Such huge number of biological processes discourage the process model construction.(2)
Uncertainties. The main goal of a WWTP is to maximize the reduction of wastewater while meeting the environmental requirement from government. However, uncertainties such as influent flow rate, water quality, and contamination loads add more complexity for process monitor model building.
(a)
Uncertainty about the Initial Conditions. The large variety and complexity of all present biological species in WWTPs make it identify the role of each one in a process. Moreover, even if it is possible to determine a total initial biomass concentration, this does not generally correspond to the concentration of the effectively active species.(b)
Uncertainty in Process Kinetics. Again, the large variety complexity of biological species presented in WWTPs makes it impossible to determine all kinetic parameters for all biological reactions. Moreover, bacteria can be sensitive to environmental changes and even if the kinetics could be globally described, their parameters vary with time as a function of the time-varying distribution of the biological population.(c)
Uncertainty in the Input Contractions. In contrast with food and pharmaceutical biological processes [26, 27], where feed concentrations are usually well known, the inputs of a WWTP have to be considered as disturbances. Input concentrations depend on several random factors. For instance, in municipal WWTPs, the organic load varies during the day according to the level of human activity. Also, they themselves are strongly affected by weather conditions and seasonal change. Such fluctuations often result in the degradation of the performance or even plant failure.


## 2.2. BSM1 (Benchmark Simulation Model 1) Platform

The simulation platform used is the BSM1 developed by the second IWA Respirometry Task Group, offering an unbiased benchmarking system for comparing various control strategies (http://www.benchmarkwwtp.org/) [28]. In this study, a relatively simple plant layout, which is most commonly used for removal of organic matter, nitrification, and denitrification, is selected in the simulation benchmark (see Figure 1). The plant consists of five compartment biological tanks (5999 m3) and a secondary settler (6000 m3), designing to treat an average flow of 20,000 m3/day with an average biodegradable chemical oxygen demand (COD) concentration of 300 mg/L. The last three compartments of the bioreactor are aerated, whereas the others are not. The nonreactive secondary settler is modeled with a series of 10 layers with an area of 1500 m2 and a depth of 4 m. There are two internal recycles: the nitrate internal recycle from the last tank to the first tank and the activated sludge recycle from the underflow of the secondary settler to the front end of the plant. The feed from the final biological reactor is connected to the sixth layer from the bottom of the secondary settler. The scenario with dry weather condition is tested in this case study. The simulated BSM1 comprises 14 days of influent data recorded at 15-minute intervals from a real plant. In order to make sure that the methodology is able to translate to a real WWTP easily in the future, the BSM1 model and proposed method were run in separate computer and connected by two I/O cards.Figure 1
The schematic of wastewater plant for BSM1.


## 3. Factor Analysis and Mixture Factor Analysis Model Based on Variational Bayesian Learning

### 3.1. Factor Analysis

Factor analysis is a statistical method with the capability to describe variability among observed correlated variables in terms of a potentially lower number of unobserved variables called factors. The original observed measurement variablesx∈Rm are treated as linear combination of factors t plus the noise e∈Rm and mean vector μ∈Rm; that is, (1)x=Λt+μ+e,where Λ∈R(m×k) is the linear transformation known as the factor loading matrix and e is the zero-mean Gaussian noise with diagonal covariance matrix Ψ. Given the observed variables x∈Rm, the parameters θ={μ,Λ,Ψ} in the model are needed to be identified. One of common ways is to use the expectation and maximization (EM) algorithm, also termed as the iterative likelihood maximization algorithm.


### 3.2. Mixture Factor Analysis Model Based on Variational Bayesian Learning (VBMFA)

It is a common sense that high dimensional data in fact lie on a low dimensional manifold generally [29]. In factor analysis, each factor dictates the amount of each linear transformation on the observed variables. However, such transformation can only well explain a small region of the manifold in which it is locally linear, even though the manifold is globally nonlinear.One way to deal with this problem is to use mixture models to tile the data manifold. A mixture of factor analyzers models the density for a data pointxi as a weighted average of factor analyzer densities as follows:(2)pxi∣π,Λ,μ,Ψ=∑siSpsi∣πpxi∣si,Λ,μ,Ψ,where S is the number of mixture components in the model, π is the vector of mixing proportions, si is the discrete indicator variable for the mixture component selected for modeling the data point i, and Λ={Λs}s=1S is a set of factor loadings with the factor matrix Λs for analyzer s.Due to the complexity of integrating out all those parameters, it is necessary to choose priors that the likelihood terms greatly simplify inference and interpretability. Therefore, the Dirichlet prior for mixing proportionπ, with strength α∗, is chosen as follows:(3)pπ∣α∗m∗=Dirπ∣α∗m∗with m∗=[1/S,…,1/S]. Thus, only a single hyperparameter is needed regardless of the dimensionality of π. The column of each factor loading matrix has a Gaussian prior with zero mean and a different precision parameter (drawn from a gamma distribution with fixed hyperparameters):(4)pΛ∣v=∏s=1SpΛs∣vs=∏s=1S∏j=1kspΛ.js∣vjs=∏s=1S∏j=1kspΛ.js∣0,Ivjs,where Λ.js denotes the vector of entries in the jth column of the sth analyzer in the mixture and vjs is the same scalar precision for each entry in the corresponding column.Since the number of hyperparameters inv={{vjs}j=1ks} increases with the number of analyzers, a hyperprior on every element of each vs precision vector is as follows:(5)pv∣a∗,b∗=∏s=1S∏j=1ksGavjs∣a∗,b∗,where a∗ and b∗ are shape and inverse-scale hyperparameters for a gamma distribution. Finally, the means of each analyzer in the mixture need to be integrated out. A Gaussian prior with mean μ∗ and axis-aligned precision diag⁡(v∗) is placed on each mean μ∗:(6)pμ∣μ∗,v∗=∏s=1SNμs∣μ∗,diag⁡v∗-1⁡.The directed acyclic graph for the generative model for VBMFA is shown graphically in Figure 2.Figure 2
The formulation for VBMFA.The aim of VBMFA is to optimize(7)L=ln⁡x=ln∫dπpπ∣α∗m∗∫dvpv∣a∗,b∗∫ΛpΛ∣v∫μpμ∣μ∗,v∗∏i=1N∏si=1Spsi∣π∫dtiptipxi∣si,Λ,μ,Ψ.An arbitrary distribution q(π,v,Λ,μ) is introduced to lower bound. By using the Kullback-Leibler divergence [30], L can be reformulated as follows: (8)L⩾∫dπln⁡pπ∣α∗m∗qπ+∑s=1S∫dvsqvsln⁡pvs∣a∗,b∗qvs+∫dΛ~sqΛ~sln⁡pΛ~s∣a∗,b∗qΛ~s+∑i=1N∑siSqsi∫dπln⁡psi∣πqsi+∫dtiqti∣siln⁡ptiqti∣si+∫Λ~sqΛ~s∫dtiqti∣siln⁡pxi∣si,ti,Λ~,Ψ=Fqθ,qs,t,Θ,where Θ=(α∗m∗,a∗,b∗,μ∗,v∗,Ψ). Thus, the lower bound is a sum of a foundational of variational posterior distribution over the parameters, denoted by qθ, and a foundational of variational posterior distribution over the hidden variables of each data point.In order to optimize the lower bound, functional derivatives with respect to eachq· distribution are taken and are equal to zero to find the distributions. By doing so, the corresponding q· can be obtained as follows:(9)qπ=Dirπ∣αm,where each element of the variational parameter αm is given by (10)αms=α∗ms∗+∑i=1Nqsiwith α=α∗+n, ms∗=1/S, and ∑s=1Sms=1:(11)qvls=Gavls∣a∗+p2,b∗+12∑j=1pΛjlsqΛs,qΛ~s=∏j=1pqΛ~js=∏j=1pNΛ~js∣Λ~-js,Γ~js,where Λ~js denotes the column vector associated with the jth row of Λ~s, which has ks+1 dimensions. These variational poster parameters are given by(12)Γ~js=ΣΛΛj,s-1ΣΛμj,s-1ΣμΛj,s-1Σμμj,s-1-1,Λ~-js=Λ-jsμ-jswith(13)ΣΛΛj,s-1=diagvsqvs+Ψjj-1∑i=1NqsititiTqti∣si,Σμμj,s-1=vj∗+Ψjj-1∑i=1Nqsi,ΣΛμj,s-1=Ψjj-1∑i=1Nqsitiqti∣si=ΣμΛj,s-1,Λ-js=Γ~jsΛΛΨjj-1∑i=1Nqsixijtiqti∣si,μ-js=Γ~jsμμΨjj-1∑i=1Nqsixij+vj∗μj∗,qti∣s=Nti∣t-is,Σswith Σs-1=I+ΛsTΨ-1Λsq(Λ-s), t-is=ΣsΛsTΨ-1(xi-μs)q(Λ-s).By using a result of Dirichlet distributions in [21], qsi is computed by(14)qsi=1Ziexp⁡Ψαmsi-Ψα+12ln⁡Σs-12trΨ-1xi-Λ-siti1xi-Λ-siti1TqΛ-siqti∣si,where Zi is a normalization constant for each data point, such that ∑si=1Sqsi=1:(15)Ψ-1=diag⁡1N∑i-1Nxi-Λ-siti1xi-Λ-siti1TqΛ-sqsiqti∣si,ψa∗=ln⁡b∗+1Sk∑s=1S∑j=1kln⁡vlsqvls,b∗-1=1α∗Sk∑s=1S∑j=1kvlsqvls,ψα∗-ψα∗S=1S∑s=1Sψα-ψαms.The solution for μ∗ and v∗ can be found by maximizing (8):(16)μ∗=1S∑s=1Sμsqμs,v∗=v1∗,…,vp∗with(17)vj∗=1S∑s=1Sμjs-μj∗μjs-μj∗qμs,where the update for v∗ uses the already updated μ∗.


## 3.1. Factor Analysis

Factor analysis is a statistical method with the capability to describe variability among observed correlated variables in terms of a potentially lower number of unobserved variables called factors. The original observed measurement variablesx∈Rm are treated as linear combination of factors t plus the noise e∈Rm and mean vector μ∈Rm; that is, (1)x=Λt+μ+e,where Λ∈R(m×k) is the linear transformation known as the factor loading matrix and e is the zero-mean Gaussian noise with diagonal covariance matrix Ψ. Given the observed variables x∈Rm, the parameters θ={μ,Λ,Ψ} in the model are needed to be identified. One of common ways is to use the expectation and maximization (EM) algorithm, also termed as the iterative likelihood maximization algorithm.


## 3.2. Mixture Factor Analysis Model Based on Variational Bayesian Learning (VBMFA)

It is a common sense that high dimensional data in fact lie on a low dimensional manifold generally [29]. In factor analysis, each factor dictates the amount of each linear transformation on the observed variables. However, such transformation can only well explain a small region of the manifold in which it is locally linear, even though the manifold is globally nonlinear.One way to deal with this problem is to use mixture models to tile the data manifold. A mixture of factor analyzers models the density for a data pointxi as a weighted average of factor analyzer densities as follows:(2)pxi∣π,Λ,μ,Ψ=∑siSpsi∣πpxi∣si,Λ,μ,Ψ,where S is the number of mixture components in the model, π is the vector of mixing proportions, si is the discrete indicator variable for the mixture component selected for modeling the data point i, and Λ={Λs}s=1S is a set of factor loadings with the factor matrix Λs for analyzer s.Due to the complexity of integrating out all those parameters, it is necessary to choose priors that the likelihood terms greatly simplify inference and interpretability. Therefore, the Dirichlet prior for mixing proportionπ, with strength α∗, is chosen as follows:(3)pπ∣α∗m∗=Dirπ∣α∗m∗with m∗=[1/S,…,1/S]. Thus, only a single hyperparameter is needed regardless of the dimensionality of π. The column of each factor loading matrix has a Gaussian prior with zero mean and a different precision parameter (drawn from a gamma distribution with fixed hyperparameters):(4)pΛ∣v=∏s=1SpΛs∣vs=∏s=1S∏j=1kspΛ.js∣vjs=∏s=1S∏j=1kspΛ.js∣0,Ivjs,where Λ.js denotes the vector of entries in the jth column of the sth analyzer in the mixture and vjs is the same scalar precision for each entry in the corresponding column.Since the number of hyperparameters inv={{vjs}j=1ks} increases with the number of analyzers, a hyperprior on every element of each vs precision vector is as follows:(5)pv∣a∗,b∗=∏s=1S∏j=1ksGavjs∣a∗,b∗,where a∗ and b∗ are shape and inverse-scale hyperparameters for a gamma distribution. Finally, the means of each analyzer in the mixture need to be integrated out. A Gaussian prior with mean μ∗ and axis-aligned precision diag⁡(v∗) is placed on each mean μ∗:(6)pμ∣μ∗,v∗=∏s=1SNμs∣μ∗,diag⁡v∗-1⁡.The directed acyclic graph for the generative model for VBMFA is shown graphically in Figure 2.Figure 2
The formulation for VBMFA.The aim of VBMFA is to optimize(7)L=ln⁡x=ln∫dπpπ∣α∗m∗∫dvpv∣a∗,b∗∫ΛpΛ∣v∫μpμ∣μ∗,v∗∏i=1N∏si=1Spsi∣π∫dtiptipxi∣si,Λ,μ,Ψ.An arbitrary distribution q(π,v,Λ,μ) is introduced to lower bound. By using the Kullback-Leibler divergence [30], L can be reformulated as follows: (8)L⩾∫dπln⁡pπ∣α∗m∗qπ+∑s=1S∫dvsqvsln⁡pvs∣a∗,b∗qvs+∫dΛ~sqΛ~sln⁡pΛ~s∣a∗,b∗qΛ~s+∑i=1N∑siSqsi∫dπln⁡psi∣πqsi+∫dtiqti∣siln⁡ptiqti∣si+∫Λ~sqΛ~s∫dtiqti∣siln⁡pxi∣si,ti,Λ~,Ψ=Fqθ,qs,t,Θ,where Θ=(α∗m∗,a∗,b∗,μ∗,v∗,Ψ). Thus, the lower bound is a sum of a foundational of variational posterior distribution over the parameters, denoted by qθ, and a foundational of variational posterior distribution over the hidden variables of each data point.In order to optimize the lower bound, functional derivatives with respect to eachq· distribution are taken and are equal to zero to find the distributions. By doing so, the corresponding q· can be obtained as follows:(9)qπ=Dirπ∣αm,where each element of the variational parameter αm is given by (10)αms=α∗ms∗+∑i=1Nqsiwith α=α∗+n, ms∗=1/S, and ∑s=1Sms=1:(11)qvls=Gavls∣a∗+p2,b∗+12∑j=1pΛjlsqΛs,qΛ~s=∏j=1pqΛ~js=∏j=1pNΛ~js∣Λ~-js,Γ~js,where Λ~js denotes the column vector associated with the jth row of Λ~s, which has ks+1 dimensions. These variational poster parameters are given by(12)Γ~js=ΣΛΛj,s-1ΣΛμj,s-1ΣμΛj,s-1Σμμj,s-1-1,Λ~-js=Λ-jsμ-jswith(13)ΣΛΛj,s-1=diagvsqvs+Ψjj-1∑i=1NqsititiTqti∣si,Σμμj,s-1=vj∗+Ψjj-1∑i=1Nqsi,ΣΛμj,s-1=Ψjj-1∑i=1Nqsitiqti∣si=ΣμΛj,s-1,Λ-js=Γ~jsΛΛΨjj-1∑i=1Nqsixijtiqti∣si,μ-js=Γ~jsμμΨjj-1∑i=1Nqsixij+vj∗μj∗,qti∣s=Nti∣t-is,Σswith Σs-1=I+ΛsTΨ-1Λsq(Λ-s), t-is=ΣsΛsTΨ-1(xi-μs)q(Λ-s).By using a result of Dirichlet distributions in [21], qsi is computed by(14)qsi=1Ziexp⁡Ψαmsi-Ψα+12ln⁡Σs-12trΨ-1xi-Λ-siti1xi-Λ-siti1TqΛ-siqti∣si,where Zi is a normalization constant for each data point, such that ∑si=1Sqsi=1:(15)Ψ-1=diag⁡1N∑i-1Nxi-Λ-siti1xi-Λ-siti1TqΛ-sqsiqti∣si,ψa∗=ln⁡b∗+1Sk∑s=1S∑j=1kln⁡vlsqvls,b∗-1=1α∗Sk∑s=1S∑j=1kvlsqvls,ψα∗-ψα∗S=1S∑s=1Sψα-ψαms.The solution for μ∗ and v∗ can be found by maximizing (8):(16)μ∗=1S∑s=1Sμsqμs,v∗=v1∗,…,vp∗with(17)vj∗=1S∑s=1Sμjs-μj∗μjs-μj∗qμs,where the update for v∗ uses the already updated μ∗.


## 4. Units Process Monitoring Based on VBMFA

### 4.1. Fault Detection

In order to perform process monitoring, monitoring statisticT2 and SPE is constructed first of all. Similar to PCA method, corresponding monitoring statistics for each local FA can be represented as follows:(18)Ts2=t^sTvar⁡t^s-1t^swith(19)t^s=Ets∣x,s,θ=ΣsΛsTΨ-1x-μsqΛ-s,var⁡t^s=Σs=I+ΛsTΨ-1ΛsqΛ-s-1.Similarly, the associated SPE statistic can be given as follows:(20)e^s=x-μs-Λst^s,SPEs=e^sTΨ-1e^s.The control limits of the T2 and SPE statistic can be both determined by the χ2 distribution; thus,(21)Ts2≤Tlim2=χα2ks,SPEs≤SPElim=χα2p,where s=1,2,…,S. k is the number of retained factors in each local FA model, p is the number of process variables, and α is the significance level of both monitoring statists. As for a new coming sample data xnew, both of statists are updated as follows:(22)t^news=Etnew∣xnew,s,θ=ΣsΛsΨ-1xnew-μsqΛ-s,Ts,new2=t^newsTvar⁡t^s-1t^news,e^news=xnew-μs-Λst^news,SPEs,new=e^newsTΨ-1e^news.However, it is cumbersome to monitor the process with all of the local monitoring charts, which also potentially frustrates the decision making for fault diagnosis with such many monitoring charts. Therefore, coordination of all the local monitoring charts is imperative. Different from the hard assignment method employed by [31], an adaptive method is implemented in this section. As the monitoring results of the new sample xnew in each local subspace have been identified, we can combine them through the estimated posterior probabilities, which are given as follows:(23)qsnew=1Znewexp⁡Ψαms-Ψα+12ln⁡Σs-12trΨ-1xnew-Λ-sts1xnew-Λ-sts1TqΛ-sqt∣swhich represents the responsibility of the components for the new coming data point xnew. Both of synthetic monitor statists T2 and SPE can be obtained as follows:(24)STnew2=∑s=1SqsnewTs,new2,SSPEnew=∑s=1SqsnewSPEs,new.Since the control limit of the local T2 is correlated with the number of factors ks in each FA model, there is a need to combine the control limit as well. Thus, qsnew is used to weight associated T2 control limit as follows:(25)T2≤Tlim2=∑s=1Sqsnewχα2ks.Note that, due to constant control limit of each local SPE χα2(p), it is not necessary to calculate the combined SPE control limit any more.Even though the structure for each local FA model is unchanged after model training, our process monitoring model can be adjusted by updating mixture proportions in an adaptive way. Hence, process behaviors can be judged by simply checking both of statistics accordingly if they exceed predefined control limits.


### 4.2. Fault Identification

Once a fault is detected, more information is necessarily obtained to find the root cause of the corresponding faults. A contribution plot is a common way to fill this gap, which is able to show the contribution of each process variable to the statistic calculated. A high contribution of a process variable usually indicates a problem with this specific variable. Plots of the contribution to the associated statistic are similar to standard plots of squared residuals obtained, but the contribution plots presented here also have control limits. These control limits are used for comparing the residuals of the new values for the present to the residuals of the fault-free data. If in the fault-free data a certain process variable had high residuals, it can also be expected to have high residuals in the present. However, if there are high residuals for a certain process variable that had low residuals in the fault-free data; this probably is due to a special event in the new process. If the contributions of a large group of process variables are studied, it is usually found that several process variables have high contributions. Using the control limits, it is easier to find those process variables that are really different, compared to the fault-free data. In this paper, a control limit for the contribution for each process variable is calculated as the mean of the contributions plus three times the standard deviation of the contributions for each process variable at each time period. That means contribution of each variable follows Gaussian distribution in a long period, but the matching mean and variance will be different. Control limits for contributions are studied by many researchers [32, 33]. Here, the control limits are determined according to the 3σ method [32]. The reason why we used the 3σ method is because of its efficiency and simplicity.


## 4.1. Fault Detection

In order to perform process monitoring, monitoring statisticT2 and SPE is constructed first of all. Similar to PCA method, corresponding monitoring statistics for each local FA can be represented as follows:(18)Ts2=t^sTvar⁡t^s-1t^swith(19)t^s=Ets∣x,s,θ=ΣsΛsTΨ-1x-μsqΛ-s,var⁡t^s=Σs=I+ΛsTΨ-1ΛsqΛ-s-1.Similarly, the associated SPE statistic can be given as follows:(20)e^s=x-μs-Λst^s,SPEs=e^sTΨ-1e^s.The control limits of the T2 and SPE statistic can be both determined by the χ2 distribution; thus,(21)Ts2≤Tlim2=χα2ks,SPEs≤SPElim=χα2p,where s=1,2,…,S. k is the number of retained factors in each local FA model, p is the number of process variables, and α is the significance level of both monitoring statists. As for a new coming sample data xnew, both of statists are updated as follows:(22)t^news=Etnew∣xnew,s,θ=ΣsΛsΨ-1xnew-μsqΛ-s,Ts,new2=t^newsTvar⁡t^s-1t^news,e^news=xnew-μs-Λst^news,SPEs,new=e^newsTΨ-1e^news.However, it is cumbersome to monitor the process with all of the local monitoring charts, which also potentially frustrates the decision making for fault diagnosis with such many monitoring charts. Therefore, coordination of all the local monitoring charts is imperative. Different from the hard assignment method employed by [31], an adaptive method is implemented in this section. As the monitoring results of the new sample xnew in each local subspace have been identified, we can combine them through the estimated posterior probabilities, which are given as follows:(23)qsnew=1Znewexp⁡Ψαms-Ψα+12ln⁡Σs-12trΨ-1xnew-Λ-sts1xnew-Λ-sts1TqΛ-sqt∣swhich represents the responsibility of the components for the new coming data point xnew. Both of synthetic monitor statists T2 and SPE can be obtained as follows:(24)STnew2=∑s=1SqsnewTs,new2,SSPEnew=∑s=1SqsnewSPEs,new.Since the control limit of the local T2 is correlated with the number of factors ks in each FA model, there is a need to combine the control limit as well. Thus, qsnew is used to weight associated T2 control limit as follows:(25)T2≤Tlim2=∑s=1Sqsnewχα2ks.Note that, due to constant control limit of each local SPE χα2(p), it is not necessary to calculate the combined SPE control limit any more.Even though the structure for each local FA model is unchanged after model training, our process monitoring model can be adjusted by updating mixture proportions in an adaptive way. Hence, process behaviors can be judged by simply checking both of statistics accordingly if they exceed predefined control limits.


## 4.2. Fault Identification

Once a fault is detected, more information is necessarily obtained to find the root cause of the corresponding faults. A contribution plot is a common way to fill this gap, which is able to show the contribution of each process variable to the statistic calculated. A high contribution of a process variable usually indicates a problem with this specific variable. Plots of the contribution to the associated statistic are similar to standard plots of squared residuals obtained, but the contribution plots presented here also have control limits. These control limits are used for comparing the residuals of the new values for the present to the residuals of the fault-free data. If in the fault-free data a certain process variable had high residuals, it can also be expected to have high residuals in the present. However, if there are high residuals for a certain process variable that had low residuals in the fault-free data; this probably is due to a special event in the new process. If the contributions of a large group of process variables are studied, it is usually found that several process variables have high contributions. Using the control limits, it is easier to find those process variables that are really different, compared to the fault-free data. In this paper, a control limit for the contribution for each process variable is calculated as the mean of the contributions plus three times the standard deviation of the contributions for each process variable at each time period. That means contribution of each variable follows Gaussian distribution in a long period, but the matching mean and variance will be different. Control limits for contributions are studied by many researchers [32, 33]. Here, the control limits are determined according to the 3σ method [32]. The reason why we used the 3σ method is because of its efficiency and simplicity.


## 5. Experimental Study

### 5.1. Data Pretreatment

In this case, 16 important variables shown in Table1 were selected for monitoring, covering almost all the wastewater plant qualities required by engineers and the government.Table 1
Selected sensors for process monitoring.
Variables
Comments
Symbol
1
Total nitrogen concentration in effluent (mg N/L)
TN
2
Total chemical oxygen demand in the effluent (mg COD/L)
TCOD
3
NH4+ + NH3 nitrogen in effluent (mg N/L)
S
N
H
4
Suspended solids concentration in effluent (mg SS/L)
SS
5
Biological oxygen demand concentration in effluent [g COD m−3]
BOD
6
Suspended solids concentration in influent (mg SS/L)
SS
7
NH4+ + NH3 nitrogen in influent (mg N/L)
S
N
H
8
Oxygen in the third reactor g (−COD)⋅m−3;
S
O
9
Nitrate and nitrite nitrogen in the third reactor (mg N/L)
S
N
O
10
NH4+ + NH3 nitrogen in the third reactor (mg N/L)
S
N
H
11
Oxygen in the fifth reactor g (COD)/m3;
S
O
12
Nitrate and nitrite nitrogen in the fifth reactor (mg N/L)
S
N
O
13
NH4+ + NH3 nitrogen in the fifth reactor (mg N/L)
S
N
H
14
Flowrate in the third reactor m3/d
Flowrate
15
Oxygen transfer coefficient in the third reactor
KLa
16
Flowrate in the fifth reactor m3/d
FlowrateAs depicted in Figures3(a) and 3(b), the historical data sets reveal that, during dry weather periods, daily SS and SNH concentration follow a very similar pattern, which is also observed in other variables. In this respect, the typical data profiles for each variable under dry weather conditions were generated. This was done by stacking all the historical data in one set (from 0 to 1 day), and from it calculating the average flow at each time step. It is worth noting that daily data from Monday to Friday are very similar, whereas a slightly different pattern (shift in amplitude) is observed on the weekend profiles. Therefore, two typical data profiles were generated for each variable: one for weekdays and one for weekends. These patterns are consistent with the water consumptions by inhabitants during daily life. Results about SS and SNH concentrations are presented in Figures 3(c) and 3(d). By subtracting typical data from historical data, dynamic behaviors during the wastewater plant dataset are capable of being broken up, thereby providing better data for fault detection, identification, and even prediction. Moreover, to make sure all data from different variables are consistent mutually, they were centered and scaled accordingly. We denote all these procedures as “data whitening.” Finally, the whitening data were further used for ensuing VBMFA model building.Figure 3
Typical variables for workdays and weekdays.
(a)
(b)
(c)
(d)


### 5.2. Fault Diagnosis

#### 5.2.1. Fault Detection for the Fault-Free Scenario

1344 samples in two weeks under normal operation condition were collected to assess the proposed methodology. 900 data points were used for VBMFA model training; the remaining was for testing. As depicted in Figures4(a) and 4(c), the FA model performed well with most of points under the 95% control limit in terms of T2 and SPE statistics. However, it is worthy to emphasize that unacceptable numbers of false alarm occurred during this process. This is due to the fact that, when the measurements reach the peck or the valley, FA model is incapable of capturing such strong nonlinearity properly, thus leading to too many false alarms frequently. On the contrary, better performance was achieved by VBMFA with almost all of test points under 95% control limit. It is obvious that the control limit changed along with new data points coming, resulting in a more compact control limit and less probability of false alarms.Figure 4
Fault diagnosis under normal condition using FA and VBMFA.
(a)
(b)
(c)
(d)Table2 depicts the performance of four different monitor strategies, including PCA, FA, VBPCA, and VBMFA, under fault-free scenario in dry weather conditions. The frequency of false alarm for PCA and FA was very similar in terms of T2, but the SPE for both of strategies were quite different. The FA relying on different variance for different variables achieved less false alarms in comparison to the PCA. Furthermore, variational Bayesian learning is incorporated into PCA to learn its parameters. The number of false alarms was reduced as expected. Since variational Bayesian learning is the sensitive posterior probability mass rather than the posterior probability density, thus making it more resistant against overfitting compared to other estimation methods (least-squares, point estimation). Furthermore, because uncertainty information for unknown quantities can be obtained by this variational Bayesian learning, it is useful to detect unreliable results. However, due to strong nonlinearity and dynamics, VBPCA failed to capture all the characteristics properly. Obviously, VBMFA achieved the best performance. This is an important improvement as the peaks or the valleys that frequently occurred in the concentrations of different variables were described accurately, thus leading to better monitor effectiveness. It is necessary to emphasize that all model parameters excluding VBMFA were obtained somehow with the help of manual operation. This would potentially frustrate the use of these methods in practice. Conversely, the number of factors for each factor analyzer in VBMFA was determined automatically.Table 2
Selected sensors for process monitoring.
Statistics\methods
PCA
FA
VBPCA
VBMFA
T
2
26/444
24/444
16/444
0/444
SPE
21/444
15/444
13/444
1/444
Parameters
PC: 4
Factors: 4
PC: 3
Maximum factors: 14Local factors: 2, 3, 1


#### 5.2.2. Fault Detection for an Abrupt Change and a Drift Fault

Abrupt change and drift fault are two types of common faults in the process control. Thus, to access the fault detection ability of VBMFA, two cases with different faults were generated during the testing data on the second sensor (TCOD) and the fourth sensor (SS), respectively, tabulated as Table3. Due to cost saving demand and stricter regulations from governments, TCOD and SS sensors are fairly important for WWTP monitoring.Table 3
Faults generated on two important sensors.
Fault scenarios
Comments on the fault type
Case1 (the 2nd sensor)
30% positive bias on the second sensor (an abrupt change fault)
Case2 (the 4th  sensor)
35% positive bias of standard deviation (drift fault)Therefore, it is monitored not only by the engineers from WWTP, but also by the government using a remote monitoring system. In the first type of fault, normal data was corrupted by a 30% abrupt change from day from the sample 200 to 300. Results from Figures5(a) and 5(c) suggest that ST2 and SSPE index can recognize the fault accurately. Even if a valley occurred in Figure 5(a) between data points 200 and 300, which potentially leads to a real fault being failed to be identified, VBMFA was able to adjust its control limit to adapt to such variation. To further illustrate the efficiency of VBMFA, the second scenario with a drift fault exerting on the fourth sensor was also simulated, showing in Figures 5(b) and 5(d) that it achieved a good performance with accurate detection when the drift fault was happening. It is important to notice that the control limit of ST2 decreased to lower values once it recognized significant changes in the updating data, thereby providing more potentials to identify the drift fault earlier.Figure 5
Fault detection for an abrupt change and a drift fault using VBMFA.
(a)
(b)
(c)
(d)


### 5.3. Fault Identification

In order to find the root cause of the corresponding fault, the contribution plot was employed in this paper. The contribution plot of SSPE for the abrupt change can discriminate between the faulty sensor and the normal sensors (Figure6(a)). It is also obvious in Figure 6(a) that the faulty sensor (TCOD) deviated from the major pattern and exceeded its control limit due to the negative influence of abrupt changes. Similarly, the contribution plot for the drift fault was shown in Figure 6(b). The result clearly shows that the contribution plot of the VBMFA strategy successfully separated the out of control data from normal operation ones.Figure 6
Fault identification for an abrupt change and a drift fault using VBMFA.
(a)
(b)


### 5.4. Parameters Analysis

When optimizing (8), some qsi could be zero. This can be attributed to the fact that local data is insufficient to support the dimensional complexity on the factor loading matrices. Such redundant components should be removed to increase F, because they are no longer needed. This is the component death.Component birth does not happen spontaneously; thus, a heuristic is introduced. A parent-component stochastically has a birth with probability proportional toe-βF, and we attempt to split it into two. The parameter distributions of the two Gaussians created from the split are initialized by partitioning the responsibilities for the data, qsi. This usually causes F to decrease, so by monitoring the progress of F, this attempted birth can be rejected if F does not recover.Furthermore, in the VBMFA model, the dimensionality of each factor analyzer is determined automatically by automatic relevance determination (ARD) with the potentials to improve VBMFA model performance dramatically [21, 34]. As such, the corresponding weight of each factor analyzer can be obtained in an automatic manner accordingly without necessarily resorting to trial and error method. The ST2 and SSPE index can be weighted and updated adaptively.As for the maximum factors for each factor analyzer, it represents the freedom of parameters the ARD method is able to learn. Thus, this parameter is just set top-1 or p-2.Each of the mixture components takes for a particular data pointxnew and on the parameters θ of the mixture components, not on the values of other data points. Also, due to the parameters θ having been identified in the training stage already, the update of qsi is only related to the data point xnew in fact.


### 5.5. Discussions

The present work describes the development of a VBMFA model to monitor a WWTP. This is an attempt to provide water utilities with a tool to improve the process monitor ability and formulate a guideline for system operation.The VBMFA model has been successfully applied to a simulated WWTP with strong nonlinearity and dynamics in the dry weather. Results showed that false alarms can be avoided significantly under fault-free scenario. Also, accurate fault detection and identification can be achieved if an abrupt change or a drift fault occurred. Due to the mixture of factor analyzers for modeling, VBMFA is able to describe a small region of process in a locally linear way and the entire process in a globally nonlinear manner and therefore facilitate catching the characteristic of process even in the peaks or the valleys in terms of results in the previous sections.In this paper, we show the benefits of VBMFA model for a WWTP process monitor. The results clearly demonstrate the superior performance of VBMFA model compared with other methods. However, future work should update not only the weights of each factor analyzer and corresponding monitor indexes, but also the associated variance and mean. Also, future work should aim to test and validate this methodology in a real WWTP system, which incurs us to implement our work in separate computer as shown in Figure2. Additionally, due to the demands of cost saving in the water communities, it is imperative to assimilate economic index into the process monitoring in our future work.


## 5.1. Data Pretreatment

In this case, 16 important variables shown in Table1 were selected for monitoring, covering almost all the wastewater plant qualities required by engineers and the government.Table 1
Selected sensors for process monitoring.
Variables
Comments
Symbol
1
Total nitrogen concentration in effluent (mg N/L)
TN
2
Total chemical oxygen demand in the effluent (mg COD/L)
TCOD
3
NH4+ + NH3 nitrogen in effluent (mg N/L)
S
N
H
4
Suspended solids concentration in effluent (mg SS/L)
SS
5
Biological oxygen demand concentration in effluent [g COD m−3]
BOD
6
Suspended solids concentration in influent (mg SS/L)
SS
7
NH4+ + NH3 nitrogen in influent (mg N/L)
S
N
H
8
Oxygen in the third reactor g (−COD)⋅m−3;
S
O
9
Nitrate and nitrite nitrogen in the third reactor (mg N/L)
S
N
O
10
NH4+ + NH3 nitrogen in the third reactor (mg N/L)
S
N
H
11
Oxygen in the fifth reactor g (COD)/m3;
S
O
12
Nitrate and nitrite nitrogen in the fifth reactor (mg N/L)
S
N
O
13
NH4+ + NH3 nitrogen in the fifth reactor (mg N/L)
S
N
H
14
Flowrate in the third reactor m3/d
Flowrate
15
Oxygen transfer coefficient in the third reactor
KLa
16
Flowrate in the fifth reactor m3/d
FlowrateAs depicted in Figures3(a) and 3(b), the historical data sets reveal that, during dry weather periods, daily SS and SNH concentration follow a very similar pattern, which is also observed in other variables. In this respect, the typical data profiles for each variable under dry weather conditions were generated. This was done by stacking all the historical data in one set (from 0 to 1 day), and from it calculating the average flow at each time step. It is worth noting that daily data from Monday to Friday are very similar, whereas a slightly different pattern (shift in amplitude) is observed on the weekend profiles. Therefore, two typical data profiles were generated for each variable: one for weekdays and one for weekends. These patterns are consistent with the water consumptions by inhabitants during daily life. Results about SS and SNH concentrations are presented in Figures 3(c) and 3(d). By subtracting typical data from historical data, dynamic behaviors during the wastewater plant dataset are capable of being broken up, thereby providing better data for fault detection, identification, and even prediction. Moreover, to make sure all data from different variables are consistent mutually, they were centered and scaled accordingly. We denote all these procedures as “data whitening.” Finally, the whitening data were further used for ensuing VBMFA model building.Figure 3
Typical variables for workdays and weekdays.
(a)
(b)
(c)
(d)


## 5.2. Fault Diagnosis

### 5.2.1. Fault Detection for the Fault-Free Scenario

1344 samples in two weeks under normal operation condition were collected to assess the proposed methodology. 900 data points were used for VBMFA model training; the remaining was for testing. As depicted in Figures4(a) and 4(c), the FA model performed well with most of points under the 95% control limit in terms of T2 and SPE statistics. However, it is worthy to emphasize that unacceptable numbers of false alarm occurred during this process. This is due to the fact that, when the measurements reach the peck or the valley, FA model is incapable of capturing such strong nonlinearity properly, thus leading to too many false alarms frequently. On the contrary, better performance was achieved by VBMFA with almost all of test points under 95% control limit. It is obvious that the control limit changed along with new data points coming, resulting in a more compact control limit and less probability of false alarms.Figure 4
Fault diagnosis under normal condition using FA and VBMFA.
(a)
(b)
(c)
(d)Table2 depicts the performance of four different monitor strategies, including PCA, FA, VBPCA, and VBMFA, under fault-free scenario in dry weather conditions. The frequency of false alarm for PCA and FA was very similar in terms of T2, but the SPE for both of strategies were quite different. The FA relying on different variance for different variables achieved less false alarms in comparison to the PCA. Furthermore, variational Bayesian learning is incorporated into PCA to learn its parameters. The number of false alarms was reduced as expected. Since variational Bayesian learning is the sensitive posterior probability mass rather than the posterior probability density, thus making it more resistant against overfitting compared to other estimation methods (least-squares, point estimation). Furthermore, because uncertainty information for unknown quantities can be obtained by this variational Bayesian learning, it is useful to detect unreliable results. However, due to strong nonlinearity and dynamics, VBPCA failed to capture all the characteristics properly. Obviously, VBMFA achieved the best performance. This is an important improvement as the peaks or the valleys that frequently occurred in the concentrations of different variables were described accurately, thus leading to better monitor effectiveness. It is necessary to emphasize that all model parameters excluding VBMFA were obtained somehow with the help of manual operation. This would potentially frustrate the use of these methods in practice. Conversely, the number of factors for each factor analyzer in VBMFA was determined automatically.Table 2
Selected sensors for process monitoring.
Statistics\methods
PCA
FA
VBPCA
VBMFA
T
2
26/444
24/444
16/444
0/444
SPE
21/444
15/444
13/444
1/444
Parameters
PC: 4
Factors: 4
PC: 3
Maximum factors: 14Local factors: 2, 3, 1


### 5.2.2. Fault Detection for an Abrupt Change and a Drift Fault

Abrupt change and drift fault are two types of common faults in the process control. Thus, to access the fault detection ability of VBMFA, two cases with different faults were generated during the testing data on the second sensor (TCOD) and the fourth sensor (SS), respectively, tabulated as Table3. Due to cost saving demand and stricter regulations from governments, TCOD and SS sensors are fairly important for WWTP monitoring.Table 3
Faults generated on two important sensors.
Fault scenarios
Comments on the fault type
Case1 (the 2nd sensor)
30% positive bias on the second sensor (an abrupt change fault)
Case2 (the 4th  sensor)
35% positive bias of standard deviation (drift fault)Therefore, it is monitored not only by the engineers from WWTP, but also by the government using a remote monitoring system. In the first type of fault, normal data was corrupted by a 30% abrupt change from day from the sample 200 to 300. Results from Figures5(a) and 5(c) suggest that ST2 and SSPE index can recognize the fault accurately. Even if a valley occurred in Figure 5(a) between data points 200 and 300, which potentially leads to a real fault being failed to be identified, VBMFA was able to adjust its control limit to adapt to such variation. To further illustrate the efficiency of VBMFA, the second scenario with a drift fault exerting on the fourth sensor was also simulated, showing in Figures 5(b) and 5(d) that it achieved a good performance with accurate detection when the drift fault was happening. It is important to notice that the control limit of ST2 decreased to lower values once it recognized significant changes in the updating data, thereby providing more potentials to identify the drift fault earlier.Figure 5
Fault detection for an abrupt change and a drift fault using VBMFA.
(a)
(b)
(c)
(d)


## 5.2.1. Fault Detection for the Fault-Free Scenario

1344 samples in two weeks under normal operation condition were collected to assess the proposed methodology. 900 data points were used for VBMFA model training; the remaining was for testing. As depicted in Figures4(a) and 4(c), the FA model performed well with most of points under the 95% control limit in terms of T2 and SPE statistics. However, it is worthy to emphasize that unacceptable numbers of false alarm occurred during this process. This is due to the fact that, when the measurements reach the peck or the valley, FA model is incapable of capturing such strong nonlinearity properly, thus leading to too many false alarms frequently. On the contrary, better performance was achieved by VBMFA with almost all of test points under 95% control limit. It is obvious that the control limit changed along with new data points coming, resulting in a more compact control limit and less probability of false alarms.Figure 4
Fault diagnosis under normal condition using FA and VBMFA.
(a)
(b)
(c)
(d)Table2 depicts the performance of four different monitor strategies, including PCA, FA, VBPCA, and VBMFA, under fault-free scenario in dry weather conditions. The frequency of false alarm for PCA and FA was very similar in terms of T2, but the SPE for both of strategies were quite different. The FA relying on different variance for different variables achieved less false alarms in comparison to the PCA. Furthermore, variational Bayesian learning is incorporated into PCA to learn its parameters. The number of false alarms was reduced as expected. Since variational Bayesian learning is the sensitive posterior probability mass rather than the posterior probability density, thus making it more resistant against overfitting compared to other estimation methods (least-squares, point estimation). Furthermore, because uncertainty information for unknown quantities can be obtained by this variational Bayesian learning, it is useful to detect unreliable results. However, due to strong nonlinearity and dynamics, VBPCA failed to capture all the characteristics properly. Obviously, VBMFA achieved the best performance. This is an important improvement as the peaks or the valleys that frequently occurred in the concentrations of different variables were described accurately, thus leading to better monitor effectiveness. It is necessary to emphasize that all model parameters excluding VBMFA were obtained somehow with the help of manual operation. This would potentially frustrate the use of these methods in practice. Conversely, the number of factors for each factor analyzer in VBMFA was determined automatically.Table 2
Selected sensors for process monitoring.
Statistics\methods
PCA
FA
VBPCA
VBMFA
T
2
26/444
24/444
16/444
0/444
SPE
21/444
15/444
13/444
1/444
Parameters
PC: 4
Factors: 4
PC: 3
Maximum factors: 14Local factors: 2, 3, 1


## 5.2.2. Fault Detection for an Abrupt Change and a Drift Fault

Abrupt change and drift fault are two types of common faults in the process control. Thus, to access the fault detection ability of VBMFA, two cases with different faults were generated during the testing data on the second sensor (TCOD) and the fourth sensor (SS), respectively, tabulated as Table3. Due to cost saving demand and stricter regulations from governments, TCOD and SS sensors are fairly important for WWTP monitoring.Table 3
Faults generated on two important sensors.
Fault scenarios
Comments on the fault type
Case1 (the 2nd sensor)
30% positive bias on the second sensor (an abrupt change fault)
Case2 (the 4th  sensor)
35% positive bias of standard deviation (drift fault)Therefore, it is monitored not only by the engineers from WWTP, but also by the government using a remote monitoring system. In the first type of fault, normal data was corrupted by a 30% abrupt change from day from the sample 200 to 300. Results from Figures5(a) and 5(c) suggest that ST2 and SSPE index can recognize the fault accurately. Even if a valley occurred in Figure 5(a) between data points 200 and 300, which potentially leads to a real fault being failed to be identified, VBMFA was able to adjust its control limit to adapt to such variation. To further illustrate the efficiency of VBMFA, the second scenario with a drift fault exerting on the fourth sensor was also simulated, showing in Figures 5(b) and 5(d) that it achieved a good performance with accurate detection when the drift fault was happening. It is important to notice that the control limit of ST2 decreased to lower values once it recognized significant changes in the updating data, thereby providing more potentials to identify the drift fault earlier.Figure 5
Fault detection for an abrupt change and a drift fault using VBMFA.
(a)
(b)
(c)
(d)


## 5.3. Fault Identification

In order to find the root cause of the corresponding fault, the contribution plot was employed in this paper. The contribution plot of SSPE for the abrupt change can discriminate between the faulty sensor and the normal sensors (Figure6(a)). It is also obvious in Figure 6(a) that the faulty sensor (TCOD) deviated from the major pattern and exceeded its control limit due to the negative influence of abrupt changes. Similarly, the contribution plot for the drift fault was shown in Figure 6(b). The result clearly shows that the contribution plot of the VBMFA strategy successfully separated the out of control data from normal operation ones.Figure 6
Fault identification for an abrupt change and a drift fault using VBMFA.
(a)
(b)


## 5.4. Parameters Analysis

When optimizing (8), some qsi could be zero. This can be attributed to the fact that local data is insufficient to support the dimensional complexity on the factor loading matrices. Such redundant components should be removed to increase F, because they are no longer needed. This is the component death.Component birth does not happen spontaneously; thus, a heuristic is introduced. A parent-component stochastically has a birth with probability proportional toe-βF, and we attempt to split it into two. The parameter distributions of the two Gaussians created from the split are initialized by partitioning the responsibilities for the data, qsi. This usually causes F to decrease, so by monitoring the progress of F, this attempted birth can be rejected if F does not recover.Furthermore, in the VBMFA model, the dimensionality of each factor analyzer is determined automatically by automatic relevance determination (ARD) with the potentials to improve VBMFA model performance dramatically [21, 34]. As such, the corresponding weight of each factor analyzer can be obtained in an automatic manner accordingly without necessarily resorting to trial and error method. The ST2 and SSPE index can be weighted and updated adaptively.As for the maximum factors for each factor analyzer, it represents the freedom of parameters the ARD method is able to learn. Thus, this parameter is just set top-1 or p-2.Each of the mixture components takes for a particular data pointxnew and on the parameters θ of the mixture components, not on the values of other data points. Also, due to the parameters θ having been identified in the training stage already, the update of qsi is only related to the data point xnew in fact.


## 5.5. Discussions

The present work describes the development of a VBMFA model to monitor a WWTP. This is an attempt to provide water utilities with a tool to improve the process monitor ability and formulate a guideline for system operation.The VBMFA model has been successfully applied to a simulated WWTP with strong nonlinearity and dynamics in the dry weather. Results showed that false alarms can be avoided significantly under fault-free scenario. Also, accurate fault detection and identification can be achieved if an abrupt change or a drift fault occurred. Due to the mixture of factor analyzers for modeling, VBMFA is able to describe a small region of process in a locally linear way and the entire process in a globally nonlinear manner and therefore facilitate catching the characteristic of process even in the peaks or the valleys in terms of results in the previous sections.In this paper, we show the benefits of VBMFA model for a WWTP process monitor. The results clearly demonstrate the superior performance of VBMFA model compared with other methods. However, future work should update not only the weights of each factor analyzer and corresponding monitor indexes, but also the associated variance and mean. Also, future work should aim to test and validate this methodology in a real WWTP system, which incurs us to implement our work in separate computer as shown in Figure2. Additionally, due to the demands of cost saving in the water communities, it is imperative to assimilate economic index into the process monitoring in our future work.



---

*Source: 2034826-2016-03-21.xml*