python - 使用 GEKKO 的 ODE 系统中的参数估计

Question

我有一个非线性常微分方程系统，它有 3 个参数作为输入：mu、target 和 omega。我正在尝试通过最小化我的观察值和 ODE 系统的输出 (x5) 之间的最小二乘误差来估计这些参数。我尝试了不同的目标函数以及不同的 IMODES，但解决方案从未收敛。此外，由于我总是没有找到解决方案的错误，我看不到最接近的最佳参数是什么。你对我的优化问题有什么建议吗？谢谢！

这是我的yData

我的代码是：

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
import pickle 

with open('y_observed','rb') as f:
    yData = pickle.load(f)

m = GEKKO()
dt = 0.001
m.time = np.arange(0,1.5,dt)

# Parameters
mu = m.FV(value=0.2, lb=-0.1, ub=1)
omega = m.FV(value=70, lb=0, ub=100)
target = m.FV(value=0.1, lb=-1, ub=3)

# Variables
x1 = m.Var(value=0)
x2 = m.Var(value=0.05)
x3 = m.Var(value=0)
x4 = m.Var(value=0)
x5 = m.Var(value=0.05)
x6 = m.Var(value=0)
y = m.Param(value = yData)
m.Minimize((y-x5)**2)

# Constants 
A = 75
B = 70
C = 45

# ODEs
m.Equation(x1.dt()== 1-x1)
m.Equation(x2.dt()== x3)
m.Equation(x3.dt()== x1*70*(x1*B*0.25*(target-x2) - x3))
m.Equation(x4.dt()== C*(mu - x4))
m.Equation(x5.dt()== A*1.0/np.abs(mu)*(x4 - ((x5-x2)**2 + x6**2))*(x5-x2) - omega*x6)
m.Equation(x6.dt()== A*1.0/np.abs(mu)*(x4 - ((x5-x2)**2 + x6**2))*x6 + omega*(x5-x2))


#Application options
m.options.SOLVER = 3 #APOPT solver
m.options.IMODE = 5 #Dynamic Simultaneous - estimation
m.options.EV_TYPE = 2 #absolute error
#m.options.NODES = 3 #collocation nodes (2,5)

if True:
    mu.STATUS=1
    target.STATUS=1
    omega.STATUS=1

m.solve(disp=True)

print('Objective: ' + str((x5-yData)**2))

print('Solution')
print('mu = ' + str(mu.value[0]))
print('omega = ' + str(omega.value[0]))
print('target = ' + str(target.value[0]))

plt.figure(1)

plt.subplot(2,1,1)
plt.plot(m.time,x5.value,'k:',LineWidth=2,label=r'$x_5$')
plt.plot(m.time,yData,'b-',LineWidth=2,label=r'$yData$')
plt.ylabel('Value')
plt.legend(loc='best')

plt.subplot(2,1,2)
plt.plot(m.time,mu.value,'r-',LineWidth=2,label=r'$mu$')
plt.plot(m.time,target.value,'r-',LineWidth=2,label=r'$target$')
plt.plot(m.time,omega.value,'r-',LineWidth=2,label=r'$omega$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')

plt.show()

标签： pythonoptimizationodeleast-squaresgekko

Answer 1

对绝对值使用 Gekko 函数（m.abs2()或）而不是 Numpy函数。m.abs3()np.abs()

amu = m.abs3(mu)

这需要切换到 APPT 求解器m.options.SOLVER=1以确保满足模型的混合整数部分。此外，它通过将您的方程重新公式化为避免被零除：

m.Equation(amu*x5.dt()== A*(x4 - ((x5-x2)**2 + x6**2))*(x5-x2) - omega*x6)
m.Equation(amu*x6.dt()== A*(x4 - ((x5-x2)**2 + x6**2))*x6 + omega*(x5-x2))

完整的问题有很多变量：

 Number of state variables:          47971
 Number of total equations: -        44970
 Number of slack variables: -         5996
 ---------------------------------------
 Degrees of freedom       :          -2995

您可以忽略 ，Warning: DOF <= 0因为并非所有不等式约束都处于活动状态。这给出了一个可能是本地解决方案的成功解决方案。它可能需要更好的参数初始猜测来获得全局解决方案，例如如果您需要它来跟踪循环数据。

Objective: 260.68493685
Solution
mu = 0.0
omega = 70.0
target = 0.13906027337

这是完整的源代码，其中包含使用您的 PKL 文件检索数据的结果。

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
import pickle 

with open('y_observed.pkl','rb') as f:
    yData = pickle.load(f)

m = GEKKO()
dt = 0.001
m.time = np.arange(0,1.5,dt)

# Parameters
mu = m.FV(value=0.2, lb=-0.1, ub=1)
omega = m.FV(value=70, lb=0, ub=100)
target = m.FV(value=0.1, lb=-1, ub=3)

# Variables
x1 = m.Var(value=0)
x2 = m.Var(value=0.05)
x3 = m.Var(value=0)
x4 = m.Var(value=0)
x5 = m.Var(value=0.05)
x6 = m.Var(value=0)
y = m.Param(value = yData)
m.Minimize((y-x5)**2)

# Constants 
A = 75
B = 70
C = 45

# abs(mu)
amu = m.abs3(mu)

# ODEs
m.Equation(x1.dt()== 1-x1)
m.Equation(x2.dt()== x3)
m.Equation(x3.dt()== x1*70*(x1*B*0.25*(target-x2) - x3))
m.Equation(x4.dt()== C*(mu - x4))
m.Equation(amu*x5.dt()== A*(x4 - ((x5-x2)**2 + x6**2))*(x5-x2) - omega*x6)
m.Equation(amu*x6.dt()== A*(x4 - ((x5-x2)**2 + x6**2))*x6 + omega*(x5-x2))

#Application options
m.options.SOLVER = 1  # APOPT MINLP solver
m.options.IMODE = 5   # Dynamic Simultaneous - estimation
m.options.EV_TYPE = 2 # squared error
m.options.NODES = 3   #collocation nodes (2,5)

if True:
    mu.STATUS=1
    target.STATUS=1
    omega.STATUS=1

m.solve(disp=True)

print('Objective: ' + str(m.options.OBJFCNVAL))

print('Solution')
print('mu = ' + str(mu.value[0]))
print('omega = ' + str(omega.value[0]))
print('target = ' + str(target.value[0]))

plt.figure(1)

plt.subplot(2,1,1)
plt.plot(m.time,x5.value,'k:',lw=2,label=r'$x_5$')
plt.plot(m.time,yData,'b-',lw=2,label=r'$yData$')
plt.ylabel('Value')
plt.legend(loc='best')

plt.subplot(2,1,2)
plt.plot(m.time,mu.value,'b-',lw=2,label=r'$mu$')
plt.plot(m.time,target.value,'r--',lw=2,label=r'$target$')
plt.plot(m.time,omega.value,'g:',lw=2,label=r'$omega$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')

plt.show()

关于这个问题的一个复杂的事情是允许mu消极。如果给它一个较小的数字（例如0.01）的下限以避免解决方案，则更容易mu=0解决。

# Parameters
mu = m.FV(value=0.2, lb=0.01, ub=1)

Objective: 260.68477646
Solution
mu = 0.042062519923
omega = 69.992948651
target = 0.13871376807