Dynamic fault detection and diagnosis for alkaline water electrolyzer with variational Bayesian Sparse principal component analysis
Qi Zhang Weihua Xu Lei Xie Hongye Su




	
	
	
	label1]Qi Zhang 
	label1]Weihua Xu
	label1]Lei Xie cor1
	
	label1]Hongye Su
	
	
	
	[label1]
	State Key Laboratory of Industrial Control Technology, Zhejiang University, 310027 Hangzhou, China
	[cor1]Corresponding author




Electrolytic hydrogen production serves as not only a vital source of green hydrogen but also a key strategy for addressing renewable energy consumption challenges. 
For the safe production of hydrogen through alkaline water electrolyzer (AWE), dependable process monitoring technology is essential. 
However, random noise can easily contaminate the AWE process data collected in industrial settings, presenting new challenges for monitoring methods. 
In this study, we develop the variational Bayesian sparse principal component analysis (VBSPCA) method for process monitoring. 
VBSPCA methods based on Gaussian prior and Laplace prior are derived to obtain the sparsity of the projection matrix, which corresponds to ℓ_2 regularization and ℓ_1 regularization, respectively.
The correlation of dynamic latent variables is then analyzed by sparse autoregression and fault variables are diagnosed by fault reconstruction.
The effectiveness of the method is verified by an industrial hydrogen production process, and the test results demonstrated that both Gaussian prior and Laplace prior based VBSPCA can effectively detect and diagnose critical faults in AWEs.





 Alkaline water electrolyzer, sparse principle component analysis, fault diagnosis, sparse Bayesain learning, variational Bayesian 







§ INTRODUCTION







As the global energy transition accelerates,  Alkaline water electrolyzer (AWE) for hydrogen production has become the most widely used hydrogen production technology.  
As industrial hydrogen production grows in scale, the safety of hydrogen production has attracted more attention.
The process of producing hydrogen involves multiple scientific disciplines, and the system is particularly complex and difficult to control.
With the widespread adoption of Distributed control systems (DCS) and sensor technology in industrial settings, it is possible to collect and store large amounts of data from the AWE process.
Moreover, data mining and database technologies offer strong technical support for developing and applying data-driven modeling approaches in industrial AWE processes<cit.>. 
Hence, by extracting and analyzing information from process data, process monitoring models can be built for industrial applications to ensure safe operation of industrial water electrolyzer processes<cit.>.

Latent variable models such as Principal component analysis (PCA) and Partial least squares (PLS) enable the extraction of key information from the process to build efficient monitoring systems by projecting high-dimensional process data into a low-dimensional space<cit.>.
However, both PCA and PLS under the assumption that the data involved in the process exhibit no temporal correlation.
 In fact, there is significant dynamic correlation in the process data. 
Therefore, neither of these methods exploits the dynamic correlation between the data. 
There are three types of Dynamic latent variable (DLV) methods that have been most extensively studied<cit.>. 
The most conventional approach is to construct time lag matrices, such as the dynamic PCA method, but this approach significantly reduces the interpretability of faults<cit.>.
In addition, dynamic relationships can be modelled explicitly by state space models, such as variational Bayesian state space models<cit.>. 
There are also methods such as dynamic inner PCA <cit.> and dynamic inner canonical correlation analysis (DiCCA) <cit.> which focus more on the dynamic relationships between latent variables.
However, the disadvantage of the latent variable models is the high requirement for data quality<cit.>. 
The latent variable models look for the orthogonal direction with the largest data variance to separate noise and data as much as possible<cit.>. 
Given that industrial process data are easily contaminated by noise, latent variable models are sensitive to noise and outliers, which reduces the reliability of monitoring.

Sparse methods have a strong noise reduction capability that removes noise from the data while retaining critical information<cit.>. 
This is especially important for fault detection in noisy industrial environments.
Traditionally, sparsity is obtained by utilizing regularization in the loss function. ℓ_1 regularization and ℓ_2 regularization are the most common methods, corresponding to Lasso regression and Ridge regression, respectively<cit.>.
The core of Lasso regression is that it uses the sum of the absolute values of the weight coefficients as the penalty term, which makes the loss function of Lasso regression discontinuously derivable<cit.>.
The selection of the regularization factor is usually difficult, and it is empirically selected, resulting in non-robustness to changes in the algorithmic environment.
Bayesian sparsity treats regularization as a Bayesian framework for prior knowledge<cit.>.
In this approach, the parameter estimates are updated with the posterior distributions, treating the regularisation terms as prior distributions of the model parameters<cit.>.
In addition, sparse Bayesian learning (SBL) has become a hot topic of research because it improves the interpretability and generalisability of the model by taking parameter uncertainty into account<cit.>.
Expectation maximisation (EM), which reduces computational cost and handles missing values, is a common approach to implementing sparse Bayesian methods<cit.>.
However, EM can only provide point estimates of parameters that may fall into local optima or overfitting, and requires strict assumptions about the posterior distributions of the hidden variables. Therefore, the full Bayesian uncertainty of the model cannot be captured by the EM method<cit.>.
A more advantageous approach involves utilizing variational Bayesian methods, which approximate the complete posterior distribution of parameters and hidden variables, thereby circumventing the constraints associated with the EM algorithm<cit.>.



