Dynamic Calibration Method of Sensor Drift Fault in HVAC System Based on Bayesian Inference

Sensor drift fault calibration is essential to maintain the operation of heating, ventilation and air conditioning systems (HVAC) in buildings. Bayesian inference (BI) is becoming more and more popular as a commonly used sensor fault calibration method. However, this method focused mainly on sensor bias fault, and it could be difficult to calibrate drift fault that changes with time. Therefore, a dynamic calibration method for sensor drift fault of HVAC systems based on BI is developed. Taking the drift fault calibration of the chilled water supply temperature sensor of the chiller as an example, the performance of the proposed dynamic calibration method is evaluated. Results show that the combination of the Exponentially Weighted Moving-Average (EWMA) method with high detection accuracy and the proposed BI dynamic calibration method can effectively improve the calibration accuracy of drift fault, and the Mean Absolute Percentage Error (MAPE) value between the calibrated and normal data is less than 5%.


Introduction
Building energy consumption is huge, which is close to 40% of the global energy consumption. In the next 20 years, this proportion will increase by more than 1.5% every year [1]. As the largest energy consumer in the construction field, the energy consumption of HVAC systems accounts for more than 40% of total building energy consumption [2]. Therefore, improving the energy efficiency of HVAC systems is one of the key ways to reduce building energy consumption. At present, energy-saving technologies such as system operation diagnosis, building automatic control and building energy management [3] have been widely used to improve the energy efficiency of HVAC systems. The application of these technologies is based on the measurement data of various sensors in the system [4]. However, with the passage of time, the sensors used to detect and record system operation data may inevitably suffer various types of faults in their whole life cycles. Sensor faults not only affect the energy-saving technology, increase the energy consumption of HVAC systems [5,6], but also reduce indoor thermal comfort directly [7,8]. Among them, as one of the common sensor fault types, drift fault is very harmful to the system because its fault offset changes with time. Hence, the sensor drift fault calibration of HVAC system sensors is very important.
In practical application, the field calibration of sensors after a fault usually includes two steps: fault detection and fault calibration. Sensor fault detection methods are mainly divided into three categories: quantitative model-based methods, qualitative model-based methods and data-driven methods. In recent years, due to the accelerated upgrading of computer software and hardware and the rapid development of data storage and transmission technology, data-driven fault detection methods have attracted more and time-independent fixed value, which may be very difficult to solve the problem that the sensor offset changes with time. So This means the original BI may be not suitable for directly application on sensor drift fault calibration for HVAC systems.
Based on the aforementioned literature review, this study presents a dynamic calibration method of sensor drift fault in an HVAC system based on Bayesian inference. In this method, the time characteristic variable is introduced to deal with the dynamic characteristics of the drift, and then the dynamic calibration of the BI based sensor drift fault is realized using both simulation and practical cases, combined with the data-driven fault detection algorithm. Taking the drift fault calibration of the chilled water supply temperature sensor of the chiller as an example, the performance of the proposed dynamic calibration method is evaluated by introducing different drift fault amplitudes, and the influence of time length on the calibration accuracy is further studied. The main contributions of this paper could be: (1) introducing time characteristic variables into Bayesian inference to deal with the dynamic characteristics of drift so as to realize the dynamic calibration of sensor drift fault; (2) combining the proposed dynamic calibration method with EWMA method with high detection accuracy, the time from fault to fault calibration is effectively shortened, and the calibration accuracy of sensor drift fault is improved; (3) the relationship between calibration time length and calibration accuracy is given.
The paper is organized as follows: Section 2 introduces three detection algorithms and the basic principle of the dynamic calibration method based on Bayesian inference. Section 3 introduces the research framework of this paper. Section 4 gives a description of the system used for the case study. Section 5 presents the calibration effect of the proposed dynamic calibration method in case application. Section 6 compares the dynamic calibration method based on Bayesian inference with the principal component analysis method. Section 6 compares this method with the least square method. Section 7 is the conclusion.

Principles of Three Sensor Fault Detection Methods
This study adopts the Laida criterion, Box-plot and EMWA control chart detection algorithm. These are implemented on the python platform.

