Kalman Filter Estimation and Its Implementation

Abstract

In this chapter, we use the Kalman filter to estimate the future state of a system. We present the theory, design, simulation, and implementation of the Kalman filter. We use as a case example the estimation of temperature using a Resistance Temperature Detector (RTD), which has not been reported before. After a brief literature review, the theoretical analysis of a Kalman filter is presented along with that of the RTD. The dynamics of the RTD system are analytically derived and identified using Matlab. Then, the design of a time-varying Kalman filter using Matlab is presented. The solution to the Riccati equation is used to estimate the future state. Then, we implement the design using C-code for a microprocessor ATMega328. We show under what conditions the system may be simplified. In our case, we reduced the order of the system to that of a system having a 1st order response, that of an RC system, giving us satisfactory results. Furthermore, we can find two first order systems whose response defines two boundaries inside which the evolution of a second order system remains.

1. Introduction

A deterministic system is a system whose governing physical laws are specified so that if the state of the system at some time is known, then one can precisely predict the state at a later time. Nondeterministic systems are divided into two categories: stochastic and random. A stochastic system has governing physical laws that even if the state at some point in time is known precisely, it is impossible to determine the state of the system at a later time precisely. It is possible to determine the probability of a state, rather than the state itself. A random system is one which has no apparent governing physical laws. Practically, we treat all unpredictable systems, stochastic or random as stochastic systems, since we employ the same methods to study them. While we are unable to predict the state of a random process, we can evolve a strategy to deal with such processes. Such a strategy is based on a branch of mathematics dealing with unpredictable systems, called statistics.

Estimation is the process of extracting information from data which can be used to predict the behavior of state variables in a system. The estimation uses statistical criteria to infer the actual value of unknown variables. Estimation models are used to process noisy measurements, filter them, and detect inaccuracies. When random signals are passed through a deterministic system, their statistical properties are modified. A deterministic system to which random signals are input, so that the output is a random signal with desired statistical properties is called a filter. Filters can be linear or nonlinear, time-invariant or time varying. However, for simplicity we will usually consider linear, time-invariant filters. Linear, time-invariant filters are commonly employed in control systems to reduce the effect of measurement noise on the control system. In such systems, the output is usually a superposition of a deterministic signal and a random measurement noise.

The output of a filter not only has a frequency content different from the input signal, but also certain other characteristics of the filter, such as a phase-shift or a change in magnitude. In other words, the signal passing through a filter is also distorted by the filter, which is undesirable. A filter would produce an output signal based upon its characteristics, described by the transfer-function, frequency or impulse response, or a state-space representation of the filter. However, a filter can be designed to achieve a desired set of performance objectives, i.e. the numerator and denominator polynomials of the filter’s transfer function, or coefficient matrices of the filter’s state-space model, can be selected by a design process to achieve the conflicting requirements of maximum noise attenuation and minimum signal distortion.

There are several prediction models to infer the system state, although, it can be shown that of all estimation tools Kalman Filter (KF) is the one that minimizes the variance of the estimation error which enables accurate estimation of the process.

1.1 Literature review

The first application of state estimation was in the aerospace field to solve problems related to the prediction of position in aerospace vehicles. Nowadays, estimation has been applied in several fields of engineering and control systems. One common application is in data acquisition, to solve the problem of predicting the state of a system that cannot be measured directly due to the characteristics and complexity of the environment.

KF is an estimator proposed by Rudolph E. Kalman in 1960. It is an algorithm to estimate the evolution of a dynamic system, especially when data has a lot of noise. The principle of the filter is to find the probability of the hypothesis of predicted state and using the data from the measurement to correct it and improve the future estimation at each time. It is a suitable algorithm to apply in dynamic systems, linking real-time measurements and predicting the state of system parameters through time approaches. KF has been implemented in several fields, such as in navigation systems [1, 2, 3, 4], financial models [5, 6, 7], tracking vehicles [8, 9] and image processing [10, 11, 12]; only to mention some of them. Nevertheless, this statistical tool is useful for two main purposes: estimation and performance analysis of estimators.

In the field of IC technology, it has been implemented for thermal estimation. Multicore processors use a dynamic thermal management mechanism that use embedded thermal sensors for monitoring the real-time thermal behavior of the processor, this kind of sensors are susceptible to a variety of source of noise and this causes the discrepancies between actual temperatures observed by on-chip thermal sensors. Therefore, to fix the discrepancies in sensing, Kalman’s prediction is used to estimate real values from noisy sensor readings [13]. Another novel application of KF is in the electric vehicle industry, the estimation of the charge state of lithium-ion battery is an important parameter in order to guarantee a safe operation of them. The battery performance is influenced by aging; this fact makes difficult to predict the battery state, to overcome this issue the application of KF in combination with other methods is a suitable methodology [14, 15, 16, 17].

Recently, KF has been applied in several industrial applications. With the development of manufacturing process, welding automation emerges as one important tool to speed up the production rate in the assembly line in stronger and high-quality welds. Nevertheless, there are several factors that could influence the welding quality and the most important is the arc length, which could be influenced by the irregular surface of the workpiece and the loss of the tungsten electrode. To enhance the quality during the Gas-Tungsten Arc Welding (GTAW) process, KF is applied in order to keep the arc length stable and minimize the external noise [18]. In the field of sensorless control, KF have been used in intelligence electrical drives. To control induction motor drives without mechanical speed sensor at the motor shaft allows reduced hardware complexity, and low costs. Additionally, the use of induction motors without position sensor is useful for applications with abrasive and hard surface. Thereby, the application of an estimation method it’s necessary in order to predict the position and velocity of the shaft [19, 20, 21].

In applications related with radio astronomy, KF has been applied for the analysis of Very-Long-Baseline Interferometry (VLBI) data, in order to analyze parameters such as base line lengths, earth orientation parameters, radio source coordinates and tropospheric delays. Nowadays, modern antennas are being constructed and equipped with highly accurate broadband receiving systems. Besides the accurate observations gotten by astronomic instruments, it is necessary to implement estimation methods in order to optimize the models applied in data analysis [22, 23]. In power systems, one of the main difficulties is power quality due to total harmonics distortion (THD) that is mainly caused by nonlinear loads. THD effects are strongly correlated with issues as device heating, break down electronic components, network interference, etc. Several filters have been performed to decrease the effect of harmonics; nevertheless, the application of KF has shown an important reduction in the effect of harmonics [24, 25, 26]. In the field of biomedicine KF is widely used over other estimation methodologies to overcome the different sources of noise. Specifically, KF has been used to smooth and predict signals from Electroencephalogram and Electrocardiogram signals [27, 28]. Recently in the literature there are reports on a new methodology to protect the confidentiality of the transmitted data based on a Kalman filter. This strategy proposes the implementation of encrypted algorithm using KF, and is suggested to be used in Industrial cyber–physical systems (ICPSs) to protected data privacy [29, 30].

