Mass-Spring-Damper Model Estimation

System Description

This is a really simple model: A 1-DOF mass-spring-damper system.

In [1]:
from IPython.display import Image

Image(url='https://github.com/stuckeyr/msd/raw/master/mass_spring_damper.png')
Out[1]:

See here for a more detailed explanation.

The equations of motion can be written:

$$ \sum F_x = F(t) - b \dot{x} - kx = m \ddot{x} $$

where $x$ is the displacement, $m$ is the mass, $k$ is the spring constant, $b$ is the damping constant and $F$ is the input force, expressed here as a function of time.

Reframing the system in state-space form, using the state vector:

$$ \mathbf{x} = \begin{bmatrix} x \\\ \dot{x} \end{bmatrix} $$

the equations are written:

$$ m \mathbf{\dot{x}} = \begin{bmatrix} \dot{x} \\\ \ddot{x} \end{bmatrix} = \begin{bmatrix} 0 & m \\\ -k & -b \end{bmatrix} \begin{bmatrix} x \\\ \dot{x} \end{bmatrix} + \begin{bmatrix} 0 \\\ 1 \end{bmatrix} F(t) $$

Simulation in Python

First, reset the workspace. Import all necessary libraries (the following assumes you have built both Boost and Cython models).

In [2]:
%matplotlib notebook
# %matplotlib inline
In [3]:
%reset -f
In [4]:
import numpy as np
import numpy.matlib as ml
from scipy import interpolate, integrate

Create the model. Below is the pure Python version, which I have included for information.

In [5]:
# ------------------------------------------------------------------------------
# MSD class
# ------------------------------------------------------------------------------
class MSD(object):
    """
    The MSD class represents a Mass-Spring-Damper system.
    """
    # System parameters
    m = 30.48

    def __init__(self, name, **kwargs):
        """
        Initialise the msd object.

        :param: name  = system name
        """
        self.name = name

        # Pass through any other keyword arguments
        for key in kwargs:
            self.__dict__[key] = kwargs[key]

        self.c_idx = [ 'k', 'b', 'd' ]

        # Model coefficients
        self.C = { 'k': -50.0, 'b': -10.0, 'd': 1.0, 'z': 0.0 }

        # External force function
        self.d_func = None

        # State noise function
        self.w_func = None

        self.init()

    def __str__(self):
        return self.name

    def init(self):
        """
        Construct the force and moment matrices.
        """
        # Rigid body mass matrix
        self.C_M_I = 1.0/self.m

        self.C_SD = ml.zeros((2, 2))
        self.M_EF = ml.zeros((2, 1))

    def get_coeffs(self):
        """
        Get the model coefficients.
        """
        return self.C

    def set_coeffs(self, C):
        """
        Set the model coefficients.
        """
        for ck in C.keys():
            self.C[ck] = C[ck]

    def set_external_forces(self, T_S, D_S, interp_kind):
        """
        Set the external force interpolant points.
        """
        self.d_func = interpolate.interp1d(T_S, D_S, kind=interp_kind, axis=0, bounds_error=False)

    def add_state_noise(self, T_S, W_S):
        """
        Set the state noise interpolant points.
        """
        self.w_func = interpolate.interp1d(T_S, W_S, kind='linear', axis=0, bounds_error=False)

    def rates(self, x, t):
        """
        Calculate the system state-rate for the current state x.

        :param: x = current system state [ xp, xpd ]
        :param: t = current time

        :returns: xdot = system state-rate
        """
        # Spring-damper forces
        # C_SD = np.mat([[    0.0,      self.m     ],
        #                [ self.C['k'], self.C['b']]])
        self.C_SD[0, 0] = 0.0
        self.C_SD[0, 1] = self.m
        self.C_SD[1, 0] = self.C['k']
        self.C_SD[1, 1] = self.C['b']

        M_SD = self.C_SD*x.reshape((-1, 1))

        d = np.nan_to_num(self.d_func(t))

        # External force
        # M_EF = np.mat([[           0.0 ],
        #                [ self.C['d']*f ]])
        self.M_EF[0, 0] = 0.0
        self.M_EF[1, 0] = self.C['d']*d

        xdot = np.ravel(self.C_M_I*(M_SD + self.M_EF))

        if (self.w_func is not None):
            xdot[1] += np.nan_to_num(self.w_func(t))

        return xdot

    def rrates(self, t, x):
        """
        Rates method with arguments reversed.
        """
        return self.rates(x, t, self.d_func)

    def forces(self, xdot, x):
        """
        Calculate the forces from recorded state data.

        :param: xdot = system state rate
        :param: x = system state [ xp, xpd ]

        :returns: f = state forces
        """
        xpddot = xdot[1]

        f = self.m*xpddot

        return f

    def integrate(self, x0, T):
        """
        Integrate the differential equations and calculate the resulting rates and forces.

        :param: x0 = initial system state
        :param: T = sequence of time points for which to solve for x

        :returns: X = system state array
        :returns: Xdot = state rates array
        :returns: F = state force array
        """
        N = T.shape[0]

        # Initialise the model
        self.init()

        # Perform the integration
        X = integrate.odeint(self.rates, x0, T, rtol=1.0e-6, atol=1.0e-6)

        Xdot = np.zeros((N, len(x0)))
        for n in range(N):
            Xdot[n] = self.rates(X[n], T[n])

        # Force and moment matrix
        F = np.zeros((N, 1))
        for n in range(N):
            F[n] = self.forces(Xdot[n], X[n])

        return X, Xdot, F

