background-shape
feature-image
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('default')
import logging
logging.getLogger('matplotlib.font_manager').disabled = True
from scipy.signal import savgol_filter
import seaborn
import typing
!pip install pygad
import pygad
from skimage.transform import resize
import warnings
warnings.filterwarnings('ignore')

Profil optimal de puissance en fonction du dénivelé et du vent

Construire un modèle permettant de prédire une courbe de puissance optimale minimisant le temps de parcours sur un tracé GPS est un problème complexe qui peut être séparé en deux parties :

  1. Établir une relation puissance-vitesse (et inversement).
  2. Minimiser la puissance en fonction de paramètres dépendant du parcours étant donné un profil de puissance pour l’athlète.

On fera ici la démonstration d’un modèle le plus simple possible, des modèles plus complexe peuvent être trouvé dans la littérature. L’implémentation de ce modèle sera faite en Python.

1. Relation puissance-vitesse

Nous allons utiliser un modèle simplifié qui estime la puissance à fournir égale à la puissance nécessaire pour surmonter la somme des résistances que rencontre le cycliste. Un modèle plus complexe peut être dérivé de Martin, James C., et al. “Validation of a mathematical model for road cycling power.” Journal of applied biomechanics 14.3 (1998): 276-291. On peut bien sûr complexifier ce modèle en y ajoutant toutes sortes de résistances (rendement de la transmission du vélo, résistance des roulements etc…).

Dans ce modèle, on considèrera 3 résistances:

  • La force aérodynamique
  • La résistance au roulement
  • La gravité
  • L’efficience générale du vélo

Force aérodynamique \( F_a \)

C’est la force que le cycliste doit produire pour se mouvoir dans l’air et qui peut être estimée à l’aide de 4 paramètres :

  • velocity : vitesse du cycliste \(v\)
  • air_density : densité de l’air (\(\rho = ~1,292 kg.m^{-3}\))
  • cda : la trainée aérodynamique \(CdA\) (CLM : 0.25 ; route position aérodynamique : 0.30 ; route position détendue : 0.33)
  • wind_speed : la vitesse du vent \(v_{wind}\) (positive pour un vent de face)

$$ F_a = \frac{1}{2} .\rho .CdA. (v + v_{wind})^2 $$

def aero_force(velocity, air_density, cda, wind_speed):
  return 0.5*cda*air_density*(velocity + wind_speed)**2

Résistance au roulement \(F_r\)

La résistance au roulement peut être estimée à l’aide de 3 paramètres et 1 constante :

  • rolling_friction le coefficient de friction \(Crr\) (piste en bois : 0.0020 ; asphalte bon : 0.0040 ; asphalte usé : 0.0045)
  • masse poids total \(m\) (cycliste + vélo + équipement)
  • slope pente de la route \(\alpha\) exprimée en pourcentage
  • \(g\) l’intensité de la pesanteur égale à \(9.8067 m.s^{-2}\)

$$ F_r = Crr . m . g .\cos(\arctan(0.01\alpha)) $$

def rolling_resistance(rolling_friction, masse, slope):
  return rolling_friction*masse*9.8067*np.cos(np.arctan(slope*0.01))

Force de gravité \(F_g\)

La force de gravité peut être estimée à l’aide de 3 paramètres et 1 constante :

  • velocity : vitesse du cycliste \(v\)
  • masse poids total \(m\) (cycliste + vélo + équipement)
  • slope pente de la route α exprimée en pourcentage
  • g la constante de gravitation égale à 9.8067

$$ F_g = m . g .\sin(\arctan(0.01\alpha)) $$

def gravity(masse, slope):
  return masse*9.8067*np.sin(np.arctan(slope*0.01))

Puissance de propulsion

La seule force pour contrer l’effet de l’air, la friction sur la route et la gravité est la force de propulsion que produit le cycliste. Cette force doit être majorée un coefficient (2 à 3%) à cause de la perte de puissance liée à la dissipation d’énergie par la transmission, les roulements et autres (voir Jeukendrup, Asker E, and James Martin. 2001. “Improving Cycling Performance.” Sports Medicine 31 (7): 559–69.).
La puissance de propulsion \(P\) peut être estimée à partir de 8 paramètres :

  • velocity : vitesse du cycliste \(v\)
  • air_density : densité de l’air ρ (~1,292 kg. m -3)
  • cda : la trainée aérodynamique \(CdA\) (CLM : 0.25 ; route position aérodynamique : 0.30 ; route position détendue : 0.33)
  • wind_speed : la vitesse du vent \(v_{wind}\) (positive pour un vent de face)
  • rolling_friction la friction \(Crr\) (piste en bois : 0.0020 ; asphalte bon : 0.0040 ; asphalte usé : 0.0045)
  • masse poids total \(m\) (cycliste + vélo + équipement)
  • slope pente de la route α exprimé en pourcentage
  • efficiency l’efficience \(Eff\) (0.97 à 0.98)