As it has been mentioned above, KF has been used in diverse fields of science and technology to predict specific parameters of interest according to the application. Temperature evolution is an important parameter to measure and predict, in order to study or control the temperature in an environment [31, 32], device [13, 33, 34] and process [18]. It is well known that RTDs are commercial devices very useful to monitor the temperature due their stability and accuracy. However, RTDs are self-heating causing noisy readings making the RTD a suitable example to implement KF for temperature estimation. Importantly, we searched in the literature and found no evidence of previous work reporting the use of a KF to filter the noise and predict the temperature behavior from RTD readings.

2. Theoretical analysis of a Kalman Filter

The final objective of this study is to obtain the specification of a linear dynamic system (Wiener filter [35]) which accomplishes the prediction, separation, or detection of a random signal [36]. With the state-transition method, a single derivation covers a large variety of problems: growing and infinite memory filters, stationary and non-stationary statistics, etc. Having guessed the “state” of the estimation (i.e., filtering or prediction) problem correctly, one is led to a nonlinear difference (or differential) equation for the covariance matrix of the optimal estimation error. From the solution of the equation for the covariance matrix we obtain the coefficients characterizing the optimal linear filter [36]. The following is a simplified derivation described previously in the references [37, 38].

2.1 Defining statistical quantities of use

The initial state, x(0), of a stochastic system is insufficient to determine its future state, x(t). Thus, based upon a statistical analysis of similar systems, and taking the average of their future states at a given time, t, we can calculate the mean state-vector as follows:

$${\mathbf{x}}_{\mathrm{m}}\left(\mathrm{t}\right)=\left(1/\mathrm{N}\right)\sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right)$$E1

Thus ${\mathit{x}}_m\left(t\right)$ is the expected state vector after studying N systems. It is also called the expected value of the state-vector, ${\mathit{x}}_m\left(t\right)=E\left[\mathit{x}\left(t\right)\right]$. Another statistical quantity of use is the correlation matrix of the state-vector:

$${\mathbf{P}}_x\left(\mathrm{t},\mathrm{\tau }\right)=\left(\mathrm{l}/N\right)\sum _{i=1}^N{x}_i\left(t\right)x_i^T\left(\tau \right)$$E2

The correlation matrix, ${\mathbf{P}}_x\left(\mathrm{t},\mathrm{\tau }\right),$is a measure of correlation, a statistical property among the different state variables, and between the same state variable at two different times. Two scalar variables, $x_1\left(t\right)$ and $x_2(t$), are said to be uncorrelated if the expected value of $x_1\left(t\right)x_2(\tau$), i.e. $E\left[x_1\left(t\right)x_2\left(\tau \right)\right]=0$, where τ is different from t.

The correlation matrix is the expected value of the matrix ${\mathit{x}}_i\left(t\right){\mathit{x}}_i^T\left(\tau \right)$, or ${\mathbf{P}}_x\left(\mathrm{t},\mathrm{\tau }\right)=\mathrm{E}\left[x_i\left(t\right)x_i^T\left(\tau \right)\right]$. When t = τ, the correlation matrix ${\mathbf{P}}_x\left(\mathrm{t},\mathrm{t}\right)=\mathrm{E}\left[x_i\left(t\right)x_i^T\left(t\right)\right]$, is called the covariance matrix. The covariance matrix, ${\mathbf{P}}_x\left(\mathrm{t},\mathrm{t}\right)$, is symmetric. If ${\mathbf{P}}_x\left(\mathrm{t},\mathrm{\tau }\right)$ is a diagonal matrix i.e. $E\left[x_i\left(t\right)x_j\left(\tau \right)\right]=0$, where $\mathrm{i}\ne \mathrm{j}$, it implies that all the state variables are uncorrelated.

2.2 Defining the filter in state space - discrete domain

Consider a plant which we cannot model accurately using only a deterministic model, because of the presence of uncertainties called process noise and measurement noise:

$${\mathit{x}}_{k+1}=\mathit{A}{\mathit{x}}_{\mathit{k}}+{\mathit{w}}_{\mathit{k}}$$E3

$${\mathit{y}}_{\mathit{k}}=\mathit{C}{\mathit{x}}_{\mathit{k}}+{\mathit{v}}_{\mathit{k}}$$E4

In the linear, time-varying state-space representation above, w is the process noise vector which may arise due to modeling errors such as neglecting nonlinear dynamics, and v is the measurement noise vector. The random noises, w and v, are assumed to be stationary white noises. The covariance matrices of stationary white noises, w and v, can be expressed as follows:

$$\mathit{Q}=\mathrm{E}\left[w_k{w}_k^T\right]$$E5

$$\mathbf{R}=\mathrm{E}\left[v_k{v}_k^T\right]$$E6

Since we cannot predict the state-vector, x of a stochastic plant, an observer is required for estimating the state-vector, based upon a measurement of the output, y and a known input, u. We need an observer that calculates the estimated state-vector, $\widehat{\mathbf{x}}$, optimally, based upon statistical description of the vector output and plant state. Such an observer is the Kalman Filter, which minimizes a statistical measure of the estimation error, ${\mathbf{e}}_k={\mathbf{x}}_{\mathit{k}}-{\widehat{\mathbf{x}}}_k$. This statistical measure is the covariance of the estimation error:

$${\mathbf{P}}_{\mathit{k}}=\mathrm{E}\left[{\mathit{e}}_k{\mathit{e}}_k^T\right]=E\left[\left({\mathbf{x}}_{\mathit{k}}-{\widehat{\mathbf{x}}}_k\right){\left({\mathbf{x}}_{\mathit{k}}-{\widehat{\mathbf{x}}}_k\right)}^T\right]$$E7

Since the state-vector, x, is a random vector and the estimated state $\widehat{\mathbf{x}}$, is based on the measurement of the output, y, for a finite time, say T, where $\mathrm{t}\ge \mathrm{T}$ then a true statistical average of x would require measuring the output for an infinite time.

If T < t, this is a data-smoothing (interpolation) problem. If T = t, this is called filtering. If T > t, we have a prediction problem. Since the original treatment is general enough, the collective term estimation is used [36].

Hence, the best estimate to obtain for x is not the true mean, but a conditional mean, ${\mathit{x}}_{\mathit{m}}$, based on only a finite time record of the output, y:

$${\mathit{x}}_m=E\left[\mathit{x}:\mathit{y},T\le t\right]$$E8

Taking in consideration the deviation of the estimated state- vector, $\widehat{\mathbf{x}}$, from the conditional mean, ${\mathit{x}}_{\mathit{m}}$, we can write the estimated state- vector as:

$$\widehat{\mathbf{x}}={\mathit{x}}_{\mathit{m}}+∆\mathit{x}$$E9