The main objective of this study is to develop a process monitoring technique for the AWE hydrogen production process that can effectively handle noise.
To achieve this, we have derived variational Bayesian principal component analysis (VBSPCA) methods for process monitoring based on Gaussian and Laplace priors, respectively, inspired by Bayesian sparsity.
The introduction of ℓ_2 regularization is equivalent to introducing a Gaussian prior, while the introduction of ℓ_1 regularization is equivalent to introducing a Laplace prior.
By employing matrices with a shared sparse structure and independent Gaussian priors for low-rank decomposition, Gaussian prior-based VBSPCA is capable of accurately modeling and estimating noise, thereby mitigating its effects. 
The Laplace prior-based VBSPCA designs a two-level hierarchical decomposition of the Laplace distribution to construct the conjugate structure for Bayesian posterior inference.
In addition, the parameters of VBSPCA are estimated by variational Bayesian inference, which not only quantifies the model uncertainty, but also makes the algorithm easy to handle.
Finally, to make full use of the structurally relevant information of the variables, we incorporate the vector autoregression (VAR) with ℓ_1 regularization to better extract dynamic features from the latent variables.
The latent variables processed by VBSPCA can reduce the influence of noise. Dynamic analysis of these latent variables can improve the reliability of fault detection and diagnosis.
The effectiveness of the proposed method is verified by an industrial hydrogen production process, and the test results show that both Gaussian prior and Laplace prior-based VBSPCA can effectively detect and diagnose critical faults in AWEs.





The main contributions of this work are summarized as follows:

	
	
	
	
  1) 
	This study develops VBSPCA techniques for AWE process diagnosis and fault detection, based on Gaussian prior and Laplace prior, respectively.
	
  2) 
	The developed method constructs a sparse prior that naturally reduces the impact of noise and quantifies uncertainty through variational inference.
	
  3) 
    ℓ_1 regularized vector autoregression extracts dynamic features from latent variables, which can exploit structural information and reduce noise disturbances.



The remainder of the paper is organized as follows. The related works are described in Section <ref>. The variational Bayesian Sparse PCA method is described in Section <ref>. The variational inference for these two VBSPCA techniques is presented in Section <ref>,
while dynamic fault detection and diagnosis approach are described in Section <ref>.
Simulation results are provided in Section <ref>, followed by concluding remarks in Section <ref>.







§ RELATED WORK



 §.§ PCA and Probabilistic PCA
 


Let x denote a sample vector of m sensors. Assuming that there are n samples for each sensor, a data matrix X=[x_1,x_2,…,x_n ] ∈ℝ^m × n is composed with each column representing a sample x_i∈ℝ^m.
With the projection matrix P∈ℝ^m × r,(r<m) satisfies P^⊤P=I_r, the PCA model can be obtained


    x_i= Pt_i +x̅+ξ_i ,    i = 1,2,…,n


where t_i denotes the latent variable with r dimensions,  x̅=1/n∑_i=1^nx_i∈ℝ^m is the mean vector,
ξ_i∼𝒩(0,σ^2I) represents additive process noise.

The PCA model can be reformulated by creating matrices T=(t_1,t_2,...,t_n ) ∈ℝ^r × n, X̅=(x̅_1,x̅_2,...,x̅_n ) ∈ℝ^m × n, and Ξ=(ξ_1,ξ_2,...,ξ_n ) ∈ℝ^m × n, PCA can be rewritten as the following

    X= PT +X̅+Ξ.




In order to maximize the variance of the projected data,  PCA can obtain the optimal P and T by optimizing the following objective function


    min_P,T | | X - X̅ - PT | | _  F  ^  2   s.t.P ^ ⊤P = I_r


where ||·||_F is the Frobenious norm. 


	In general, the optimal P can be obtained by stacking the r singular vectors corresponding to the first r largest singular values of the sample covariance matrix S=∑_i=1^n(x_i-x̅)(x_i-x̅)^⊤/(n-1). As a result, latent variable t_i=(x_i-x̅)P can be obtained<cit.>.


The disadvantage of PCA is that it is not possible to obtain probability distributions with respect to the projection matrix P and the latent variable t. 
 Probabilistic PCA (PPCA) can be obtained from a Gaussian latent variable model, and its parameters are estimated iteratively using the EM method<cit.>.
The conditional probability distribution of x_i can be expressed as


    p(x_i|t_i)=𝒩(x_i|Pt_i+x̅,σ^2I).


The prior distribution of the latent variable t_i is assumed to be a standard Gaussian distribution, the marginal distribution of x_i can thus be calculated by integrating out the latent variables


    p(x_i)=∫p(x_i|t_i)p(t_i)dt_i=𝒩(x_i|x̅ ,PP^⊤+σ^2I).


According to the Maximum likelihood estimation (MLE), the log-likelihood function is given by

    ℒ ( P , x̅ , σ ^  2 I| X ) =    ∑ _  i = 1  ^  n ln p ( x _  i  ) 
    	= ∑ _  i = 1  ^  n ln𝒩( x_i|x̅ ,PP^⊤+σ^2I) .


Maximizing ℒ ( P , x̅ , σ ^  2 I| X ) yields the optimal P, x̅, and σ^2. 




	Finding the projection matrix P, where r<m accomplishes the reduced dimensionality and feature extraction, is the goal of both PCA and PPCA. 
	In this way, the variable x_i can be mapped to the low-dimensional space and the latent variable (extracted feature) can be obtained.
	This is actually the process of feature extraction.



	Typically, the matrix P derived from PCA and PPCA is not sparse, leading to poor interpretability. If P is sparse, then the latent variable t can be expressed as a weighted combination of a small number of original features, thereby improving the interpretability of the model.






 §.§ Sparsity analysis

