Last year for my cycling team’s secret santa I organised a time-delay chase race up “The Grade”, a virtual climb on Zwift. This is a 10-20 minute effort for most, so I used the same methodology as I did for the race in Innsbruck to set time delays based on people’s 30 day best power.

This year I’ve decided unilaterally that we’re tackling Alpe du Zwift instead. This is a approximately 1 hour effort for most normal people, so I may need to refine my methodology. People might not necessarily have done a full gas 60+ minute effort in the last 90 days, and a critical power model won’t necessarily be valid over this duration. So Let’s explore some approaches we can take.

*12.2km at an average gradient of 8%. Graphic from VeloViewer.*

Use best previous efforts?

The most straightforward way would be to look at each rider’s best effort up the Alpe:

Show code

fastest_efforts %>% 
  select(elapsed, wkg, date) %>% 
  mutate(
    elapsed = seconds_to_period(round(elapsed, 0)),
    date = as_datetime(date)
    )

# A tibble: 12 × 3
   elapsed     wkg date               
   <Period>  <dbl> <dttm>             
 1 1H 0M 16S   2.9 2024-11-30 07:08:00
 2 50M 3S      3.9 2024-10-13 06:48:36
 3 1H 4M 29S   3   2024-03-16 12:12:18
 4 1H 1M 45S   2.8 2022-12-29 07:14:57
 5 52M 20S     3.6 2024-10-13 06:50:53
 6 43M 15S     4.4 2025-01-02 19:39:28
 7 47M 47S     3.9 2022-07-15 21:55:55
 8 56M 6S      3.4 2021-12-30 18:49:06
 9 50M 4S      3.7 2022-07-14 12:27:44
10 56M 39S     3.4 2024-01-06 12:14:55
11 58M 44S     3.3 2023-11-18 09:01:32
12 51M 28S     3.8 2024-03-16 09:53:26

The problem here is the most recent of these efforts was over a year ago, and the latest half a decade ago. And one of the contestents has never attempted it. So we’re going to have to use a model to account for people’s changes in fitness.

Estimate plausible times from a model?

Modelling time to complete the climb for a given w/kg

We can start by building a model of time to climb the segment. Similar to previous approaches, I can pull the best efforts for my team:

Show code

data_for_plot <- segment_data_tv %>% 
  select(zwid, power, weight, elapsed) %>% 
  mutate(
    wkg = power/weight,
    elapsed_h = elapsed / 3600)

data_for_plot %>% ggplot() +
  aes(x = wkg, y = elapsed_h) +
  geom_point() +
  scale_x_continuous(name = "Average W/kg", 
                     breaks = seq(0.5, 5, by = 0.5), 
                     minor_breaks = seq(0.5, 5, by = 0.1)) +
  scale_y_continuous(name = "Elapsed time (hours)") +
  theme_bw()

This once again looks quite power law so we can try to fit \(t = \alpha (\text{w/kg}) ^ \beta\), given that on Zwift the effects of weather changes are non-existence and rolling resistance shouldn’t change (unless someone tried their best effort on a mountain bike):

Show code

library(tidymodels)

data_for_analysis <- data_for_plot %>% 
  filter(wkg >= 2.0 & wkg <= 4.5) %>%
  mutate(
    t_s    = as.numeric(elapsed),
    ln_t   = log(as.numeric(elapsed)),
    ln_wkg = log(as.numeric(wkg))
  )

lm_spec <- linear_reg() %>%
  set_mode("regression") %>%
  set_engine("lm")

lm_fit <- lm_spec %>% 
  fit(ln_t ~ ln_wkg, data = data_for_analysis)

data_for_analysis %>% 
  ggplot() +
  aes(x = ln_wkg, y = ln_t) +
  geom_point() +
  scale_x_continuous(labels = function(x) exp(x),
                     breaks = function(x) log(pretty(exp(x)))) +
  scale_y_continuous(labels = function(y) exp(y),
                     breaks = function(y) log(pretty(exp(y)))) +
  labs(x = "Average W/kg", y = "Elapsed time (seconds)") +
  geom_smooth(method = "lm") +
  theme_bw()

This doesn’t look too bad. We could try and improve on this with a generalised additive model but to be honest the fit (in blue) isn’t much better, and the power-law model (red) remains biologically plausible:

Show code

library(mgcv)

gam_fit <- gam(
  t_s ~ s(wkg, k = 8),
  data = data_for_analysis,
  method = "REML"
)

grid <- data_for_analysis %>%
  summarise(wkg_min = min(wkg, na.rm = TRUE),
            wkg_max = max(wkg, na.rm = TRUE)) %>%
  { tibble(wkg = seq(.$wkg_min, .$wkg_max, length.out = 300)) } %>% 
  mutate(ln_wkg = log(wkg))

