Research on automatic diagnosis and optimal control of faults in housing heating systems

Recurrent neural networks (LSTM) are a class of neural networks with short-term memory that are very suitable for sequence data analysis and modelling. In order to process sequence data and utilise its historical information, the hidden layer of a recurrent neural network employs neurons with self-feedback, which can receive not only the output information of other neurons, but also their output information, thus forming a network structure with loops. Figure 3 shows the structure of a simple recurrent neural network, which contains an input layer, a hidden layer and an output layer. In the figure, vector xR,t∈ℝNR$${x_{R,t}} \in {\mathbb{R}^{{N_R}}}$$ represents the input at the moment of t, variable N_R represents the number of input variables, vector hR,t∈ℝDR$${h_{R,t}} \in {\mathbb{R}^{{D_R}}}$$ represents the state of the hidden layer at the moment of t (referred to as the hidden state), variable D_R represents the number of neurons in the hidden layer, and variable y^t∈ℝ$${\hat y_t} \in \mathbb{R}$$ represents the output at the moment of t. The expression of the simple recurrent neural network is given by. (1)hR,t=fact(WRxR,t+URhR,t−1+bR)$${h_{R,t}} = {f_{act}}\left( {{W_R}{x_{R,t}} + {U_R}{h_{R,t - 1}} + {b_R}} \right)$$ (2)y^t=VRThR,t$${\hat y_t} = V_R^T{h_{R,t}}$$

Figure 3.

Architecture of recurrent neural network

Where, fact(⋅)$${f_{act}}\left( \cdot \right)$$ - nonlinear activation function, usually Logistic or Tanh function.

WR∈ℝDR×NR$${W_R} \in {\mathbb{R}^{{D_R} \times {N_R}}}$$ - State-input weight matrix of the recurrent neural network, learnable parameters.

UR∈ℝDR×DR$${U_R} \in {\mathbb{R}^{{D_R} \times {D_R}}}$$ - State-state weight matrix of the recurrent neural network, learnable parameters.

bR∈ℝDR$${b_R} \in {\mathbb{R}^{{D_R}}}$$ - Bias vector of the recurrent neural network, learnable parameters.

VR∈ℝDR$${V_R} \in {\mathbb{R}^{{D_R}}}$$ - Linear transformation weight vector of the recurrent neural network, learnable parameter.

A recurrent neural network (LSTM) based fault prediction model is trained using the small batch gradient descent algorithm and the Adam optimisation algorithm to minimise the loss function of the model. The loss function of the model L(⋅)$$\mathcal{L}\left( \cdot \right)$$ is the mean square error between the predicted and true values of the target time series, calculated as: (3)L(θ)=1K∑i=tt+K−1(yi−y^i)2$$\mathcal{L}\left( \theta \right) = \frac{1}{K}\sum\limits_{i = t}^{t + K - 1} {{{\left( {{y_i} - {{\hat y}_i}} \right)}^2}}$$

Where: θ - all learnable parameters in the predictive model.

K - batch size, i.e., the number of training samples used to calculate the gradient and update the model parameters in one iteration of training.

