Gekko：Spring-Mass 系统上的参数识别

Question

我想在弹簧质量系统上进行参数估计

采用直接搭配法。参数 k 应根据响应确定。

我通过

模拟了这个系统

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]

然后我有这些代码来识别参数 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

通过调整k的初始值，我发现如果k的初始值为25到65，k将收敛到实际值50。否则结果将是-0.39，这不好。我很困惑，因为这个系统是线性的，应该很容易解决。所以我的安静：如何微调上面的代码，使 k 收敛到 50 并具有任意初始值？

Answer 1

-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

如果约束未知或存在多个局部解决方案，则可以使用多启动方法来搜索全局最优值，如全局优化上的设计优化课程所示。下面是对初始猜测空间的全局搜索。

hyperopt

包使用贝叶斯优化来寻找全局解决方案。

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

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