Optimizing cycling pacing strategy with Gekko, Converged to a point of local infeasibility

Question

I am trying to produce a model of the motion of a cyclist which optimises the time taken to complete a 12km hill time trial, subject to various constraints such as total anaerobic energy as max power. I am attempting to reproduce the results in "Road Cycling Climbs Made Speedier by Personalized Pacing Strategies" by Wolf et al. The model so far is described by the code below: where Pow is the cyclist power, tf is the time taken to complete the course, E_kin is the kinetic energy of the cyclist, E_an is the anaerobic energy reserve (initially E_an_init) which decreases when the cyclist outputs a power higher than P_c, s is the gradient of the course and mu is the rolling resistance constant, eta is the mechanical efficiency of the bike and C_p is the aerodynamic drag coefficient. I am currently using IMODE 6 but I've also tried 9. When currently running the code I get an error saying the solution could not be found and I am not sure why this is occurring, I've tried checking the infeasibilities file but I can't make sense. I would appreciate some pointers as to where I might be going wrong. This is my first time using Gekko so I assume it may be something related to my implementation of the model.

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")

The results of the failed attempt look okay when plotted but it says that that the 12km are covered in 48.6 seconds which can't be correct when the velocity of the cyclist is around 20km/hr for most of the 12km. Seems like the results are a function of nt, which makes no sense. Moreover, the velocity of the cyclist derived from dx/dt doesn't match that derived from the kinetic energy equation.

I have used this solution as a reference so far: https://apmonitor.com/do/index.php/Main/MinimizeFinalTime

Answer 1

You need to reform the final equation for position. Otherwise, it is 0==12000 for all but the final point. The original form gives an infeasible solution.

m.Equation(xpos*final == x_f)

Here is the correct form of that equation.

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

In this form, it is 0==0 everywhere except at the end where it is 1*(xpos-12000)==0. It also helps to give an objective function term to guide the final condition with:

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

This guides the solver towards the feasible solution to meet the end point.

The inequality constraint

m.Equation(E_an*final >= 0)

isn't needed if E_an has a lower bound of zero with E_an = m.Var(lb=0) or by setting E_an.LOWER=0. It should always be positive throughout the entire time horizon, not just at the final solution. It is easier for the solver to have a constraint on the variable (e.g. lower bound) than to add an additional inequality constraint (e.g. with 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()

The solution doesn't look quite right for a competitive cyclist to cover 12 km. It should be much faster so there may be a problem with the current equations. The final position constraint is met and the final energy is zero but the time should be faster.

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