But we are going to use the Cython implementation here, as the Python model is very slow!. Simulate and add (Gaussian) noise to the output.

In [6]:
# Simulation model
MODEL = 'pyublas' # ['python', 'cython', 'pyublas', 'numba', 'boost']

if (MODEL == 'python'):
    # Pure Python
    from msd import MSD
elif (MODEL == 'cython'):
    # Cython
    from msd.msdc import MSD_CYTHON
elif (MODEL == 'pyublas'):
    # PyUblas extension
    from msd.msdu import MSD_PYUBLAS
elif (MODEL == 'numba'):
    # Numba JIT
    from msd.msdn import MSD_NUMBA
elif (MODEL == 'numba_jc'):
    # Numba JIT
    from msd.msdn import MSD_NUMBA_JC, msd_integrate
elif (MODEL == 'boost'):
    # Boost extension
    from msd.msdb import MSD_BOOST
In [7]:
# Measurement and state noise standard deviations
NOISE_SD = [ 0.001 for _ in range(3) ]
STATE_NOISE_SD = 0.0

VERBOSE = False

# Zero the RNG seed
np.random.seed(1)

# Initial system state and external force input
x0 = np.zeros((2, ))
d0 = 0.0

# Sample period
dt = 0.01

# Start and end time
t0 = 0.0
tN = 15.0

# Create the time vector
T = np.arange(t0, tN, dt)
N = T.shape[0]

# Create the predefined external force vector
T_S0 = np.hstack((t0, np.arange(t0 + 1.0, tN + 1.0, 1.0), tN))
D_S0 = np.hstack((d0, np.array([ d0 + ((_ % 2)*2 - 1) * 1.0 for _ in range(T_S0.shape[0] - 2) ]), d0))
interpfun = interpolate.interp1d(T_S0, D_S0, kind='zero', axis=0, bounds_error=False)
D0 = np.array([ [ interpfun(t) ] for t in T ])
T_S = T_S0.copy()
D_S = D_S0.copy()
D = D0.copy()

# Create the simulation model
if (MODEL == 'python'):
    # Pure Python
    msd = MSD("Mass-Spring-Damper (Python)")
    msd.set_external_forces(T_S, D_S, 'zero')