Laida Criterion
The Laida criterion [12] takes the standard deviation δ of three times the normal value as the selection standard. The detection algorithm has the advantages of a simple model, small computation and good fault detection performance. Equation (1) for standard deviation δ is as follows, it is the parameter calculated after a large number of repeated observations and its calculation Equation (1) is as follows: In Equation (1), n is the number of observed values in the HVAC system, Y n represents normal data and Y is the arithmetic mean of normal data.
The calculation of its threshold is shown in Equation (2): where UCL Laida and LCL Laida are the upper and lower limits of the Laida criterion, respectively. If the HVAC system test data is greater than UCL Laida or less than LCL Laida , the sensor has failed.

Box-Plot
Box-plot [13] is an effective graph for processing data. Its algorithm is efficient. This method uses all kinds of quantiles (upper quartile, median and lower quartile) of the  (3) and (4).
where, Y 3 and Y 1 are the upper quartile and lower quartile of normal data of the HVAC and R systems in ascending order, respectively. If the test data of the HVAC system is greater than UCL Box−plot or less than LCL Box−plot , it indicates that the sensor has failed.

EWMA Control Diagram
EWMA control chart [15] is widely used to monitor the operation process of industrial equipment because of its simple detection model, easy construction and small amount of calculations. For the system operation process data, the EWMA control chart algorithm allocates the corresponding weight factors on the basis of fully considering the importance of data at the current time and the previous time and constructs the corresponding upper or lower limit of control limit to effectively detect abnormal small fluctuations. Its calculated value represents the time weight mean of all data in the monitoring process. EWMA calculated values are usually defined as: In Equation (5) Z m is the m-th EWMA calculated value, Z m−1 is the m-1st EWMA calculated value, y n is the nth data, and N represents the number of iterations, λ Is the weight factor (0 ≤ λ ≤ 1). Usually, a smaller λ will lead to smaller displacement and a faster detection process [29]. The best λ is usually between 0.2 and 0.3 [30]; in this study, λ = 0.2, and the number of iterations is 6305. The control limit definition of the EWMA control chart is shown in Equation (6): In Equation (6), factor L represents the width of the control limit, L = 3. δ is the standard deviation of data, and n represents the number of data. When the observed value y is greater than UCL EW MA or less than LCL EW MA , it indicates that the numeric is a fault numeric.

Dynamic Calibration Method Based on Bayesian Inference
Bayesian inference [31][32][33] is a statistical inference method that updates the probability of specific assumptions according to new information, which is widely used in many fields. In the building system, by comparing the difference between the data of the sensor under normal conditions and the data under fault conditions, it obtains the required calibration results through a series of formula derivations.
After the sensor has a bias fault, the offset between the normal value and the fault numeric is a fixed value, which can be obtained by the Bayesian inference method. However, after the sensor has drift fault, the offset between the normal value and the fault numeric will change with time, while the traditional Bayesian method cannot get the dynamic offset and cannot solve the drift fault. Therefore, this study introduces a time characteristic based on Bayesian inference and proposes a dynamic calibration method of drift fault based on Bayesian inference. The details are as follows: The error is between normal value and fault numeric of sensor x = (x 1 , x 2 ..., x i ), where x i represents the offset at time i, and π (x) represents the prior function of error x; a priori is the first guess of x probability. The system model g (x) = (x i , Y) is constructed, and the likelihood Equation (7) in the posterior distribution function is obtained. The normalized  [34] algorithm so as to obtain the posterior distribution Equation (10) and the maximum posterior distribution value x m . The time t m at which x m is located is obtained by Equation (10). The drift slope k is obtained by Equation (11), and then the drift fault calibration value is obtained by Equation (12): where, t 0 is the time when the fault is detected.
where T i0 is the calibration value at time i and T i1 is the fault numeric at time i.

Lag Time and Calibration Time
The trend of drift fault is that zero increases gradually. Any actual fault detection method cannot be detected part of the time, which leads to a lag between fault detection and fault occurrence. The calibration fault needs fault data, so the calibration process will also lag behind the fault detection. Figure 1 shows the drift fault results of 0.12 • C/h, 0.18 • C/h and 0.24 • C/h under the EWMA control chart. In Figure 1, τ 1 is the time point when the fault occurs, at which time the data begins to drift. τ 2 is the time point when the fault of 0.24 • C/h is detected, and the period between τ 1 and τ 2 is the time when the fault occurs but is not detected, which is called the lag time t 1 . τ 3 is the time to start calibration. Fault calibration occurs by modeling the fault data of τ 2 and τ 3 , which is called calibration time t 2 . ; a priori is the first guess of probability. The system model ( ) = ( , ) is constructed, and the likelihood Equation (7) in the posterior distribution function is obtained. The normalized function ( ) of the fault numeric in Equation (8) is sampled by the the Markov chain Monte Carlo (MCMC) [34] algorithm so as to obtain the posterior distribution Equation (10) and the maximum posterior distribution value . The time at which is located is obtained by Equation (10). The drift slope k is obtained by Equation (11), and then the drift fault calibration value is obtained by Equation (12): = ( ) (10) where, is the time when the fault is detected.
where is the calibration value at time and is the fault numeric at time .