$$ P = v . Eff . (F_a + F_r + F_g) $$

On peut inverser cette formule pour calculer la puissance étant donné une vitesse. On obtient alors un polynôme de degré 3 qui avec les contraintes de notre système a toujours une racine réelle et positive.

$$ \frac{1}{2}.Eff.\rho.CdA.v^3 + Eff.\rho.CdA.v_{wind}.v^2 + Eff.[\frac{1}{2}.\rho.CdA.v_{wind} + F_r + F_g]v - P = 0 $$

def propulsive_power(velocity, air_density, cda, wind_speed, rolling_friction, masse, slope, efficiency):
  return velocity*efficiency*(aero_force(velocity, air_density, cda, wind_speed) + rolling_resistance(rolling_friction, masse, slope) + gravity(masse, slope))


def velocity_from_power(power, air_density, cda, wind_speed, rolling_friction, masse, slope, efficiency):

  # Allows for arrays
  if not isinstance(power, typing.Iterable):
    power = [power]
  if not isinstance(slope, typing.Iterable):
    slope = [slope]
  if not isinstance(wind_speed, typing.Iterable):
    wind_speed = [wind_speed]
  
  if len(power) != len(wind_speed) != len(slope):
    return np.nan

  res = []
  for i, j, k  in zip(power, slope, wind_speed):
    order_3 = efficiency*0.5*air_density*cda
    order_2 = efficiency*air_density*cda*k
    order_1 = efficiency*(0.5*air_density*cda*k+rolling_resistance(rolling_friction, masse, j) + gravity(masse, j))
    order_0 = -i
    roots = np.roots([order_3, order_2, order_1, order_0])
    roots = roots[~np.iscomplex(roots) & (roots>0)]
    res.append(np.max(roots).real)
  return np.asarray(res)

Exemple

Pour tester notre modèle, on peut prendre l’exemple d’un cycliste de \(70kg\) sur un vélo de contre-la-montre de \(8kg\) roulant sur une piste en bois à \(30km.h^{-1}\) en l’absence de vent :

  • velocity = 8.3 m/s
  • air_density = 1,292 kg/m3
  • cda = 0.29
  • wind_speed = 0
  • rolling_friction = 0.0020
  • masse = 8kg
  • slope = 0
  • efficiency = 0.8
print("Ce cycliste devra produire {} W.".format(propulsive_power(velocity=8.3, air_density=1.292, cda=0.29, wind_speed=0, rolling_friction=0.002, masse=78, slope=0, efficiency=0.98)))
print("On peut test la formule inverse et y injectant la puissance précédemment calculée et voir qu'on retrouve bien {} m/s".format(velocity_from_power(propulsive_power(velocity=8.3, air_density=1.292, cda=0.29, wind_speed=0, rolling_friction=0.002, masse=78, slope=0, efficiency=0.98), air_density=1.292, cda=0.29, wind_speed=0, rolling_friction=0.002, masse=78, slope=0, efficiency=0.98)[0]))
Ce cycliste devra produire 117.41996590520002 W.
On peut test la formule inverse et y injectant la puissance précédemment calculée et voir qu'on retrouve bien 8.300000000000004 m/s

Puissance en fonction de la pente

On prend maintenant ce même cycliste voulant maintenir \(15 km.h^{-1}\) dans une pente à pourcentage croissant. On remarque que la puissance à produire est bien linéaire avec la pente de la route.

slope = np.linspace(0, 25, 50)
power = propulsive_power(velocity=4.15, air_density=1.292, cda=0.29, wind_speed=0, rolling_friction=0.002, masse=78, slope=slope, efficiency=0.98)
with plt.xkcd():
  plt.plot(slope, power, label="Pour une vitesse de 15 km/h")
  plt.xlabel("Pente en pourcents")
  plt.ylabel("Puissance en watts")
  plt.legend()