To obtain sparse solutions for P, Lasso regression (ℓ_1 regularization) and Ridge regression (ℓ_2 regularization) are commonly used. 
The optimization objective of Lasso regression is expressed as

    P^∗ = min_P1/n∑_i=1^n ||  x_i - Pt_i   ||+ λ ||P ||_1 .


The optimization objective of Ridge regression is expressed as

    P^∗ = min_P1/n∑_i=1^n ||  x_i -Pt_i   ||+ λ ||P ||_2^2 .



	Using the ℓ_1 norm as a regularization term allows for feature selection, resulting in a sparse model where some regression coefficients can be compressed to zero; while using the ℓ_2 norm as a regularization term can help prevent overfitting.



Following the conditional probability distribution of x_i in (<ref>), sparsity analysis will be reformulated as maximizing the posterior probability

    P^∗   = max_Pln( ∏ _  i = 1  ^  n   p(x_i|t_i, P ) p(P)   ) 
       = max_P∑_i=1^n(    ln p(x_i|t_i, P )+ ln∑_i=1^n p(p_i)               )

where ln p(x_i|t_i, P ) ∝( x_i - Pt_i )^2. 

It can be seen that the target loss part of (<ref>) is solved here by a probability  distribution.
The significance of p(P) is the prior assumption about the probability distribution of P. The Ridge Regression will be obtained by assuming that the prior distribution of P is a zero-mean Gaussian distribution P ( p _  i  ) = γ _  i /√( 2π)exp ( - γ _  i  | |  p _  i  || ^  2  ).
Correspondingly, (<ref>) can be calculated as

    P^∗ =max_P∑_i=1^n(    ln p(x_i|t_i, P )+   γ∑_i=1^n ||p_i ||^2      ) .



The Lasso Regression will be obtained by assuming that the prior distribution of P is a zero-mean Laplace distribution P ( p _  i  ) = φ_i/ 2 exp ( - φ | p_i | ),
Correspondingly, (<ref>) can be calculated as

    P^∗ =max_P∑_i=1^n(    ln p(x_i|t_i, P )+   φ∑_i=1^n |p_i |     ) .



Consequently, it can be shown that in Bayesian inference, the Laplace prior is in fact equivalent to ℓ_1 regularization, and the Gaussian prior is equivalent to ℓ_2 regularization. It can usually be solved by EM methods, however, as a point estimation method, the EM method can only estimate the maximum posterior of the single most probable value for each parameter. This parameter learning method may lead to local optima and overfitting, thus limiting the performance of the model. 



§ VARIATIONAL BAYESIAN SPARSE PCA
 






 §.§ Gaussian prior-based VBSPCA
 



Suppose  X can be decomposed into a low-rank matrix H and the noise Ξ


    X   = H + ξ
    		=∑_i=1^kp_it_i^⊤+ Ξ


where H=PT^⊤ such that 𝖱𝖺𝗇𝗄(H) = r ≪ min (n,m).



According to the principle of ℓ_2 regularized sparsity,  a sparsity inducing Gaussian priors is used on each coefficient of the unknown vector based on Automatic relevance determination (ARD).
The columns of P and T are associated with Gaussian priors of precisions γ to ensure that the majority of the columns in each are set to zero, so γ is closely related to the sparsity of the solution<cit.>


    p ( P | γ ) = ∏ _  i = 1  ^  k 𝒩 ( p_i | 0 , γ _  i  ^  - 1 I _  m  ), 
        p ( T | γ ) = ∏ _  i = 1  ^  k 𝒩 ( t_i | 0 , γ _  i  ^  - 1 I _  n  )


where I_m is the m × m identity matrix, I_n is the n × n identity matrix. Thus, the columns of P and T have the same sparsity profile enforced by the common precisions γ_i. 

 

In this inference process, many of the precisions γ_ i  will assume very large values, which effectively removes the corresponding outer-products from H  and thus reduces the rank of the estimate.

Gamma priors  p ( γ _  i  ) can be placed on the precision parameters of the above Gaussian distributions of P and T. The form of conjugate Gamma hyperprior is expressed as


    p ( γ _  i  ) = Gam ( a_0 ,  1 / b_0  ) ∝γ _  i  ^  a_0 - 1 exp ( - b_0 γ _  i  )





	To ensure that the hyperpriors are sufficiently vague and diffuse, we treat the parameters a_0 and b_0 as fixed constants and assign them very small values, such as 10^-5, instead of estimating them from the data.




	The noise in the measurements is assumed to be i.i.d. drawn from a Gaussian distribution, which has an effect on each of the measurements. Noting that this construction is only within the prior, and that an approximation to the posterior ξ can be inferred for perturbations to the prior.




Mathematically, the noise in (<ref>) is modeled as


    p(Ξ|β)=∏_j=1^m∏_i=1^n𝒩(Ξ_ji|0,β^-1).


The parameter selection of the prior α can largely affect the parameter estimation of the posterior distribution.
 To reduce the effect of the parameter selection of the original distribution, the noise precision β is set to the non-informative Jeffrey’s prior


    p ( β ) = β^-1 .