Due to the dynamic characteristics of the heating system operating conditions, it is difficult for the prediction results of the prediction model to match the real change curves exactly for some normal fluctuations of larger magnitude of the observed values in the time series. Therefore, the spike fluctuations in the error vector are first filtered out by error smoothing. Specifically, the exponentially weighted moving average (EWMA) algorithm is used to smooth the error, and the following formula can express the smoothing process: (4)eS,t={ett=1wcet+(1−wc)eS,t−1t>1$${e_{S,t}} = \left\{ {\begin{array}{*{20}{l}} {{e_t}}&{t = 1} \\ {{w_c}{e_t} + \left( {1 - {w_c}} \right){e_{S,t - 1}}}&{t > 1} \end{array}} \right.$$

Where, e_{S, t} - smoothing error at moment t.

w_e - weight decay coefficient, calculated as: (5)we=2/(WEWMA+1)$${w_e} = {2 /{\left( {{W_{EWMA}} + 1} \right)}}$$

Where, W_EWMA - the window size of the exponentially weighted moving average.

At moment t, the smoothed error values of the current moment and the past N_e moments are selected to form an error vector eS=(eS,t−Ne,eS,t−Nc+1,…,eS,t)T$${e_S} = {\left( {{e_{S,t - {N_e}}},{e_{S,t - {N_c} + 1}}, \ldots ,{e_{S,t}}} \right)^T}$$, and the inter-anomaly values are determined based on the error vector.

Given the scalar vector z=(2.5,3,3.5,…,10)T$$z = {\left( {2.5,3,3.5, \ldots ,10} \right)^T}$$, calculate the candidate anomaly threshold as: (6)ε=μ(es)+zσsd(eS)$$\varepsilon = \mu \left( {{e_s}} \right) + z{\sigma _{sd}}\left( {{e_S}} \right)$$

Where: ε - the set of all candidate anomaly thresholds.

μ(⋅)$$\mu \left( \cdot \right)$$ - A function to calculate the mean.

σsd(⋅)$${\sigma _{sd}}\left( \cdot \right)$$ - A function to calculate the standard deviation.

The optimal anomaly threshold is determined by maximising the objective function equation: (7)ε*=argmaxε∈ε[μ(eS)−μ(eN)]/μ(eS)+[σsd(eS)−σsd(eN)]/σsd(eS)|eA|+|eseq|2{}$$\begin{array}{l} {\varepsilon ^*} = \mathop {\arg \max }\limits_{\varepsilon \in \varepsilon } \\ \frac{{{{\left[ {\mu \left( {{e_S}} \right) - \mu \left( {{e_N}} \right)} \right]} / {\mu \left( {{e_S}} \right)}} + {{\left[ {{\sigma _{sd}}\left( {{e_S}} \right) - {\sigma _{sd}}\left( {{e_N}} \right)} \right]} / {{\sigma _{sd}}\left( {{e_S}} \right)}}}}{{\left| {{e_A}} \right| + {{\left| {{e_{seq}}} \right|}^2}}}\left\{ {} \right\} \\ \end{array} $$

Where, e_N - set of normal error values eN={eS∈eS|eS<ε}$${e_N} = \left\{ {{e_S} \in {e_S}\left| {{e_S}} \right. < \varepsilon } \right\}$$.

e_A - set of candidate abnormal error values eA={eS∈eS|eS≥ε}$${e_A} = \left\{ {{e_S} \in {e_S}\left| {{e_S}} \right. \ge \varepsilon } \right\}$$.

e_seq - a set of continuous error vectors eseq={eseq}$${e_{seq}} = \left\{ {{e_{seq}}} \right\}$$, and a continuous error vector e_seq is a vector of candidate anomaly error values e_A ∈ e_A at successive time steps. In particular, an isolated candidate anomaly error value is also considered as a continuous error vector.

|eA|$$\left| {{e_A}} \right|$$ - A total number of candidate anomaly error values.

|eseq|$$\left| {{e_{seq}}} \right|$$ - The total number of continuous error vectors.

If the smoothed error value eS,t≥ε*$${e_{S,t}} \ge {\varepsilon ^*}$$ at moment t, then the error value at that moment is a candidate anomaly error value and the corresponding sensor observation is a candidate anomaly value.

By anomaly pruning, candidate anomalies that are unreliable and less suspicious are excluded. This step aims to reduce the false alarm rate of system fault diagnosis.

The maximum values in the set of normal error values e_N and the maximum values of each successive error vector e_seq in the set of successive error vectors e_seq are selected and arranged in descending order to form a new error vector e_max. Obviously, the length of this vector is |eseq|+1$$\left| {{e_{seq}}} \right| + 1$$. Iterate over the error vector e_max and compute the rate of descent of each element in the vector with respect to the latter element Δe: (8)Δei=emax,i−emax,i+1emax,ii=1,2,…,|eseq|$$\Delta {e_i} = \frac{{{e_{\max ,i}} - {e_{\max ,i + 1}}}}{{{e_{\max ,i}}}}i = 1,2, \ldots ,\left| {{e_{seq}}} \right|$$

Given a minimum descent rate Δe^*, if the i rd element e_max,i ∈ e_max of vector e_max has a descent rate Δe_i ≥ Δe^*. Then all the continuous error vectors e_seq corresponding to the first i elements {emax ,j∈emax |j≤i}$$\left\{ {{e_{{\text{max }},j}} \in {{\mathbf{e}}_{{\text{max }}}}|j \le i} \right\}$$ in vector e_max are still diagnosed as anomalies.

After pruning, the smoothing error values in those continuous error vectors that are still diagnosed as anomalies are the final anomalies. If the smoothed error value e_S,t at moment t is finally detected as an anomaly, it means that the heating system is in a faulty condition at moment t, and the algorithm will immediately send a fault alert to the operation manager.

The adjustment of the system parameters is mainly to consider the speed of response, the accuracy value in steady state conditions, the size of the overshooting amount, and whether the system control has stability. The specific rules of adjustment need to be considered in conjunction with the system output response curve. PID three parameter adjustment needs to achieve the fuzzy PID control under different circumstances of the input to the control parameters of the requirements of the specific circumstances are as follows: 1)

When the absolute value of error e is large, in order to avoid the phenomenon of calculus saturation in the heating system control, in order to make the system have good tracking performance, the value of proportional parameter K_p is taken as a larger value, the value of differential parameter K_d is taken as a smaller value, and the value of integral parameter K_i is taken as 0.

Fig. 4 shows the relationship between K_p, K_i, and K_d of the temperature control of the heating system and the error (e) and the rate of change of the error following the change of the error (ec). The specific analyses are as follows: 1)