$∆\mathit{x}$ is the deviation from the conditional mean. The conditional covariance matrix of the estimation error based on a finite record of the output is then:

$${\mathbf{P}}_{\mathit{k}}=\mathrm{E}\left[{\mathit{e}}_k{\mathit{e}}_k^T:\mathit{y},T\le t\right]=E\left[\mathit{x}{\mathit{x}}^T\right]-{\mathit{x}}_m{\mathit{x}}_m^T+∆\mathit{x}∆{\mathit{x}}^T$$E10

The best estimate of state-vector happens if $∆\mathit{x}=0,\text{or}\widehat{\mathbf{x}}={\mathit{x}}_{\mathit{m}}$, and would result in a minimization of the conditional covariance matrix, or error covariance matrix, ${\mathbf{P}}_{\mathit{k}}$. In other words, minimization of ${\mathbf{P}}_{\mathit{k}}$ yields the optimal observer, which is the Kalman filter.

2.3 Defining the Kalman gain

The state-equation of the Kalman filter is that of a time-varying observer, and can be written as follows:

$${\widehat{\mathbf{x}}}_{\mathit{k}+\mathbf{1}}=\mathit{A}{\widehat{\mathbf{x}}}_{\mathit{k}}+\mathit{B}{\mathit{u}}_{\mathit{k}}+{\mathit{K}}_{\mathit{k}}\left[{\mathit{y}}_{\mathit{k}}-\mathit{C}{\widehat{\mathbf{x}}}_{\mathit{k}}\right]$$E11

${\mathit{K}}_{\mathit{k}}$ is the gain matrix of the Kalman filter. Assuming the prior estimate of ${\widehat{\mathbf{x}}}_{\mathit{k}}$ is called ${\widehat{\mathbf{x}}}_k^{\prime }$, gained by knowledge of the system. We write an update equation for the new estimate, combing the old estimate with measurement data, ${\widehat{\mathbf{x}}}_k^{\prime }\mathit{=}\mathit{A}{\widehat{\mathbf{x}}}_{\mathit{k}}+\mathit{B}{\mathit{u}}_{\mathit{k}}$ and:

$${\widehat{\mathbf{x}}}_{\mathit{k}}={\widehat{\mathbf{x}}}_k^{\prime }+{\mathit{K}}_{\mathit{k}}\left[{\mathit{y}}_{\mathit{k}}-\mathit{C}{\widehat{\mathbf{x}}}_k^{\prime }\right]$$E12

If we substitute Eq. (4) into Eq. (12) we get:

$${\widehat{\mathbf{x}}}_{\mathit{k}}={\widehat{\mathbf{x}}}_k^{\prime }+{\mathit{K}}_{\mathit{k}}\left[\mathit{C}{\mathit{x}}_{\mathit{k}}+{\mathit{v}}_{\mathit{k}}-\mathit{C}{\widehat{\mathbf{x}}}_k^{\prime }\right]$$E13

Substituting Eq. (13) into Eq. (7)

$${\mathbf{P}}_{\mathit{k}}=E\left[\left(\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}\mathit{C}\right)\left({\mathbf{x}}_{\mathit{k}}-{\widehat{\mathbf{x}}}_k^{\prime }\right)-{\mathit{K}}_{\mathit{k}}{\mathit{v}}_{\mathit{k}}\right){\left(\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}\mathit{C}\right)\left({\mathbf{x}}_{\mathit{k}}-{\widehat{\mathbf{x}}}_k^{\prime }\right)-{\mathit{K}}_{\mathit{k}}{\mathit{v}}_{\mathit{k}}\right)}^T\right]$$E14

Here $\left({\mathbf{x}}_{\mathit{k}}-{\widehat{\mathbf{x}}}_k^{\prime }\right)$ is the error of the prior estimate. Since there is no correlation among the input, process noise and measurement noise, then the expectation may be re-written as;

$${\mathbf{P}}_{\mathit{k}}=\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}\mathit{C}\right)E\left[\left({\mathbf{x}}_{\mathit{k}}-{\widehat{\mathbf{x}}}_k^{\prime }\right){\left({\mathbf{x}}_{\mathit{k}}-{\widehat{\mathbf{x}}}_k^{\prime }\right)}^T\right]{\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}\mathit{C}\right)}^T+{\mathit{K}}_{\mathit{k}}\mathit{E}\left[{\mathit{v}}_{\mathit{k}}{{\mathit{v}}_{\mathit{k}}}^{\mathit{T}}\right]{\mathit{K}}_{\mathit{k}}^{\mathit{T}}$$E15

Using Eqs. (6) and (7), we obtain:

$${\mathbf{P}}_{\mathit{k}}=\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}\mathit{C}\right){\mathit{P}}_{\mathit{k}}^{\prime }{\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}\mathit{C}\right)}^T+{\mathit{K}}_{\mathit{k}}\mathit{R}{\mathit{K}}_{\mathit{k}}^{\mathit{T}}$$E16

Eq. (16) is the error covariance update equation, where ${\mathit{P}}_{\mathit{k}}^{\prime }$ is the prior estimate of ${\mathbf{P}}_{\mathit{k}}$.

The trace of the error covariance matrix is the sum of the mean squared errors. The mean squared error may be reduced by minimizing the trace of ${\mathbf{P}}_{\mathit{k}}$. This requires to differentiate the trace of ${\mathbf{P}}_{\mathit{k}}$ with respect to ${\mathit{K}}_{\mathit{k}}$, then the result set to zero to find ${\mathit{K}}_{\mathit{k}}$ that minimizes the trace of ${\mathbf{P}}_{\mathit{k}}$.

We rewrite Eq. (16);

$${\mathbf{P}}_{\mathit{k}}={\mathit{P}}_{\mathit{k}}^{\prime }-{\mathit{P}}_{\mathit{k}}^{\prime }{\mathit{C}}^{\mathit{T}}{\mathit{K}}_{\mathit{k}}^{\mathit{T}}-{\mathit{K}}_{\mathit{k}}\mathit{C}{\mathit{P}}_{\mathit{k}}^{\prime }+{\mathit{K}}_{\mathit{k}}\mathit{C}{\mathit{P}}_{\mathit{k}}^{\prime }{\mathit{C}}^{\mathit{T}}{\mathit{K}}_{\mathit{k}}^{\mathit{T}}+{\mathit{K}}_{\mathit{k}}\mathit{R}{\mathit{K}}_{\mathit{k}}^{\mathit{T}}$$E17

Taking the trace of this expression gives:

$$\mathbf{T}\left[{\mathbf{P}}_{\mathit{k}}\right]=\mathrm{T}\left[{\mathit{P}}_{\mathit{k}}^{\prime }\right]-2\mathrm{T}\left[{\mathit{K}}_{\mathit{k}}\mathit{C}{\mathit{P}}_{\mathit{k}}^{\prime }\right]+T\left({\mathit{K}}_{\mathit{k}}\left(\mathit{C}{\mathit{P}}_{\mathit{k}}^{\prime }{\mathit{C}}^{\mathit{T}}+\mathit{R}\right){\mathit{K}}_{\mathit{k}}^{\mathit{T}}\right)$$E18