In accordance with the form of PCA in (<ref>), the mean X̅ in (<ref>) is modeled as

    p(X̅|α)=∏_j=1^m∏_i=1^n𝒩(X̅_ji|0,α_ji^-1).



	In classical PCA fault detection method, a preliminary step includes zero-centering. Hence, the decomposition of X doesn't consider the mean values.
	For the rigor of the algorithm, we give here the form of the mean X̅, which can also be interpreted as the sparse error.


Similarly, the precision α is set to the non-informative Jeffrey’s prior


    p ( α_ji ) = α_ji^-1 .


The conditional probability distribution for X can be expressed as follows


    p(X|P,T,β)=∏_j=1^m∏_i=1^n𝒩(y_ji|x_ji,β^-1) .


Finally, with the conditional probability and priors, the joint distribution is expressed as

    p ( X , P , T ,X̅ , γ , β ) =    p ( X | P , T ,X̅, β ) p ( P | γ ) p ( T | γ ) p(X̅|α) p(α) p (  γ )  p ( β ) .













 §.§ Laplace prior-based VBSPCA
 


The Laplace distribution has coarser tails than the Gaussian distribution, hence the name heavy-tailed distribution. It means that the Laplace distribution can fit heavy noise and outliers better than the Gaussian distribution.


	It is necessary to employ another prior, a two-level hierarchical decomposition of the Laplace distribution, as the Laplace distribution is not a conjugate prior to the Gaussian distribution and does not support posterior inference within a Bayesian framework.



According to the principle of ℓ_2 regularized sparsity, the first level of P introduces a Gaussian prior p(P_ji| η _  j i ) ∼𝒩(0 , η_ji), while the second level of P is  an exponential distribution hyperprior to the variance η _  j i, p ( η _  j i  ) =  1 /φexp ( - η _  j i /φ )

    p(P_ji|η _  j i )    = ∫ p ( P _  j i  | η _  j i  ) p ( η _  j i  ) d η _  j i 
       =  1 / 2 √( 2 /φ)exp ( - √( 2 /φ) | P _  j i  | )



Assuming that the noise ξ  follows a Gaussian distribution ξ∼𝒩(0,ϑ^-1), a prior over the inverse noise variance ϑ is introduced

    p(ϑ)=Gam(ϑ|c_0,d_0)  .



In accordance with the form of PCA in (<ref>), it is assumed that both the mean x̅ and the latent variable t have Gaussian distributions

    p(x̅)    ∼𝒩(0, ς^-1I )
    
    		p(t_i)    ∼𝒩(0,I )


  
Finally, with the conditional probability and priors, the joint distribution is expressed as

    p ( X , P,t, x̅, ϑ ) = 
    		∏_i=1^n p( y | P,t, x̅, ϑ, η
    		)P(t)p(P|η)p(η)p(x̅)p(ϑ)









§ VARIATIONAL BAYESIAN INFERENCE



 §.§ A review of variational inference

A variational Bayesian approach is implemented in this section to approximate the posterior distribution of the parameters.
The variational Bayesian of learning p(Θ|y) avoids intractable integration by learning an approximating distribution Q(Θ).
As a result, the marginal log-likelihood function ln p(y) can be separated into two components: the Evidence Lower Bound (ELBO) ℒ(Q(Θ)) and Kullback-Leibler (KL) divergence 𝕂𝕃(Q(Θ)||p(Θ|y))<cit.>


    ln p(y)   =ℒ(Q(Θ))+𝕂𝕃( Q(Θ)||p(Θ|y))  
       = ln∫_Θ Q(Θ)  {p(Θ,y)/Q(Θ)} d Θ - ∫_ΘQ(Θ) ln{p(Θ,y)/Q(Θ)} d Θ


where Q(Θ) approximates as closely as possible the true distribution p(Θ|y) through KL divergence.









	 Since KL(Q∥ p)≥ 0, ℒ(Q(Θ)) is a rigorous lower bound for ln p( y).
	Maximizing ℒ(Q(Θ)) is equivalent to minimizing KL( Q∥ p ), and thus the posterior distribution p(Θ |y) can be approximated by Q( Θ).




With the Mean field variational inference (MFVI), the free joint distribution is chosen to be factorable <cit.>

    Q(Θ)= ∏_i=1^Iq_i(Θ_i).


In order to find an approximate posterior distribution for each latent variable, the expected values of all parameters in the joint distribution are relative to their most recent distribution.
The posterior approximation q_i(Θ_i) for each latent variable Θ_i is optimized by


    q_j(Θ_j)=exp𝔼_i≠ j(ln p(y,Θ)) /∫_Θ_jexp(  𝔼_i≠ j(ln p(y,Θ))  )  .

where 𝔼(·) is the expectation.




 §.§ Variational inference for Gaussian prior-based VBSPCA


The latent variable and parameters of the Gaussian prior-based VBSPCA are Θ = (  P,T, X̅, γ, α, β). Furthermore, following (<ref>), Θ can be written as 

    Q(Θ) =q( P) q(T) q(X̅) q(γ)    q(α)    q(β).