pred_lm  <- predict(lm_fit,  new_data = grid)
pred_gam <- predict.gam(gam_fit, newdata = grid) %>% as_tibble()

grid_pred <- grid %>% mutate(
  fit_lm  = exp(pred_lm$.pred),
  fit_gam = pred_gam$value
  )

data_for_analysis %>% 
  ggplot() +
  aes(x = wkg, y = t_s) +
  geom_point() +
  geom_line(data = grid_pred, aes(x = wkg, y = fit_lm), color = "red") +
  geom_line(data = grid_pred, aes(x = wkg, y = fit_gam), color = "blue") +
  labs(x = "Average W/kg", y = "Elapsed time (seconds)") +
  theme_bw()

Turning power-curves into plausible power-duration models

Next step is to pull everyone’s power curves from ZwiftPower:

Show code

power_curves %>% 
  ggplot() + 
  aes(x = x, y = y_90) +
  geom_line() +
  scale_x_log10() +
  labs(x = "Duration (seconds)", y = "Maximal power (W)") +
  facet_wrap(~id, nrow = 4) +
  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.text = element_blank()
  )

The problem here is that it’s clear some people haven’t done maximal efforts beyond 40-50 minutes in the last 90 days, as shown by dramatic steps in what should be a smooth curve. This means that modelling time to climb a ~60 minute effort will potentially be a huge under-estimate, and their time delay will be inappapriately huge, leading to a fairly boring chase race.

Here’s one particular offender. Names are removed to protect the innocent, but despite clearly exceptional sprint power he hasn’t done a maximal effort beyond 40 minutes (dotted line) in the last three months:

Show code

mmp_offender %>% 
  ggplot() +
  aes(x = x, y = y_90) +
  geom_point() +
  scale_x_log10() +
  geom_vline(xintercept = 40*60, linetype = "dotted") +
  labs(x = "Duration (seconds)", y = "Maximal power (W)") +
  theme_bw()

Option 1: Critical power model

We could try to fit a CP model to this using their 1 to 15 minute power:

Show code

mmp_cp <- mmp_offender %>% filter(x >= 1*60 & x <= 15*60)

mmp_fit <- nls(
  y_90 ~ CP + Wp / x,
  data = mmp_cp,
  start = list(CP = 100,
               Wp = 1000),
  algorithm = "port",
  lower = c(CP = 0, Wp = 0)
)

pred_cp <- tibble(x = seq(60*1, 60*60*2, by = 10)) %>%
  mutate(pred = predict(mmp_fit, newdata = .))

ggplot() +
  geom_point(data = mmp_offender %>% filter(x >= 30 & x <= 7200), aes(x, y_90)) +
  geom_line(data = pred_cp, aes(x, pred), colour = "darkred") +
  scale_x_log10("Duration (s)") +
  scale_y_continuous("Power (W)") +
  theme_bw()

This might work, but CP tends to over-estimate the true power that can be held for longer duration efforts - you can already see the model diverging around the 30 minute mark.

Option 2: Power-law tail model

Instead we need a “credible” modelled power-duration function, using the observed maximal power over a reasonable endurance-ish timeframe (say 10 to 30 minutes) and extrapolating a tail from this using a power law. We also need a function to detect if they’ve been sandbagging even within this window, and cut it off at any obvious step.

Show code

# Clean up the power curve
prep_mmp <- function(mmp) {
  mmp %>%
    transmute(t = as.numeric(x), P = as.numeric(y_90)) %>%
    filter(is.finite(t), is.finite(P), t > 0, P > 0) %>%
    group_by(t) %>% summarise(P = max(P), .groups = "drop") %>%
    arrange(t)
}

# Fit a power-law tail on 10-30 minute efforts
fit_tail_powerlaw <- function(df,
                              t_min = 10*60,
                              t_max = 30*60) {
  d <- df %>%
    filter(t >= t_min, t <= t_max) %>%
    mutate(logt = log(t), logP = log(P))

  fit <- lm(logP ~ logt, data = d)

  list(
    a = exp(unname(coef(fit)[["(Intercept)"]])),
    b = unname(coef(fit)[["logt"]]),
    fit = fit
  )
}

# Detect significant steps in the power curve
tail_missing_flag <- function(df, 
                              tail,
                              check_s = 3600,
                              ratio_cut = 0.90) {
  P_obs  <- approxfun(df$t, df$P, rule = 2)
  P_pred <- tail$a * (check_s ^ tail$b)
  P_obs(check_s) < ratio_cut * P_pred
}

