estimate_cda_crr.rb

#!/usr/bin/env ruby
# frozen_string_literal: true

require "csv"
require "optparse"

options = {
  mass: 70.0,
  efficiency: 0.97,
  wind: 0.0,
  window: 30.0,
  step: nil,
  csv: File.join(__dir__, "fujihilllog.csv"),
  include_ke: true,
  fixed_crr_values: [0.0025, 0.0030, 0.0035, 0.0040, 0.0045, 0.0050, 0.0055, 0.0060],
  fixed_cda_values: [0.20, 0.25, 0.30, 0.35, 0.40]
}

OptionParser.new do |opts|
  opts.banner = "Usage: ruby estimate_cda_crr.rb [options] [csv]"

  opts.on("--mass KG", Float, "Total rider+bike mass, default: #{options[:mass]}") { |v| options[:mass] = v }
  opts.on("--efficiency N", Float, "Drivetrain efficiency, default: #{options[:efficiency]}") { |v| options[:efficiency] = v }
  opts.on("--tailwind-kph KPH", Float, "Tailwind speed along route, default: 0") { |v| options[:wind] = v / 3.6 }
  opts.on("--headwind-kph KPH", Float, "Headwind speed along route, default: 0") { |v| options[:wind] = -v / 3.6 }
  opts.on("--window SECONDS", Float, "Fit window size, default: #{options[:window]}") { |v| options[:window] = v }
  opts.on("--step SECONDS", Float, "Window step, default: same as --window") { |v| options[:step] = v }
  opts.on("--no-ke", "Ignore kinetic-energy changes in each window") { options[:include_ke] = false }
end.parse!

options[:csv] = ARGV.shift if ARGV[0]
options[:step] ||= options[:window]

G = 9.80665

def density_at_altitude(meters)
  # International Standard Atmosphere, troposphere. Good enough for this estimate.
  1.225 * [(1.0 - 2.25577e-5 * meters), 0.0].max**4.25588
end

def f(row, name)
  value = row[name]
  raise "missing #{name}" if value.nil?

  value.to_f
end

rows = []
File.open(options[:csv], "r:bom|utf-8") do |file|
  CSV.new(file, headers: true).each do |row|
    rows << {
      t: f(row, "Timestamp"),
      power: f(row, "Power (watt)"),
      elevation: f(row, "Elevation (meter)"),
      distance: f(row, "Distance (meter)"),
      speed: f(row, "Speed (m/s)")
    }
  end
end

rows.sort_by! { |row| row[:t] }
rows = rows.chunk { |row| row[:t] }.map { |_t, same_time| same_time.first }

intervals = rows.each_cons(2).filter_map do |a, b|
  dt = b[:t] - a[:t]
  ds = b[:distance] - a[:distance]
  next unless dt.between?(0.5, 2.0) && ds.positive?

  v = ds / dt
  next unless v.between?(0.5, 20.0)

  elevation = (a[:elevation] + b[:elevation]) / 2.0
  rho = density_at_altitude(elevation)
  power = (a[:power] + b[:power]) / 2.0
  air_speed = [v - options[:wind], 0.0].max

  {
    dt: dt,
    ds: ds,
    dh: b[:elevation] - a[:elevation],
    v0: a[:speed],
    v1: b[:speed],
    input_energy: options[:efficiency] * power * dt,
    gravity_work: options[:mass] * G * (b[:elevation] - a[:elevation]),
    rolling_x: options[:mass] * G * ds,
    aero_x: 0.5 * rho * air_speed**2 * v * dt
  }
end

def prefix_sums(intervals, key)
  sums = [0.0]
  intervals.each { |interval| sums << sums[-1] + interval[key] }
  sums
end

prefix = %i[dt ds input_energy gravity_work rolling_x aero_x].to_h do |key|
  [key, prefix_sums(intervals, key)]
end

window_count = options[:window].round
step_count = options[:step].round
windows = []

window_start = 0
while window_start + window_count <= intervals.length
  window_finish = window_start + window_count
  dt = prefix[:dt][window_finish] - prefix[:dt][window_start]
  ds = prefix[:ds][window_finish] - prefix[:ds][window_start]
  current_start = window_start
  window_start += step_count
  next unless dt >= options[:window] * 0.8 && ds.positive?

  gravity_work = prefix[:gravity_work][window_finish] - prefix[:gravity_work][current_start]
  grade = gravity_work / (options[:mass] * G * ds)
  next unless grade.between?(-0.05, 0.20)

  kinetic_energy = if options[:include_ke]
                     0.5 * options[:mass] * (intervals[window_finish - 1][:v1]**2 - intervals[current_start][:v0]**2)
                   else
                     0.0
                   end

  input_energy = prefix[:input_energy][window_finish] - prefix[:input_energy][current_start]
  rolling_x = prefix[:rolling_x][window_finish] - prefix[:rolling_x][current_start]
  aero_x = prefix[:aero_x][window_finish] - prefix[:aero_x][current_start]

  windows << {
    y: input_energy - gravity_work - kinetic_energy,
    rolling_x: rolling_x,
    aero_x: aero_x,
    dt: dt
  }