Optimizing the approximate distribution Q(Θ)  based on the joint distribution of (<ref>) gives the form of these parameters and latent variables

	
    q( p_i)    = 𝒩 ( p _  i   , | μ^p _  i   , Σ^p_  i   )  
     q( t_j)    = 𝒩 ( t _  j   , | μ^t _  j   , Σ^t_  j   )  
     q( X̅_ji)    =𝒩 ( X̅_j i , | μ^X̅_j i  , Σ^X̅_j i )  
     q(γ)     = Gam(γ_i| a,b)




Following (<ref>) and (<ref>), α and β introduce a non-informative Jeffrey’s prior with the following expectations 

    𝔼(α_ji)    = 𝔼(β)+𝔼(α_ji) / 1+ 𝔼(X̅_ji) (𝔼(β)+𝔼(α_ji))  
    𝔼(β)    =  m n /𝔼(  | | X - PT ^ ⊤ -X̅ | | _  F  ^  2 )


The approximate distribution q( p_i) in (<ref>) is given by

    μ^p _  i      =   𝔼(β) Σ^p_  i  𝔼(𝕋)(y_i-𝔼(x̅_i))
    		
    Σ^p_  i      = (  𝔼(β)  𝔼 (T ^ ⊤T)  + 𝔼(γ)I)  ^  - 1   .


The approximate distribution q( t_i) in (<ref>) is given by

    μ^t _  j      =   𝔼(β) Σ^t_  j  𝔼(ℙ) (y_j-𝔼(x̅_j))
    		
    Σ^t_  j      = (  𝔼(β)  𝔼 (P ^ ⊤P)  + 𝔼(γ)I)  ^  - 1   .



The approximate distribution q( X̅_ji) in (<ref>) is given by

    μ^X̅_j i    =  𝔼(β)/𝔼(β)+𝔼(α_ji) (    X_ji -𝔼(p_i )  𝔼(t_i )^⊤  )
    		
    Σ^X̅_j i    =1/𝔼(β)+𝔼(α_ji) .


The approximate distribution q(γ) in (<ref>) is given by

    a      =     a_0 +m + n /2
    
    		b     =   b_0 + 𝔼( p_i^⊤p_i ) + 𝔼( t_i^⊤t_i )/2  .


The required moments are easily evaluated using basic properties of probability distributions

    𝔼 (P ^ ⊤P)    = 𝔼 (P)^⊤𝔼 (P)+mΣ^p_  i  , 
    𝔼 (T ^ ⊤T)    = 𝔼 (T)^⊤𝔼 (T)+nΣ^t_  i  , 
    𝔼( p_i^⊤p_i )    = 𝔼 (p_i)^⊤𝔼 (p_i)+m(Σ^p_  i  )_ii,
    𝔼( t_i^⊤t_i )    = 𝔼 (t_i)^⊤𝔼 (t_i)+n(Σ^t_  i  )_ii,
    𝔼(  | | X - PT ^ ⊤ -X̅ | | _  F  ^  2 )    = | | X - 𝔼(P) 𝔼 (T) ^ ⊤ - 𝔼( X̅_ji) | | _  F  ^  2  
    		+ t r (  n 𝔼 (P)^⊤𝔼 (P) Σ ^ t_i )  
       +  t r (  m 𝔼 (T)^⊤𝔼 (T) Σ ^ p_i ) 
    		+ t r (  mn Σ ^ p_i Σ ^ t_i ) +∑_i=1^n∑_j=1^mΣ^X̅_j i .




 §.§ Variational inference for VBSPCA with Laplace prior



Following (<ref>), the latent variables and parameters of the Laplace distribution-based VBSPCA are Θ = (  P,t, x̅, ϑ, η). Furthermore, Θ can be written as 

    Q(Θ) =q( P) q(t)  q(x̅)    q(ϑ)    q(η).





Optimizing the approximate distribution Q(Θ)  based on the joint distribution of (<ref>) gives the form of these parameters and latent variables

	
    q( P)    =∏ _  i = 1  ^  r 𝒩 ( P _  i   , | μ^P _  i   , Σ^P_  i   )  
    	q( t)    =∏ _  i = 1  ^  n 𝒩 ( t _  i   , | μ^t _  i   , Σ^t_  i   )  
    	q( x̅)    =𝒩 ( x̅ , | μ^x̅  , Σ^x̅ )  
    	q( ϑ)    = Gam(ϑ| c,d)  
    	q(η)    =  1 /√(πη _  i  j φ)exp(  -  1 / 2 η _  ji  ( P _  j i  ) ^ ⊤ ( P _  j i  ) - η_j i/φ+ √(2/φ) | P_j i |    )



The approximate distribution q( P) in(<ref>) is given by

    μ^P _  i      =   𝔼(ϑ) Σ^P_  i  ∑_i=1^n𝔼(x_i)
    		(y_i-𝔼(x̅))
    		
    Σ^P_  i      = (  𝔼(ϑ)  ∑_i=1^n𝔼(x_ix_i^⊤) dg(𝔼(η_i)) + I)  ^  - 1 dg(𝔼(η_i)) .

where dg converts the vector to a diagonal matrix.

The approximate distribution q( t) in (<ref>) is given by

    μ^t _  i      =   𝔼(ϑ) 	Σ^t_  i  𝔼(ℙ)(y_i-𝔼(x̅))
    		
    Σ^t_  i      = (  𝔼(ϑ)  𝔼 (P ^ ⊤P)  + I)  ^  - 1   .