# Predict the power curve
make_power_predictor <- function(mmp,
                                 switch_s = 1800) {
  df <- prep_mmp(mmp)
  tail <- fit_tail_powerlaw(df, t_min = 600, t_max = 1800)

  missing_tail <- tail_missing_flag(df, tail, check_s = 3600, ratio_cut = 0.90)
  
  # Cut at 20 min if tail missing
  switch_eff <- if (missing_tail) 20*60 else switch_s   

  P_obs  <- approxfun(df$t, df$P, rule = 2)
  P_tail <- function(t) tail$a * (t ^ tail$b)
  P_hat <- function(t) ifelse(t <= switch_eff, P_obs(t), P_tail(t))

  list(
    P_hat = P_hat,
    tail_a = tail$a, tail_b = tail$b,
    switch_s = switch_eff,
    tail_missing = missing_tail,
    df = df
  )
}

Let’s try fitting this to our offending rider:

Show code

pm <- make_power_predictor(mmp_offender, switch_s = 1800)

df <- pm$df %>% filter(t >= 30, t <= 7200)

grid <- tibble(t = exp(seq(log(min(df$t)),
                           log(max(df$t)),
                           length.out = 400))) %>%
  mutate(P_hat = pm$P_hat(t),
         source = if_else(t <= pm$switch_s, "Observed (interp)", "Model tail"))

ggplot() +
  geom_point(data = df, aes(t, P)) +
  geom_line(data = grid, aes(t, P_hat), color = "darkgreen") +
  scale_x_log10("Duration (seconds)") +
  scale_y_continuous("Power (W)") +
  theme_bw()

Now we need to do the same for everyone, which looks a lot more reasonable:

Show code

# Fit models for all riders
riders_mod <- riders %>%
  mutate(power_model = map(mmp, make_power_predictor))

# Build predicted curves (log-spaced duration grid)
curves <- riders_mod %>%
  transmute(
    zwid, name,
    df = map(power_model, "df"),
    P_hat_fn = map(power_model, "P_hat"),
    switch_s = map_dbl(power_model, "switch_s"),
    tail_missing = map_lgl(power_model, "tail_missing")
  ) %>%
  mutate(
    grid = pmap(
      list(df, P_hat_fn),
      \(df, P_hat_fn) {
        tmin <- max(30, min(df$t, na.rm = TRUE))
        tmax <- min(7200, max(df$t, na.rm = TRUE))
        t_grid <- exp(seq(log(tmin), log(tmax), length.out = 400))

        tibble(
          t = t_grid,
          P_hat = P_hat_fn(t_grid)
        )
      }
    )
  ) %>%
  select(zwid, name, switch_s, tail_missing, grid) %>%
  unnest(grid)

# Keep the observed power curves for plotting
obs <- riders_mod %>%
  transmute(
    zwid, name,
    df = map(power_model, "df")
  ) %>%
  unnest(df) %>%
  filter(t >= 30, t <= 7200) %>% 
  arrange(zwid)

# Plot the result
ggplot() +
  geom_line(data = obs, aes(t, P)) +
  geom_line(data = curves, aes(t, P_hat), color = "darkgreen") +
  scale_x_log10("Duration (seconds)") +
  ylab("Power (W)") +
  facet_wrap(~ zwid, nrow = 4) +
  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.text = element_blank()
  )

Predict climbing time from the model

Given we now have more reliable looking estimates of power-duration, we can solve for time up the Alpe in a similar way to before. Plotting these predictions on the previous attempts looks reassuring:

Show code

# Helper function to turn wkg into a predicted time
time_from_wkg <- function(wkg, wf_fit) {
  stopifnot(all(is.finite(wkg)), all(wkg > 0))

  new_data <- tibble(ln_wkg = log(wkg))

  ln_t_hat <- predict(wf_fit, new_data = new_data) %>%
    pull(.pred)

  as.numeric(exp(ln_t_hat))
}

# Function to predict time up AdZ which we can itterate over riders
predict_adz_time <- function(weight_kg, power_model, time_model = lm_fit,
                             t_lower = 20*60, t_upper = 120*60,
                             use_mean = FALSE) {

  time_fun <- if (use_mean) time_from_wkg_mean else time_from_wkg

  f <- function(T) {
    wkg <- power_model$P_hat(T) / weight_kg
    time_fun(wkg, time_model) - T
  }

  out <- tryCatch(uniroot(f, lower = t_lower, upper = t_upper),
                  error = function(e) NULL)

  if (is.null(out)) return(NA_real_)
  out$root
}

# Do the prediction
riders_pred <- riders %>%
  mutate(
    power_model = map(mmp, make_power_predictor),
    adz_time_s = map2_dbl(
      weight_kg, power_model,
      ~ predict_adz_time(weight_kg = .x, power_model = .y, time_model = lm_fit,
                         t_lower = 10*60, t_upper = 150*60, use_mean = FALSE)
    ),
    wkg_at_pred = map2_dbl(power_model, adz_time_s, \(pm, t) pm$P_hat(t)) / weight_kg,
    tail_missing = map_lgl(power_model, "tail_missing"),
    switch_min = map_dbl(power_model, "switch_s")/60
  )