Fig. 4(A) shows the relationship between the scale parameter K_p and e and ec. It can be seen that K_p increases with e and the two are positively correlated. In the interval ec ∈ [1,-1], K_p achieves the minimum value when ec=0.

Figure 4.

The relationship between PID parameter setting and error and error change rate

Under the environment of housing heating system, the three control methods of PID, fuzzy PID, and fuzzy self-tuning PID (Smith-fuzzy PID) proposed in this paper are simulated and compared and analysed. When the housing heating system is at primary load, the simulation comparison results of the three temperature control methods are shown in Fig. 5. It can be seen that the maximum value of PID is 26.9°C, and the time to reach stability from the coordinates reaches stability (23°C) after 0.813 × 10⁴ s, with an oscillation of 16.96%. The fuzzy PID has a maximum value of 24.6°C and reaches a steady state in 0.572 × 10⁴ s with an oscillation of 6.96%. The Smith-fuzzy PID has a maximum temperature value of 23.6, a minimum time required to reach a steady state of 0.332 × 10⁴ s and a minimum oscillation of 2.61%. In the whole temperature control process, Smith-fuzzy PID is obviously better than PID and fuzzy PID. PID has higher overshooting and a longer time to reach a steady state. Fuzzy PID has the second-best effect, with less overshooting and a relatively smaller stabilisation time. While Smith-Fuzzy PID has a high adaptive ability, low overshooting, short time to reach stabilisation, less time consuming, better stability and robustness. Compared with the other two control methods, it is more indicative that the selection of the Smith-Fuzzy PID control method is more suitable to be applied in this system.

Figure 5.

Simulation results of three temperature control methods

In order to verify whether the Smith-Fuzzy PID control method is better than the other two control methods when the heating system is in fault conditions. Under the same conditions, the effect of factors such as opening and closing the door back and forth on the temperature is considered, and the sensor fault signal is added, and the simulation results are shown in Fig. 6. It can be seen that under the condition of adding the fault signal for 2 × 10⁴ s, from the perspective of oscillation amplitude, stabilisation time and the speed of achieving stability after the disturbance, by comparing and analysing the effects of the three control methods, it can be concluded that the Smith-fuzzy PID is able to restore stability in the case of a fault disturbance in a shorter period of time of only 0.240 × 10⁴ s. This shows that the optimized control strategy proposed in this paper is able to restore stability faster in the case of a system fault compared to the other two control methods. Two control methods, it is able to carry out coordinated control faster, which further indicates that this method is the preferred method for the control of the housing heating system.