Then, we differentiate with respect to ${\mathit{K}}_{\mathit{k}}$;

$$\frac{\mathrm{d}\mathbf{T}\left[{\mathbf{P}}_{\mathit{k}}\right]}{\mathrm{d}{\mathit{K}}_{\mathit{k}}}=-2\mathrm{T}\left[\mathit{C}{\mathit{P}}_{\mathit{k}}^{\prime }\right]+2T\left({\mathit{K}}_{\mathit{k}}\left(\mathit{C}{\mathit{P}}_{\mathit{k}}^{\prime }{\mathit{C}}^{\mathit{T}}+\mathit{R}\right)\right)$$E19

Setting to zero and solving for ${\mathit{K}}_{\mathit{k}}$ we obtain the Kalman gain equation:

$${\mathit{K}}_{\mathit{k}}={\mathit{P}}_{\mathit{k}}^{\prime }{\mathit{C}}^{\mathit{T}}{\left(\mathit{C}{\mathit{P}}_{\mathit{k}}^{\prime }{\mathit{C}}^{\mathit{T}}+\mathit{R}\right)}^{-\mathbf{1}}$$E20

Substitution of Eq. (20) into [17], gives:

$${\mathbf{P}}_{\mathit{k}}=\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}\mathit{C}\right){\mathit{P}}_{\mathit{k}}^{\prime }$$E21

Eq. (21) is the update equation for the error covariance matrix with optimal gain.

State projection is derived using;

$${\widehat{\mathbf{x}}}_{\mathrm{k}+1}^{\prime }=\mathbf{A}{\widehat{\mathbf{x}}}_{\mathrm{k}}^{\prime }+{\mathbf{w}}_{\mathbf{k}}$$E22

To project the error covariance matrix into the next time interval, k + 1 we first find an expression for the error based on the prior error;

$$\begin{array}{cc}{\mathbf{e}}_{\mathbf{k}+\mathbf{1}}^{\prime }& ={\mathbf{x}}_{\mathrm{k}+1}-{\widehat{\mathbf{x}}}_{\mathrm{k}+1}^{\prime }\\ & =\mathbf{A}{\mathbf{x}}_{\mathbf{k}}+{\mathbf{w}}_{\mathbf{k}}-\mathbf{A}{\widehat{\mathbf{x}}}_{\mathrm{k}}^{\prime }\\ & =\mathbf{A}{\mathbf{e}}_{\mathbf{k}}+{\mathbf{w}}_{\mathbf{k}}\end{array}$$E23

Eq. (7) in time k + 1 is;

$${\mathbf{P}}_{\mathit{k}+\mathbf{1}}^{\prime }=\mathrm{E}\left[{\mathit{e}}_{\mathit{k}+\mathbf{1}}^{\prime }{{\mathit{e}}_{\mathit{k}+\mathbf{1}}^{\prime }}^{\mathit{T}}\right]=E\left[\left(\mathit{A}{\mathit{e}}_{\mathit{k}}+{\mathit{w}}_{\mathit{k}}\right){\left(\mathit{A}{\mathit{e}}_{\mathit{k}}+{\mathit{w}}_{\mathit{k}}\right)}^T\right]$$E24

Assuming that ${\mathit{e}}_{\mathit{k}}$ and $w_k$ have zero cross-correlation.

$$\begin{array}{cc}{\mathbf{P}}_{\mathit{k}+\mathbf{1}}^{\prime }& =\mathrm{E}\left[{\mathit{e}}_{\mathit{k}+\mathbf{1}}^{\prime }{{\mathit{e}}_{\mathit{k}+\mathbf{1}}^{\prime }}^{\mathit{T}}\right]\\ & =E\left[\left(\mathit{A}{\mathit{e}}_{\mathit{k}}{{\mathit{e}}_{\mathit{k}}}^{\mathit{T}}{\mathit{A}}^{\mathit{T}}+{\mathit{w}}_{\mathit{k}}{{\mathit{w}}_{\mathit{k}}}^{\mathit{T}}\right)\right]\\ & =\mathbf{A}{\mathbf{P}}_{\mathbf{k}}{\mathbf{A}}^{\mathbf{T}}+\mathbf{Q}\end{array}$$E25

This completes the description of the filter.

2.4 Algorithm loop

An algorithm loop is required to make the program in MATLAB and in C-code for the microprocessor. The loop is summarized in the Figure 1.

Figure 1.

Recursive algorithm for the Kalman filter.

The KF assumes that the system model is linear and known, the system and measurement noises are white, and the states have initial conditions with known means and variances. The power spectral densities used can be treated as tuning parameters to design an observer with excellent performance and robustness. The linear Kalman filter can also be used to design observers for nonlinear plants, by treating nonlinearities as process noise with appropriate power spectral density matrix.

2.5 Derivation of the Riccati equation

Since the Kalman filter is an optimal observer the appearance of matrix Riccati equation is not surprising. We are interested in a steady Kalman filter, i.e. the Kalman filter for which the covariance matrix converges to a constant in the limit $t\to \infty$. This happens when the plant is time invariant. The derivation goes as follows [39, 40]:

From the projections into $\infty$ we get:

$${\widehat{\mathbf{x}}}_{\infty +1}^{\prime }=\mathbf{A}{\widehat{\mathbf{x}}}_{\infty }$$E26

$${\mathbf{P}}_{\infty +\mathbf{1}}=\mathbf{A}{\mathbf{P}}_{\infty }{\mathbf{A}}^{\mathbf{T}}+\mathbf{Q}$$E27

$${\mathbf{P}}_{\infty }=\left(\mathit{I}-{\mathit{K}}_{\infty }\mathit{C}\right){\mathit{P}}_{\infty }^{\prime }$$E28

$${\widehat{\mathbf{x}}}_{\mathbf{\infty }}={\widehat{\mathbf{x}}}_{\infty }^{\prime }+{\mathit{K}}_{\mathbf{\infty }}\left[{\mathit{y}}_{\mathbf{\infty }}-\mathit{C}{\widehat{\mathbf{x}}}_{\infty }^{\prime }\right]$$E29

$${\mathit{K}}_{\infty }={\mathit{P}}_{\infty }^{\prime }{\mathit{C}}^{\mathit{T}}{\left(\mathit{C}{\mathit{P}}_{\infty }^{\prime }{\mathit{C}}^{\mathit{T}}+\mathit{R}\right)}^{-\mathbf{1}}$$E30

Using the Eqs. 26–29 we get:

$${\widehat{\mathit{x}}}_{\infty +1}^{\prime }=\mathit{A}{\widehat{\mathbf{x}}}_{\infty }^{\prime }+{\mathit{AK}}_{\mathbf{\infty }}\left[{\mathit{y}}_{\mathbf{\infty }}-\mathit{C}{\widehat{\mathbf{x}}}_{\infty }^{\prime }\right]$$E31

$${\mathit{P}}_{\infty +\mathbf{1}}=\mathit{A}\left(\mathit{I}-{\mathit{K}}_{\infty }\mathit{C}\right){\mathit{P}}_{\infty }^{\prime }{\mathit{A}}^{\mathit{T}}+\mathit{Q}$$E32

Using Eq. 30 in Eqs. 31 and 32 we get:

$${\widehat{\mathit{x}}}_{\infty +1}^{\prime }=\mathit{A}{\widehat{\mathbf{x}}}_{\infty }^{\prime }+\mathit{A}{\mathit{P}}_{\infty }^{\prime }{\mathit{C}}^{\mathit{T}}{\left(\mathit{C}{\mathit{P}}_{\infty }^{\prime }{\mathit{C}}^{\mathit{T}}+\mathit{R}\right)}^{-\mathbf{1}}\left[{\mathit{y}}_{\infty }-\mathit{C}{\widehat{\mathbf{x}}}_{\infty }^{\prime }\right]$$E33

$${\mathit{P}}_{\infty +\mathbf{1}}=\mathit{A}\left(\mathit{I}-{\mathit{P}}_{\infty }^{\prime }{\mathit{C}}^{\mathit{T}}{\left(\mathit{C}{\mathit{P}}_{\infty }^{\prime }{\mathit{C}}^{\mathit{T}}+\mathit{R}\right)}^{-\mathbf{1}}\mathit{C}\right){\mathit{P}}_{\infty }^{\prime }{\mathit{A}}^{\mathit{T}}+\mathit{Q}$$E34

Rewriting Eq. 34 we get:

$${\mathit{P}}_{\infty +\mathbf{1}}=\mathit{A}{\mathit{P}}_{\infty }^{\prime }{\mathit{A}}^{\mathit{T}}-\mathit{A}{\mathit{P}}_{\infty }^{\prime }{\mathit{C}}^{\mathit{T}}{\left(\mathit{C}{\mathit{P}}_{\infty }^{\prime }{\mathit{C}}^{\mathit{T}}+\mathit{R}\right)}^{-\mathbf{1}}\mathit{C}{\mathit{P}}_{\infty }^{\prime }{\mathit{A}}^{\mathit{T}}+\mathit{Q}$$E35

When in steady state:

$${\mathit{P}}_{\infty +\mathbf{1}}={\mathit{P}}_{\infty }^{\prime }={\mathit{P}}_{\infty }$$E36

Then we arrive at the Riccati equation:

$$\mathit{A}{\mathit{P}}_{\infty }{\mathit{A}}^{\mathit{T}}-\mathit{A}{\mathit{P}}_{\infty }{\mathit{C}}^{\mathit{T}}{\left(\mathit{C}{\mathit{P}}_{\infty }{\mathit{C}}^{\mathit{T}}+\mathit{R}\right)}^{-\mathbf{1}}\mathit{C}{\mathit{P}}_{\infty }{\mathit{A}}^{\mathit{T}}-{\mathit{P}}_{\infty }+\mathit{Q}=\mathbf{0}$$E37

The iterative solution of the Riccati equation is not required in real time. The observer gain is calculated off-line for predictive control applications [40]. Riccati equations are mainly used to control large scale systems, estimation, and, detection processes.

2.6 Solution to the Riccati equation using MATLAB

In this work the discrete-time algebraic Riccati equation (DARE) was solved to obtain the covariance matrix P of the Kalman gain. The discrete-time algebraic Riccati equation is represented by the next form [41]:

$$X=A^T\mathit{XA}+Q-\left(A^T\mathrm{XB}\right){\left(R+B^T\mathit{XB}\right)}^{-1}\left(B^T\mathit{XA}\right)$$E38

Where$A,X,Q=A^T\in {\mathbb{R}}^{n\times n}$, $B\in {\mathbb{R}}^{n\times m}$, $R\in {\mathbb{R}}^{m\times m}\left(m\le n\right)$, and $R=B^T>0$.

Eq. (38) can be written in the short form:

$${\mathrm{A}}^{\mathrm{T}}\mathrm{X}{\left(\mathrm{I}+\mathrm{SX}\right)}^{-1}\mathrm{A}-\mathrm{X}+\mathrm{Q}=0$$E39

Where:

$$S=B{R}^{-1}{B}^T$$E40

The application of the Kalman filter implies solving the DARE, which can be solved by several solution methods. Computational methods to solve Riccati equations can be categorized into three classes: invariant subspace methods, deflating subspace methods, and Newton’s methods. The generalized Schur method that is classified as a deflating subspace method is used to solve DARE. The generalized Schur algorithm is a strong algebraic tool that allows computing classical decompositions of matrices, such as the QR and LU factorizations [42]. The next algorithm was used to solve DARE [43]:

Input arguments:

$$\begin{array}{l}\mathrm{A}-\text{An}\mathrm{n}\times \mathrm{n}\text{matrix}\\ \mathrm{B}-\text{An}\mathrm{n}\times \mathrm{m}\text{matrix}\\ \mathrm{Q}-\text{An}\mathrm{n}\times \mathrm{n}\text{symetric matrix}\\ \mathrm{R}-\text{An}\mathrm{m}\times \mathrm{m}\text{symetrix matrix}\end{array}$$

Output arguments: $X-\mathit{\text{DARE solution}}$

1. Form the pencil $P_{\mathit{\text{DARE}}}-\lambda N_{\mathit{\text{DARE}}}$, where

   $${\mathrm{P}}_{\text{DARE}}=\left(\begin{array}{cc}\mathrm{A}& 0\\ -\mathrm{Q}& \mathrm{I}\end{array}\right),$$E41
   $${\mathrm{N}}_{\text{DARE}}=\left(\begin{array}{cc}\mathrm{I}& \mathrm{S}\\ 0& {\mathrm{A}}^{\mathrm{T}}\end{array}\right)$$E42

2. Transform the pencil $P_{\mathit{\text{DARE}}}-\lambda N_{\mathit{\text{DARE}}}$ to the generalized real Schur form apply QZ algorithm, that is, find orthogonal matrices Q₁ and Z₁ such that:

   $$Q_1{P}_{\mathit{\text{DARE}}}{Z}_1=P_1=\left(\begin{array}{cc}{P}_{11}& P_{12}\\ 0& P_{22}\end{array}\right),$$E43
   $${\mathrm{Q}}_1{\mathrm{N}}_{\text{DARE}}{\mathrm{Z}}_1={\mathrm{N}}_1=\left(\begin{array}{cc}{\mathrm{N}}_{11}& {\mathrm{N}}_{12}\\ 0& {\mathrm{N}}_{22}\end{array}\right)$$E44