elif (MODEL == 'cython'):
    # Cython
    msd = MSD_CYTHON("Mass-Spring-Damper (Cython)")
    msd.set_external_forces(T_S, D_S, 'zero')
elif (MODEL == 'pyublas'):
    # PyUblas extension
    msd = MSD_PYUBLAS("Mass-Spring-Damper (PyUblas)", N)
    msd.set_external_forces(T_S, D_S, 'zero')
elif (MODEL == 'numba'):
    # Numba JIT
    msd = MSD_NUMBA("Mass-Spring-Damper (Numba)", N)
    msd.set_external_forces(T_S, D_S, 'zero')
elif (MODEL == 'numba_jc'):
    # Numba JIT
    msd = MSD_NUMBA_JC(N)
    msd.set_external_forces(T_S, D_S, 0)
elif (MODEL == 'boost'):
    # Boost extension
    msd = MSD_BOOST("Mass-Spring-Damper (Boost)", N)
    msd.set_external_forces(T_S, D_S, 'zero')

# Identification keys
c_idx = ['k', 'b', 'd']

# True parameter set
CT = [ msd.get_coeffs()[ck] for ck in c_idx ]

# Initial parameter set
C0 = [ 0.5*msd.get_coeffs()[ck] for ck in c_idx ]

# Add any state noise
if (STATE_NOISE_SD > 0.0):
    sdw = STATE_NOISE_SD
    W = np.random.randn(N, 1)*sdw
    msd.add_state_noise(T, W)

# Compute the response
if (MODEL in ['python', 'cython', 'pyublas', 'numba', 'boost']):
    X, Xdot, F = msd.integrate(x0, T)
elif (MODEL == 'numba_jc'):
    X, Xdot, F = msd_integrate(msd, x0, T)

# State noise standard deviation vector
sdz = np.zeros((len(x0),))
if any(w > 0.0 for w in NOISE_SD[:2]):
    sdz = np.array(NOISE_SD[:2])

# Measured state matrix
Z = X + np.random.randn(N, len(x0))*sdz

# Set the initial measured state equal to the initial true state
z0 = x0

Nu = Z[:,1]
sdnu = sdz[1]

# State rate noise standard deviation vector
sdzdot = np.zeros((len(x0),))
if any(w > 0.0 for w in NOISE_SD[:2]):
    sdzdot = np.array(NOISE_SD[:2])

# Measured state rate matrix
Zdot = Xdot + np.random.randn(N, len(x0))*sdzdot

# External force noise standard deviation vector
sde = 0.0
if (NOISE_SD[2] > 0.0):
    sde = NOISE_SD[2]

# Measured external force matrix
E = D + np.random.randn(N, 1)*sde

# Set the initial measured external force equal to the initial true external force
e0 = d0

# Compute the inertial force and noise standard deviation
G = F.copy()
sdg = 0.0
if any(w > 0.0 for w in NOISE_SD):
    # Forces are calculated from (measured) accelerations, not measured directly
    for n in range(N):
        G[n] = msd.forces(Zdot[n], Z[n])
    sdg = np.std(F - G)

Plot the measured response and control input.

In [8]:
import matplotlib.pyplot as pp
In [9]:
x_str = [ r'$x (m)$', r'$\dot x (m/s)$' ]
x_fac = [ 1.0, 1.0 ]
xd_str = [ r'$\dot x (m/s)$', r'$\ddot x (m/s^2)$' ]
xd_fac = [ 1.0, 1.0 ]
d_str = [ r'$\delta_F (N)$' ]
d_fac = [ 1.0 ]
f_str = [ r'$X (N)$' ]
f_fac = [ 1.0 ]
ff_str = [ r'ln ${\Delta X}^2 (N^2)$' ]

nc = 2 # number of columns
nr = 3 # number of rows

fig, AxesArr = pp.subplots(nr, nc, figsize=(10.0, 6.0))

Axes = np.ravel(AxesArr)