The approximate distribution q( x̅) in (<ref>) is given by

    μ^x̅    =   𝔼(ϑ) 	Σ^x̅∑_i=1^n(  y_i - 𝔼(P) 𝔼(x_i)     ) 
    		
    Σ^x̅    = (  n𝔼(ϑ)  + ς)  ^  - 1 I .


The approximate distribution q( ϑ) in (<ref>) is given by

    c       =   c_0+nm/2
    
    		d     = d_0+ 1/2∑_i=1^n|| x̅ ||^2 +tr( 𝔼(P^⊤P) 𝔼(x_i x_i)^⊤ ) 
       +2𝔼(x̅)^⊤𝔼(P)𝔼(x_i) -2y_i^⊤𝔼(P)𝔼(x_i) - 2y_i^⊤𝔼(x̅)
    		+|| y_i ||^2.


The approximate distribution q( η) (<ref>) is given by

    μ_η_ji = 1/2 (φ+√(2φ) | P_ji| ).











§ DYNAMIC PROCESS MONITORING
  


 §.§ Fault detection


In this section, fault detection is accomplished and the dynamics of the latent variables are examined after the projection loading matrix P is calculated by variational inference using the VBSPCA approach, which is based on Gaussian prior and Laplace prior. 


	Fault detection and diagnosis depends more on the dynamics of the latent variables t^∗_k of the newly collected samples, although it also depends on the projection matrix P.
	Therefore, it is more interpretable to analyze the dynamics of the latent variables directly.



A new collected sample y_k∈ℝ^m is obtained at time k with time interval L, the latent variable  t^∗_k of the new collected sample y_k is obtained by a linear transformation of P

    t^∗_k = P^⊤y_k

where t^∗_k is a potential time series, which means that t^∗_k at the current moment is influenced by the previous moment.


As the dynamic relationship of the latent variable t^∗_k is not explained, it is necessary to build a dynamic model to characterize the auto correlation.
As a result, t^∗_k can be expressed as the following VAR model


    t^∗_k=d=1τ∑ ω_dt^∗_k-d+ϵ_k, k = τ+1,⋯,L


where τ is the order of dynamics,  ϵ_k is zero mean Gaussian noise, ω_i∈ℝ^r × r, d = 1,2, ⋯, τ represents the coefficient matrix of the VAR model.

The objective is to learn the coefficient matrix ω from the latent variable, and the new linear expression for ω and χ_k is rewritten as

    ω   = [ω_1,⋯,ω_τ] ^⊤∈ℝ^(r ·τ)× r
    χ_k   = [t^∗_k-1^⊤,⋯, t^∗_k-τ^⊤]  ∈ℝ^(r ·τ)





The VAR model for the latent variable t can be rewritten as


    t_k^∗≈d=1τ∑ ω_kt_k-d^∗=ω^⊤χ_k , k = τ+1,⋯,L



In high dimensional scaling, where the transition matrix driving the temporal dynamics displays a more complex structure. In order to train the sparse projection in the regression model, ℓ_1 regularization is therefore required


    ω̂=min_ω1/2Z-Qω_2^2+λω_1



where Z and Q are defined as Z=[ t_τ+1^⊤,...,t_L^⊤]^⊤∈ℝ^(F-τ)× r, Q =[ χ_τ+1^⊤,...,χ_L^⊤]^⊤∈ℝ^(F-τ)×(r·τ), respectively.


	A sparse VAR about ω is built such that ω becomes a sparse matrix for feature selection on the latent variable t^∗_k of the newly collected sample y_k. This further reduces the noise perturbation of the newly collected sample y_k and improves the model interpretability.



The prediction of t̂_k from the sparse VAR model is 


    t̂_k = ∑ _d=1^τω̂_d	t^∗_k-d = ∑ _d=1^τω̂ _i P^⊤ y_k-d , k = τ+1,⋯,L   .



From the predicted latent variable t̂_k, the T_2 statistics and SPE statistics can be obtained, which are the essential indicators for performing fault detection.


The T_2 statistic can be constructed as

    T_2 = t̂_k Λt̂_k ^⊤.

where Λ is a diagonal matrix.

The SPE statistic can be constructed as

    SPE = ( I-PP^⊤)t̂_k ^2.


In this work, the kernel density estimation (KDE) method is used to calculate the control limit, which is a non-parametric method to build an unknown density function


    f̂ _h( e )=1/Li=1L∑ K_h( e-e_i)=1/Lhi=1L∑ K_h( e-e_i/h).

where K_h is a non-negative kernel function and h is a relative smoothing parameter called bandwidth.  e is the reconstruction error, to construct the density estimate accordingly. The confidence limit can be obtained with a confidence level α = 0.95.
When the statistic exceeds the confidence limit, the test sample y_k is considered faulty.







 §.§ Fault diagnosis




A detectable fault should have an impact on the measurement vector that deviates from the normal scenario. 
The measurement vector for the fault-free part is denoted by ψ and will be unavailable when a fault occurs<cit.>. 
In the presence of a process fault, y is represented as

    y=ψ+fϕ

where ϕ is orthonormal so that ||f|| represents the magnitude of the fault.
ψ is the fault-free part of the measurement.


The objective of fault reconstruction is to estimate the fault-free part ψ by removing the effects of faults. 
A reconstructed sample vector ψ_k ∈ψ is calculated by correcting the effect of a fault on the process data y

    ψ_k=y-f_k ϕ_k