Lag Time and Calibration Time
The trend of drift fault is that zero increases gradually. Any actual fault detection method cannot be detected part of the time, which leads to a lag between fault detection and fault occurrence. The calibration fault needs fault data, so the calibration process will also lag behind the fault detection. Figure 1 shows the drift fault results of 0.12 °C/h, 0.18 °C/h and 0.24 °C/h under the EWMA control chart. In Figure 1, is the time point when the fault occurs, at which time the data begins to drift.
is the time point when the fault of 0.24 °C/h is detected, and the period between and is the time when the fault occurs but is not detected, which is called the lag time .
is the time to start calibration. Fault calibration occurs by modeling the fault data of and , which is called calibration time .  Figure 2 shows the main framework of this study. The research framework of this study is divided into four steps: data preparation, fault detection, fault calibration and discussion of influencing factors.  Figure 2 shows the main framework of this study. The research framework of this study is divided into four steps: data preparation, fault detection, fault calibration and discussion of influencing factors. •

Data preparation
The chilled water supply temperature sensor in the chiller is selected as the target sensor under both simulation and actual conditions. The detection algorithms are the Laida criterion, Box-plot chart and the EWMA control chart, and some or all of the data in May, June and July are selected for modeling.

Fault detection
The time point of fault is detected by three detection methods under different modeling data, and the detection algorithm with the shortest lag time and modeling data are selected to detect the fault. •

Fault calibration
From a practical point of view, the mode of detection before calibration is adopted. After detection by the detection algorithm obtained in the previous section, the fault is calibrated by the dynamic calibration method to obtain the calibration value. The calibration slope and the MAPE value between the normal value and calibration value are calculated. •

Influence factor
The factors influencing the calibration accuracy are discussed as follows: (1) different actual working conditions; (2) different calibration times; (3) different fault amplitudes; (4) different sampling intervals.
The calculation formula of MAPE is shown in Equations (13) and (14): •

Data preparation
The chilled water supply temperature sensor in the chiller is selected as the target sensor under both simulation and actual conditions. The detection algorithms are the Laida criterion, Box-plot chart and the EWMA control chart, and some or all of the data in May, June and July are selected for modeling.

• Fault detection
The time point of fault is detected by three detection methods under different modeling data, and the detection algorithm with the shortest lag time and modeling data are selected to detect the fault.

• Fault calibration
From a practical point of view, the mode of detection before calibration is adopted. After detection by the detection algorithm obtained in the previous section, the fault is calibrated by the dynamic calibration method to obtain the calibration value. The calibration slope and the MAPE value between the normal value and calibration value are calculated.
The calculation formula of MAPE is shown in Equations (13) and (14): where, Y 0 is fault data, Y cal is calibration data and t 2 is calibration time.

Case 1 Simulated Chiller System
The target building is a typical small office building. Figure 3 is the stereoscopic view of the target building. The single-story building has five hot areas, one core hot area in the middle and four hot areas around. The whole building has a total length of 30.5 m, a total width of 15.2 m, a total height of 3.0 m and a building area of 463.6 m 2 .
where, is fault data, is calibration data and is calibration time.

Case 1 Simulated Chiller System
The target building is a typical small office building. Figure 3 is the stereoscopic view of the target building. The single-story building has five hot areas, one core hot area in the middle and four hot areas around. The whole building has a total length of 30.5 m, a total width of 15.2 m, a total height of 3.0 m and a building area of 463.6 m . The cold source of the building's central air conditioning is a chiller. Figure 4 shows the system diagram of the chiller. As shown in the figure, the target sensor is the chilled water supply temperature sensor. In order to better meet the actual working conditions, the working time of the simulation system is set as 7:00-18:00 of the working day, and the meteorological parameters are taken from Wuhan city. We selected part of the data of May, June and July as the normal data. The data on August 1 was taken as the test data. The fault occurred at 12:00 noon on August 1. The drift slopes were 0.12 °C/h (1.5% of the average value of the sensor under normal working conditions), 0.18 °C/h (2.25%) and 0.24 °C/h (3%), respectively. The sampling interval was once every 5 min.