3. Using an orthogonal transformation and reorder the generalized real Schur form. So that all the pencil $P_{11}-\lambda N_{11}$ has all the eigenvalues with moduli less than 1. Find Q_z and Z₂ orthogonal matrices, such that:

   $$Q_2{Q}_1{P}_{\mathit{DARE}}{Z}_1{Z}_2=\mathit{\text{quasi}}-\mathit{\text{upper triangular}}$$E45
   $$Q_2{Q}_1{N}_{\mathit{\text{DARE}}}{Z}_1{Z}_2=\mathit{\text{upper triangular}}$$E46

4. Form the matrix:

   $$\mathrm{Z}={\mathrm{Z}}_1{\mathrm{Z}}_2=\left(\begin{array}{cc}{\mathrm{Z}}_{11}& {\mathrm{Z}}_{12}\\ {\mathrm{Z}}_{21}& {\mathrm{Z}}_{22}\end{array}\right)$$E47

5. $\mathit{\text{Compute}}X=Z_{21}{Z}_{11}^{-1}$

3. Application example: resistive temperature detectors (RTD)

Resistive temperature detectors (RTD) have attracted attention to be employed as thermal health monitors. As clinical thermometers they are stable and reliable presenting high accuracy and resolution [44]. One of the most widely used RTD is the emerging thin-film resistor which has minimal impact on complex circuits due to its small size and due to their negligible mass.

The basic function of the sensor is determined by a proportional increment of resistance when temperature is applied. RTDs can be employed on a rigid or flexible substrate [45, 46, 47], the metal combination with a flexible o rigid substrate can cover conformal applications. RTD fabrication can be done by metals like Pt [48, 49, 50], Cu [51], Ag [52], and Ni [53], among other materials. Nickel presents a suitable option for RTD fabrication due to its wide temperature linear range of operation and its relatively low price.

Clinical thermometers require a high definition and reliability because less than 1°C difference can indicate a health problem. The thermometer signal can be amplified by electronic means, but it is desirable to filter such readings. This work is focused to the filtering and prediction of an highly sensitive Nickel based thin film RTD (range, 273–325 K), to be incorporated to complex circuits [54], we present the theoretical analysis about the relation sensibility-resistance that matches with experimental results.

3.1 Theoretical analysis of an RTD

All metals produce an increase in its resistance to an increase in specific temperature, which means that resistance is linearly proportional to temperature change. This dependence between electrical resistance and temperature is the principle of operation used by a resistance temperature detector (RTD). The relation between temperature-resistance for Pt wire (RTD), is described by the equation known as the Calendar-Van Dusen, Eq. 41) [50].

$$R_{\left(T\right)}=R_{0°C}\left(1+\mathit{\alpha T}+{\mathit{\beta T}}^2\right)$$E48

Where R_0°C is the resistance at 0°C, α and β are temperature coefficients and T is temperature, the temperature coefficients depend only on material properties. In addition, the RTD resistance depends on its geometrical design, according to Eq. 42.

$$\mathrm{R}=\frac{\mathrm{\sigma }\ast \mathrm{L}}{\mathrm{A}}=\frac{\mathrm{\sigma }\ast \mathrm{L}}{\mathrm{w}\ast \mathrm{t}}$$E49

Where “σ” is the resistivity, “L” length, “A” lateral area, “w” channel width, and “t” channel height. Only by increasing the length “L” or decreasing the area “A” that means reducing the “t” film thickness or the “w” channel wide, the resistance can increase.

3.1.1 State-space description of an RTD

The estimation of the thermal system is represented by the linear stochastic state-space description $x_k=A{x}_{k-1}+{\mathit{Bu}}_{k-1}+w_{k-1}$ and $y_k=C{x}_k+v_k$. Where A is an nxn state transition matrix applied to the previous state vector $x_{k-1}$, B is the control-input matrix applied to the control vector $u_{k-1}$, and $w_{k-1}$ is the process noise vector. The linear combination of the measurement noise and the signal value is represented by $y_k$, where $C$ is the measurement matrix, and $v_k$ is the measurement noise vector with covariances matrices represented by Q and R. The covariances are assumed to be independent and are given by $Q=E\left[w_k{w}_k^T\right]$ and $=E\left[v_k{v}_k^T\right]$ .

Generally, the RTD system is modeled as an RLC circuit, which consists of a resistor a capacitor and an inductor in series with an input voltage. The output that we analyzed is the voltage across the resistor which is related to temperature change. The RLC circuit is represented by a second-order differential equation $L\frac{d^{2}i\left(t\right)}{{\mathit{dt}}^2}+R\frac{\mathit{di}\left(t\right)}{\mathit{dt}}+\frac{1}{C}i\left(t\right)=0$. To solve the above equation we implement the next matrix system:

$$\mathrm{A}=\left(\begin{array}{cc}0& 1\\ -1/\mathrm{L}·\mathrm{C}& -\mathrm{R}/\mathrm{L}\end{array}\right)\mathrm{B}=\left(\begin{array}{c}0\\ 1/\mathrm{L}·\mathrm{C}\end{array}\right)$$E50

$$\mathrm{C}=\left(10\right)\mathrm{D}=0$$

Also, we may simplify the response of the system to that of a first-order RC circuit. This implies to solve a first-order ordinary differential equation: $\mathit{RC}\frac{\mathit{dq}}{\mathit{dt}}+q=\mathit{VC}$. The dynamic model is defined by the following system:

$$\mathrm{A}=\left(-\frac{1}{\mathrm{R}·\mathrm{C}}\right)\mathrm{B}=\left(\frac{1}{\mathrm{R}·\mathrm{C}}\right)$$E51

$$\mathrm{C}=1\mathrm{D}=0$$

4. Design and simulations

4.1 Kalman filter in resistance thermal detectors (RTD)

In this work, a Kalman Filter is proposed to decrease the time response to improve the speed feedback and filtering of the perturbations by signal noise from physical signals as thermal detectors. In some instances, a reduced model is advisable to use in an embedded system due to easy implementation and low computational complexity [2].

Kalman filter can be embedded in a temperature system made by Resistance Thermal Detectors (RTD).RTD’s are robust elements that require relatively easy measurement, as a consequence are a useful thermal sensor for industry and medical applications. Nevertheless, these devices are exposing to vibration, electrical noise, and measurement errors generated by the thermoelectric effect caused by the temperature difference between electrical contacts, which affects the response time of the sensor. The implementation of the Kalman filter in a temperature system produces an optimal estimative of thermal behavior and decreases the uncertainties about the prediction of the temperature.

In order to describe the system in the state space, it is necessary to apply system identification methods using MATLAB. Then, after obtaining the system’s state space model we are able to use the Kalman filter algorithm to estimate the future output of the system.

