用Gekko优化循环起搏策略，收敛到局部不可行点

Question

我正在尝试制作一个自行车运动员的运动模型，该模型优化了完成12公里坡道计时试验所需的时间，并受到各种限制，例如总无氧能量作为最大功率。我试图重现Wolf等人的“通过个性化步调策略使公路自行车爬坡更快”中的结果。到目前为止，该模型由以下代码描述：Pow是自行车运动员的力量，tf是完成课程所需的时间，E_kin是自行车运动员的动能，E_an是无氧能量储备（最初为E_an_init）当骑车人输出的功率高于P_c时，它会减小，s是路线的梯度，mu是滚动阻力常数，eta是自行车的机械效率，c_P是空气动力阻力系数。我目前正在使用IMODE 6，但我也尝试了9。当前运行代码时，我得到一个错误，表示解决方案n找不到，我也不确定为什么会出现这种情况，我已经尝试检查不可行性文件，但我没有任何意义。我希望有一些关于我可能会出错的地方的指针。这是我第一次使用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, 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 = m.FV(value=1,lb = 0, name='tf')
tf.STATUS = 1
# control changes every time period
Pow = m.MV(value=0,lb=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得出的骑车人的速度不是m由动能方程导出的atch

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

Answer 1

你需要修改位置的最终等式。否则，除了最后一点之外，它是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,lb=0,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,1)
plt.plot(t,xpos.value,'r-',label='xpos')
plt.legend()
plt.subplot(4,1,2)
plt.plot(t,E_kin.value,'b-',label='E_kin')
plt.legend()
plt.subplot(4,1,3)
plt.plot(t,E_an.value,'k-',label='E_an')
plt.legend()
plt.subplot(4,1,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