Case 2 Actual Chiller System
The data of case 2 is the data of an actual chiller system collected in an industrial plant along the east coast of China from May to August 2018. Some data from May, June and July 2018 were selected as normal data. The data on 6 August 2018 was used as test data. The target sensor was also a chilled water supply sensor. The working hours of the system were 9:00-18:00. The fault occurred at 12:00 noon on 1 August 2018. The drift slopes were 0.12 °C/h, 0.18 °C/h and 0.24 °C/h, respectively. Data were recorded every 5 min. The cold source of the building's central air conditioning is a chiller. Figure 4 shows the system diagram of the chiller. As shown in the figure, the target sensor is the chilled water supply temperature sensor. In order to better meet the actual working conditions, the working time of the simulation system is set as 7:00-18:00 of the working day, and the meteorological parameters are taken from Wuhan city. We selected part of the data of May, June and July as the normal data. The data on August 1 was taken as the test data. The fault occurred at 12:00 noon on August 1. The drift slopes were 0.12 • C/h (1.5% of the average value of the sensor under normal working conditions), 0.18 • C/h (2.25%) and 0.24 • C/h (3%), respectively. The sampling interval was once every 5 min.
where, is fault data, is calibration data and is calibration time.

Case 1 Simulated Chiller System
The target building is a typical small office building. Figure 3 is the stereoscopic view of the target building. The single-story building has five hot areas, one core hot area in the middle and four hot areas around. The whole building has a total length of 30.5 m, a total width of 15.2 m, a total height of 3.0 m and a building area of 463.6 m . The cold source of the building's central air conditioning is a chiller. Figure 4 shows the system diagram of the chiller. As shown in the figure, the target sensor is the chilled water supply temperature sensor. In order to better meet the actual working conditions, the working time of the simulation system is set as 7:00-18:00 of the working day, and the meteorological parameters are taken from Wuhan city. We selected part of the data of May, June and July as the normal data. The data on August 1 was taken as the test data. The fault occurred at 12:00 noon on August 1. The drift slopes were 0.12 °C/h (1.5% of the average value of the sensor under normal working conditions), 0.18 °C/h (2.25%) and 0.24 °C/h (3%), respectively. The sampling interval was once every 5 min.

Case 2 Actual Chiller System
The data of case 2 is the data of an actual chiller system collected in an industrial plant along the east coast of China from May to August 2018. Some data from May, June and July 2018 were selected as normal data. The data on 6 August 2018 was used as test data. The target sensor was also a chilled water supply sensor. The working hours of the system were 9:00-18:00. The fault occurred at 12:00 noon on 1 August 2018. The drift slopes were 0.12 °C/h, 0.18 °C/h and 0.24 °C/h, respectively. Data were recorded every 5 min.

Case 2 Actual Chiller System
The data of case 2 is the data of an actual chiller system collected in an industrial plant along the east coast of China from May to August 2018. Some data from May, June and July 2018 were selected as normal data. The data on 6 August 2018 was used as test data. The target sensor was also a chilled water supply sensor. The working hours of the system were 9:00-18:00. The fault occurred at 12:00 noon on 1 August 2018. The drift slopes were 0.12 • C/h, 0.18 • C/h and 0.24 • C/h, respectively. Data were recorded every 5 min. Data-driven detection algorithm: EWMA control chart, Laida criterion method and Box-plot method are widely used in fault detection of HVAC systems. Therefore, these three methods were used to detect sensor data in this study in order to reduce the influence of time lag in modeling and detection. Figure 5 shows the thresholds obtained from different modeling data and different detection algorithms. The minimum threshold was obtained by the EWMA control chart algorithm through one-month modeling, and the maximum threshold was obtained by the Laida criterion algorithm through three-month modeling.

Discussion on Detection Methods
Data-driven detection algorithm: EWMA control chart, Laida criterion method and Box-plot method are widely used in fault detection of HVAC systems. Therefore, these three methods were used to detect sensor data in this study in order to reduce the influence of time lag in modeling and detection. Figure 5 shows the thresholds obtained from different modeling data and different detection algorithms. The minimum threshold was obtained by the EWMA control chart algorithm through one-month modeling, and the maximum threshold was obtained by the Laida criterion algorithm through three-month modeling.