To study the dynamics of our system, we used MATLAB functions etfe and spa to firstly estimate the empirical transfer functions and then estimate the frequency response with fixed frequency resolution using spectral analysis. The continuous time-identified transfer function obtained is:

$$\frac{2.278\mathrm{s}+0.1711}{{\mathrm{s}}^2+2.488\mathrm{s}+0.1695}$$E52

Using MATLAB we are able to acquire the Discrete-time identified state-space model:

$$\mathrm{x}\left(\mathrm{t}+\mathrm{Ts}\right)=\mathrm{A}\mathrm{x}\left(\mathrm{t}\right)+\mathrm{B}\mathrm{u}\left(\mathrm{t}\right)+\mathrm{K}\mathrm{e}\left(\mathrm{t}\right)$$E53

$$\mathrm{y}\left(\mathrm{t}\right)=\mathrm{C}\mathrm{x}\left(\mathrm{t}\right)+\mathrm{D}\mathrm{u}\left(\mathrm{t}\right)+\mathrm{e}\left(\mathrm{t}\right)$$

with:

$$A=\left[\begin{array}{cc}0.8342& -0.08908\\ 0.0942& -0.9716\end{array}\right],B=\left[\begin{array}{c}0.01966\\ -0.03341\end{array}\right],$$E54

$$C=\left[\begin{array}{cc}7.966& 0.3005\end{array}\right],D=0,K=\left[\begin{array}{c}0.006289\\ -0.2834\end{array}\right]$$

Estimated using N4SID on time domain data. Fit to estimation data: 90.27% (prediction focus) with FPE: 0.4532 and MSE: 0.2292. Figure 2 shows the Input–output model for which the input was set to a constant value of 38°C. The output, the step response, is that of the second order system as can be seen in the bode plot shown in Figure 3. Figure 4 shows the evolution of the measured versus the modeled step responses.

Figure 2.

System ID using MATLAB. Input output model for a step response defined problem.

Figure 3.

Bode diagram indicating the system is a second order system as described by the system transfer function.

Figure 4.

Systems model and measured evolutions in time. Fit to estimation data: 90.27%.

4.2 Kalman filter

We modify the MATLAB example for the time-varying case found in [55] and we code our own function to solve the Discrete Algebraic Riccati Equation. MATLAB functions like predict or forecast were found useful to understand the problem at hand, however they were not used in the code we present here.

4.2.1 Time-varying Kalman filter using MATLAB

w(1:n) = sqrt(Q)*randn(n,1);

v(1:n) = sqrt(R)*randn(n,1);

systv = ss(A,B,C,0,Ts);

ytv(1:n) = lsim(systv,U(1:n) + w(1:n)).

yvtv(1:n) = ytv(1:n) + v(1:n);

Ptv(:,:) = B(:,:)*Q*B(:,:)’;   % Initial error covariance.

x = zeros(order,1);  % Initial condition on the state.

order = 2;

yetv(1:n) = zeros(n,1);

ycov(1:n) = zeros(n,1);

for i = 1:n.

% Measurement update.

Mn(:,:) = Ptv(:,:)*C(:,:)’/(C(:,:)*Ptv(:,:)*C(:,:)’ + R);

x = x + Mn(:,:)*(yvtv(i)-C(:,:)*x); % x[n|n].

Ptv(:,:) = (eye(order)-Mn(:,:)*C(:,:))*Ptv(:,:);  % P[n|n].

yetv(i) = C(:,:)*x;

errcov(i) = C(:,:)*Ptv(:,:)*C(:,:)’;

% Time update.

x = A(:,:)*x + B(:,:)*U(i);  % x[n + 1|n].

Ptv(:,:) = A(:,:)*Ptv(:,:)*A(:,:)’ + B(:,:)*Q*B(:,:)’; P[n + 1|n].

end

%% DARE. We coded our own dare function [X,L,G] = sdare(A,B,Q,R).

[P_inf,L,M_inf] = sdare(atv,ctv’,Q,R);

for i = 1:p

% Measurement update.

x = x + M_inf’*(yvtv(i)-ctv*x); % x[n|n].

yetv_inf(i) = ctv*x;

errcov_inf(i) = ctv*P_inf*ctv’;

% Time update.

x = atv*x + btv*U(i); % x[n + 1|n].

P_inf = atv*P_inf*atv’ + btv*Q*btv’; % P[n + 1|n].

end

function [SD] = sdare(A,B,Q,R).

At = transpose(A);

Bt = transpose(B);

S1 = size(A);

E = eye(S1);

Z = zeros(S1);

Ri = inv.(R);

S = B*Ri*Bt;

Pdare = [A Z; −Q E];

Ndare = [E S; Z At];

[AA,BB,L,Z] = qz(Pdare,Ndare);

[AAS1,BBS1,QS1,ZS1] = ordqz(AA,BB,L,Z,‘udi’);

O = ZS1(1:2,1:2);

P = ZS1(3:4,1:2);

H = inv.(O);

SD = P*H;

end

5. Simulations

Matlab was used to simulate the response of an RTD modelled as a second order system. In Figure 5(A) we show the plot of the true response y (cyan line) and the filtered response (red line). In Figure 5(B) the plot compares the measurement error with the estimation error. As can be seen in Figure 5(C) the time-varying filter also estimates the covariance errcov of the estimation error at each sample which shows when the filter reached steady state. As it can be seen, we have the possibility to predict the state after approximately 8 seconds. Also, we show the evolution of the estimated temperature response showing an error of −0.0948°C in the best of the cases and less than 1°C in the worst of the cases after 45 seconds.

Figure 5.

(A) Evolution of the estimated temperature response showing an error of −0.0948°C in the best of the cases and less than 1°C in the worst of the cases. (B) Evolution of the measurement and estimation errors. (C) Evolution of the covariance of the error showing the possibility to predict the state after approximately 8 seconds.

6. Implementation

The unit step response depends on the roots of the characteristic equation. If both roots are real-valued, the second-order system behaves like a chain of two first-order systems, and the step response has two exponential components. If the roots are complex, the step response is a harmonic oscillation with an exponentially decaying amplitude [56]. In our case, the roots of the characteristic polynomial: ${\mathrm{s}}^2+2.488\mathrm{s}+0.1695$ are −2.4179 and − 0.0701. Thus our system behaves like two first order systems in series.

The state description for an RC system is described above. From there, we know that the dynamics are dependent only on the RC constant. In addition, there is an amplificator in the system electronics that has a gain of 260. To solve for the RC constant of the system we use the least-squares method (Chi square minimization). The system has a solution of the form $y=e^{\mathit{Bx}}$ and we take n data points to form the vectors $x_i$ and $Y_i$. The problem is to minimize the error function, $\mathit{err}={\sum }_{i=1}^n{\left(Y_i-A{e}^{B{x}_i}\right)}^2$. The trick on the algorithm goes as follows:

