Data-Driven Fault Detection for HVAC Control Systems in Pharmaceutical Manufacturing Workshops

1. Introduction

In pharmaceutical manufacturing, the environmental conditions of production workshops—including temperature, humidity, cleanliness, and pressure differentials—directly affect the quality and safety of pharmaceutical products [1,2]. Different types of pharmaceuticals impose distinct requirements on these environmental parameters. To ensure that production environments meet stringent regulatory standards, it is essential to employ reliable HVAC (heating, ventilation, and air conditioning) control systems to regulate the operation of all critical equipment [3]. Consequently, the adoption of effective monitoring and fault diagnosis techniques is of great significance for maintaining high product quality and production yield.

The HVAC control system in pharmaceutical workshops typically adjusts chilled water valves, hot water valves, humidification valves, and other equipment to maintain air parameters (e.g., temperature, humidity, and pressure) at their setpoints [4]. Compared to physical-model-based or empirical-model-based methods, data-driven models offer the advantages of not requiring detailed knowledge of system structure and exhibiting superior generalization capability [5,6]. Modern HVAC systems are often equipped with SCADA (Supervisory Control and Data Acquisition) platforms, which integrate extensive sensor networks to collect and archive large volumes of operational data. Proper utilization of these data enables the monitoring of equipment states, and data-driven state variable models constructed from historical data under normal operating conditions can facilitate fault detection by analyzing the residuals between predicted and measured values [7].

Artificial neural networks, particularly those employing backpropagation algorithms, have been widely applied in the fault detection and diagnosis of HVAC systems and related equipment [8]. Some studies have integrated backpropagation neural networks with decision tree methods to enhance the detection and classification of both known and unknown faults [9,10]. The introduction of algorithms such as Round Robin Dithering (RDP) has further simplified neural network models for fault diagnosis, improving diagnostic efficiency without sacrificing accuracy [11]. While neural networks possess powerful learning capabilities [8], their diagnostic performance heavily depends on the quality and quantity of training data, and the resulting models often lack interpretability and extensibility, limiting their practical application.

Nonlinear State Estimation Technology (NSET), originally proposed by Singer et al. [9], is a data-driven, nonparametric, and nonlinear empirical modeling method. It has been successfully applied to early fault warning in wind turbine generators [10], and further improvements have been achieved by integrating additional fault detection algorithms to enhance prediction accuracy [11]. However, the effectiveness of NSET relies on the selection of a representative process memory matrix constructed from historical data. This approach has notable limitations: (i) expert-selected samples may lack representativeness, and (ii) selecting too many historical samples can result in an excessively large memory matrix [12].

To address these challenges, this work proposes an integrated approach combining Principal Component Analysis (PCA) with NSET [13]. The main idea is to utilize PCA for data cleaning and dimensionality reduction of SCADA-acquired datasets, thereby constructing a temperature prediction model for the HVAC system under normal operation [14]. This enables real-time and accurate fault diagnosis, facilitating the timely identification of potential failures and the implementation of preventive maintenance strategies to ensure the reliable and safe operation of the system.

2. Materials and Methods

2.1. Data Preprocessing

When the number of measurements is sufficiently large, it can be assumed that the sample data follow a normal distribution [15]. Assuming the data contain only random errors, the standard deviation of the dataset can be calculated to define a confidence interval. Any measurement error that falls outside this interval is considered a gross error, which can be identified as an outlier or abnormal data point [16].

Let a set of measurements be denoted as$x_i\left(i=1,2,\dots ,n\right)$, The mean $\overline{x}$ and deviation $v_i=x_i-\overline{x}$ of the dataset are calculated. If, for a particular data point $x_k$, the deviation $v_k$ $\left(1\le k\le n\right)$ satisfies the following condition:

$$\left|v_k\right|=\left|x_k-\overline{x}\right|>3\sigma$$

(1)

then the data point$x_k$is considered to be an outlier.

Prior to further processing, the raw data should be normalized to eliminate the influence of different units and scales [17]. Z-score normalization is a commonly used method for standardizing data. This approach standardizes the original data based on its mean and standard deviation, which alters the scale of the data but does not change the type of its distribution [18]. Z-score normalization is particularly suitable for situations where the maximum and minimum values of the dataset are unknown.

