Vertual Elevation Method

パワーとスピードのログデータから高度変化をグラフにして描き、CdAやCrrを推定する。

ロードバイクのパワーとスピードのデータから高度変化をグラフにして描き、CdAやCrrを推定する。インプットデータとして、GoldenCheetah( http://www.goldencheetah.org/ )からエクスポートしたJSONファイル(.json)を利用する。

library(jsonlite)

## Warning: package 'jsonlite' was built under R version 3.2.5

d <- fromJSON("2016-05-19sample.json") #バンクの中盤を回っているデータ

## Warning: JSON string contains (illegal) UTF8 byte-order-mark!

d <- d$RIDE$SAMPLES
head(d)

##   SECS    KM WATTS CAD    KPH ALT      LAT      LON SLOPE TEMP LRBALANCE
## 1    0 0.000   163  64 20.772 5.2 31.20088 130.7948  10.9   28        44
## 2    1 0.006   177  64 20.592 5.2 31.20093 130.7948   7.6   28        48
## 3    2 0.012   154  64 20.988 5.2 31.20103 130.7948   7.6   28        52
## 4    3 0.018    69  76 20.988 6.6 31.20107 130.7948   7.6   28        63
## 5    4 0.023    74  74 21.168 6.6 31.20112 130.7948   9.1   28        59
## 6    5 0.029   116  65 20.988 6.6 31.20116 130.7948   9.1   28        54
##    LTE RTE  LPS  RPS
## 1 75.0   0 23.0  0.0
## 2 91.5   0 34.5  0.0
## 3 83.0  80 26.5 26.5
## 4 77.5   0 24.0  0.0
## 5  0.0   0  0.0  0.0
## 6  0.0   0  0.0  0.0

幾つかのパラメタを設定する。参考資料( https://github.com/GoldenCheetah/GoldenCheetah/blob/master/doc/contrib/ChungVE.pdf )

m = 72 #バイクと人間の総重量
g = 9.81 #重力加速度
Crr = 0.004 #推定路面転がり抵抗　参考：アスファルト*ロードタイヤ 0.005
#コンクリート*ロードタイヤ 0.004
#木製バンク*ピストタイヤ 0.001
rho = 1.20 #大気密度　気圧・気温から推定
CdA = 0.219 #推定CdA
eta = 0.98 #ドライブトレインの効率。この場合、クランクのパワーメーターで計測された98%がホイールに伝達されているものとする

各時点のスロープ（タンジェント）を推定し、各時点の速さとかけて高度変化を計算し、足し合わせる

library(ggplot2) #グラフ作成パッケージ
library(dplyr) #データ操作パッケージ

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

d %>% mutate(KPS = KPH*1000/3600, #秒速になおす
             a = KPS - lag(KPS)) -> d #各時点の加速度を計算する
d %>% mutate(s = WATTS*eta/m/g/KPS - Crr - a/g - rho * CdA * KPS^2/2/m/g) -> d #スロープsを計算する
d %>% mutate(delta = s * KPS) -> d #s各時点の標高変化deltaを計算
d$delta[-1] %>% cumsum -> v #標高変化の各時点までの累積和ベクトルvを生成
v[166:813] -> v #ポジションを一定にしていた部分を切り出し
qplot(seq_along(v), v, geom="line") #グラフを生成

500秒あたりで大きくポジションを変えてしまっているが、それ以外の部分でほぼ周期的かつ水平にグラフが推移していることから、そこそこよいパラメタの数値である

これはバンクの中盤を回っているので、本来一周ごとではなく、半周ごとに標高変化の周期は現れるはずだが、半周ごとでは標高が合わない ->風の影響

緯度、経度の推移から進行方向を決め（北が0度、東側がプラス、西側がマイナス、南が+-180度）当時のデータから東側7.5km/h（90度、7.5km/h）の風を対空速度に補正する

library(geosphere) #緯度経度まわりの計算用パッケージ

## Loading required package: sp

## Warning: package 'sp' was built under R version 3.2.5

direc <- bearing(cbind(d$LON,d$LAT),cbind(lead(d$LON),lead(d$LAT))) #各時点1秒前の緯度経度と比べて進行方向を割り出す
d$KPH_w <- d$KPH + cospi((direc-90)/180)*7.5 #90度、すなわち東風7.5m/hの補正を対空速度に加える。
d %>% mutate(KPS_w = KPH_w/3600*1000) -> d #秒速に変換
d %>% mutate(s_w = WATTS*eta/m/g/KPS - Crr - a/g - rho * CdA * KPS_w^2/2/m/g) -> d #対空速度を利用してスロープを算出
d %>% mutate(delta_w = s_w * KPS) -> d #標高変化に変換
d$delta_w[-1] %>% cumsum -> v_w #各時点までの累積和を計算
v_w[166:813] -> v_w #ラップの切り出し
qplot(seq_along(v_w), v_w, geom="line") #グラフ化

周期が半周ごとになり、うまく補正できていることがわかる。

Vertual Elevation Method

パワーとスピードのログデータから高度変化をグラフにして描き、CdAやCrrを推定する。

課題