Predicting process behavior using FFT and FIR models

Predicting the future behavior of a system based on past data is a cool idea and can be used in decision making. Model predictive control (MPC) uses this technique to make feedforward actions anticipating the future behavior of the system. Traditionally the systems are represented by transfer functions, but when we use them for large scale control appellations, time to solve the optimization problem may exceed the sampling frequency (Typically 60 seconds) of controller. So most of the commercially scalable MPC adopt a time domain representation of transfer function (Eg: Honeywell Profit Controller, Aspentech DMCPlus/DMC3).

Transfer function is the Laplace transform of impulse response of a linear time invariant (LTI) system when initial conditions set to zero.

Finite Impulse Response (FIR)

Finite impulse response models are simple 1D floating point array but can represent almost any dynamic behavior of a system. A sample numpy array representing a process dynamics (A distillation column) of 120 minutes is given below.

array([ 5.33662890e-03,  1.08608828e-03, -8.45890854e-04, -2.34296534e-03
       -3.24092270e-03, -3.52775132e-03, -3.28877642e-03, -2.66050355e-03,
       -1.79460808e-03, -8.32199356e-04,  1.12431318e-04,  9.59946480e-04,
        1.66764266e-03,  2.22640709e-03,  2.65459311e-03,  2.99025901e-03,
        3.28295687e-03,  3.58597613e-03,  3.94966951e-03,  4.41623635e-03,
        5.01612323e-03,  5.76602926e-03,  6.66837611e-03,  7.71201759e-03,
        8.87391581e-03,  1.01214949e-02,  1.14153920e-02,  1.27123514e-02,
        1.39680472e-02,  1.51396623e-02,  1.61881020e-02,  1.70797607e-02,
        1.77878058e-02,  1.82929741e-02,  1.85839060e-02,  1.86570644e-02,
        1.85162996e-02,  1.81721288e-02,  1.76408033e-02,  1.69432356e-02,
        1.61038525e-02,  1.51494338e-02,  1.41079886e-02,  1.30077094e-02,
        1.18760366e-02,  1.07388550e-02,  9.61983580e-03,  8.53992938e-03,
        7.51700810e-03,  6.56565165e-03,  5.69706410e-03,  4.91910845e-03,
        4.23644192e-03,  3.65073468e-03,  3.16095387e-03,  2.76369560e-03,
        2.45354810e-03,  2.22347082e-03,  2.06517590e-03,  1.96950012e-03,
        1.92675746e-03,  1.92706443e-03,  1.96063211e-03,  2.01802064e-03,
        2.09035375e-03,  2.16949217e-03,  2.24816629e-03,  2.32006935e-03,
        2.37991361e-03,  2.42345250e-03,  2.44747221e-03,  2.44975694e-03,
        2.42903169e-03,  2.38488708e-03,  2.31769017e-03,  2.22848552e-03,
        2.11889006e-03,  1.99098545e-03,  1.84721085e-03,  1.69025895e-03,
        1.52297753e-03,  1.34827833e-03,  1.16905482e-03,  9.88109825e-04,
        8.08093705e-04,  6.31453421e-04,  4.60392477e-04,  2.96841475e-04,
        1.42438786e-04, -1.47938193e-06, -1.33880350e-04, -2.54028744e-04,
       -3.61473716e-04, -4.56030691e-04, -5.37758726e-04, -6.06934542e-04,
       -6.64024262e-04, -7.09653811e-04, -7.44578842e-04, -7.69654993e-04,
       -7.85809151e-04, -7.94012295e-04, -7.95254403e-04, -7.90521779e-04,
       -7.80777061e-04, -7.66942090e-04, -7.49883697e-04, -7.30402423e-04,
       -7.09224083e-04, -6.86994051e-04, -6.64274068e-04, -6.41541366e-04,
       -6.19189834e-04, -5.97532967e-04, -5.76808313e-04, -5.57183116e-04,
       -5.38760903e-04, -5.21588709e-04, -5.05664719e-04, -4.90946069e-04]),        

So what can we infer from this curve? not rally much! What about step response? We can convert impulse response to step response just by calling numpy cusum method on the same array.

step_response = np.cumsum(impulse_response)        

What about bode plot?

Frequency response of the system in the Fourier transform of impulse response!

I hope, now you might have understood how fundamental is impulse response of a system. From impulse response we can generate all other representations of LTI system.

System Simulation

Linear system simulations are an important tool to study what-if analyses and estimate model quality by comparing measurement and perdition plots. In time domain system response can be simulated using convolution of the input signal and the system impulse response.

This a computationally expensive operation, the same operation in frequency domain is a simple multiplication. In terms of number of operations convolution is Quadratic whereas frequency domain transform, then multiplication, finally inverse transform altogether is N log N.

Quadratic is worse than N log N

Do we need to worry about implementation?

Not really! Scipy.signal has a neat implementation of this (Scipy.signal.convolve). By default, convolve use method='auto', which calls choose_conv_method to choose the fastest method using pre-computed values (choose_conv_method can also measure real-world timing with a keyword argument).

A simulate Python function

def simulate_fir(fir_model, data)
    """
    Returns a pandas dataframe of dependant variables predictions.
        
        Parameters
        ----------
        fir_model (dict(dict(numpy.array))): nested dictionary of numpy array containig FIR coeficiants
        data (pandas.DataFrame): pandas data frame containing independant and dependant variables data.


        Returns
        -------
        predictions (pandas.DataFrame): pandas data frame containing dependant variables predictions.
    """
    N = len(data)
    deps = list(fir_model.keys())
    inds = list(fir_model[deps[0]].keys())
    predictions = data[deps].copy(deep=True)
    for dep in deps:
        predictions[dep].values[:] = 0.0
        for ind in inds:
            predictions[dep] += signal.convolve(data[ind], fir_model[dep][ind], mode='full')[:N]
        tss = len(fir_model[dep][inds[0]])
        predictions[dep][:tss] = predictions[dep].values[tss+1]
        predictions[dep] = predictions[dep] - predictions[dep].mean() + data[dep].mean()
    return predictions:        

Conclusion

Scipy has an efficient implementation for array convolution and it is computationally efficient. This may be used for MPC applications as perdition tool.

Notebooks

References

DSP Tutorial PyCon 2017