Fault Calibration
According to the previous section, we found that the EWMA control chart method based on one-month modeling had the smallest threshold, so the fault detection time was the shortest compared with other methods. In order to improve the timeliness and accuracy of calibration, the detection method of the modeling data is combined with the dynamic calibration method based on Bayesian inference, and the EWMA control chart is used for calibration after detection. Figure 6 shows the calibration results under three drift slopes. Table 1 shows the drift slope calculated by the dynamic calibration method and the absolute percentage error (MAPE) between the calibrated data and the average. As can be seen from Figure 6, the calibrated data gradually approaches the normal value with increasing time. Finally, it remains the same as the normal data. The larger the drift slope, the smaller the lag time and the shorter the calibration time. Table 1 shows that the calculated slopes under the three drift slopes are about 0.6 °C/h, and the MAPE value increases with the increase of the drift slope.

Fault Calibration
According to the previous section, we found that the EWMA control chart method based on one-month modeling had the smallest threshold, so the fault detection time was the shortest compared with other methods. In order to improve the timeliness and accuracy of calibration, the detection method of the modeling data is combined with the dynamic calibration method based on Bayesian inference, and the EWMA control chart is used for calibration after detection. Figure 6 shows the calibration results under three drift slopes. Table 1 shows the drift slope calculated by the dynamic calibration method and the absolute percentage error (MAPE) between the calibrated data and the average. As can be seen from Figure 6, the calibrated data gradually approaches the normal value with increasing time. Finally, it remains the same as the normal data. The larger the drift slope, the smaller the lag time and the shorter the calibration time. Table 1 shows that the calculated slopes under the three drift slopes are about 0.6 • C/h, and the MAPE value increases with the increase of the drift slope.

Discussion on Detection Methods
Data-driven detection algorithm: EWMA control chart, Laida criterion method and Box-plot method are widely used in fault detection of HVAC systems. Therefore, these three methods were used to detect sensor data in this study in order to reduce the influence of time lag in modeling and detection. Figure 5 shows the thresholds obtained from different modeling data and different detection algorithms. The minimum threshold was obtained by the EWMA control chart algorithm through one-month modeling, and the maximum threshold was obtained by the Laida criterion algorithm through three-month modeling.

Fault Calibration
According to the previous section, we found that the EWMA control chart method based on one-month modeling had the smallest threshold, so the fault detection time was the shortest compared with other methods. In order to improve the timeliness and accuracy of calibration, the detection method of the modeling data is combined with the dynamic calibration method based on Bayesian inference, and the EWMA control chart is used for calibration after detection. Figure 6 shows the calibration results under three drift slopes. Table 1 shows the drift slope calculated by the dynamic calibration method and the absolute percentage error (MAPE) between the calibrated data and the average. As can be seen from Figure 6, the calibrated data gradually approaches the normal value with increasing time. Finally, it remains the same as the normal data. The larger the drift slope, the smaller the lag time and the shorter the calibration time. Table 1 shows that the calculated slopes under the three drift slopes are about 0.6 °C/h, and the MAPE value increases with the increase of the drift slope.   Figure 7 shows the changes in k value and MAPE value under three drift faults, and the black dotted line represents the set slope. It can be seen that the three calibration slopes and MAPE values tend to be stable with the increase in calibration time. The calibration slope and MAPE value fluctuate greatly when the calibration time is short. Furthermore, the change curves of the three drift slopes are not consistent. The calibration slope at 0.12 • C/h decreases gradually with the increasing calibration time. The change curve at 0.18 • C/h rises first and then decreases, while it rises gradually at 0.24 • C/h. All three calibration slopes are above the set slope. The trend of the MAPE curve and drift slope curve under the three fault amplitudes is opposite.  Figure 7 shows the changes in value and MAPE value under three drift faults, and the black dotted line represents the set slope. It can be seen that the three calibration slopes and MAPE values tend to be stable with the increase in calibration time. The calibration slope and MAPE value fluctuate greatly when the calibration time is short. Furthermore, the change curves of the three drift slopes are not consistent. The calibration slope at 0.12 °C/h decreases gradually with the increasing calibration time. The change curve at 0.18 °C/h rises first and then decreases, while it rises gradually at 0.24 °C/h. All three calibration slopes are above the set slope. The trend of the MAPE curve and drift slope curve under the three fault amplitudes is opposite.