ax = Axes[0]
ax.grid(color='lightgrey', linestyle=':')
ax.plot(T, G*f_fac[0], color='#DFD4F4', linestyle='-', linewidth=1.5)
ax.set_xlim(T[0], T[N - 1])
ax.autoscale(enable=False)
ax.set_ylabel(f_str[0])

for j in range(2):
    ax = Axes[j + 1]
    ax.grid(color='lightgrey', linestyle=':')
    ax.plot(T, Z[:,j]*x_fac[j], color='#BCE8E6', linestyle='-', linewidth=1.5)
    ax.set_xlim(T[0], T[N - 1])
    ax.autoscale(enable=False)
    ax.set_ylabel(x_str[j])

ax = Axes[3]
ax.grid(color='lightgrey', linestyle=':')
ax.plot(T, E*d_fac[0], color='#BDDCBD', linestyle='-', linewidth=1.5)
ax.set_xlim(T[0], T[N - 1])
ax.autoscale(enable=False)
ax.set_ylabel(d_str[0])

for j in range(2):
    ax = Axes[j + 4]
    ax.grid(color='lightgrey', linestyle=':')
    ax.plot(T, Zdot[:,j]*xd_fac[j], color='#BCE8E6', linestyle='-', linewidth=1.5)
    ax.set_xlim(T[0], T[N - 1])
    ax.autoscale(enable=False)
    ax.set_ylabel(xd_str[j])

pp.subplots_adjust(left=0.1, wspace=0.3)

Least-Squares Estimation

Perform a Least-Squares estimation on the measured system.

In [10]:
from scipy import linalg, stats
In [11]:
FF = ml.repmat(None, 50, 1)

print("LINEAR REGRESSION:")

A = np.c_[Z[:,0],Zdot[:,0],E[:,0]]

# Least squares solution
C, resid, rank, sigma = linalg.lstsq(A, G)

C_LS = C.flatten().tolist()

sigma2 = np.sum((G - np.dot(A, C))**2.0)/(N - len(c_idx)) # RMSE

cov = sigma2*np.linalg.inv(np.dot(A.T, A)) # covariance matrix
se = np.sqrt(np.diag(cov)) # standard error

alpha = 0.05
conf = 100.0*(1.0 - alpha) # confidence level

sT = stats.distributions.t.ppf(1.0 - alpha/2.0, N - len(c_idx)) # student T multiplier
CI = sT*se

print("            TRUE       EST    {:2.0f}% CONF".format(conf))
for i in range(len(c_idx)):
    ck = c_idx[i]
    print("{:5s}: {:10.4f} {:10.4f} +/-{:-.4f}".format(c_idx[i], msd.get_coeffs()[ck], C_LS[i], CI[i]))

SS_tot = np.sum((G - np.mean(G))**2.0)
SS_err = np.sum((np.dot(A, C) - G)**2)

#  http://en.wikipedia.org/wiki/Coefficient_of_determination
Rsq = 1.0 - SS_err/SS_tot

print("R^2 = {:.4f}".format(Rsq))

# Estimated force matrix
H = np.dot(A, C_LS)
LINEAR REGRESSION:
            TRUE       EST    95% CONF
k    :   -50.0000   -49.8210 +/-1.9412
b    :   -10.0000    -7.3768 +/-0.8025
d    :     1.0000     0.9765 +/-0.0120
R^2 = 0.9708
In [12]:
# Create the simulation model
if (MODEL == 'python'):
    # Pure Python
    msd_est = MSD("Mass-Spring-Damper_REG_EST")
    msd_est.set_external_forces(T, E, 'linear_unifom')
elif (MODEL == 'cython'):
    # Cython
    msd_est = MSD_CYTHON("Mass-Spring-Damper_REG_EST (Cython)")
    msd_est.set_external_forces(T, E, 'linear_uniform')