png

Puissance en fonction de la force du vent

On prend maintenant ce même cycliste voulant maintenir \(30 km.h^{-1}\) avec un vent variable. On remarque que la puissance nécessaire pour lutter avec le vent est bien quadratique.

wind = np.linspace(-20, 50, 100)
power = propulsive_power(velocity=8.30, air_density=1.292, cda=0.29, wind_speed=wind*(1e3/3600), rolling_friction=0.002, masse=78, slope=0, efficiency=0.98)
with plt.xkcd():
  plt.plot(wind, power, label="Pour une vitesse de 30 km/h")
  plt.xlabel("Dos       | Vent en km/h |      Face")
  plt.ylabel("Puissance en watts")
  plt.legend()

png

Influence du CdA

On peut maintenant comparer l’influence d’un vélo de contre-la-montre par rapport à un vélo classique. On remarque que comme prévu, le vélo de contre-la-monte est plus efficace que le vélo classique. Cette différence est plus marquée quand la vitesse relative par rapport au vent augmente.

velocity = np.linspace(0, 70, 100)
power_tt = propulsive_power(velocity=velocity*(1e3/3600), air_density=1.292, cda=0.29, wind_speed=0, rolling_friction=0.002, masse=78, slope=0, efficiency=0.98)
power_road = propulsive_power(velocity=velocity*(1e3/3600), air_density=1.292, cda=0.35, wind_speed=0, rolling_friction=0.002, masse=78, slope=0, efficiency=0.98)
power_tt_w = propulsive_power(velocity=velocity*(1e3/3600), air_density=1.292, cda=0.29, wind_speed=50*(1e3/3600), rolling_friction=0.002, masse=78, slope=0, efficiency=0.98)
power_road_w = propulsive_power(velocity=velocity*(1e3/3600), air_density=1.292, cda=0.35, wind_speed=50*(1e3/3600), rolling_friction=0.002, masse=78, slope=0, efficiency=0.98)
with plt.xkcd():
  plt.plot(velocity, power_tt, label="CLM, sans vent")
  plt.plot(velocity, power_road, label="Route, sans vent")
  plt.plot(velocity, power_tt_w, label="CLM, 50km/h vent de face")
  plt.plot(velocity, power_road_w, label="Route, 50km/h vent de face")
  plt.legend()
  plt.xlabel("Vitesse km/h")
  plt.ylabel("Puissance en watts")

png

2. Profil optimal de puissance

Trouver le profil optimal de puissance sur un parcours revient à minimiser le temps de parcours, ie. maximiser la vitesse moyenne, en minimisant la puissance requise. La solution optimale à ce problème dans le cas général où la pente et le vent varient sur le parcours peut s’avérer difficile à obtenir. Une solution du simplifié sans vent du problème en utilisant le principe de Pontryagin peut être trouvé dans De Jong, Jenny, et al. “The individual time trial as an optimal control problem.” Proceedings of the Institution of Mechanical Engineers, Part P: Journal of Sports Engineering and Technology 231.3 (2017): 200-206.

Une autre approche du problème est possible en utilisant un algorithme génétique pour se rapprocher d’une solution optimale par évolution.
L’algorithme génétique est une méthode permettant de résoudre des problèmes d’optimisation basée sur le principe de sélection naturelle. À chaque génération, l’algorithme modifie une population de solution (ici le profil de puissance) et utilise cette population pour produire d’autres solutions par mutations et recombinaisons. Au fil des générations, la population évolue vers une solution optimale. Pour cela, on doit définir une fonction de “fitness” qui va évaluer la chance de survie de la solution durant le processus.

Dans notre cas, cette fonction doit privilégier les solutions proches de la vitesse moyenne de référence à laquelle on veut parcourir le parcours. Elle doit aussi favoriser la solution qui dépense le moins d’énergie et est contrainte par une puissance minimale et maximale.

$$ F = max[1, \frac{1}{abs(\overline{v}-v_{target})}] + \frac{p_{eq}}{\overline{p}}$$ Avec \(\overline{v}\) la vitesse moyenne sur le parcours, \(v_{target}\) la vitesse cible, \(p_{eq}\) l’équivalent en puissance pour \(v_{target}\) et \(\overline{p}\) la puissance moyenne calculée comme la moyenne des puissances sur les intervalles de distance.