Discussion of Sensor Sampling Interval
This section discusses the effect of sensor data sampling interval on calibration accuracy. Under the drift fault of 0.24 °C/h, the sampling interval of the sensor is 5 min, 10 min and 15 min. Figure 8 shows the changes in value and MAPE value at three different time intervals. With the increase in calibration time, the calibration slope and MAPE value under the three time intervals tended to be the same. Under the same calibration time, the fluctuation of the and MAPE values with a five-minute sampling interval is the smallest, followed by the sampling interval of ten minutes. The most volatile is when the sampling interval is 15 min. The shorter the calibration time, the greater the impact of sampling interval. The longer the calibration time, the smaller the impact of sampling interval on calibration.  Figure 9 shows the detection results under three drift slopes. It can be seen that the detection threshold of the EWMA control chart method is smaller than that of the other two methods. At the same time, under the same fault condition, EWMA is the first to detect the fault, with the shortest lag time, followed by the Laida criterion method, and finally the Box-plot method. Under the same detection algorithm, the lag time decreases with the increase of drift slope.  Figure 9 shows the detection results under three drift slopes. It can be seen that the detection threshold of the EWMA control chart method is smaller than that of the other two methods. At the same time, under the same fault condition, EWMA is the first to detect the fault, with the shortest lag time, followed by the Laida criterion method, and finally the Box-plot method. Under the same detection algorithm, the lag time decreases with the increase of drift slope.  Figure 9 shows the detection results under three drift slopes. It can be seen that the detection threshold of the EWMA control chart method is smaller than that of the other two methods. At the same time, under the same fault condition, EWMA is the first to detect the fault, with the shortest lag time, followed by the Laida criterion method, and finally the Box-plot method. Under the same detection algorithm, the lag time decreases with the increase of drift slope.

Influence of Modeling Data on Test Results
According to the discussion in Section 5.2.1, compared with the other two methods, the EWMA control chart has the minimum threshold and the shortest lag time. Therefore, the EWMA control chart algorithm is selected for fault detection. Due to the modeling characteristics, this section selects the amount of modeling data for the EMWA control chart. The data volume of one month in July, two months in June and July and three months in May, June and July are used for modeling. Three different temperature thresholds were obtained. As shown in Figure 10, the detection results of three different modeling data quantities are shown. We found that the threshold obtained by modeling with one month of data is the smallest, followed by the threshold obtained by modeling with two months of data and the threshold obtained by modeling with three months is the largest. The first fault detected is the threshold value modeled from one month's data.

Influence of Modeling Data on Test Results
According to the discussion in Section 5.2.1, compared with the other two methods, the EWMA control chart has the minimum threshold and the shortest lag time. Therefore, the EWMA control chart algorithm is selected for fault detection. Due to the modeling characteristics, this section selects the amount of modeling data for the EMWA control chart. The data volume of one month in July, two months in June and July and three months in May, June and July are used for modeling. Three different temperature thresholds were obtained. As shown in Figure 10, the detection results of three different modeling data quantities are shown. We found that the threshold obtained by modeling with one month of data is the smallest, followed by the threshold obtained by modeling with two months of data and the threshold obtained by modeling with three months is the largest. The first fault detected is the threshold value modeled from one month's data.

Fault Calibration
According to Sections 5.2.1 and 5.2.2, the EWMA control chart method based on one month of modeling is combined with the dynamic calibration method based on Bayesian inference. Figure 11 Table 2 shows the drift slope after calibration and the MAPE values of calibrated data and normal data. Figure 11 shows that the calibrated data are distributed near the normal data. It shows that this method is effective in calibrating drift fault. It is found in Table 2 that as the drift slope increases, the difference between the calculated slope and the set slope becomes smaller. At the same time, it can be seen that the MAPE value between normal data and calibration data is less than 5%.

Fault Calibration
According to Sections 5.2.1 and 5.2.2, the EWMA control chart method based on one month of modeling is combined with the dynamic calibration method based on Bayesian inference. Figure 11 shows the calibration results under three drift slopes, in which (a) is 0.12 • C/h drift fault, (b) is 0.18 • C/h drift fault and (c) is 0.24 • C/h drift fault. Table 2 shows the drift slope after calibration and the MAPE values of calibrated data and normal data. Figure 11 shows that the calibrated data are distributed near the normal data. It shows that this method is effective in calibrating drift fault. It is found in Table 2 that as the drift slope increases, the difference between the calculated slope and the set slope becomes smaller. At the same time, it can be seen that the MAPE value between normal data and calibration data is less than 5%.

