Giving uncertainties to interpolated values in time series: application to mass balance reconciliation

Question

My problem concerns a mass balance reconciliation in an industrial system.

I have a node with several flows in-coming and one or several out-coming, what's incoming should equals what's out-coming (mass conservation, Lavoisier 1777). Of course, when you perform direct or indirect measurements in real life, it is rarely the case (measurments errors, aparatus precision, mistakes, sampling errors, ..). Fortunately, many methods to reconciliate the data exit. In our case, it is a problem of optimisation: minimising an objective function while respecting constraints.

Here is a presentation illustrating my words:

https://www.slideshare.net/AntoineCollet9/20190804massbalancereconciliationslideshare

And here some chunks of code both in R (Rmd) and python, linked with the presentation:

https://gitlab.com/antoinecollet5/datareconciliation

Now that you get the big picture, let me introduce the issue I am facing:

In the previous example, I am treating a case where all flows have been measured and then, the uncertainties are the same i.e. it is the aparatus/method precision.

Now, lets assume that all variables are in fact day-based time series. I will have to perform this reconciliation each day, which is very simple to do from the previous code (simply by adding a for loop).

In real life of course, some measures will be missing (problems with the aparatus, wrong reporting, etc.) or removed because considered as outliers. I then interpolate all missing data (linear, cubic splines, kalman,...). The big question now is how to calculate the uncertainties associated with these interpolated values ? I was not able to find papers/examples clearly talking about that... I guess that the further from measured data, the higher the uncertainty should be, but how to quantify it from the base uncertainty? Is there a rule of thumb ? I precise that I know it will depends from the interpolation methods and that some already include uncertainties calculation. I am looking for something more general...

Here is a graph to illustrate my words:

Thank you in advance for your help :)

Answer 1

The Kalman filter seems to answer your need :)

https://en.wikipedia.org/wiki/Kalman_filter

Answer 2

As Yrieix suggested, Kalman is indeed a nice way to answer that problem. In python, there is a library available: pykalman.

There are already post on how to handle missing values with pykalman.

I took the code given by Anton there and I simply added a smoothing:

# estimate hidden states
(smoothed_state_means, smoothed_state_covariances) = kf.smooth(X)

The result is quite amazing:

Here is the full code code to obtain that. The raw data are the ones from Anton's post and are available here.

from pykalman import KalmanFilter
import numpy as np
from numpy import ma
import matplotlib.pyplot as plt


# reading data (quick and dirty)
Time=[]
X=[]

for line in open('Kalman_dataset.csv'):
    f1, f2  = line.split(';')
    Time.append(float(f1))
    X.append(float(f2))

#Masking some values
X = ma.asarray(X)
X[200:750] = ma.masked

# Filter Configuration

# time step
dt = Time[2] - Time[1]

# transition_matrix  
F = [[1,  dt,   0.5*dt*dt], 
     [0,   1,          dt],
     [0,   0,           1]]  

# observation_matrix   
H = [1, 0, 0]

# transition_covariance 
Q = [[   1,     0,     0], 
     [   0,  1e-4,     0],
     [   0,     0,  1e-6]] 

# observation_covariance 
R = [0.04] # max error = 0.6m

# initial_state_mean
X0 = [0,
      0,
      0]

# initial_state_covariance
P0 = [[ 10,    0,   0], 
      [  0,    1,   0],
      [  0,    0,   1]]

n_timesteps = len(Time)
n_dim_state = 3

filtered_state_means = np.zeros((n_timesteps, n_dim_state))
filtered_state_covariances = np.zeros((n_timesteps, n_dim_state, n_dim_state))

# Kalman-Filter initialization
kf = KalmanFilter(transition_matrices = F, 
                  observation_matrices = H, 
                  transition_covariance = Q, 
                  observation_covariance = R, 
                  initial_state_mean = X0, 
                  initial_state_covariance = P0)

# estimate hidden states
(smoothed_state_means, smoothed_state_covariances) = kf.smooth(X)

position_sigma = np.sqrt(smoothed_state_covariances[:, 0, 0])

# plot of the resulted trajectory
plt.plot(Time, smoothed_state_means[:, 0], "g-", label="Filtered position", markersize=1)
plt.plot(Time, smoothed_state_means[:, 0] + position_sigma, "r--", label="+ sigma", markersize=1)
plt.plot(Time, smoothed_state_means[:, 0] - position_sigma, "r--", label="- sigma", markersize=1)

plt.grid()
plt.legend(loc="upper left")
plt.xlabel("Time (s)")
plt.ylabel("Position (m)")
plt.savefig(fname = 'Kalman_5.png', dpi = 300,
            bbox_inches='tight')
plt.show()