使用 Gekko 优化自行车步调策略，收敛到局部不可行点

问题描述

我正在尝试制作一个骑自行车者的运动模型，该模型优化了完成 12 公里山地计时赛所需的时间，受到各种限制，例如总无氧能量作为最大功率。我试图重现 Wolf 等人在“Road Cycling climbs Made Speedier by Personalized Pacing Strategies”中的结果。到目前为止的模型由以下代码描述：其中 Pow 是骑车人的功率，tf 是完成课程所需的时间，E_kin 是骑车人的动能，E_an 是无氧能量储备（最初是 E_an_init），当骑车人输出的功率高于 P_c，s 是赛道的坡度，mu 是滚动阻力常数，eta 是自行车的机械效率，C_p 是空气阻力系数。我目前正在使用 IMODE 6，但我也尝试过 9。当当前运行代码时，我收到一条错误消息，说找不到解决方案，我不确定为什么会发生这种情况，我已经尝试检查不可行性文件，但我没有意义。我会很感激一些关于我可能出错的地方的指示。这是我第一次使用 Gekko，所以我认为这可能与我的模型实现有关。

from gekko import GEKKO
import numpy as np

m = GEKKO(remote=False) # initialize gekko
nt = 501
m.time = np.linspace(0,1,nt)

#Parameters,taken from paper,idepdendant of the path
#mass = 80
g = 9.81
mu = 0.004
r_w = 0.335
I_s = 0.658
M = 80#mass + (I_s/r_w**2)
C_p = 0.7
rho = 1.2
A = 0.4
eta = 0.95
P_c = 150
P_max = 1000

#Path dependant
x_f = 12000
s = 0.081

#Define anaerobic energy reserve,dependant on the cyclist
E_an_init = 20000

#Variables
xpos = m.Var(value = 0,lb = 0,ub = x_f,name='xpos')
E_kin = m.Var(value = 0,lb=0,name='E_kin')#start at 1m/s
E_an =  m.Var(value = E_an_init,ub=E_an_init,name='E_an')

p = np.zeros(nt) # mark final time point
p[-1] = 1.0
final = m.Param(value=p)

# optimize final time
tf = m.FV(value=1,name='tf')
tf.STATUS = 1
# control changes every time period
Pow = m.MV(value=0,ub=P_max,name='Pow')
Pow.STATUS = 1

# Equations
m.Equation(xpos.dt()==((2*E_kin/M)**0.5)*tf)
m.Equation(E_kin.dt()==(eta*Pow-((2*E_kin/M)**0.5)*(M*g*(s+mu*m.cos(m.atan(s))) + C_P*rho*A*E_kin/M))*tf)
m.Equation(E_an.dt()==(P_c-Pow)*tf)
m.Equation(xpos*final == x_f)
m.Equation(E_an*final >= 0)

#m.options.MAX_ITER = 1000
m.Minimize(tf*final) # Objective function
m.options.IMODE = 6 # optimal control mode
m.open_folder() # open folder if remote=False to see infeasibilities.txt
m.solve(disp=True) # solve

print('Final Time: ' + str(tf.value[0]) + "s")

失败尝试的结果在绘制时看起来不错，但它表示 12 公里在 48.6 秒内完成，这在 12 公里的大部分时间里，当骑车人的速度约为 20 公里/小时时，这是不正确的。结果好像是nt的函数，没意义。此外，由 dx/dt 导出的骑自行车者的速度与由动能方程导出的速度不匹配。

到目前为止，我已将此解决方案用作参考： https://apmonitor.com/do/index.php/Main/MinimizeFinalTime

解决方法

您需要修改位置的最终方程。否则，除了最后一点之外的所有点都是 0==12000。原始形式给出了一个不可行的解决方案。

m.Equation(xpos*final == x_f)

这是该等式的正确形式。

m.Equation(final*(xpos-x_f)==0)

在这种形式中，除了最后是 0==0 外，其他地方都是 1*(xpos-12000)==0。它还有助于给出一个目标函数项来指导最终条件：

m.Minimize(final*(xpos-x_f)**2)

这会引导求解器找到可行的解决方案以达到终点。

不等式约束

m.Equation(E_an*final >= 0)

如果 E_an 的下限为 E_an = m.Var(lb=0) 或通过设置 E_an.LOWER=0，则不需要

。它应该在整个时间范围内始终为正，而不仅仅是在最终解决方案中。求解器对变量施加约束（例如下限）比添加额外的不等式约束（例如使用 m.Equation()）更容易。

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

m = GEKKO() # initialize gekko
nt = 101
m.time = np.linspace(0,1,nt)

#Parameters,taken from paper,idepdendant of the path
#mass = 80
g = 9.81
mu = 0.004
r_w = 0.335
I_s = 0.658
M = 80 #mass + (I_s/r_w**2)
C_p = 0.7
rho = 1.2
A = 0.4
eta = 0.95
P_c = 150
P_max = 1000

#Path dependant
x_f = 12000
s = 0.081

#Define anaerobic energy reserve,dependant on the cyclist
E_an_init = 20000

#Variables
xpos = m.Var(value = 0,lb = 0,ub = x_f,name='xpos')
E_kin = m.Var(value = 0,lb=1e-3,name='E_kin')#start at 1m/s
E_an =  m.Var(value = E_an_init,lb=0,ub=E_an_init,name='E_an')

p = np.zeros(nt) # mark final time point
p[-1] = 1.0
final = m.Param(value=p)

# optimize final time
tf_guess = 900
tf = m.FV(value=tf_guess,lb=300,ub=10000,name='tf')
tf.STATUS = 1
# control changes every time period
Pow = m.MV(value=100,ub=P_max,name='Pow')
Pow.STATUS = 1

# Equations
Esr = m.Intermediate((2*E_kin/M)**0.5)
m.Equation(xpos.dt()==Esr*tf)
m.Equation(E_kin.dt()==(eta*Pow-Esr*(M*g*(s+mu*m.cos(m.atan(s))) \
                          + C_p*rho*A*E_kin/M))*tf)
m.Equation(E_an.dt()==(P_c-Pow)*tf)

m.Minimize(final*(xpos-x_f)**2)
m.Equation(final*(xpos-x_f)==0)
m.Minimize(final*(tf**2))

m.options.IMODE = 6 # optimal control mode
m.solve(disp=True) # solve

print('Final Time: ' + str(tf.value[0]) + "s")
print('Final Position: ' + str(xpos.value[-1]))
print('Final Energy: ' + str(E_an.value[-1]))

t = m.time*tf.value[0]/60.0
plt.figure(figsize=(9,6))
plt.subplot(4,1)
plt.plot(t,xpos.value,'r-',label='xpos')
plt.legend()
plt.subplot(4,2)
plt.plot(t,E_kin.value,'b-',label='E_kin')
plt.legend()
plt.subplot(4,3)
plt.plot(t,E_an.value,'k-',label='E_an')
plt.legend()
plt.subplot(4,4)
plt.plot(t,Pow.value,'-',color='orange',label='Pow')
plt.legend()
plt.xlabel('Time (min)')
plt.show()

对于一个有竞争力的自行车手来说，完成 12 公里的解决方案看起来不太合适。它应该快得多，所以当前的方程可能有问题。满足最终位置约束，最终能量为零，但时间应该更快。

Final Time: 5550.6961268s
Final Position: 12000.0
Final Energy: 0.0