data_for_plot %>% ggplot() +
  aes(x = wkg, y = elapsed_h) +
  geom_point() +
  geom_point(data = riders_pred, aes(x = wkg_at_pred, y = adz_time_s/(60*60)), color = "red") +
  scale_x_continuous(name = "Average W/kg", 
                     breaks = seq(0.5, 5, by = 0.5), 
                     minor_breaks = seq(0.5, 5, by = 0.1)) +
  scale_y_continuous(name = "Elapsed time (hours)") +
  theme_bw()

Finally, we can turn this into time delays:

Show code

t_slowest = tail(riders_pred %>% arrange(adz_time_s), 1)$adz_time_s

prediction_table <- riders_pred %>% 
  mutate(t_delay = round(t_slowest - adz_time_s, 0))

prediction_table %>%
  left_join(fastest_efforts, by = "name") %>% 
  arrange(adz_time_s) %>% 
  mutate(
    t_pred    = seconds_to_period(round(adz_time_s, 0)),
    wkg_pred  = round(wkg_at_pred, 1),
    t_best    = seconds_to_period(round(elapsed, 0)),
    wkg_best  = round(wkg, 1),
    date_best = as_datetime(date),
    t_delay   = seconds_to_period(round(t_delay, 0))
  ) %>% 
  select(t_pred, wkg_pred, t_best, wkg_best, date_best, t_delay)

# A tibble: 12 × 6
   t_pred    wkg_pred t_best    wkg_best date_best           t_delay 
   <Period>     <dbl> <Period>     <dbl> <dttm>              <Period>
 1 45M 26S        4.3 43M 15S        4.4 2025-01-02 19:39:28 16M 17S 
 2 49M 20S        3.9 58M 44S        3.3 2023-11-18 09:01:32 12M 23S 
 3 49M 44S        3.9 1H 4M 29S      3   2024-03-16 12:12:18 12M 0S  
 4 50M 48S        3.8 NA            NA   NA                  10M 55S 
 5 51M 14S        3.7 50M 3S         3.9 2024-10-13 06:48:36 10M 29S 
 6 52M 7S         3.7 52M 20S        3.6 2024-10-13 06:50:53 9M 37S  
 7 53M 11S        3.6 1H 0M 16S      2.9 2024-11-30 07:08:00 8M 33S  
 8 54M 41S        3.5 1H 1M 45S      2.8 2022-12-29 07:14:57 7M 3S   
 9 55M 43S        3.4 47M 47S        3.9 2022-07-15 21:55:55 6M 1S   
10 56M 33S        3.3 51M 28S        3.8 2024-03-16 09:53:26 5M 10S  
11 56M 52S        3.3 56M 39S        3.4 2024-01-06 12:14:55 4M 51S  
12 1H 1M 43S      3   56M 6S         3.4 2021-12-30 18:49:06 0S

For people doubting this, we can also show where the model gets the prediction from:

Show code

adz_pts <- riders_pred %>%
  arrange(zwid) %>% 
  transmute(
    zwid, name, weight_kg, adz_time_s,
    P_hat_fn = map(power_model, "P_hat")
  ) %>%
  mutate(
    P_at_adz = map2_dbl(P_hat_fn, adz_time_s, \(f, t) ifelse(is.na(t), NA_real_, f(t))),
    wkg_at_adz = P_at_adz / weight_kg
  ) %>%
  select(-P_hat_fn)

ggplot() +
  geom_line(data = obs, aes(t, P)) +
  geom_line(data = curves, aes(t, P_hat), color = "darkgreen") +
  geom_vline(data = adz_pts, aes(xintercept = adz_time_s),
             linetype = "dotted", linewidth = 0.6, alpha = 0.8) +
  geom_hline(data = adz_pts, aes(yintercept = P_at_adz),
             linetype = "dotted", linewidth = 0.6, alpha = 0.8) +
  geom_point(data = adz_pts, aes(adz_time_s, P_at_adz), size = 1.6, color = "red") +
  scale_x_log10("Duration (seconds)") +
  ylab("Power (W)") +
  facet_wrap(~ zwid, ncol = 3) +
  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.text = element_blank()
  )

Finally, it goes without saying that this is very much a back of the envelope calculation to allow me to set vaguely competitive time delays. There are a load of assumptions involved, and the potential error in the power-law model will get wider the further it gets from the window used to derive it, so I’ll be interested to see how the results stack up.