$$x_n{\left(m\right)}^{\prime }={\displaystyle \frac{x_n\left(m\right)-\mu }{\sigma }}$$

(2)

In Equation (2),$x_n{\left(m\right)}^{\prime }$denotes the standardized data,$\mu$represents the mean of the sample data, and$\sigma$is the standard deviation of the sample data. After standardization, the processed data follow a distribution with a mean of 0 and a standard deviation of 1.

2.2. Feature Selection and Principal Component Analysis

The object of this study is the HVAC system of a pharmaceutical manufacturing facility located in Zhejiang Province, China. This system controls three separate rooms: the weighing room (w-room) for measuring raw materials, the preprocessing room (p-room) for preliminary handling of ingredients, and the mixing room (m-room) for final product manufacturing. To develop the NSET model for temperature prediction in the mixing room, it is necessary to identify the parameters in the SCADA system that are most closely related to the mixing room temperature [18,19]. The parameters recorded by the SCADA system, along with their Pearson correlation coefficients with the mixing room temperature and their respective operating ranges, are summarized in Table 1.

As shown in Table 1, the mixing room temperature exhibits a strong positive correlation with both the weighing room temperature and the preprocessing room temperature. This is primarily because all three rooms are regulated by the same HVAC system. However, it should be noted that the Pearson correlation coefficient has certain limitations: it is only suitable for describing linear relationships, is sensitive to outliers, and does not provide information about causal relationships [20]. Therefore, it is necessary to combine other methods for more robust and comprehensive variable selection in the modeling process.

Principal Component Analysis (PCA) is a widely used dimensionality reduction algorithm [21,22]. The fundamental idea of PCA is to transform a set of correlated variables into a smaller number of uncorrelated composite variables, known as principal components, through linear combination. These principal components are designed to capture as much of the useful information from the original data as possible [23]. The specific steps for determining the principal components are as follows:

Suppose there are$n$observed variables, denoted as$X_i=\left(x_{1i},x_{2i},\dots ,x_{Ni}\right)$,$i=1,2,\dots ,n$. Let the correlation coefficient between variables$X_s$and$X_t$$\left(s,t=1,2,\dots ,n\right)$be$r_{st}$. The principal steps for performing PCA are as follows:

1. Data preprocessing. Prepare the raw data for analysis, typically involving normalization or standardization.

2. Calculation of correlation coefficients. Compute the pairwise correlation coefficients among all variables to obtain the correlation coefficient matrix$R=\left(r_{st}\right),s,t=1,2,\dots ,n$.

3. Eigenvalue and eigenvector computation. Calculate the eigenvalues${\lambda }_i$ and corresponding eigenvectors$e_i=\left(e_{i1},e_{i2},\dots ,e_{ip}\right)$of the correlation coefficient matrix$R$. The eigenvectors represent the directions of the principal components.

4. Calculation of variance contribution rates. Determine the variance contribution rate$Vcr$ of each principal component, as well as the cumulative variance contribution rate$Cvcr$of the first$l$ principal components.

$$Vcr={\displaystyle \frac{{\lambda }_i}{{\displaystyle \sum _{i=1}^n{\lambda }_i}}},Cvcr={\displaystyle \frac{{\displaystyle \sum _{i=1}^l{\lambda }_i}}{{\displaystyle \sum _{i=1}^n{\lambda }_i}}}$$

(3)

Principal component analysis was performed on the parameters recorded by the SCADA system using Equation (3). The contribution rates and cumulative contribution rates of all principal component vectors were calculated and then ranked in descending order of their contribution rates, as illustrated in Figure 1.

The variance contribution rates and cumulative variance contribution rates of the principal components for all operating parameters are systematically summarized in Table 2. These values were calculated based on the results of the PCA conducted on the dataset acquired from the SCADA system [25]. As shown, each operating parameter exhibits a distinct contribution to the overall variance, and the cumulative contribution rates demonstrate how much of the total information from the original data is retained as additional components are included. This comprehensive summary provides an essential foundation for subsequent variable selection and model construction in this study.

