Python Gekko：SpringMass系统的参数识别

-0.39是优化问题的局部最小值.由于最初的猜测与解决方案相距甚远，因此它会找到不同的本地解决方案.为了防止非物理解决方案，请添加约束，让求解器仅在界限内搜索.这可以在初始化期间通过以下方式完成:

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

或初始化后:

k.LOWER = 10
k.UPPER = 100

这是一个完整的脚本，无论最初的猜测如何，它都会返回正确的解决方案.

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}')

输出显示，无论最初的猜测如何，都找到了正确的解决方案.

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

如果没有约束，就会找到局部解，但目标函数较高，这表明它不是很好的匹配.

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

如果约束未知或存在多个局部解，则可以使用多起点方法来搜索全局最优值，如Design Optimization course上Global Optimization所示.以下是对初始猜测空间的全球搜索.hyperopt包使用Bayesian优化来寻找全局解决方案.

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']}")

解决方案如下:

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

虽然这个问题不需要，但有逻辑可以检测不良的初始猜测何时产生失败的解决方案并消除该初始猜测.