<figure><img src="image_1.png" /></figure>

HAL
open science

# Heat sources failure instant estimation: comparison of Kalman smoothing and supervised learning

Sylvain Verron, Mohamed Salim Bidou, Laetitia Perez, Laurent Autrique

► To cite this version:

Sylvain Verron, Mohamed Salim Bidou, Laetitia Perez, Laurent Autrique. Heat sources failure instant estimation: comparison of Kalman smoothing and supervised learning. MQDS 2023 - Congress on Measurement, Quality and Data Science, Jun 2023, Bordeaux, France. hal-04320806

HAL Id: hal-04320806

https://univ-angers.hal.science/hal-04320806v1

Submitted on 4 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Heat sources failure instant estimation: comparison of Kalman smoothing and supervised learning

Sylvain Verron¹, Mohamed Salim Bidou¹, Laetitia Perez¹, Laurent Autrique¹

¹LARIS, Polytech, University of Angers, 62, av. Notre Dame du Lac, 49000, Angers, France, sylvain.verron@univ-angers.fr

**ABSTRACT:** Fault detection has become a critical aspect in industrial systems due to the increased complexity and safety concerns. This paper focuses on the potential failure of a fixed heat source located on a thin aluminum plate, with the added difficulty of another non-failing heat source being present on the same plate. Only three fixed sensors are integrated to the plate, which is considered a 2D problem. The objective of this study is to compare two approaches for estimating the failure instant of the heat source. The first approach uses a Kalman filter, followed by a smoothing phase and an exploitation phase to better estimate the failure instant. The second approach is a supervised learning method using multiple regression, where the relationship between the failure instant and the observed sensor data is established through a training dataset generated using the model developed for the Kalman filter. The proposed application considers a 1m x 1m aluminum plate with 3mm thickness, 3 sensors and 2 heat sources, with an observation time of 3600 seconds. The effect of measurement noise on the identification of defects is analyzed for different signal-to-noise ratios through Monte-Carlo simulations.

**Keywords:** Fault Detection, Kalman filter, multiple regression, partial differential equations.

Received January 1, 2023; In final form March 31, 2023; Published June 2023

**Copyright:** This is an open-access article distributed under the terms of the Creative Comm seons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Corresponding author:** Sylvain Verron, e-mail: sylvain.verron@univ-angers.fr

## 1. INTRODUCTION

Fault detection has become a crucial area of research over the past two decades due to the increasing complexity of industrial systems and the growing concern for their safety. The detection of faults in complex physical systems is typically achieved by knowledge models based on partial differential equations (PDEs) [1]. In the field of thermic, for instance, heat exchange is described by Fourier's law, which results in a system of PDEs. Conventional thermal studies usually consider fixed or mobile heat sources in a finite dimensional space [2].

In this paper, we focus on the potential failure of a fixed heat source located on a thin aluminum plate, which can be considered a 2D problem. A challenging aspect of this study is the presence of another heat source, considered to be non-failing, on the plate. Additionally, we make the assumption that only 3 fixed sensors (with defined and fixed positions) are integrated into the plate [3].

The objective of this paper is to compare two methods for estimating the instant of failure of the heat source. The first approach we develop is based on a Kalman filter, followed by a smoothing phase and an exploitation phase of the latter to better estimate the failure instant a posteriori [4]. To model

the problem as a Kalman filter (linear dynamical system), a modelling step is necessary. For this, we will use the finite difference method [5].

The second proposed method is a supervised learning method, a multiple regression, exploiting the heat evolution model developed for the Kalman filter. For this, a training dataset is generated using the defined model, in which we simulate different possible failure times. The regression will then establish the link between the failure instants and the observed sensor data [6].

The proposed application is an aluminium plate 1 meter by 1 meter, 3mm thick, with 3 sensors and 2 supposedly known heat sources. The observation time is 3600 seconds, after which the proposed methods must announce the estimated instant of failure of the source of interest. An interesting point is that the frequency of sensors reading is low: every 9s. Particular attention is paid to the effect of measurement noise on defect identification to highlight the relevance of the two approaches. Several signal-to-noise ratios are considered through an analysis based on Monte-Carlo simulations.