As shown in Table 2, the cumulative contribution rate of the first six principal components reaches 98.4%. Therefore, selecting the first six principal components ensures that a sufficient amount of the original data’s information is retained while also maintaining high computational efficiency [26].

Based on the above analysis, the variables selected for model construction are as follows: mixing room temperature, preprocessing room temperature, weighing room temperature, supply air duct humidity, mixing room humidity, and fan frequency [25].

The temperature of the air conditioning system is affected by seasonal variations; therefore, the modeling period should not be excessively long. In this study, the model was developed using operational data collected over a three-month period, specifically from 00:50 on December 22, 2023, to 16:20 on March 22, 2024. The dataset consists of 13,198 records, with a sampling interval of 10 minutes.

The variation in the mixing room temperature during this period is illustrated in Figure 2.

As shown in Figure 2, the mixing room temperature is noticeably lower during the period corresponding to data points 6420 to 7676 compared to other times. This is because the air conditioning unit was not in operation during this interval, resulting in the room temperature reflecting the ambient (uncontrolled) conditions. This observation is further supported by the variation in the supply fan frequency of the air conditioning unit, as depicted in Figure 3.

In addition to the non-operational period between data points 6420 and 7676, there are also isolated instances at other times where the supply fan frequency is zero. These data points are considered abnormal and are treated as outliers. Therefore, thorough data preprocessing is required before constructing the temperature prediction model for the air conditioning system [26].

2.3. Nonlinear State Estimation Modeling

In this study, a nonlinear state estimation (NSET) method is employed to develop a temperature prediction model for the air conditioning system under normal operating conditions [27,28]. The established model is then used to predict the system’s output, and the residuals—defined as the differences between the measured and predicted values—are analyzed. The magnitude, range, and variation of these residuals provide important information for assessing the operational status of the HVAC unit [29]. When the system experiences an abnormal condition, its dynamic characteristics deviate from those observed during normal operation, resulting in increased residuals. This forms the basis for early fault detection and warning for the air conditioning unit or related equipment.

At a specific time$n$, the air conditioning unit collects$i$interrelated variables, which are represented as an observation vector [30].

$$X=\left[x_1,x_2,\dots ,x_n\right]$$

(4)

The system observation matrix$P_{n\times m}$can thus be expressed as follows:

$$P_{n\times m}=\left[X\left(1\right),X\left(2\right),\dots ,X\left(m\right)\right]$$

(5)

In Equation (5),$m$denotes the number of observation vectors, and$n$represents the number of variables in each observation vector.

The process memory matrix$D$serves as a memory and representation of the system’s normal operating conditions [31]. It is constructed by selecting$k$observation vectors under different normal operating scenarios.

$$D=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \dots & x_{1k}\\ x_{21}& x_{22}& \dots & x_{2k}\\ ⋮& ⋮& ⋮& ⋮\\ x_{n1}& x_{n2}& \dots & x_{nk}\end{array}\right]$$

(6)

Each column of the process memory matrix$D$represents an observation vector corresponding to a specific normal operating condition of the system. Therefore, it is essential to carefully select the data before constructing the matrix, ensuring that the$k$historical observation vectors used in$D$adequately represent the full range of normal operating scenarios [32]. Ideally, the process memory matrix$D$should encompass all typical normal operating states of the air conditioning system.

The input to the NSET model is an observation vector$X_{obs}$at a given time, and the output is the corresponding predicted vector$X_{est}$. For any input observation vector$X_{obs}$, the NSET model generates a$k$-dimensional weight vector $W$[33]. The vector $W$is a column vector whose dimension is equal to the number of components in the input observation vector$X$.

$$W={\left[\begin{array}{cccc}{\omega }_1& {\omega }_2& \dots & {\omega }_k\end{array}\right]}^{{\rm T}}$$

(7)

The predicted vector$X_{est}$can be expressed as follows:

$$X_{est}=D\cdot W={\omega }_1\cdot X\left(1\right)+\dots +{\omega }_k\cdot X\left(k\right)$$

