Search code examples
pythongekko

Gekko: parameter identification on a Spring-Mass system

I want to do parameter estimation on a Spring-Mass system

enter image description here

with direct collocation method. The parameter k should be determined from response.

I simulated this system by

from scipy.integrate import odeint
import numpy as np
def dy(y, t):
    x, xdot = y
    return [xdot, -50*x]

t = np.linspace(0, 1, 40)
sol = odeint(dy, [2.0, 1.0], t)
sol_x = sol[:, 0]
sol_xdot = sol[:, 1]

Then I have these code to identify parameter k:

from gekko import GEKKO

m = GEKKO(remote=False)
m.time = t
x = m.CV(value=sol_x); x.FSTATUS = 1  # fit to measurement
xdot = m.CV(value=sol_xdot); xdot.FSTATUS = 1
k = m.FV(value = 40.0); k.STATUS = 1  # change initial value of k here

m.Equation(x.dt() == xdot)            # differential equation
m.Equation(xdot.dt() == -k*x)

m.options.IMODE = 5   # dynamic estimation
m.options.NODES = 40   # collocation nodes
m.options.EV_TYPE = 2

m.solve(disp=False)   # display solver output

By playing around initial value of k, I found k will converge to real value 50 if its initial value is 25 to 65. Otherwise the result will be -0.39 which is not good. I'm quite confusing because this system is linear and should be easy to be solved. So my quietion: how to fine tune the above code so that k converge to 50 with arbitry initial value?

Solution

  • The -0.39 is a local minimum to the optimization problem. As the initial guess is further from the solution, it finds a different local solution. To prevent non-physical solutions, add a constraint for the solver to search only within bounds. This can be done during initialization with:

    k = m.FV(value=ki,lb=10,ub=100)

    or after initialization with:

    k.LOWER = 10
k.UPPER = 100

    Here is a complete script that returns the correct solution regardless of the initial guess.

    from gekko import GEKKO
from scipy.integrate import odeint
import numpy as np

# generate data
def dy(y, t):
    x, xdot = y
    return [xdot, -50*x]

t = np.linspace(0, 1, 40)
sol = odeint(dy, [2.0, 1.0], t)
sol_x = sol[:, 0]
sol_xdot = sol[:, 1]

# regression with range of initial guess values
kx = np.linspace(0,100,11)
for i,ki in enumerate(kx):
    m = GEKKO(remote=False)
    m.time = t
    x = m.CV(value=sol_x); x.FSTATUS = 1
    xdot = m.CV(value=sol_xdot); xdot.FSTATUS = 1
    k = m.FV(value=ki,lb=10,ub=100); k.STATUS = 1
    m.Equation(x.dt() == xdot)
    m.Equation(xdot.dt() == -k*x)

    m.options.IMODE = 5   # dynamic estimation
    m.options.NODES = 6   # collocation nodes (2-6)
    m.options.EV_TYPE = 2

    m.solve(disp=False)
    print(f'Initial Guess: {ki} ' +
          f'Solution: {k.value[0]} ' +
          f'ObjFcn: {m.options.OBJFCNVAL}')

    The output shows that the correct solution is found, regardless of the initial guess.

    Initial Guess: 0.0 Solution: 50.000000284 ObjFcn: 1.8231507179e-10
Initial Guess: 10.0 Solution: 50.000000284 ObjFcn: 1.8229910145e-10
Initial Guess: 20.0 Solution: 50.000000284 ObjFcn: 1.8231237768e-10
Initial Guess: 30.0 Solution: 50.000000284 ObjFcn: 1.8229796958e-10
Initial Guess: 40.0 Solution: 50.000000284 ObjFcn: 1.8229906845e-10
Initial Guess: 50.0 Solution: 50.000000284 ObjFcn: 1.8229906842e-10
Initial Guess: 60.0 Solution: 50.000000284 ObjFcn: 1.8229906933e-10
Initial Guess: 70.0 Solution: 50.000000284 ObjFcn: 1.8229906849e-10
Initial Guess: 80.0 Solution: 50.000000284 ObjFcn: 1.8230214413e-10
Initial Guess: 90.0 Solution: 50.000000284 ObjFcn: 1.8229908716e-10
Initial Guess: 100.0 Solution: 50.000000284 ObjFcn: 1.8230004454e-10

    Without the constraints, a local solution is found but the objective function is higher that indicates it is not a good fit.

    Initial Guess: 0.0 Solution: -0.39654383084 ObjFcn: 88263.987254
Initial Guess: 10.0 Solution: -0.39654383084 ObjFcn: 88263.987254
Initial Guess: 20.0 Solution: -0.39654383084 ObjFcn: 88263.987254
Initial Guess: 30.0 Solution: 50.000000284 ObjFcn: 1.8229906872e-10
Initial Guess: 40.0 Solution: 50.000000284 ObjFcn: 1.8229906866e-10
Initial Guess: 50.0 Solution: 50.000000284 ObjFcn: 1.8229906856e-10
Initial Guess: 60.0 Solution: 50.000000284 ObjFcn: 1.8229906861e-10
Initial Guess: 70.0 Solution: -0.39654383084 ObjFcn: 88263.987254
Initial Guess: 80.0 Solution: -0.39654383084 ObjFcn: 88263.987254
Initial Guess: 90.0 Solution: -0.39654383084 ObjFcn: 88263.987254
Initial Guess: 100.0 Solution: -0.39654383085 ObjFcn: 88263.987254

    If constraints are unknown or there are multiple local solutions then a multi-start method can be used to search for a global optimum as shown in the Design Optimization course on Global Optimization. Below is a Global search over an initial guess space. The hyperopt package uses Bayesian optimization to find the global solution.

    from gekko import GEKKO
from scipy.integrate import odeint
import numpy as np
from hyperopt import fmin, tpe, hp
from hyperopt import STATUS_OK, STATUS_FAIL

# generate data
def dy(y, t):
    x, xdot = y
    return [xdot, -50*x]

t = np.linspace(0, 1, 40)
sol = odeint(dy, [2.0, 1.0], t)
sol_x = sol[:, 0]
sol_xdot = sol[:, 1]

def objective(params):
    ki = params['kx']
    m = GEKKO(remote=False)
    m.time = t
    x = m.CV(value=sol_x); x.FSTATUS = 1
    xdot = m.CV(value=sol_xdot); xdot.FSTATUS = 1
    k = m.FV(value=ki); k.STATUS = 1
    m.Equation(x.dt() == xdot)
    m.Equation(xdot.dt() == -k*x)

    m.options.IMODE = 5   # dynamic estimation
    m.options.NODES = 6   # collocation nodes (2-6)
    m.options.EV_TYPE = 2

    m.solve(disp=False)
    obj = m.options.OBJFCNVAL
    if m.options.APPSTATUS==1:
        s=STATUS_OK
    else:
        s=STATUS_FAIL
    m.cleanup()    
    return {'loss':obj, 'status': s, 'k':k.value[0]}
          
# Define the search space for the hyperparameters
space = {'kx': hp.quniform('kx', 0, 100, 10),}

best = fmin(objective, space, algo=tpe.suggest, max_evals=10)
sol = objective(best)
print(f"Solution Status: {sol['status']}")
print(f"Objective: {sol['loss']:.2f}")
print(f"Initial Guess: {best['kx']}")
print(f"Solution: {sol['k']}")

    Here is the solution:

    Solution Status: ok
Objective: 0.00
Initial Guess: 30.0
Solution: 50.000000284

    While not needed for this problem, there is logic to detect when the poor initial guess produces a failed solution and eliminates that initial guess.