Figure 1 show the position of the 2 sources and the positions of the three sensors. On figure 2, one can view the heat flux of each source, for normal operating conditions (no failure).
<figure><img src="image_3.png" /><figcaption>Figure 1. Sources and Sensors Positions</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. Heat flux of sources 1 and 2</figcaption></figure>

## 2. METHODS

This methodology section will be subdivided into several subsections. We will first approach the method of modelling the phenomenon of heat diffusion by the finite difference approach. Indeed, this modelling in the form of a linear system will be exploited by the two types of methods proposed to estimate the instant of failure of the heat source. The first proposed method will exploit the power of the Kalman filter and smoother, followed by an optimization phase. The third subsection will present a classic machine learning approach, namely multiple regression.

### 2.1. Finite difference modelling

The finite difference method is a numerical approach to solving differential equations. The principle is to cut the plate into a finite number of

small segments, or elements, and to use mathematical approximations to represent the conditions at each edge of these elements.

In our application, temperature evolution is described [2] by:

$$ \rho C \frac{\partial T}{\partial t} - \lambda \frac{\partial^2 T}{\partial x^2} - \lambda \frac{\partial^2 T}{\partial y^2} = \frac{G}{e} - \frac{2h(T-T_0)}{e} \quad (1) $$

Where T is temperature ($T_0$ is the initial one), t is time, $\rho C$ is the heat density, $\lambda$ is the thermal conductivity (considered as constant in this study), e is the thickness of the plate, h is the convective transfer coefficient, and finally, G is the heating flux.

To model a metal plate with the finite difference method, we discretize the plate to adapt it to a numerical model. For that, we define temperature in each discretised point {x=i, y=j}, at instant k by $T = T_{i,j}^k$. So, we can apply the finite difference method in order to approximate the first and second derivatives:

$$ \frac{\partial T_{i,j}^k}{\partial t} \approx \frac{T_{i,j}^{k+1} - T_{i,j}^k}{\Delta t} \quad (2) $$

$$ \frac{\partial^2 T}{\partial x^2} \approx \frac{T_{i-1,j}^k - 2T_{i,j}^k + T_{i+1,j}^k}{\Delta x^2} \quad (3) $$

$$ \frac{\partial^2 T}{\partial y^2} \approx \frac{T_{i,j-1}^k - 2T_{i,j}^k + T_{i,j+1}^k}{\Delta y^2} \quad (4) $$

Where $\Delta t$ is the time step, and $\Delta x$ and $\Delta y$ are the spatial steps on respectively the x dimension and the y dimension. Same discretization is done for flux heat $G = G_{i,j}^k$.

Using these formulas, we can construct a system of linear equations representing the temperature conditions at each edge of each element. This system can be solved to determine the temperature at each location and time, giving us a numerical approximation of the solution to the original differential equation.

$$ T_{i,j}^{k+1} = \left(1 - \frac{2\Delta t \cdot \lambda}{\rho C \cdot \Delta x^2} - \frac{2\Delta t \cdot \lambda}{\rho C \cdot \Delta y^2} - \frac{2\Delta t \cdot h}{\rho C \cdot e}\right) T_{i,j}^k + \left(\frac{\Delta t \cdot \lambda}{\rho C \cdot \Delta x^2}\right) T_{i-1,j}^k + \left(\frac{\Delta t \cdot \lambda}{\rho C \cdot \Delta x^2}\right) T_{i+1,j}^k + \left(\frac{\Delta t \cdot \lambda}{\rho C \cdot \Delta y^2}\right) T_{i,j-1}^k + \left(\frac{\Delta t \cdot \lambda}{\rho C \cdot \Delta y^2}\right) T_{i,j+1}^k + \left(\frac{\Delta t}{\rho C \cdot e}\right) G_{i,j}^k + \left(\frac{\Delta t \cdot 2h \cdot \theta_0}{\rho C \cdot e}\right) \quad (5) $$