end

def fit_two_parameter(windows)
  s_rr = windows.sum { |w| w[:rolling_x] * w[:rolling_x] }
  s_ra = windows.sum { |w| w[:rolling_x] * w[:aero_x] }
  s_aa = windows.sum { |w| w[:aero_x] * w[:aero_x] }
  s_ry = windows.sum { |w| w[:rolling_x] * w[:y] }
  s_ay = windows.sum { |w| w[:aero_x] * w[:y] }
  det = s_rr * s_aa - s_ra * s_ra
  raise "fit is singular" if det.abs < 1e-9

  crr = (s_ry * s_aa - s_ay * s_ra) / det
  cda = (s_rr * s_ay - s_ra * s_ry) / det
  [crr, cda]
end

def rmse_watts(windows, crr, cda)
  mse = windows.sum do |w|
    err = crr * w[:rolling_x] + cda * w[:aero_x] - w[:y]
    err * err
  end / windows.length
  avg_dt = windows.sum { |w| w[:dt] } / windows.length
  Math.sqrt(mse) / avg_dt
end

def fit_fixed_crr(windows, crr)
  numerator = windows.sum { |w| w[:aero_x] * (w[:y] - crr * w[:rolling_x]) }
  denominator = windows.sum { |w| w[:aero_x] * w[:aero_x] }
  numerator / denominator
end

def fit_fixed_cda(windows, cda)
  numerator = windows.sum { |w| w[:rolling_x] * (w[:y] - cda * w[:aero_x]) }
  denominator = windows.sum { |w| w[:rolling_x] * w[:rolling_x] }
  numerator / denominator
end

duration = intervals.sum { |i| i[:dt] }
distance = intervals.sum { |i| i[:ds] }
gain = rows[-1][:elevation] - rows[0][:elevation]
avg_power = intervals.sum { |i| i[:input_energy] } / duration
gravity_power = intervals.sum { |i| i[:gravity_work] } / duration

crr, cda = fit_two_parameter(windows)

puts "Input"
puts "  csv: #{options[:csv]}"
puts "  mass: #{format('%.1f', options[:mass])} kg"
puts "  drivetrain efficiency: #{format('%.3f', options[:efficiency])}"
puts "  tailwind: #{format('%.3f', options[:wind])} m/s (#{format('%.1f', options[:wind] * 3.6)} kph)"
puts "  air density: ISA by logged altitude"
puts
puts "Ride"
puts "  samples: #{rows.length}"
puts "  duration: #{format('%.0f', duration)} s"
puts "  distance: #{format('%.1f', distance)} m"
puts "  elevation gain from endpoints: #{format('%.1f', gain)} m"
puts "  avg logged/wheel power used: #{format('%.1f', avg_power)} W"
puts "  avg gravity power: #{format('%.1f', gravity_power)} W"
puts
puts "Fit, #{options[:window].round}s windows"
puts "  windows: #{windows.length}"
puts "  unconstrained Crr: #{format('%.5f', crr)}"
puts "  unconstrained CdA: #{format('%.3f', cda)} m^2"
puts "  RMSE: #{format('%.1f', rmse_watts(windows, crr, cda))} W"
puts
puts "Fixed Crr sweep"
options[:fixed_crr_values].each do |fixed_crr|
  fixed_cda = fit_fixed_crr(windows, fixed_crr)
  puts "  Crr #{format('%.4f', fixed_crr)} -> CdA #{format('%.3f', fixed_cda)} m^2, RMSE #{format('%.1f', rmse_watts(windows, fixed_crr, fixed_cda))} W"
end
puts
puts "Fixed CdA sweep"
options[:fixed_cda_values].each do |fixed_cda|
  fixed_crr = fit_fixed_cda(windows, fixed_cda)
  puts "  CdA #{format('%.2f', fixed_cda)} m^2 -> Crr #{format('%.5f', fixed_crr)}, RMSE #{format('%.1f', rmse_watts(windows, fixed_crr, fixed_cda))} W"
end