The fault detection index for reconstructing variables can be written as


    Index( ψ_k)=y_k^⊤Φy_k=y_k_Φ^2=ψ-f_kϕ_k_Φ^2 .


The purpose of reconstruction is to find the value that minimizes f_k. 
	Taking the derivation of (<ref>), 
f_k=( ϕ_k^⊤Φϕ_k)^-1ϕ_k^⊤Φy can be obtained. 
The reconstruction-based contribution RBC_k^Index along the fault direction can be expressed as


    RBC_k^Index=ϕ_kf_k_y^2
    	=ψ^⊤Φϕ_k( ϕ_k^⊤Φϕ_k)^-1ϕ_k^⊤Φy .


The reconstructed index Index(ψ_k ) in (<ref>) can be expressed as

    Index(ψ_k )= y^⊤Φ[ I -  ϕ_k( ϕ_k^⊤Φϕ_k) ^-1ϕ_k^⊤Φ]    y
    	=y^⊤Φy-RBC_k^Index


As a result, Index(x)=Index(ψ_k )+RBC_k^Index.

For T_2 indexes, RBC_k^T_2 can be expressed as

    RBC_k^T_2= (  ϕ_k^⊤PΛ^-1P^⊤y)^2 /( PΛ^-1P^⊤) _kk

where ( PΛ^-1P^⊤) _kk is the ith diagonal element of ( PΛ^-1P^⊤).

For SPE indexes, RBC_k^SPE can be expressed as

    RBC_k^SPE= (  ϕ_k^⊤( I-P^⊤P)  y)^2 /( I-P^⊤P)_kk

where ( I-P^⊤P)_kk is the ith diagonal element of ( I-P^⊤P).
















§ INDUSTRIAL AWE HYDROGEN PRODUCTION STUDIES
 


	The AWE hydrogen production industrial process consists of three main parts, electrolysis module, gas-liquid separation module, and circulation module. As shown in Figure <ref>, they include alkaline electrolyzer, heat exchanger, gas-liquid separator, condenser and dryer.
		The direct current supply powers the alkaline water electrolyzer, and the electrolyte is then decomposed to generate hydrogen and oxygen.
	Hydrogen and oxygen are dissolved in the electrolyte. When the pressure rises to the rated value, the mixtures is sent to the gas-liquid separator through the back pressure valve, where gas and electrolyte are separated.
	After separation, the gases and liquid go through different routes.
	The filtered electrolyte is cooled to the working temperature by the condenser, and then the electrolyte is returned to the electrolytic cell through the circulating pump, which can reduce the temperature of the electrolyzer.
	The separated gas is cooled by the condenser, and the cooled gas enters the gas dryer for dehydration. In order to improve the purity of hydrogen, hydrogen is purified and stored by H_2 purification equipment.
	The liquid flows through the circulation pump and filter and then returns to the electrolysis cell from the bottom of the separators. The hydrogen is tested for purity and then collected and stored.
	During the electrolysis process, the circulation module cools the system while circulating the electrolyte. Pure water is consumed by electrolysis, so it is necessary to supplement water to get electrolyte. The generated electrolyte and the separated electrolyte are circulated into the electrolyzer to maintain the set concentration.
	
		
	

	The data of hydrogen production process used in this study comes from the real data of a large-scale hydrogen production plant, including normal data and fault data. 
	According to privacy requirements of the plant, the raw data is desensitized in order to ensure the privacy and security of the information.
 The variables of the electrolytic water hydrogen production process can be divided into control variables and process variables, and the regulation and monitoring of the variables can be realized by the automated control system to ensure the stability, efficiency and safety of the electrolytic water hydrogen production process. The plant collects 32 variables, including 10 control variables and 22 process variables. Table <ref> lists the control variables and Table <ref> lists the process variables.
	The main faults that may occur in a AWE hydrogen production system are summarized in the Table <ref>.