$${\mathit{Yn}}_i=\mathit{ln}{Y}_i=\mathit{ln}\left(A{e}^{B{x}_i}\right)=\mathit{lnA}+ln\left(e^{B{x}_i}\right)=C+{\mathit{Bx}}_i$$E55

Which is a linear equation. Using a linear fitting program:

$$\mathrm{B}=\frac{\mathrm{n}\sum {\mathrm{x}}_{\mathrm{i}}{\mathrm{y}}_{\mathrm{i}}-\sum {\mathrm{x}}_{\mathrm{i}}\sum {\mathrm{y}}_{\mathrm{i}}}{\mathrm{n}\sum {\mathrm{x}}_{\mathrm{i}}^2-{\left(\sum {\mathrm{x}}_{\mathrm{i}}\right)}^2}\text{and}\mathrm{c}=\frac{\sum {{\mathrm{x}}_{\mathrm{i}}}^2\sum {\mathrm{y}}_{\mathrm{i}}-\sum {\mathrm{x}}_{\mathrm{i}}\sum {\mathrm{x}}_{\mathrm{i}}{\mathrm{y}}_{\mathrm{i}}}{\mathrm{n}\sum {\mathrm{x}}_{\mathrm{i}}^2-{\left(\sum {\mathrm{x}}_{\mathrm{i}}\right)}^2}$$E56

We obtain $B$ and $A=exp\left(c\right)$. We coded two Kalman filters with different calibration parameters (Q and R) as is written below. The KF was implemented in an ATMEGA328 microprocessor. Code for an Arduino was generated using C language which is available in the following section. Only the relevant portion is written.

6.1 Arduino code

readings[readIndex] = analogRead(inputPin); // read from the sensor.

total = total + readings[readIndex]; // add the reading to the total.

readIndex = readIndex +1; // advance to the next position in the array.

time_equis_readings[time_equis_readIndex] = time_equis_readIndex;

time_equis_readIndex = time_equis_readIndex +1;

if (readIndex > = numReadings) // if we’re at the end of the array.

{

  for(i = 0;i < =numReadings-1;i++).

   {

   Y[i] = log(readings[i]);

   time1[i] = time_equis_readings[i];

   sumx = (sumx +time_equis_readings[i]);

   sumx2 = (sumx2 + time_equis_readings[i]*time_equis_readings[i]);

   sumy = (sumy +Y[i]);

   sumxy = (sumxy +time_equis_readings[i]*Y[i]);

   }

 den = (numReadings*sumx2-sumx*sumx);

 a = (sumx2*sumy -sumx*sumxy)/den;

 Bc = (n*sumxy-sumx*sumy)/den;

// State description.

 A = -Bc;B=Bc;C = 260;D = 0;

//wrap around to the beginning:

 readIndex = 0;time_equis_readIndex = 0;

 }

// KALMAN.

 errcov = C*P*C;

 for(i = 0;i < =numReadings-1;i++).

 {

  Mn = P*C/((C*P*C + R)); // initial estimate.

  X = X + Mn*(readings[i]-C*X); // update estimate Average_readings[i].

  P = (1-Mn*C)*P; // update covariance.

y_e[i] = C*X;

  errcov = C*P*C;

X = A*X + B*U; // project into k + 1.

  P = A*P*A + B*Q*B; // project into k + 1.

}

 timer0_millis = millis();

 // Solution to the Riccati equation.

 F = -Bc;H = 260;

 SQ = sqrt((H*H*Q*R) + (F*F*R*R));

SR = F * R;

P_inf = (SQ + SR)/(H*H);

M_inf = P_inf*C/(C*P_inf*C + R);

for(i = 0;i < =numReadings-1;i++).

 {

  // Measurement update.

  // M_inf;

  x_inf = x_inf + M_inf*(readings[i]-C*x_inf); // % x[n|n].

//P_inf; % P[n|n].

  y_e_inf = C*x_inf;

errcov_inf = C*P_inf*C;

  // Time update.

  x_inf = A*x_inf + B*U;

  P_inf = A*P_inf*A + B*Q*B;

 }

}

7. Results

The experiments were performed in a thermal bath giving step responses to the desired setup temperature. Figure 6 depicts the upward and downward evolution of the temperature, the Kalman filter and the two predictors (using two different Q and R settings). As can be seen the predictors follow the Temperature of the sensor closely, especially for the upward way, while the Kalman filter lags behind.

Figure 6.

Implemented system. Step response for the upwards and downwards evolution. Two different Kalman filters were used to predict (by solving the DARE equation) the evolution of the future state with different Q and R to calibrate the desired response. In blue the evolution of the RTD sensor analog input, in Yellow and red the two Kalman predictors and the Kalman estimation in cyan color.

8. Boundary layer

Sliding control [57] is an additional tool to predict the behavior of a second order system basically smoothing the system by boundary layers. The prediction of the system state trajectory is given using an uncertain model of the system. The subspace which represents the quantity of uncertainties in the prediction process, forces the estimate state trajectory to switching gain to converge the estimates to within a boundary of the real state values. To predict the state trajectory of our RLC system it’s possible to switch its gain by the subspace represented by a first-order RC model. The estimated state trajectory is forced to keep a switch back and forth within the boundary layer represented in our case by a RC model. By creating a boundary layer, the system is further constrained to have a solution existing in between two RC model solutions.

In Figure 7 it can be clearly seen that the use of two estimators may help predict the behavior of the RTD in a much better way. The system needs to be calibrated first in order to have the two Kalman filters enveloping the required solution. As can be seen in the upward direction, both predictors (yellow and red) envelope the desired response (blue), that of the RTD sensor improving the response of the Kalman filter without boundaries (cyan). Unfortunately, this is not the case in the downward evolution. From the nonlinear control systems point of view these two evolutions demark a region where the RTD stands thus making possible to program a better estimator. It is left as an outlook to program a third estimator using this boundary layer in order to have a better predictor, especially for the downward evolution.

Figure 7.

Ascending and descending step responses of the Kalman filter and two Predictors which function in real time. In blue the RTD sensor response, in cyan the estimator response, in yellow and in red the two differently calibrated Kalman predictors.

9. Conclusions

As it can be shown the implementation of the Kalman filter brings the opportunity to estimate the forecast in real time of a second order system using first, MATLAB and second that of two first order systems using a simple RC system coded in C-language for a microprocessor. It has been shown that the program is able to predict the evolution of temperature for a RTD system. Even if the system is implemented using a first order system we can find evolving solutions for our estimation and prediction to be good enough. We predict the state after approximately 8 seconds showing an error of −0.0948°C in the best of the cases. In addition, a boundary layer may be programmed using two first order Kalman predictors which may be tuned by setting Q and R properly. We believe this is the first report on the use of a Kalman filter to predict the evolution of temperature from a RTD.