Figure 6.

Comparison of simulation results of heating system with fault signal

In order to validate the automatic diagnosis of heating system faults (LSTM-DT) method based on recurrent neural network (LSTM) and dynamic threshold algorithm (DT) proposed in this paper, intelligent classification and diagnosis research is carried out for the user’s behaviour of window opening, pipe blockage, actuator failure, and measurement of coarse error conditions in the indoor heating system. The sample data obtained through simulation is large, and nine feature quantities of indoor temperature difference, return water temperature difference, and indoor/outdoor temperature difference between the user and four neighbouring users are considered for classification and diagnosis. The simulation model is used to collect the initial data, and the corresponding sample data is obtained after a brief calculation to create a sample data set that includes data from the four operating conditions.

Since the features of the data are not independent of each other, according to Occam’s razor principle, it is hoped to use as few dimensions as possible for calculation. At the same time, reducing the dimension of the computation also reduces computational efficiency and reduces the possibility of overfitting phenomena in the classification diagnostic model.

As it is desired to use a smaller amount of features for diagnostic analysis, minimising the dimensionality of the sample data but not overly eliminating the amount of information contained in the samples, correlation analysis between 9-dimensional data is performed by plotting a correlation matrix (i.e., heatmap). The correlation between two or two feature measures of the sample data can be analysed using a heat map, the metric of which is the Pearson correlation coefficient, which is the quotient of the covariance and standard deviation between two feature measures, as in Eq: (9)ρX,Y=cov(X,Y)σXσY=E(XY)−E(X)−E(Y)E(X)2−E(Y)2E(Y)2−E(X)2$${\rho _{X,Y}} = \frac{{\operatorname{cov} (X,Y)}}{{{\sigma _X}{\sigma _Y}}} = \frac{{E(XY) - E(X) - E(Y)}}{{\sqrt {E{{(X)}^2} - E{{(Y)}^2}} \sqrt {E{{(Y)}^2} - E{{(X)}^2}} }}$$

Where: ρ_X,Y - Pearson’s correlation coefficient.

cov(X, Y) - Covariance of variables X and Y.

σ_X, σ_Y - Standard deviation of variables X and Y.

E() - Mathematical expectation of a parameter.

Fig. 7 shows the heat map of the correlation analysis of the eigenvalues of the sample data, and the values in the cells indicate the correlation coefficients between the two eigenvalues. Among them, A1~A5 is the temperature difference between indoor and outdoor, with the upstairs household, with the downstairs household, with the left-neighbouring household and with the right-neighbouring household, respectively. B1~B4 is the return water temperature difference with the upstairs household, with the downstairs household, with the left-neighbouring household and with the right-neighbouring household, respectively. As can be seen from the figure, the correlation coefficients between the four-room temperature difference parameters of the heating system and the four return water temperature difference parameters are all between 0.82 and 0.95, which is an extremely strong correlation, so there is no need for all of them to be used as the data feature quantity. Combined with the actual after comprehensive consideration, this paper selects the user indoor and outdoor temperature difference, upstairs room temperature difference, and downstairs room temperature difference, and upstairs return water temperature difference, and downstairs return water temperature difference five parameters composed of characteristic parameters, composed of indoor heating system working condition diagnostic data.

Figure 7.

Sample data characteristic quantity correlation analysis heat map

The eigenparameters can be expressed as: (10)K=(tn−tw,tn(i,j)−tn(i+1,j),tn(i,j)−tn(i−1,j),th(i,j)−th(i+1,j),th(i,j)−tn(i−1,j))$$\begin{array}{l} K = ({t_n} - {t_w},{t_{n(i,j)}} - {t_{n(i + 1,j)}},{t_{n(i,j)}} - {t_{n(i - 1,j)}}, \\ {t_{h(i,j)}} - {t_{h(i + 1,j)}}, {t_{h(i,j)}} - {t_{n(i - 1,j)}}) \\ \end{array}$$