Fault Calibration
According to Sections 5.2.1 and 5.2.2, the EWMA control chart method based on one month of modeling is combined with the dynamic calibration method based on Bayesian inference. Figure 11 shows the calibration results under three drift slopes, in which (a) is 0.12 °C/h drift fault, (b) is 0.18 °C/h drift fault and (c) is 0.24 °C/h drift fault. Table 2 shows the drift slope after calibration and the MAPE values of calibrated data and normal data. Figure 11 shows that the calibrated data are distributed near the normal data. It shows that this method is effective in calibrating drift fault. It is found in Table 2 that as the drift slope increases, the difference between the calculated slope and the set slope becomes smaller. At the same time, it can be seen that the MAPE value between normal data and calibration data is less than 5%.

Relationship between Calibration Accuracy and Calibration Time
In order to consider whether the accuracy of calibration is related to the time required for calibration, the curve of drift slope and MAPE changing with calibration time is made. Figure 12 shows the changes in value and MAPE with calibration time under three drift slopes. As can be seen from Figure 12:

Relationship between Calibration Accuracy and Calibration Time
In order to consider whether the accuracy of calibration is related to the time required for calibration, the curve of drift slope k and MAPE changing with calibration time is made. Figure 12 shows the changes in k value and MAPE with calibration time under three drift slopes. As can be seen from Figure 12: • Under the drift slope of 0.12 • C/h, with the increase in calibration time, the k value and MAPE showed a trend of gradually increasing at first and then tending to be stable. The slope is stable at about 0.14 • C/h, and the MAPE value is stable at about 2%.

•
The k value of 0.18 • C/h increases with time, and the value gradually drops to the setting slope. The MAPE value gradually rises to around 4%.

•
The k value of 0.24 • C/h rises first, then decreases, and finally tends to be near the set slope. MAPE first decreased, then increased, and finally stabilized at around 4.5%. It can be seen from the figure that when the calibration time is short, the value and MAPE value obtained fluctuate greatly, and the value is unstable. If high calibration accuracy is required, the calibration will take a long time.
At the same time, in order to explore the change in the calibration accuracy of the same sensor with the calibration data under different working conditions, the same type of sensor drift fault of a metro chiller is selected for verification, and the sampling interval is once every 5 min. Figure 13 shows the changes in value and MAPE value under a subway. It can be found that under the drift slopes of 0.12 °C/h and 0.18 °C/h, the slope after cooling tower data calibration is closer to the set slope than that after subway data calibration in the same calibration time. Under the drift slope of 0.24 °C/h, the slope after Metro data calibration is closer to the set slope. At the same time, the calibration data of Metro is smaller than the calibration data MAPE of the cooling tower in the three cases, indicating that after calibration, the data of Metro is closer to the normal data than the data of the cooling tower.  It can be seen from the figure that when the calibration time is short, the k value and MAPE value obtained fluctuate greatly, and the value is unstable. If high calibration accuracy is required, the calibration will take a long time.
At the same time, in order to explore the change in the calibration accuracy of the same sensor with the calibration data under different working conditions, the same type of sensor drift fault of a metro chiller is selected for verification, and the sampling interval is once every 5 min. Figure 13 shows the changes in k value and MAPE value under a subway. It can be found that under the drift slopes of 0.12 • C/h and 0.18 • C/h, the slope after cooling tower data calibration is closer to the set slope than that after subway data calibration in the same calibration time. Under the drift slope of 0.24 • C/h, the slope after Metro data calibration is closer to the set slope. At the same time, the calibration data of Metro is smaller than the calibration data MAPE of the cooling tower in the three cases, indicating that after calibration, the data of Metro is closer to the normal data than the data of the cooling tower.
is once every 5 min. Figure 13 shows the changes in value and MAPE value under a subway. It can be found that under the drift slopes of 0.12 °C/h and 0.18 °C/h, the slope after cooling tower data calibration is closer to the set slope than that after subway data calibration in the same calibration time. Under the drift slope of 0.24 °C/h, the slope after Metro data calibration is closer to the set slope. At the same time, the calibration data of Metro is smaller than the calibration data MAPE of the cooling tower in the three cases, indicating that after calibration, the data of Metro is closer to the normal data than the data of the cooling tower.