Note that for a point in the domain, its temperature at time k+1 will depend on:

* Its own temperature at time k

* The temperature at time k of its neighbours (right/left/top/bottom)

* The heat flux at time k at this point

* The constant linked to h

For computational purposes, we rewrite the problem in the form of a state evolution model and
sequentially rearrange all nodes (i, j) with indices m = 1,...,M. where M =I × J. Thereby, our application could be modelled as a state representation of a discrete time invariant linear system:

$$
\left\{
\begin{array}{l}
\mathbf{T}^{k+1} = \mathbf{A} \cdot \mathbf{T}^k + \mathbf{B} \cdot \mathbf{G}^k + \mathbf{H} \\
\mathbf{T}_{\text{obs}}^k = \mathbf{C} \cdot \mathbf{T}^k
\end{array}
\right.
\qquad (6)
$$

Where $\mathbf{T}^k$ is a matrix containing temperature of all discretised points of the plate at the instant $k$, $\mathbf{G}^k$ is a matrix containing heat flux values (2 sources g1 and g2) at the instant $k$, $\mathbf{A}$ is a matrix allowing to take into account the temperature (with specific coefficient) of neighbours points to compute temperature at next instant step, $\mathbf{B}$ and $\mathbf{C}$ matrices are coding respectively, the positon of the sources on the plate and the position of the sensors on the plate. $\mathbf{H}$ is a constant vector of same value $(\frac{\Delta t \cdot 2h \cdot \theta_0}{\rho C \cdot e})$. $\mathbf{T}_{\text{obs}}^k$ is a matrix containing the value of the three sensors at instant k.

This modeling, also classically called “direct problem”, unlike the inverse problem, will make it
possible to simulate the system, to observe its behaviour, for different configurations (in particular different cases of failures). Some modeling parameters can cause the results to differ somewhat, in particular the step of discretization of the physical plate Δx and Δy, or the step of calculation time Δt.

For example, if it is assumed that a failure of
source 1 occurred at t = 1500s, the evolution of the
temperatures at the three sensors are shown in figure
3 (observations from sensors are subject to a
Gaussian noise of standard deviation of 0.5).

<figure><img src="image_5.png" /><figcaption>Figure 3. Temperatures without failure and
without noise (continuous) and with failure and
noise (points)</figcaption></figure>

2.2. First approach: Kalman Smoothing

The Kalman filter is a mathematical algorithm
used to estimate the state of a dynamic system from
noisy measurements [7]. It is a recursive algorithm
that uses a prediction step to estimate the state at the
next time step based on the current state and a
control input, followed by a correction step that
updates the state estimate based on the
measurement.

As we have modelled our application as a dynamical system (equation (6)), we can use the Kalman filter. But we can see that in this case, we need to know $\mathbf{G}^k$ for each instant. However, real values of $\mathbf{G}^k$ are not known, this is particularly what we would like to estimate: g1 and g2. So, in this context, we reformulate the dynamical system in order to integrate the $\mathbf{G}^k$ matrix in the state vector [8]. For that, we reformulate our equation (6), as presented in equation (7) below.

$$
\left\{
\begin{array}{l}
\begin{bmatrix} \mathbf{T}^{k+1} \\ \mathbf{G}^{k+1} \end{bmatrix} = \mathbf{A}' \cdot \begin{bmatrix} \mathbf{T}^k \\ \mathbf{G}^k \end{bmatrix} + \mathbf{H} + \mathbf{w}_k \\[1em]
\mathbf{T}_{\text{obs}}^k = \mathbf{C}' \cdot \begin{bmatrix} \mathbf{T}^k \\ \mathbf{G}^k \end{bmatrix} + \mathbf{v}_k
\end{array}
\right.
\qquad (7)
$$

Where: $\mathbf{A}' = \begin{bmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{0} & \mathbf{1} \end{bmatrix}$ and $\mathbf{C}' = [\mathbf{C} \quad \mathbf{0}]$, $\mathbf{w}_k$ and $\mathbf{v}_k$ are respectively the process noise and the observation noise, assumed to be Gaussian from mean zero and covariance matrices $\mathbf{Q}$ and $\mathbf{R}$.

We name $\hat{\mu}_k$ and $\hat{\Sigma}_k$ the estimate of, respectively, the mean and covariance matrix of the state vector $[\mathbf{T}^k | \mathbf{G}^k]$. With the transformed linear system of equation (7), it becomes possible to apply the Kalman filter and estimate $\hat{\mu}_k$, for each instant k, and so the temperature, but moreover the flux of the sources. The Kalman filter steps are described below:

$$
\text{Initialize the estimates: } \hat{\mu}_0 = \begin{bmatrix} \mathbf{T}^0 \\ \mathbf{G}^0 \end{bmatrix} \text{ and } \hat{\Sigma}_0 = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \end{bmatrix}
$$

For each k from 1 to 3600:

$$
\begin{align*}
\hat{\mu}_{k|k-1} &= A' \hat{\mu}_{k-1|k-1} + \mathbf{H} \\
\hat{\Sigma}_{k|k-1} &= A' \hat{\Sigma}_{k-1|k-1} A'^T + \mathbf{Q} \\
\mathbf{K}_k &= \hat{\Sigma}_{k|k-1} C'^T (C' \hat{\Sigma}_{k|k-1} C'^T + \mathbf{R})^{-1} \\
\hat{\mu}_{k|k} &= \hat{\mu}_{k|k-1} + \mathbf{K}_k (\mathbf{T}_{\text{obs}}^k - C' \hat{\mu}_{k|k-1}) \\
\hat{\Sigma}_{k|k} &= (I - \mathbf{K}_k C') \hat{\Sigma}_{k|k-1}
\end{align*}
$$

where *I* is the identity matrix and *T* denotes the
transpose.

Under the assumption that the entire time horizon
observations are available, it is possible to refine
estimation of $\hat{\mu}_k$ and $\hat{\Sigma}_k$ with the application of the Kalman smoother (also named RTS smoother) [9].

As opposed to just the current measurement used in the Kalman filter, the Kalman smoother is an algorithm that provides a more accurate estimate by using the entire sequence of measurements. The Kalman smoother can be seen as a backward pass of the Kalman filter, using the same prediction and correction steps, but using the corrected estimates from the future time steps to improve the estimates at the current time step. The Kalman smoother algorithm can be formulated as follows:

From the Kalman filter algorithm, start from the last estimates: $\hat{\mu}_{3600}$ and $\hat{\Sigma}_{3600}$, and for each k from 3599 to 1:

$$ 
\begin{aligned}
P_k &= \hat{\Sigma}_{k|k} A'^T \hat{\Sigma}_{k+1|k}^{-1} \\
\hat{\mu}_{k|1:3600} &= \hat{\mu}_{k|k} + P_k (\hat{\mu}_{k+1|1:3600} - \hat{\mu}_{k+1|k}) \\
\hat{\Sigma}_{k|1:3600} &= \hat{\Sigma}_{k|k} + P_k (\hat{\Sigma}_{k+1|1:3600} - \hat{\Sigma}_{k+1|k0}) C'^T
\end{aligned}
$$

We will name $\hat{g}$ the estimated signal composed of all $k \hat{\mu}_{k|1:3600}$. We have to identified the failure instant based on this signal, and knowing also the normal operating conditions, said theoretical g1. So, an optimization phase of the estimation of the failure instant $t_{fail}$, based on the minimization of a root mean square criteria between the $\hat{g}$ estimate signal and possible faulty signals (one can see [4] for more details).

For illustration of this methodology, we simulate an example of a failure of source 1 at $t = 1500s$. Available data is composed of the no fault g1 signal (theoretical g1), and of the data from the 3 sensors during the 3600 seconds, sampled every 9 seconds, giving 400 observed instants (with a noise sensor of $\sigma = 0.5$).

Firstly, on figure 4, we can plot the data from the 3 sensors. Secondly, from this data, we can apply the Kalman filter algorithm, followed by the Kalman smoother algorithm. Each of these algorithms, give an estimate of g1, but evidently, estimation from the smoother are the better one because this algorithm takes into account all possible measurements. The two estimations of g1 are represented on figure 5. Finally, figure 6 presents the global results after the search phase of the estimation of the failure instant $t_{fail}$. On this figure, we can observe that the failure instant estimate is $t_{fail} = 1491s$, which is very satisfying compared to the real failure instant of 1500s.

<figure><img src="image_6.png" /><figcaption>Figure 4. Available data</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. g1 estimates (filter and smoother)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. g1 best estimates</figcaption></figure>

## 2.3. Second approach: Multiple regression

Multiple regression is a machine learning algorithm used to predict one dependent variable (i.e. output or target variable) from multiple
independent variables (i.e. inputs or variables explanatory) [6].

The main objective of multiple regression is to find the best line that represents the relationships between the independent variables and the dependent variable. This line is called the regression function. The multiple regression model can be represented mathematically as follows:

$$y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots + \beta_nx_n + \varepsilon$$

Where y is the dependent variable (output), x₁, x₂, ..., xₙ are the independent variables (inputs), β₀, β₁, β₂, ..., βₙ are the regression coefficients, ε is the random error that takes into account unexplained influences on the dependent variable.

The objective of a regression learning will be to find the coefficients β₀, β₁, β₂, ..., βₙ which minimize the error between the predicted value and the actual value of the dependent variable. Optimization algorithms such as gradient descent or least squares can be used to find the optimal coefficients.

We will now explain how we use multiple regression to solve the problem of estimating the time of failure of a source. Indeed, in our problem, the instant of failure will be the variable to be predicted (y = tₕail). The descriptive variables will be the values of the 3 sensors registered over the period studied (1 hour).

In order to create a useful regression model, it is not possible to use a history since that would never cover the different possible failure times. Thus, we will use the modeling of the phenomenon (equation (6)) to generate/simulate different failure cases, and this for different sensor noise levels. The generation of this simulated history proceeds according to the following indications.

A possible instant of failure tₕail is randomly drawn over the period studied (3600s). The system is then simulated including the failure of source 1 at time tₕail. The values of the sensors C1, C2 and C3 are then recorded over the entire period. So, we obtain a number of $\frac{\text{total time} \times \text{number of sensors}}{\text{sampling time}}$ ie 3600×3/9=10800/9=1200 descriptors. Repeat this procedure for a desired number of observations, at least greater than number of descriptors.

As with all supervised learning techniques, an important step in the application of this regression will be the selection of components (descriptors), thus making it possible to make the model more robust and provide better performance. Indeed, it is

certain that among the 1200 descriptors used to estimate the instant of failure, some are much less informative than others, but still bring noise. Many component selection methods exist in the literature, generally classified into two types: wrapper and filter. Since the subject of this article is not the evaluation of this type of technique, we have chosen a classic evaluation of the importance of descriptors: univariate feature ranking for regression using F-tests. We then selected only the best descriptors.

## 3. RESULTS

In order to compare the two presented approaches (Kalman smoother and Linear regression), we propose different configurations. For each configuration, all simulations were performed on Matlab software.

The modeling parameters making it possible to solve the direct problem (finite differences) were chosen as follows: time step Δt = 1 s, and Δx = Δy = 0.1 m.

The first step is to use finite difference modeling to generate the history. The number of individuals in this history was chosen to be 5000. Once the history was created, we applied the Univariate feature ranking for regression using F-tests algorithm, allowing us to see the importance of each descriptor with respect to the instant of failure variable. The result of this evaluation can be seen in Figure 7. We thus clearly see that certain descriptors are not at all informative. We therefore made the choice to select only those that provided a lot of information.

<figure><img src="image_9.png" /><figcaption>Figure 7. univariate feature ranking using F-tests
For 30 simulations, we have run, on the same data, the two presented methods. We give, in table</figcaption></figure>
1, the results for 3 different level of sensor noise: $\sigma = 0.1$, $\sigma = 0.5$ and $\sigma = 1$. So, in order to see the impact of such a parameter, for each configuration, we give results in term of mean and standard deviation (in brackets) on 30 simulations.

<table><caption>Table 1. Monte-Carlo simulations results</caption><thead><tr><th></th><th>Kalman Smoother<br>$t_{fail}$</th><th>Linear Regression<br>$t_{fail}$</th></tr></thead><tbody><tr><td>$\sigma = 0.1$</td><td>1495,1(1,69)</td><td>1506,2 (6,53)</td></tr><tr><td>$\sigma = 0.5$</td><td>1494,9(11,09)</td><td>1508,0 (40,3)</td></tr><tr><td>$\sigma = 1$</td><td>1500,3 (14,61)</td><td>1514,9 (66,65)</td></tr></tbody></table>

As we can see in table 1, the Kalman smoother method is much more efficient than linear regression. Indeed, even if for most cases the averages of the simulations are comparable, we clearly see that the standard deviation for the linear regression method is significantly higher than for the Kalman smoother.

Another fact, absolutely prohibitive for the use of linear regression on this type of problem, is the fact that this technique cannot correctly diagnose a fault on a heat source other than source one. Indeed, if we wanted to take this possible technique to diagnose the failure of source 1 and source 2, it would amount to creating a virtual history with not 3600 possible failure instants (those of source one), but the combination of 3600 × 3600 possible failure times, then leading to a combinatorial explosion (almost 13 million cases). Moreover, the Kalman smoother technique will be able to diagnose failures but also restarts of sources 1 and 2.

## 4. CONCLUSIONS

In this article, we proposed to compare 2 approaches to solve the problem of identifying an instant of failure of a heat source on a fixed plate equipped with 3 sensors. These 2 approaches, the smoother of Kalman as well as the linear regression, required the creation of a model based on the finite differences. In the case of linear regression, this model is used to generate a simulated example database, this base then being used to correctly estimate the coefficients of a regression model.

In the case of the Kalman smoother, the model is seen as a dynamic linear system, thus allowing an estimation of the flux of each of the heat sources. In order to compare the 2 approaches, 30 simulations of the system with a failure at time 1500 seconds were studied. Although on average the 2 techniques are very close to a diagnosis at 1500 seconds, the Kalman smoother is much more efficient because the standard deviation over the 30 simulations is much lower than for linear regression.

## REFERENCES

[1] R. Isermann, Fault-Diagnosis Systems From Fault Detection to Fault Tolerance. Springer. (2006)

[2] H. Massard, H. Orlande, and O. Fudym, Estimation of position dependent transient heat source with the kalman filter, In Inverse Problems in Science and Engineering, vol. 20, no. 7, pp. 1079–1099 (2012).

[3] M.S. Bidou, L. Perez, S. Verron and L. Autrique, Identification of failure times for a system governed by a non-linear parabolic partial differential equation. In 1st IFAC Workshop on Control of Complex System (2022).

[4] M.S. Bidou, S. Verron, L. Perez and L. Autrique, Kalman smoother for detection of heat sources defects. In conference ICCAD (2022).

[5] C. Grossmann; H. Roos; M Stynes, Numerical Treatment of Partial Differential Equations. Springer-Verlag Berlin. (2007).

[6] R.O. Duda, P.E. Hart, D.G. Stork, Pattern Classification, 2nd Edition. (2000).

[7] R.E. Kalman, A new approach to linear filtering and prediction problems. (1960).

[8] J.C. Berg, A.K. Miller, Force Estimation via Kalman Filtering for Wind Turbine Blade Control. In: Structural Dynamics and Renewable Energy, Volume 1. Conference Proceedings of the Society for Experimental Mechanics Series. Springer. (2011)

[9] K. P. Murphy, Machine learning: a probabilistic perspective. MIT press. (2012).