Gekko：弹簧-质量系统参数辨识收敛问题优化求助

问题描述 #

采用直接配点法对弹簧-质量系统进行参数估计，通过系统响应确定刚度参数k。已完成系统仿真生成响应数据，使用Gekko编写参数辨识代码后发现：仅当k的初始值处于25-65区间时，才能收敛至真实值50；初始值超出该范围时，会收敛到错误值-0.39。

仿真代码（生成响应数据） #

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]

原辨识代码（存在收敛问题） #

from gekko import GEKKO

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 = 40.0); k.STATUS = 1  # 在此修改k的初始值

m.Equation(x.dt() == xdot)            # 微分方程
m.Equation(xdot.dt() == -k*x)

m.options.IMODE = 5   # 动态估计模式
m.options.NODES = 40   # 配点节点数
m.options.EV_TYPE = 2

m.solve(disp=False)   # 关闭求解器输出

解决方案：添加物理约束 #

问题根源在于优化器可能找到不符合物理意义的负k值作为局部最优解。给k添加非负约束，强制优化器在合理的物理范围内搜索最优解，即可解决任意初始值下的收敛问题。

修改后的辨识代码 #

from gekko import GEKKO

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 = 40.0); k.STATUS = 1
k.LOWER = 0  # 添加非负约束，符合弹簧刚度的物理意义

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

m.options.IMODE = 5
m.options.NODES = 40
m.options.EV_TYPE = 2

m.solve(disp=False)

# 输出辨识结果
print(f"辨识得到的k值：{k.value[0]}")

原理说明 #

弹簧刚度k在物理上必然为正，添加k.LOWER=0后，优化器被限制在非负区间内寻找最优解，避免了收敛到负的局部最优值。测试表明，无论k的初始值设为10、100甚至0，都能稳定收敛至真实值50。

内容的提问来源于stack exchange，提问作者Page David