@rlabbe Thank you so much for this resource. It is invaluable for gaining intuition with Kalman Filters. After reading the first 4 chapters I have had a go at implementing a KF for my own toy model.

Ultimately, I want to improve the output of my model simulations by updating estimates with new measurements when they come in. Great it sounds like a Kalman Filter can help!

Note

All the code for the below problem is available HERE

The hydrological Model

This is the ABC Model - see Appendix for simulating discharge, a super simple hydrological model!

The ABC Model for simulating discharge (Q)

The equations are relatively simple

Q[t] = Qf[t] +Qs[t]
Qs[t] = S[t] * c
Qf[t] = P[t] * (1 - a - b)
S[t] = S[t - 1] * (1 - c) + (P[t] * a)
L[t] = P[t] * b

Where

Parameters
a: Proportion of P[t] passing into S[t]
b: Proportion of P[t] lost (e.g. to evaporation)
c: Proportion of S[t] passing into Qs[t] (slow flow)

Values
P: Precipitation Inputs
L: Losses 
S: Storage (groundwater / soil)
Qf: Fast Flow
Qs: Slow Flow
Q: River Discharge

The model isn't perfect but I calibrated some of the parameters for a small catchment in the United Kingdom where I have rainfall and precipitation (rainfall). Here is a run of the model with parameters fit from historical data.

ABC Model against Observed Discharge

1D Kalman Filter

I am using the filterpy.kalman implementation of predict and update.

Data Simulation

I set up a mini-experiment before trying to filter the model and observations from the real data. Instead I simulated data from the ABC model.

I used the observed Precipitation (without added noise) to drive the ABC Model, calculating a simulated discharge (Qsim) for station 39034. I use this as the ‘unobserved truth’.

I then add gaussian noise to the ‘unobserved truth’ with a variance of “measurement error”. This means that we have data as below. (Note that the 'true' data is now the same as the orange line in the above plot, I am pretending a world where we know the actual discharge completely).

Simulated Data from the ABC Model

Predict Step

The way that I have understood Kalman Filters is that we want to use this process model to model u (dx as in your earlier chapters).

In the predict step of the KF, I want to use the ABC model to make a prior prediction of the discharge.

One source of confusion is that the variable we are measuring and is of interest (discharge - Q) is not fed back into the model. The only obvious state variable is S - the storage component, but since we are not observing this I don't know how to use it in the kalman filter.

Another source of confusion is that the ABC model produces a simulation of Q, not of the change of Q over time. I altered this by calculating dQ/dt by hand:

dx = Qsim - X

So the overall predict step looks like this:

    for ix, precip in enumerate(data["precip"]):
        # predict step
        Qsim, storage = abcmodel(storage, precip)   # process model (simulate Q)
        dx = Qsim - X                               # convert into a change (dx)
        X, P = kf.predict(x=X, P=P, u=dx, Q=Q)
        predicted_values.append(X)

Update Step

In the update step I use the simulated measurement of discharge to update the estimate.

        # update step
        z = data['qsim_noise'][ix]                  # measurement
        X, P, y, K, S, log_likelihood = kf.update(x=X, P=P, z=z, R=R, return_all=True)

Results

Now this more or less works

a) We have the filter between the modelled and observed outputs.

Filter, Predict, Measured

b) But, because of the experimental setup, the model is perfectly reproducing the observations (maybe I should add noise to the input precipitaiton instead).

The filter is also very noisy as you can see!

True vs. Filtered Results

Questions

Firstly, I wanted to know what you think about using the 1D kalman filter for this problem?

It isn't clear to me how to translate this back to the original problem shown in the first plot (ABC Model against Observed Discharge) where the model isn't great and my observations are the thing that i'm evaluating the model on as 'truth'.
Secondly, do I need the ABC model to take filtered discharge (X) from the filter to feed back into it's predictions for this method to be useful?

Thankyou!

If you do get a chance to read this thank you. Thank you even if you don't read this for making this book and the filterpy package. I will continue working on it and am happy to provide an update if it interests you!

Best Wishes Tommy

rlabbe / Kalman-and-Bayesian-Filters-in-Python

Help with Kalman Filter for Hydrological Model #350