Where, t_n - user indoor temperature, ℃.

t_w - Outdoor temperature, ℃.

t_h -Heating system return water temperature, ℃.

i - User’s floor.

j - User number.

It is worth noting that the sample data of users located on the bottom and top floors are incomplete, so the data of neighbouring users on the same floor are used to replace the data of upstairs or downstairs users for constructing the sample dataset.

In addition, in the actual uploaded data, there are far more normal working condition data than faulty working condition and abnormal behavioral working condition data. The possibility of a large number of abnormal data is small, so this paper randomly selects some data from the various types of data in the sample data set, including the user’s indoor and outdoor temperature difference, the indoor temperature difference with the upstairs and downstairs users, and the return temperature difference, and the experimental data set used for the verification of intelligent diagnostic methods is shown in Table 1. Among them, it includes 1892 data for normal working conditions and 900 data for fault conditions.

This subsection is designed to validate the suitability of LSTM-DT as a first-class classifier for a two-tier fault identification and localisation scheme. The experiments in this subsection focus on the variation of the accuracy and robustness of the LSTM-DT when a component failure occurs in the heating system. In this section, the effects of different data sampling window lengths, as well as noise intensity, on system fault identification are considered.

The following types of experiments are performed regarding the identification of system faults: 1)

Classification performance of the LSTM-DT model for identifying system faults, including the identification performance of LSTM-DT for different window lengths and different signal-to-noise ratios.

Figure 8 shows the recognition effect of LSTM-DT for nine types of system faults under different data window lengths and different signal-to-noise ratios. The recognition effectiveness is measured by the average F1 value, which is the average of the F1 values of the 900 heating system fault diagnosis, and is a comprehensive reflection of the recognition accuracy of the classifier. The data window lengths are 5, 10, 20, and 40 sampling points in order, and the sampling frequency is 1 sample/second. The signal-to-noise ratios are 20dB, 40dB, 60dB, 80dB, and 100dB in order, where 100dB means no noise is added. In the case of the same sample signal-to-noise ratio, it can be seen from the figure that the average F1 value is increased when the window length is changed from 5 to 40, indicating that the longer the data window length, the better the LSTM-DT recognition effect.

Figure 8.

LSTM-DT recognition effect with different window length and signal-to-noise ratio

With a constant data window length, it can be seen that the LSTM-DT remains essentially unchanged and stable for system fault identification F1 values when the signal-to-noise ratio is varied from 40dB to 100dB. When stronger noise (20dB) is added, the average F1 decreases, but both remain above 90%. It can be learnt that LSTM-DT is less affected by noise and can identify various sub-faults stably. From the overall point of view, the average F1 values of LSTM-DT are all above 90%, so it can be concluded that LSTM-DT has high accuracy and robustness for the identification of heating system faults. 2)

A comparative study of LSTM-DT with KNN, RF, and BPNN for identifying system faults, comparing the identification performance of the four methods at different window lengths and different signal-to-noise ratios.

Fig. 9 shows the recognition results of system faults under different data window lengths and different noise intensities. Where (A)~(D) are the identification results of LSTM-DT with KNN, RF and BPNN for system faults, respectively, and the measure in the matrix is the average F1 value, which is the average of the F1 values of 900 system faults. It can be seen that the flat F1 value of LSTM-DT is greater than that of KNN, RF and BPNN in the vast majority of cases. Moreover, LSTM-DT also performs more reliably when affected by noise compared to other methods. KNN also performs relatively well for the identification of system faults, while RF and BPNN perform relatively poorly.

Figure 9.

Comparison diagram of four methods to identify system faults

Overall, from the four subgraphs, the window length affects the recognition accuracy of various classification methods, so it is important to adopt the appropriate window length. In terms of method accuracy, LSTM-DT and KNN perform better, while RF and BPNN perform average. Taken together, the LSTM-DT model proposed in this paper has good accuracy and robustness for fault diagnosis in heating systems.