elif (MODEL == 'pyublas'):
    # PyUblas extension
    msd_est = MSD_PYUBLAS("Mass-Spring-Damper_REG_EST (PyUblas)", N)
    msd_est.set_external_forces(T, E, 'linear_unifom')
elif (MODEL == 'numba'):
    # Numba JIT
    msd_est = MSD_NUMBA("Mass-Spring-Damper_REG_EST (Numba)", N)
    msd_est.set_external_forces(T, E, 'linear_unifom')
elif (MODEL == 'numba_jc'):
    # Numba JIT
    msd_est = MSD_NUMBA_JC(N)
    msd_est.set_external_forces(T, E, 1)
elif (MODEL == 'boost'):
    # Boost extension
    msd_est = MSD_BOOST("Mass-Spring-Damper_REG_EST (Boost)", N)
    msd_est.set_external_forces(T, E, 'linear_uniform')

for i in range(len(c_idx)):
    msd_est.set_coeffs({ 'k': C_LS[0], 'b': C_LS[1], 'd': C_LS[2] })

# Compute the response
if (MODEL in ['python', 'cython', 'pyublas', 'numba', 'boost']):
    Xe, Xedot, Fe = msd_est.integrate(z0, T)
elif (MODEL == 'numba_jc'):
    Xe, Xedot, Fe = msd_integrate(msd_est, z0, T)

That's pretty good. Let's see how the estimated system response compares.

In [13]:
from plot import plot, addplot
In [14]:
fig, Axes, Lines, Text = plot(msd_est.name, T, E, Z, G, Xe=Xe, Fe=Fe)

Nice. Ok, now let's increase the amount of noise in the state & state rate by a factor of 5.

In [15]:
# Observation and state noise standard deviations
NOISE_SD = [ 0.005 for _ in range(3) ]
In [16]:
PLOT_SIM = False
%run -i sim.py
In [17]:
FF = ml.repmat(None, 50, 1)
fig, Axes, Lines, Text = plot(msd.name, T, E, Z, G, Xe=None, Fe=None, FF=None)

And perform a Least-Squares estimation again. The script below does exactly the same as above.

In [18]:
%run -i reg.py
LINEAR REGRESSION:
            TRUE       EST    95% CONF
k    :   -50.0000   -25.1409 +/-2.1565
b    :   -10.0000    -7.2012 +/-1.1670
d    :     1.0000     1.0512 +/-0.0173
R^2 = 0.9300

Whoa! That's terrible! How does the estimated system response look?

In [19]:
fig, Axes, Lines, Text = plot(msd_est.name, T, E, Z, G, Xe=Xe, Fe=Fe)

That is not good! This is primarily due to the fact that the error in the independant variables is introducing bias into the LS estimate.

Let's see how the same model compares with a different response dataset.

In [20]:
C_LS0 = C_LS
In [21]:
# Create the predefined external force vector
T_S1 = np.hstack((t0, np.arange(t0 + 1.5, tN + 1.5, 1.5), tN))
D_S1 = np.hstack((d0, np.array([ d0 + np.random.randint(-2, 3)/2.0*(1.0 - np.abs(d0)) for _ in range(T_S1.shape[0] - 2) ]), d0))
interpfun = interpolate.interp1d(T_S1, D_S1, kind='zero', axis=0, bounds_error=False)
D1 = np.array([ [ interpfun(t) ] for t in T ])
T_S = T_S1.copy()
D_S = D_S1.copy()
D = D1.copy()
In [22]:
%run -i sim.py

Now we'll use the new external force input to stimulate our model.

In [23]:
msd_est.set_external_forces(T, E, 'linear_uniform')

for i in range(len(c_idx)):
    msd_est.set_coeffs({ 'k': C_LS0[0], 'b': C_LS0[1], 'd': C_LS0[2] })

# Compute the response
Xe, Xedot, Fe = msd_est.integrate(z0, T)

And plot the system response.