Discussion of Sensor Sampling Interval
This section discusses the influence of sensor data sampling interval on calibration accuracy. Under the drift fault of 0.24 • C/h, the data shall be recorded every 5 min, 10 min and 15 min. Figure 14 shows the changes in k value and MAPE value at three different time intervals. It can be found that under the same calibration time, the fluctuation of k value and MAPE value obtained from the data with a sampling interval of five minutes is the smallest, followed by the case with a sampling interval of ten minutes. The most volatile is when the sampling interval is 15 min. In terms of calibration accuracy, with the increase in calibration time, the calibration results of three datasets with different time intervals are basically the same. It can be found from Figure 14 that the earlier the fault calibration is carried out, the greater the relationship between the accuracy after calibration and the sampling interval, and the longer the calibration time, the smaller the relationship between the two.  Figure 14 shows the changes in value and MAPE value at three different time intervals. It can be found that under the same calibration time, the fluctuation of value and MAPE value obtained from the data with a sampling interval of five minutes is the smallest, followed by the case with a sampling interval of ten minutes. The most volatile is when the sampling interval is 15 min. In terms of calibration accuracy, with the increase in calibration time, the calibration results of three datasets with different time intervals are basically the same. It can be found from Figure 14 that the earlier the fault calibration is carried out, the greater the relationship between the accuracy after calibration and the sampling interval, and the longer the calibration time, the smaller the relationship between the two.

Comparison with Least Square Method
The least squares method is often used to calibrate drift fault [35,36]. It finds the best fitting function of the input data set by minimizing the square of error and obtains the

Comparison with Least Square Method
The least squares method is often used to calibrate drift fault [35,36]. It finds the best fitting function of the input data set by minimizing the square of error and obtains the drift slope k.
In order to demonstrate the effectiveness of the dynamic calibration method based on Bayesian inference, this method is compared with the commonly used least squares method. Figure 15 shows the calibration results under the two calibration methods. Table 3 shows the calibration slope under the two methods and the MAPE between the calibrated data and the normal value.
The least squares method is often used to calibrate drift fault [35,36]. It finds the best fitting function of the input data set by minimizing the square of error and obtains the drift slope .
In order to demonstrate the effectiveness of the dynamic calibration method based on Bayesian inference, this method is compared with the commonly used least squares method. Figure 15 shows the calibration results under the two calibration methods. Table  3 shows the calibration slope under the two methods and the MAPE between the calibrated data and the normal value.
As can be seen from Figure 15, the difference between the calibrated data of the least squares method represented by the green line and the normal value is greater than that of the dynamic calibration method represented by the blue line. It shows that the data calibrated by the Bayesian-based dynamic calibration method is closer to the normal value. In Table 3, under the three drift amplitudes, the MAPE values obtained by the dynamic calibration method based on Bayesian inference are all smaller than those obtained by the least squares method. It shows that the dynamic calibration method based on Bayesian inference has better calibration performance than the least square method.  As can be seen from Figure 15, the difference between the calibrated data of the least squares method represented by the green line and the normal value is greater than that of the dynamic calibration method represented by the blue line. It shows that the data calibrated by the Bayesian-based dynamic calibration method is closer to the normal value. In Table 3, under the three drift amplitudes, the MAPE values obtained by the dynamic calibration method based on Bayesian inference are all smaller than those obtained by the least squares method. It shows that the dynamic calibration method based on Bayesian inference has better calibration performance than the least square method.

Conclusions
In this study, a dynamic sensor drift fault calibration method based on Bayesian inference is proposed for HVAC systems. The method has two steps. The first step is to detect the fault with a high-precision algorithm, shorten the lag time and improve the calibration accuracy. In the second step, the drift fault is calibrated by the dynamic calibration method based on Bayesian inference to complete the fault calibration from the field perspective. The proposed method is verified by using the fault dataset of a chilled water supply temperature sensor of simulation and actual chiller. Main conclusions are as follows.
(1) Through the combination of the dynamic calibration method based on Bayesian inference and the fault detection method of the EWMA control chart, the detection accuracy method is better than the other two methods-Laida criterion and Box-plotso as to shorten the lag time and improve the calibration accuracy.
(2) After the dynamic calibration method based on Bayesian inference is used to calibrate various drift faults, the MAPE value between the calibrated data and the normal data under the actual data is less than 5%. (3) When the dynamic method based on Bayesian inference is used to calibrate the drift fault, the shorter the calibration time, the lower the calibration accuracy, and the greater the influence of the sensor sampling interval on the calibration accuracy.