Fitting a Second Order Model to Input Data

Let’s imagine that we have a second order system characterized by the following characteristics:

• System gain = 2.0
• Damping coefficient = 0.5
• Natural frequency = 1.0

This system’s theoretical response to a unit step is show below:

And, let’s imagine that the data we record about this system’s input and output isn’t exactly perfect and has lots of noise (I left the real system’s output in red for comparison against the noisy data in green).

Now, let’s imagine we want to learn more about this system. That we want to learn important things that define this system’s characteristics. Mainly, its overall gain, damping coefficient and natural frequency.

That’s what this article is about. I have created a python class which receives some data and tries to fit a second order model to it.

Below are the results I get after passing that noisy data to my SecondOrderFitter.fit() method.

• System gain = 2.0314
• Damping coefficient = 0.4852
• Natural frequency = 0.9811

Pretty awesome, isn’t it? The class was able to characterize the data pretty well. The picture below shows the difference in response between the real system (without noise) and the system fitted by the SecondOrderFitter algorithm.

This is what this article is about.

Intro

It has been some time since I last wrote here. My professional commitments and personal life didn’t leave much time for me to play and have fun with python and my small custom-built simulator. But, now that I am on a business trip, and since there isn’t much for me to do at the hotel at night, I decided to build a python module that would fit a second order model to a system’s input/output data.

I decided to create this module with the goal to understand better the dynamics of a device based solely on the data measured from its input/outputs. That could come handy next time I visit a power plant and need to get a better picture of the dynamics related to a certain actuator.

I have added this module inside my simulator project, called the class that handles second order models SecondOrderFitter and its entire code can be found here. If you are curious about my small simulator project, you can read more about it on my blog (you can start with this post, for example; but more are available here).

The model

The model I want to fit to the data to is a second order model. Therefore, it’s expected that the device from which the data is collected have only two poles or at least only two dominant poles. And, for now, the model only supports SISO (single input, single output) systems. The second order system equation that defines such a system can be given by the equation below:

Where omega represents the system’s natural frequency, zeta its damping coefficient and K the system gain.

The same equation can be rewritten on the time domain as follows:

If we define the states of the model as:

We can write the following set of state equations that define the model:

This set of ODEs can then be easily implemented in Python by creating a method and passing to it, besides the time t and states x, the parameters u (input), k, e, and w (system gain, damping coefficient and natural frequency, respectively):

def model(t, x, u, p):
    # Unpack parameters
    k, e, w = p
    # Unpack states
    x0, x1 = x
    # Calculate derivatives as defined by the state equations ODEs
    dx = np.zeros(2)
    dx[0] = x1
    dx[1] = k*w**2*u - (2*e*w*x1 + w**2*x0)
    return dx

The optimization problem

With the model in hands, we can create a method that integrate those ODEs in order to calculate the system’s output y by using the solve_ivp method from the scipy.integrate library. The method below contains an attribute called data that holds the data we want to fit the model to. Therefore, lines represents the number of lines in this dataset.

def integrate(self, p):
    lines = self.data.shape[0]
    y = np.zeros(lines)
    # initial output
    y[0] = self.data[0,2]
    # Considering system is stable at initial state (y' = 0)
    x0 = [y[0], 0.0]
    for line in range(1,lines):
        x = solve_ivp(
            self.model,
            (self.data[line-1,0], self.data[line,0]),
            x0,
            args = (self.data[line,1], p)
        )
        # Update state
        x0 = x.y[:,-1]
        y[line] = x.y[0,-1]
    return y

With the calculated y, we can use scipy once again, but this time the scipy.optimize library and its minimize method to search the parameters (k, e and w) that minimizes the error between a simulated second order model and the real data passed to the class. But, before, we need to define the function that calculates the error.

def objective(self, p):
    y = self.integrate(p)
    obj = 0
    for index, value in enumerate(y):
        obj += (value-self.data[index,2])**2
    return obj

Then, we can pass this cost function to scipy's minimize method.

def fit(self, x0 = None):
    if not x0:
        x0 = self.x0
    result = minimize(self.objective, x0)
    self.k, self.e, self.w = result.x
    return = result.success

With this setup, I can pass data while initializing an instance of my SecondOrderFitter class and call its fit() method to try to find the best model’s parameters that reduce the error between the fitted model and the data passed to the class.

Next steps

• Generalize the concept and create a base class from which I can create fitters for other types of models.
• Create fitter models for first second models, n-order models, specific actuators, etc.
• Try this thing out with real field data!

Originally published at https://fabiomolinar.com on January 21, 2020.