Une expérience rapide de pensé nous montre que pour une puissance équivalente de \(300W\) le terme \(p_{eq}/\overline{p}\) va varier entre \(300/P_{max}\) et \(300/P_{min}\) donc pour une borne supérieure de l’ordre de grandeur \(\approx1\) d’où le facteur numérique de saturation dans la fonction de fitness pour éviter que \(\frac{1}{abs(\overline{v}-v_{target})}\) ne diverge et que ce terme ne prenne trop d’importance dans le résultat final.

Pour tester notre modèle, on simulera un contre-la-montre de \(50km\) avec \(25km\) de plat puis \(25km\) de côte à \(8%\) et une vitesse cible de \(8m.s^{-1}\). Tous les autres paramètres du modèle resteront constants. On calculera 20 profils optimaux que l’on moyennera pour réduire les changements de puissance abrupts irréalistes.

slope_profile = np.zeros((500))
slope_profile[250::] = 8

slope_elev = np.zeros((500))+10
slope_elev[250::] = np.linspace(10,510,250)

p_max = 500
p_min = 200

def fitness(ga_instance, solution, index):
      sol = resize(solution, (500,1))

      # Temporal mean of velocity
      v = velocity_from_power(power=sol, air_density=1.292, cda=0.29, wind_speed=np.zeros((500)), rolling_friction=0.004, masse=85, slope=slope_profile, efficiency=0.98)
      t = 1/v
      v_mean = 500/np.sum(t)

      v_target = 8
      p_eq = propulsive_power(velocity=v_target, air_density=1.292, cda=0.29, wind_speed=np.zeros((500)), rolling_friction=0.004, masse=85, slope=slope_profile, efficiency=0.98)[0] 

      return np.clip(1/(np.abs(v_mean - v_target)), 0, 1) + p_eq/np.mean(solution)


results = []
for i in range(20):
  ga_instance = pygad.GA(num_generations=20,num_parents_mating=8,
                        fitness_func=fitness,
                        sol_per_pop=50,
                        num_genes=50,
                        init_range_low=190,
                        init_range_high=210,
                        gene_type=float,
                        parent_selection_type="rank",
                        keep_parents=0,
                        crossover_type="single_point",
                        mutation_type="adaptive",
                        keep_elitism=5,
                        mutation_probability=[.10, .05],
                        gene_space={'low': p_max, 'high': p_min},
                        stop_criteria=["saturate_5"]
                        )

  ga_instance.run()
  solution, solution_fitness, solution_idx = ga_instance.best_solution()

  solution = resize(solution, (500,1))
  results.append(solution)
with plt.xkcd():
  fig, ax1 = plt.subplots()
  ax1.set_xlabel("Distance en km")

  ax3 = ax1.twinx()
  ax3.fill_between(np.linspace(0, 50, 500), slope_elev, "y", alpha=0.2, label="Elevation")


  ax2 = ax1.twinx()
  ax2.plot(np.linspace(0, 50, 500), np.mean(np.asarray(results), axis=0).flatten(), label="Puissance")
  ax1.plot(np.linspace(0, 50, 500), velocity_from_power(power=np.mean(np.asarray(results), axis=0), air_density=1.292, cda=0.29, wind_speed=np.zeros((500)), rolling_friction=0.004, masse=85, slope=slope_profile, efficiency=0.98), 'r', label="Vitesse")
  ax1.set_ylabel("Vitesse m/s")
  ax2.set_ylabel("Puissance W")
  ax2.set_ylim(0,600)
  fig.legend(loc='upper right')

png

On remarque que pour tenir une vitesse moyenne de \(8 m.s^{-1}\), le cycliste devra parcourir la partie plate à \(7.7 m.s^{-1}\) pour une puissance de \(290W\) et la côte à \(5.5 m.s^{-1}\) pour \(390W\). Cette stratégie s’apparente à la stratégie communément adoptée de dépenser plus de puissance dans les difficultés et moins dans les parties faciles de manière à minimiser le temps total de parcours.

On peut comparer notre résultat avec d’autres stratégies qui respecte la même moyenne de puissance sur les \(50km\) :

  • Stratégie constante : tout le parcours à \(340W\).
  • Stratégie inverse : \(390W\) sur le plat et \(290W\) dans la côte.
  • Stratégie aléatoire : tout le parcours à \(340W\) avec une puissance instantanée aléatoire.