Ten common faults are present in all three components of the AWE hydrogen production system. 

	






	









	In this section, two faults in the three modules of the hydrogen production process are analyzed. All faults in the test sample were introduced from the 201st sample.
	Fault 2 is the  electrolyte inter-electrode voltage is abnormal. 
	The hydrogen production rate in the alkaline water electrolyzer depends on the DC power supply output. The higher the voltage applied to the water, the faster the water molecules will split into hydrogen and oxygen gases. This is because the voltage provides enough energy to overcome the bond energy of water and drive the electrolysis reaction. 
	When the voltage increases, more electrons are transferred to the water, making the reaction more vigorous. However, too high a voltage can also generate excessive bubbles that can lower the reaction efficiency. 
	Therefore, it is important to choose an optimal voltage for water electrolyzer to achieve the best balance between reaction speed and effect.
	When the voltage is abnormal, it can cause the electrolysis reaction to not proceed properly. High voltage can cause bubbles to accumulate in the electrolyzer, posing a fire and explosion hazard. Voltage fluctuations can lead to unstable electrolysis reactions and the hydrogen production system poses an explosion risk.



	
	By capturing the dynamic relationships of multiple variables through lagged values of each variable and regression analysis, the response of a dynamic latent variable to another variable at different points in time can be calculated. 
	The predicted values of the latent variables calculated by the sparse VAR model are shown in Figure <ref>, where the five dynamic latent variables have been able to characterize the full range of changes in the system. 
	The correlations of the obtained dynamic latent variables were checked by plotting correlation plots. The correlations among the dynamic latent variables are presented in Figure <ref>. As depicted, it is evident that the correlations have been eliminated through modeling, whereby the first five dynamic latent variables are extracted.

	
		
	
	

			

	

	
	The results of the fault detection are shown in Figure <ref>, where the blue line represents the fault statistic, which presents the variation of the fault process. The red line is the confidence limit, and a statistic above the confidence limit indicates that a fault has been detected.
	In order to verify the advantages of the VBSPCA method, we compare the variational Bayesian dictionary learning (VBDL) method with the DiCCA method. 
	The VBDL obtains sparsity by constructing a dictionary matrix, but does not consider the dynamic relationships<cit.>. 
	The DiCCA method analyses the dynamic relationships of the latent variables but lacks robustness to noise perturbations. 
	 Fault 2 is heavily contaminated by noise, and three different Bayesian sparsity methods effectively detect the fault, but the DiCCA method does not detect all faults. 
	 In addition, the Laplace prior-based VBSPCA method has a better fault detection ability compared to the  Gaussian prior-based VBSPCA method.
	
	
		
	

	The fault diagnosis results based on T_2 and SPE indices are shown in Figures <ref>, respectively. The horizontal axis represents the sample index, while the vertical axis represents the variable index. The RBC for each variable and sample is represented by the color depth, where deeper colors indicate higher RBC.
	
	The VBSPCA method based on both Gaussian prior and Laplace prior can diagnose the fault.
	The fault diagnosis results indicate that variable 24 plays a crucial role in the early stage of the fault, while variables 27 and 28 are the main factors in the later stage of the fault. 
	In addition, the Laplace prior-based VBSPCA suggests that the first five variables also contribute to the fault.
	
			
	
	

	
	

	

	
	Fault 6 is a sharp rise or fall in the hydrogen or oxygen level.
	If there is a sudden change in the hydrogen level, it may result in hydrogen leakage, which can cause a hydrogen explosion accident, leading to casualties and equipment damage. 
	There are several factors that can cause changes in the hydrogen level during the production process, including blockage of the balance valve or pipeline, and blockage of the outlet of the separation and defogging device. 
	Additionally, damage to the diaphragm inside the electrolytic cell due to long-term operation can also result in improper collection of hydrogen or oxygen, leading to a sudden change in the hydrogen or oxygen liquid level.
	
	

	The predicted values of the latent variables calculated by the sparse VAR model are shown in Figure <ref>, where the five dynamic latent variables have been able to characterize the full range of changes in the system. 
	It is also evident from the fault detection results that the first dynamic latent variable and the detection results exhibit the highest degree of similarity.
	The correlations of the obtained dynamic latent variables were checked by plotting correlation plots. The correlations among the dynamic latent variables are presented in Figure <ref>. As depicted, it is evident that the correlations have been eliminated through modeling, whereby the first five dynamic latent variables are extracted.
	
	
		
	
	
			
	
		The results of the fault detection are shown in Figure <ref>.
	Fault 6 is also heavily contaminated by noise, and three different Bayesian sparsity methods effectively detect the fault, but the DiCCA method does not detect all faults. 
	VBDL fails to obtain sufficiently good detection results. This also shows the importance of dynamic analysis for fault detection.
	The Laplace prior-based VBSPCA method has a better fault detection ability compared to the  Gaussian prior-based VBSPCA method.
	This suggests that the Lasso sparsity corresponding to the Laplace prior is more favourable for feature selection and noise reduction.

		


	
	The fault diagnosis results based on T_2 and SPE indices are shown in Figures <ref>.
	The fault diagnosis results indicate that in this case, fault 2 is primarily caused by variable 8, which makes an absolute contribution. It should be noted that under different conditions, this fault may be caused by different factors. Additionally, the contribution plot based on SPE shows that variables 18 to 21 also contribute to the fault. 
	The Gaussian prior-based VBSPCA only considers variable 8 as contributing in T_2. But the Laplace prior-based VBSPCA assumes that there are many variables that contribute to the fault, although not continuously.
	
	
				
	
	

	

	
	





§ CONCLUSIONS
 


Hydrogen production from electrolytic water is a green method of hydrogen production, which is the main way to produce hydrogen from renewable energy such as photovoltaic and wind power. 
With the rapid development of renewable energy and the update and iteration of material technology, electrolytic water hydrogen production technology has a broad development prospect. 
However, the safety of the hydrogen production process is often overlooked. 
The process monitoring of electrolytic water hydrogen production refers to the real-time measurement and control of each parameter in the electrolytic water hydrogen production system to ensure the safe, stable and efficient operation of the system.
The industrial environment causes the AWE process data to be randomly contaminated by noise.
In this study, VBSPCA approach is developed to detect fault of AWE industrial process.
We derive two sparse PCA methods in the framework of variational inference. 
And it is shown that Laplace prior and Gaussian prior correspond to ℓ_1 regularisation and ℓ_2 regularisation. 
Constructing a sparse projection matrix is more conducive to feature selection and improves the robustness of the model to noise.
The proposed method can detect faults in the production environment of industrial noise and has the capability to dynamically analyze data. 
The established dynamic process monitoring scheme can effectively extract the dynamic relationship of the data and obtain reliable detection and diagnosis results.
Therefore, the proposed solution is attractive for green hydrogen production applications. 









elsarticle-num


 






