In [24]:
fig, Axes, Lines, Text = plot(msd.name, T, E, Z, G, Xe=Xe, Fe=Fe)

This really highlights the deficiency of the equation error (least squares) approach.

Nonlinear (iterative) Gradient-based Estimation (Powell's Method)

Let's try that again, this time using an iterative technique. We'll use Powell's Method here.

Going back to our original external force vector:

In [25]:
# Create the predefined external force vector
T_S = T_S0.copy()
D_S = D_S0.copy()
D = D0.copy()
In [26]:
import matplotlib.pyplot as pp
In [27]:
%run -i sim.py

Create the objective function and minimise.

In [28]:
import math
import time

from scipy import optimize
from plot import updateplot
In [29]:
FF = ml.repmat(None, 50, 1)

# Create the simulation model
if (MODEL == 'python'):
    # Pure Python
    msd_fest = MSD("Mass-Spring-Damper_FMIN_EST")
    msd_fest.set_external_forces(T, E, 'linear_unifom')
elif (MODEL == 'cython'):
    # Cython
    msd_fest = MSD_CYTHON("Mass-Spring-Damper_FMIN_EST (Cython)")
    msd_fest.set_external_forces(T, E, 'linear_uniform')
elif (MODEL == 'pyublas'):
    # PyUblas extension
    msd_fest = MSD_PYUBLAS("Mass-Spring-Damper_FMIN_EST (PyUblas)", N)
    msd_fest.set_external_forces(T, E, 'linear_unifom')
elif (MODEL == 'numba'):
    # Numba JIT
    msd_fest = MSD_NUMBA("Mass-Spring-Damper_FMIN_EST (Numba)", N)
    msd_fest.set_external_forces(T, E, 'linear_unifom')
elif (MODEL == 'numba_jc'):
    # Numba JIT
    msd_fest = MSD_NUMBA_JC(N)
    msd_fest.set_external_forces(T, E, 1)
elif (MODEL == 'boost'):
    # Boost extension
    msd_fest = MSD_BOOST("Mass-Spring-Damper_FMIN_EST (Boost)", N)
    msd_fest.set_external_forces(T, E, 'linear_uniform')

fig, Axes, Lines, Text = plot(msd_fest.name, T, E, Z, G, Xe=np.zeros(X.shape), Fe=np.zeros(F.shape), FF=FF)
fig.canvas.draw()

kws = { 'fig': fig, 'Axes': Axes, 'Lines': Lines, 'Text': Text }

c_idx = ['k', 'b', 'd']

print("POWELL'S MINIMIZATION:")

class Objfun(object):

    def __init__(self, z0, T, G, FF):
        self.z0 = z0
        self.T = T
        self.G = G
        self.FF = FF
        self.fopt_max = None
        self.it = 0;

    def __call__(self, C, fig, Axes, Lines, Text):
        msd_fest.set_coeffs({ 'k': C[0], 'b': C[1], 'd': C[2] })

        # Compute the response
        if (MODEL in ['python', 'cython', 'pyublas', 'numba', 'boost']):
            Xe, Xedot, Fe = msd_fest.integrate(z0, T)
        elif (MODEL == 'numba_jc'):
            Xe, Xedot, Fe = msd_integrate(msd_fest, z0, T)

        # For fmin, fmin_powell, fmin_bfgs, fmin_l_bfgs_b
        dF = G - Fe
        fopt_sum = np.sum(dF*dF)

        if (self.it < np.size(FF, 0)):
            self.FF[self.it, 0] = math.log(fopt_sum)
        else:
            self.FF = np.roll(self.FF, -1)
            self.FF[-1, 0] = math.log(fopt_sum)

        f_max = None
        if ((self.fopt_max is None) or (self.fopt_max < math.log(fopt_sum))):
            f_max = math.log(fopt_sum) * 1.1
            self.fopt_max = math.log(fopt_sum)
            # rescale = True
        f_txt = '{:.4f}'.format(fopt_sum)

        updateplot(fig, Axes, Lines, Text, Xe, Fe, self.FF, f_max=f_max, f_txt=f_txt, c_txt=C)

        self.it += 1

        return fopt_sum

tic = time.clock()

objfun = Objfun(z0, T, G, FF)

# Need to start with a nontrivial parameter set to avoid getting stuck in a local minima straight away...
C = optimize.fmin_powell(objfun, C0, args=( fig, Axes, Lines, Text ), maxiter=100)

toc = time.clock() - tic
print("Time elapsed: {:f} seconds".format(toc))

C_PM = C.tolist()

print()
print("            TRUE      F_EST")
for i in range(len(c_idx)):
    ck = c_idx[i]
    print("{:5s}: {:10.4f} {:10.4f}".format(ck, msd.get_coeffs()[ck], C_PM[i]))

msd_fest.set_coeffs({ 'k': C_PM[0], 'b': C_PM[1], 'd': C_PM[2] })

# Compute the response
if (MODEL in ['python', 'cython', 'pyublas', 'numba', 'boost']):
    Xe, Xedot, Fe = msd_fest.integrate(z0, T)
elif (MODEL == 'numba_jc'):
    Xe, Xedot, Fe = msd_integrate(msd_fest, z0, T)
POWELL'S MINIMIZATION:
Optimization terminated successfully.
         Current function value: 87.164648
         Iterations: 6
         Function evaluations: 434
Time elapsed: 73.060531 seconds

            TRUE      F_EST
k    :   -50.0000   -50.2010
b    :   -10.0000    -8.9262
d    :     1.0000     0.9761

Much better! Although it did take some time. And the system response is quite good too. Let's check again using a different input and response dataset.

In [30]:
C_PM0 = C_PM
In [31]:
T_S = T_S1.copy()
D_S = D_S1.copy()
D = D1.copy()
In [32]:
%run -i sim.py
In [33]:
msd_est.set_external_forces(T, E, 'linear_uniform')

for i in range(len(c_idx)):
    msd_est.set_coeffs({ 'k': C_PM0[0], 'b': C_PM0[1], 'd': C_PM0[2] })

# Compute the response
if (MODEL in ['python', 'cython', 'pyublas', 'numba', 'boost']):
    Xe, Xedot, Fe = msd_est.integrate(z0, T)
elif (MODEL == 'numba_jc'):
    Xe, Xedot, Fe = msd_integrate(msd_est, z0, T)
In [34]:
fig, Axes, Lines, Text = plot(msd_fest.name, T, E, Z, G, Xe=Xe, Fe=Fe)

Great. Other methods (from the SciPy Optimisation Library) that work just as well include BFGS, and nonlinear Least-Squares.

Let's see how a Levenberg-Marquardt optimisation performs. Using our original external force vector.

Nonlinear (iterative) Gradient-based Estimation (Levenberg-Marquardt)

In [35]:
# Create the predefined external force vector
T_S = T_S0.copy()
D_S = D_S0.copy()
D = D0.copy()
In [36]:
%run -i sim.py
In [37]:
import math
import time

import lmfit as lm
from plot import updateplot
In [38]:
FF = ml.repmat(None, 50, 1)

# Create the simulation model
if (MODEL == 'python'):
    # Pure Python
    msd_fest = MSD("Mass-Spring-Damper_FMIN_EST")
    msd_fest.set_external_forces(T, E, 'linear_unifom')
elif (MODEL == 'cython'):
    # Cython
    msd_fest = MSD_CYTHON("Mass-Spring-Damper_FMIN_EST (Cython)")
    msd_fest.set_external_forces(T, E, 'linear_uniform')
elif (MODEL == 'pyublas'):
    # PyUblas extension
    msd_fest = MSD_PYUBLAS("Mass-Spring-Damper_FMIN_EST (PyUblas)", N)
    msd_fest.set_external_forces(T, E, 'linear_unifom')
elif (MODEL == 'numba'):
    # Numba JIT
    msd_fest = MSD_NUMBA("Mass-Spring-Damper_FMIN_EST (Numba)", N)
    msd_fest.set_external_forces(T, E, 'linear_unifom')
elif (MODEL == 'numba_jc'):
    # Numba JIT
    msd_fest = MSD_NUMBA_JC(N)
    msd_fest.set_external_forces(T, E, 1)
elif (MODEL == 'boost'):
    # Boost extension
    msd_fest = MSD_BOOST("Mass-Spring-Damper_FMIN_EST (Boost)", N)
    msd_fest.set_external_forces(T, E, 'linear_uniform')

fig, Axes, Lines, Text = plot(msd_fest.name, T, E, Z, G, Xe=np.zeros(X.shape), Fe=np.zeros(F.shape), FF=FF)
fig.canvas.draw()

kws = { 'fig': fig, 'Axes': Axes, 'Lines': Lines, 'Text': Text }

c_idx = ['k', 'b', 'd']

print("LEVENBERG-MARQUARDT OPTIMIZATION:")

class Fcn2min(object):

    def __init__(self, z0, T, G, FF):
        self.z0 = z0
        self.T = T
        self.G = G
        self.FF = FF
        self.fopt_max = None
        self.it = 0;

    def __call__(self, P, **kws):
        C = [ P[c_idx[i]].value for i in range(len(c_idx)) ]
        msd_fest.set_coeffs({ 'k': C[0], 'b': C[1], 'd': C[2] })

        # Compute the response
        if (MODEL in ['python', 'cython', 'pyublas', 'numba', 'boost']):
            Xe, Xedot, Fe = msd_fest.integrate(self.z0, self.T)
        elif (MODEL == 'numba_jc'):
            Xe, Xedot, Fe = msd_integrate(msd_fest, self.z0, self.T)

        fopt = np.ravel(self.G - Fe)
        fopt_sum = np.sum(fopt*fopt)

        if (self.it < np.size(FF, 0)):
            self.FF[self.it, 0] = fopt_sum
        else:
            self.FF = np.roll(self.FF, -1)
            self.FF[-1, 0] = fopt_sum

        f_max = None
        if ((self.fopt_max is None) or (self.fopt_max < fopt_sum)):
            f_max = fopt_sum * 1.1
            self.fopt_max = fopt_sum
        f_txt = '{:.4f}'.format(fopt_sum)

        updateplot(kws['fig'], kws['Axes'], kws['Lines'], kws['Text'], Xe, Fe, self.FF, f_max=f_max, f_txt=f_txt, c_txt=C)

        self.it += 1

        return fopt

# Create a set of Parameters
P = lm.Parameters()
for i in range(len(c_idx)):
    ck = c_idx[i]
    P.add(ck, value=C0[i])

tic = time.clock()

fcn2min = Fcn2min(z0, T, G, FF)

# Do fit, here with leastsq model
res = lm.minimize(fcn2min, P, kws=kws, method='leastsq', epsfcn=0.1)

toc = time.clock() - tic
print("Time elapsed: {:f} seconds".format(toc))

C_LM = [ res.params[c_idx[i]].value for i in range(len(c_idx)) ]

print()
print("            TRUE      F_EST")
for i in range(len(c_idx)):
    ck = c_idx[i]
    print("{:5s}: {:10.4f} {:10.4f}".format(ck, msd.get_coeffs()[ck], C_LM[i]))

msd_fest.set_coeffs({ 'k': C_LM[0], 'b': C_LM[1], 'd': C_LM[2] })

# Compute the response
if (MODEL in ['python', 'cython', 'pyublas', 'numba', 'boost']):
    Xe, Xedot, Fe = msd_fest.integrate(z0, T)
elif (MODEL == 'numba_jc'):
    Xe, Xedot, Fe = msd_integrate(msd_fest, z0, T)