t_computed =  np.sum((50000/500)/velocity_from_power(power=np.mean(np.asarray(results), axis=0), air_density=1.292, cda=0.29, wind_speed=np.zeros((500)), rolling_friction=0.004, masse=85, slope=slope_profile, efficiency=0.98))
t_constante = 25000/velocity_from_power(power=340, air_density=1.292, cda=0.29, wind_speed=0, rolling_friction=0.004, masse=85, slope=0, efficiency=0.98)[0] + 25000/velocity_from_power(power=340, air_density=1.292, cda=0.29, wind_speed=0, rolling_friction=0.004, masse=85, slope=8, efficiency=0.98)[0]
t_inverted = 25000/velocity_from_power(power=390, air_density=1.292, cda=0.29, wind_speed=0, rolling_friction=0.004, masse=85, slope=0, efficiency=0.98)[0] + 25000/velocity_from_power(power=290, air_density=1.292, cda=0.29, wind_speed=0, rolling_friction=0.004, masse=85, slope=8, efficiency=0.98)[0]
t_random = np.sum((50000/500)/velocity_from_power(power=np.random.normal(340, 30, 500), air_density=1.292, cda=0.29, wind_speed=np.zeros((500)), rolling_friction=0.004, masse=85, slope=slope_profile, efficiency=0.98))

with plt.xkcd():
  plt.bar(0, t_computed/60, label="Optimale")
  plt.bar(1, t_constante/60, label="Constante")
  plt.bar(2, t_inverted/60, label="Inversee")
  plt.bar(3, t_random/60, label="Aleatoire")
  plt.ylim(100,160)
  plt.legend()
  plt.gca().axes.get_xaxis().set_visible(False)
  plt.ylabel("Temps en minutes")

png

On remarque que notre stratégie est bien la plus rapide des 4 testées. Comme prévu, la stratégie inverse de la stratégie optimale est la plus lente. La stratégie aléatoire (std \(\pm30W\)) est quant à elle équivalente à la stratégie constante.

On peut aussi simulerer un contre-la-montre de avec de vent de face puis de vent de dos pour une vitesse cible de \(11.7m.s^{-1}\). Tous les autres paramètres du modèle resteront constants. On calculera 5 profils optimaux que l’on moyennera pour réduire les changements de puissance abrupts irréalistes.

png

On remarque que la statégie optimale semble être d’augmenter la puissance vent de face et de la réduire légèrement vent de dos de manière à gagner du temps sur la partie difficile mais de ne pas trop en perdre dans la partie facile. Ce résultat semble cohérent avec Atkinson, Greg, and Adam Brunskill. “Pacing strategies during a cycling time trial with simulated headwinds and tailwinds.” Ergonomics 43.10 (2000): 1449-1460. et Anton, A. Brad. “Optimal time-trial bicycle racing with headwinds and tailwinds.” arXiv preprint arXiv:1309.1741 (2013).

Conclusion

On a dans un premier temps établit un modèle permettant de calculer la puissance nécessaire à un cycliste pour surmonter la somme des résistances au roulement. On a ensuite essayé de s’approcher de la répartition de puissance optimale, ie. celle qui consomme le moins d’énergie, pour terminer un parcours avec une vitesse moyenne donnée. Cet optimum a pu être calculé en utilisant un algorithme génétique à partir d’une fonction de fitness simplifiée.

On remarque que cette approche a tendance à donner une courbe de puissance ayant de grands changements d’amplitudes peu réalistes (ce qui se retrouve souvent dans la littérature quand cette approche est utilisée, par exemple pour le calcul de la trajectoire optimale dans un virage). Pour atténuer cet effet, on a donc pris plusieurs réalisations de trajectoires optimales et calculé leur moyenne par intervalle de distance.

Ce modèle basique pourra être amélioré pour obtenir des résultats plus réalistes :

  • On peut complexifier le modèle de relation puissance-vitesse en y ajoutant des résistances supplémentaires. Par exemple en intégrant l’inertie des éléments en rotation comme les roues.
  • Trouver une fonction de fitness plus réaliste, par exemple en ajoutant un terme lié à la répartition de puissance du cycliste. On pourra alors simuler plusieurs profils de coureur (un “rouleur” aura une répartition de puissance centrée autour de la moyenne, ie. sa FTP, alors qu’un “puncher” pourra avoir de plus grosses déviations à la moyenne).