(8)

The weight vector $W$is obtained by minimizing the residual between the predicted vector$X_{est}$and the observation vector$X_{obs}$[34]:

$$W={\left(D^{{\rm T}}\otimes D\right)}^{-1}\cdot \left(D^{{\rm T}}\otimes X_{obs}\right)$$

(9)

In Equation (9),$\otimes$denotes a nonlinear operator. There are various choices for the nonlinear operator; in this study, the Euclidean distance between$D^{{\rm T}}$and$D$is used, as follows:

$$\otimes \left(X,Y\right)=\sqrt{{\displaystyle \sum _{i=1}^n{\left(x_i-y_i\right)}^2}}$$

(10)

By substituting Equation (9) into Equation (8), the prediction vector of the NSET model can be obtained as follows:

$$X_{est}=D\cdot {\left(D^{{\rm T}}\otimes D\right)}^{-1}\left(D\otimes X_{obs}\right)$$

(11)

Using Equation (11), the predicted value of the current input data can be obtained from historical data, which ultimately provides diagnostic information regarding potential faults in the system [37].

3. Results and Discussion

3.1. Data Preprocessing Effects

The NSET model prediction residual for the mixing room temperature is given by [35]:

$$\epsilon =x_t-{\widehat{x}}_t$$

(12)

In Equation (11), $x_t$ represents the mixing room temperature component in the input observation vector for the NSET model, while ${\widehat{x}}_t$ denotes the corresponding predicted value of the mixing room temperature. In the training set, the residuals of the NSET models for the mixing room temperature are compared for both the raw (unprocessed) data and the preprocessed data. The comparison results are shown in Figure 4.

As shown in Figure 4, except for a few time points, the residuals of the mixing room temperature predicted by the model built with preprocessed data are generally smaller than those obtained with unprocessed data. This demonstrates that data preprocessing is both necessary and effective for improving the model’s prediction accuracy.

3.2. Model Comparison and Fault Detection Performance

Experimental data from one week in March were selected to evaluate the performance of the air conditioning control system models [36,37]. A comparison between the conventional NSET model and the proposed PCA-NSET model clearly demonstrates the superior accuracy of the latter, as illustrated in Figure 5. As shown, the residuals and error fluctuations of the PCA-NSET model are significantly smaller than those of the conventional NSET model, indicating that the PCA-NSET model achieves higher prediction accuracy and better overall performance.

Experimental validation was conducted using data from March 14 to March 22. During this eight-day period, the air conditioning unit operated normally for the first seven days, while a fault occurred on the eighth day. The established PCA-NSET model was applied to perform fault analysis on this dataset. The threshold for residuals can be determined either empirically by the operators or based on the Laida criterion described in Section 2.1. The predicted mixing room temperature and the corresponding temperature residuals for both the first seven days and the eighth day are shown in Figure 6 and Figure 7, respectively.

As shown in Figure 6, when the air conditioning control system is operating normally, the model’s predicted values closely match the observed values, with residuals remaining within the threshold except for a few isolated points. In contrast, Figure 7 demonstrates that when a fault occurs in the system, there is a significant deviation between the predicted and observed values. The residuals rapidly increase and consistently exceed the threshold within a short period. Therefore, by monitoring the residuals, effective fault prediction can be achieved; specifically, frequent exceedance of the threshold by the residuals can serve as an early warning indicator for system faults.

4. Conclusions

This study investigates the operational status and fault detection of air conditioning units based on SCADA monitoring data. After performing data preprocessing, a predictive model was developed using a combination of Principal Component Analysis (PCA) and Nonlinear State Estimation Technology (NSET), and its performance was compared with that of the conventional NSET model. The results demonstrate that the PCA-NSET model achieves higher accuracy than the conventional NSET model, verifying the effectiveness of the proposed approach. Further analysis of the SCADA data from the air conditioning unit shows that, under normal operating conditions, the model’s residuals remain within a reasonable range, whereas during a fault, the residuals consistently exceed the threshold. Experimental results confirm that the proposed PCA-NSET model provides excellent performance in detecting faults in air conditioning control systems.