# Fusing uncalibrated IMUs and handheld smartphone video to reconstruct knee
kinematics

J.D. Peiffer${}^{1,2,\href https://orcid.org/0000-0003-2382-8065}$, Kunal
Shah1, Shawana Anarwala1, Kayan Abdou1 and R. James Cotton${}^{1,3,\href
https://orcid.org/0000-0001-5714-1400}$ This work was generously supported by
the American Neurological Foundation, the Restore Center P2C (NIH
P2CHD101913), and the Research Accelerator Program of the Shirley Ryan
AbilityLab. JDP is supported by the National Science Foundation Graduate
Research Fellowship Program under Grant No. DGE-2234667.1 Center for Bionic
Medicine, Shirley Ryan AbilityLab, 355 E Erie St, Chicago, IL
<EMAIL_ADDRESS>Department of Biomedical Engineering, Northwestern
University, Evanston, IL 3 Department of Physical Medicine and Rehabilitation,
Northwestern University, Evanston, IL

###### Abstract

Video and wearable sensor data provide complementary information about human
movement. Video provides a holistic understanding of the entire body in the
world while wearable sensors provide high-resolution measurements of specific
body segments. A robust method to fuse these modalities and obtain
biomechanically accurate kinematics would have substantial utility for
clinical assessment and monitoring. While multiple video-sensor fusion methods
exist, most assume that a time-intensive, and often brittle, sensor-body
calibration process has already been performed. In this work, we present a
method to combine handheld smartphone video and uncalibrated wearable sensor
data at their full temporal resolution. Our monocular, video-only,
biomechanical reconstruction already performs well, with only several degrees
of error at the knee during walking compared to markerless motion capture.
Reconstructing from a fusion of video and wearable sensor data further reduces
this error. We validate this in a mixture of people with no gait impairments,
lower limb prosthesis users, and individuals with a history of stroke. We also
show that sensor data allows tracking through periods of visual occlusion.

##  I INTRODUCTION

The kinematics of human movement can provide important health information. For
instance, upper limb kinematics change significantly in the months following a
stroke [1], while lower limb kinematics can discriminate between healthy and
pathological gait and changes during rehabilitation [2]. Traditionally,
measuring kinematics was performed in motion analysis laboratories using
marker-based methods [3]. This technique detects reflective markers placed on
the surface of the skin, which are then fit with a biomechanical model to
obtain joint angle kinematic trajectories. While this method produces accurate
tracking, the marker placement requires significant domain knowledge and setup
time that limits the usability of this technique in the clinic.

Figure 1: Using video from a handheld moving camera (A), our method
disentangles camera and body orientation changes to accurately track a timed
up and go test in global space (B,C). We validated the joint kinematics
against a multi-camera system (D), finding close agreement between the two
methods.

In recent years, multiview markerless motion capture has seen considerable
improvement and expansion. This technology leverages computer vision to
estimate joint locations from 2D RGB videos captured by multiple calibrated
and synchronized cameras. Instead of physical surface markers, virtual marker
locations, known as ’keypoints’, are detected. As with marker-based systems, a
biomechanical model can be fit to follow these virtual markers. Recent work
has shown markerless systems achieve comparable accuracy to marker-based
accuracy [4] and our previous work has demonstrated its applicability to gait
analysis in clinical populations [5]. We have also shown that using
differentiable biomechanical models allows us to formulate the biomechanical
reconstruction from videos as an end-to-end optimization problem. Directly
optimizing kinematic trajectories to fit the detected keypoints outperforms
the traditional two-stage approach, which first computes marker trajectories
and then applies inverse kinematics [6].

Historically, human pose estimation (HPE) has focused on detecting 2D
keypoints in the image plane. However, contemporary methods such as MeTRAbs-
ACAE now enable precise 3D keypoint estimates directly from images [7].
Concurrently, another approach involves using established mesh body models,
specifically the Skinned Multi-Person Linear (SMPL) model, to optimize or
directly regress shape and pose parameters [8],[9]. The shaping and forward
kinematic processes of this model are differentiable and benefit from GPU
acceleration, rendering them suitable for optimization or regression
techniques. However, the kinematic tree configuration does not accurately
represent anatomical joint locations and articulations [10], thus restricting
the model’s applicability in biomechanics where specific conventions are
established [11]. In this work, we extend our previous work optimizing a
differentiable biomechanical model [6] to 3D keypoints detected from monocular
videos. We also show how this framework can be extended to fuse video data
with wearable sensor data, even if the sensor placement has not been
calibrated.

Wearable sensors, particularly Inertial Measurement Units (IMUs) have been
used with vision-based approaches to enhance human pose estimate accuracy [12,
13]. IMUs include an accelerometer, gyroscope, and magnetometer to measure
linear acceleration, angular velocity, and heading relative to magnetic north,
respectively. These measurements are fused using algorithms like the Kalman or
complementary filter to estimate orientation [14],[15]. Sensor data can be
collected at much higher rates than video, and when fused with stationary
video, they have shown significant improvements in pose estimation [16]. The
approaches are complementary: video tracks the global configuration of a
person, while IMUs provide high-resolution information that is robust to
occlusions. However, extracting biomechanics from sensor data typically
requires a cumbersome calibration process, making it impractical for routine
clinical use.

There have been a few recent approaches to creating portable motion capture
and biomechanics systems. First, OpenCap [17] uses two or more calibrated,
stationary smartphones on tripods to capture RGB images and provides a
convenient web application for users to process their data. OpenSense [18]
uses 8 IMUs to calculate the kinematics and kinetics of a biomechanical model.
Setup and calibration steps, required for both of these, limit the adoption of
motion capture and biomechanics into routine clinical use [19]. We have
previously described our portable smartphone system for capturing synchronous
sensor and video data during gait in the clinic [20]. In that work, sensor
data was only used to detect gait cycles to validate the vision algorithms.
Here we leverage the data from our Portable Biomechanics Laboratory to fuse
data from a single mobile camera and uncalibrated wearable sensors.

Our work addresses two critical gaps: (1) the absence of methods to estimate
biomechanically-grounded kinematics using a single handheld camera, and (2)
the lack of techniques for utilizing uncalibrated IMU data in kinematic
tracking, which is essential for their use in time-constrained clinical
settings. The primary contributions of this work are:

  * •

We extend our differentiable biomechanics approach [6] from multiview
markerless motion capture data to monocular video, which already produces
accurate fits.

  * •

We show that uncalibrated IMU data can be fused with smartphone video,
improving the accuracy of our knee kinematics.

  * •

We demonstrate global orientation tracking of the participant, disentangling
the camera and participant rotation, using both rotations measured from the
sensor and smartphone gyroscope (Figure 1).

##  II METHODS

###  II-A Data collection

This study was approved by the Northwestern University Institutional Review
Board. 38 participants participated in this study, of those 38, 7 participants
were lower limb prosthesis users, 11 had a history of stroke and 20 were able-
bodied controls. All recordings took place in 7.4 x 8 m room with participants
walking diagonally across the room (11 m) for one trial. Participants
completed nine walking trials, 3 at a self selected pace, 3 at a fast, and 3
at a slow pace. Gait data was collected simultaneously with both our wearable
sensor and smartphone platform as well as our multiview markerless motion
capture (MMMC) system. We analyzed 1-2 videos of each participant.

###  II-B Portable Biomechanics Laboratory

Our Portable Biomechanics Laboratory (PBL) system consists of a Samsung Galaxy
S20 smartphone and custom-made wearable IMUs [21, 14]. The smartphone collects
1080 $\times$ 1920 resolution RGB video data at 30 Hz and internal smartphone
gyroscope data at 100 Hz. It also records depth videos, which were not used in
this study. The IMUs collect accelerometer, magnetometer, and gyroscope data
at 562.5 Hz. They run a complementary filter that estimates orientation data
at 55 Hz [14]. The gyroscope and orientation data are used in our fusion. IMUs
were placed unilaterally on the participant’s thigh and lower leg using velcro
straps, typically on the lateral side. For participants with unilateral gait
impairments, sensors were placed on the impaired limb. For each trial, the
researcher followed behind the participant holding the smartphone as the
participant walked across the room (Fig 1).

###  II-C Multiview Markerless Motion Capture

The Multiview Markerless Motion Capture (MMMC) system consists of 8-12 FLIR
BlackFly S GigE cameras which acquire synchronized RGB video at 29 fps.
Cameras were arranged such that at least three cameras covered participants at
all points during a recording. A detailed description of this system can be
found in [22]. We reconstructed biomechanical fits using methods described in
[6]. This uses the same biomechanical model used in the monocular
reconstruction to provide a consistent basis for comparison.

###  II-D Video processing

We processed smartphone and MMMC videos using PosePipe [23], an open-source
tool for human pose estimation. Bounding boxes were first created using
DeepSort [24] and 3D keypoints were estimated using MeTRAbs-ACAE [7]. This
algorithm produces the superset of many keypoint sets, and we used the 87 BML-
MoVi keypoints [25], as this dense markerset is designed for biomechanics. As
recommended by the MeTRAbs-ACAE author, we computed keypoint confidence by
measuring the standard deviation of each joint location estimated from 10
different augmented versions of each video frame. This was converted to a
confidence estimate using a sigmoid function with a half maximum at 30mm and a
width of 10mm. We annotated the participant of interest if multiple
participants were present in the video.

###  II-E Notation

We follow the convention presented in [26] for rotations. In brief, rotation
$R_{ab}$ represents a rotation from reference frame $\\{a\\}$ to reference
frame $\\{b\\}$. If this rotation varies with time, as in the case of a body
segment or sensor reading, we represent that rotation as $R_{ab}(t)$. We
denote arbitrary vector quantities that exhibit temporal variation as
$\vec{V}(t)$, where $\vec{V}$ represents the vector in question and $t$
indicates its dependence on time. We represent quantities that are estimates
of our model with the $\hat{\cdot}$ symbol and measured quantities, such as
video-estimated 3D keypoints and raw sensor readings, with no hat.

###  II-F Body model

For both the multi and single camera data, we used a biomechanically grounded,
forward kinematic model implemented in MuJoCo [27],[28]. A full summary of
this model can be found in [6]. In brief, this model is based on a MuJoCo
implementation [29] of an OpenSim model [30]. We further modified it to
include a neck joint and optimized the default marker locations for the BML-
MoVi keypoints [25]. Forward kinematics of this model can be massively
accelerated and parallelized on the GPU with MuJoCo MJX. Scaling within the
model is executed via a matrix that maps a set of scale parameters to the
isotropic scaling of individual body segments. We defined eight scale
parameters: overall size, the pelvis, left thigh, left leg and foot, right
thigh, right leg and foot, the left arm, and the left leg. We also included an
$87\times 3$ matrix of marker offsets for small changes in the 87 MoVi
keypoint locations. Together, these parameters are represented as a vector
$\vec{\beta}\in\mathbb{R}^{8+87\times 3}$. The model inputs, denoted by
$\vec{\theta}\in\mathbb{R}^{40}$, include six parameters for translating and
orienting the pelvis, and 36 parameters for actuating the remaining segments
of the kinematic tree. We assigned each of the 87 MoVi keypoints to specific
body segments and positioned them using the forward kinematic process, as
defined by the following equation:

$\hat{p}_{n}(t)=\mathcal{M}_{p}(\vec{\beta},\vec{\theta}(t))$ (1)

In this equation, $\vec{\beta}$ and $\vec{\theta}$ represent the scale and
pose parameters, respectively. The term $\hat{p}_{n}\in\mathbb{R}^{87\times
3}$ denotes the 3D locations of the 87 MoVi keypoints relative to the body
model, which is oriented within a global frame ${n}$.

Additionally, the forward kinematic process calculates the orientations of
each body segment $b_{i}$ in the global coordinate system:

$\hat{R}_{nb}(t)=\mathcal{M}_{R}(\vec{\beta},\vec{\theta}(t))$ (2)

In this model, the orientation of an arbitrary body segment
$R_{nb_{i}}\in\mathbb{R}^{3\times 3}$ is represented within a matrix
$R_{nb}\in\mathbb{R}^{22\times 3\times 3}$, encapsulating the orientations of
all segments.

###  II-G Fusion Architecture

Figure 2: Our method jointly optimizes an implicit function learning the
trajectory of an individual recording and $\vec{\beta}$ the scaling parameters
of a biomechanical model. The implicit function takes time as an input and
outputs $\vec{\theta}(t)$, the joint pose parameters, and $\hat{R}_{nc}(t)$
the orientation of the camera in a global frame. The pose and scaling
parameters produce joint locations $\hat{p}_{n}(t)$ and joint orientations
$\hat{R}_{nb}(t)$ that are compared to detected keypoints and sensor readings.
Note the sensor calibrations are included in the optimziation. Changing the
keypoint reference frame using $\hat{R}_{nc}$, as in Eqs. 8 and 9, is not
explicitly shown.

A diagram of this architecture is shown in Figure 2. Based on recent work
employing an implicit function triangulating 2D markers and biomechanical
poses [5, 6], we use an implicit function $f_{\phi}$ to optimize a mapping
from recording time $t$ to poses $\vec{\theta}(t)$, implemented as a multi-
layer perceptron:

$f_{\phi}:t\rightarrow(\hat{\vec{\theta}},\hat{R}_{nc})$ (3)

The pose vector $\hat{\vec{\theta}}$ is passed through the forward kinematic
model 1,2. The implicit function also outputs $\hat{R}_{nc}(t)$ representing
the smartphone camera’s orientation relative to the global north frame
$\\{n\\}$ at time $t$.

Composing equations 1 or 2 with equation 3 yields a kinematic trajectory,
parameterized by $\phi$, that outputs estimated 3D joint locations, body
segment orientations, and camera orientations at an arbitrary timepoint in the
recording trial.

To obtain the rate of change of model outputs, such as predicted gyroscope
measurements, we differentiate through the forward process represented in
equation 2 to obtain $\hat{\dot{R}}_{nb_{i}}$ and the forward pass in equation
3 for an estimate of $\hat{\dot{R}}_{nc}$.

###  II-H Sensor Calibration Solution

IMUs produce an orientation estimate in their own reference frame
$\\{n^{\prime}\\}$. Here, we are defining $\\{n\\}$ as a horizontal reference
frame facing north, and $\\{n^{\prime}\\}$ captures any slight heading offset
from the IMU zero position due to magnetic disturbance or imperfect sensor
magnetic calibration. We denote the measured IMU orientation of sensor $i$ at
a point in time as $R_{n^{\prime}s_{i}}(t)$ where $\\{s_{i}\\}$ is the sensor
frame for a sensor attached to segment ${i}$. We transform this reading onto
its body segment $b_{i}$ in the global frame with equation 4. If our
calibration and model are optimized correctly, the calibrated sensor and body
model should agree (equation 5).

$\displaystyle
R_{nn^{\prime}}(t)R_{n^{\prime}s_{i}}(t)R_{sb_{i}}\approx\hat{R}_{nb_{i}}(t)$
(4) $\displaystyle
R_{nn^{\prime}}(t)R_{n^{\prime}s}(t)R_{sb_{i}}\approx\mathcal{M}_{R_{i}}(\hat{\vec{\beta}},\hat{\vec{\theta(t)}})$
(5)

with $R_{sb_{i}}$ being the rotation from the sensor frame to the body segment
frame and $R_{nn^{\prime}}$ a rotation from the global frame to the IMU’s home
frame. Critical to our approach, we include $R_{nn^{\prime}}$ and $R_{sb_{i}}$
as learned parameters in our optimization.

While the orientation estimation algorithm is designed to minimize drift [14],
the frame $\\{n^{\prime}\\}$ may still drift over time. We account for this by
representing $R_{nn^{\prime}}(t)$ as a piecewise spherical linear
interpolation (SLERP) [31] between three optimized quaternions
$q_{nn^{\prime}}(0)$, $q_{nn^{\prime}}(T_{mid})$ and $q_{nn^{\prime}}(T)$.

The sensor gyroscope reading $\omega_{s_{i}}$ for a sensor on segment $i$
depends on the body segment alignment. We differentiate through our implicit
function and 2 with respect to time using autograd to obtain
$\hat{\dot{R}}_{nb_{i}}$. Using the sensor calibration, $R_{sb_{i}}$, we
predict sensor $i$’s gyroscope reading as:

$R_{sb_{i}}(\hat{R}_{nb_{i}}^{-1}\hat{\dot{R}}_{nb_{i}})=[\hat{\omega}_{s_{i}}]$
(6)

where $[\hat{\omega}_{s}]$ is the skew-symmetric representation of the
gyroscope reading. For a derivation of this equation refer to [26]. Similarly,
we predict the phone gyroscope with:

$(\hat{R}_{nc}^{-1}\hat{\dot{R}}_{nc})=[\hat{\omega}_{c}].$ (7)

###  II-I Loss functions and implementation

A forward pass of our model produces estimates we denote as $\hat{p}_{n}(t)$,
$\hat{R}_{nb}(t)$, $\hat{R}_{nc}(t)$, $\hat{\omega}_{s}(t)$, and
$\hat{\omega}_{c}(t)$. We optimize the parameters of the neural network
$\phi$, the sensor calibrations $R_{sb}$, $R_{nn^{\prime}}$ and the body
scaling parameters $\vec{\beta}$ using the following losses. To reduce
notational clutter, we define our losses at a single point in time, but all of
these terms are averaged over time except for the sensor calibrations.

####  II-I1 Keypoint loss

We obtain a pure video estimate of the 3D keypoint locations $p_{c}$ using
MeTRAbs-ACAE [7]. These keypoints are in the camera reference frame. We can
rotate them into the global frame using $p_{n}=\hat{R}_{nc}p_{c}$ and define a
loss function on the euclidean distance between keypoints from the model
$\hat{p}_{n}$ with the video keypoints as

$\mathcal{L}_{keypoint}(\phi,\vec{\beta})=\frac{1}{J}\sum_{j\in
J}c(j)g(\|\hat{p}_{n}-\hat{R}_{nc}p_{c}\|_{2})$ (8)

where $c(j)\in[0,1]$ is the confidence score for joint keypoint $j$ which also
varies with time and $g(\cdot)$ is a Huber loss which is quadratic within 100
cm and linear after, which was necessary for stabilizing early training. This
loss is computed between the 3D keypoints set with their mean translation
removed.

####  II-I2 Reprojection loss

MeTRAbs-ACAE also produces 2D keypoints in the image frame
$\vec{x}\in\mathbb{R}^{87\times 2}$. The 3D keypoints produced by our model
$\hat{p}_{n}$ can be rotated into the camera frame using $\hat{R}_{nc}$ and
projected through a camera function $\Pi$ with the calibrated camera
intrinsics, to compute the error with detected 2D keypoints:

$\mathcal{L}_{reprojection}(\phi,\vec{\beta})=\frac{1}{J}\sum_{j\in
J}c(j)g(\|\Pi(\hat{R}_{nc}^{-1}\hat{p}_{n})-\vec{x}\|_{2})$ (9)

In this loss, $g(\cdot)$ is a Huber loss which is quadratic within 100 pixels
and linear after.

####  II-I3 Sensor attitude loss

The forward kinematic process produces an estimate for each body segment $i$’s
orientation in the world frame $\hat{R}_{nb_{i}}$. For each segment with a
sensor attached $B_{s}$, We transform the sensor reading $R_{n^{\prime}s_{i}}$
to this frame using equation 4 and measure the angular difference between the
two estimates with the following loss:

$\begin{split}\mathcal{L}_{att}(\phi,\vec{\beta},R_{sb_{i}},R_{n_{i}^{\prime}n})&=\\\
&\frac{1}{|B_{s}|}\sum_{i\in B_{s}}(\frac{\text{arccos
Tr}[\hat{R}_{nb_{i}}R_{nb_{i}}^{-1}]-1}{2})^{2}\end{split}$ (10)

, where $B_{s}$ represents the set of body segments with sensors attached to
them.

####  II-I4 Sensor gyroscope loss

We compare the estimate of the sensor gyroscope obtained from our model
$\hat{\omega}_{s_{i}}$ (equation 6) to the measured sensor gyroscope
measurement $\omega_{s_{i}}$ for each sensor $s_{i}$ using the mean squared
error:

$\mathcal{L}_{gyro\\_sens}(\phi,\vec{\beta},R_{sb_{i}})=\frac{1}{|B_{s}|}\sum_{i\in
B_{s}}\|\hat{\omega}_{s_{i}}-\omega_{s_{i}}\|^{2}_{2}$ (11)

####  II-I5 Camera gyroscope loss

Finally, we include a loss between the predicted camera gyroscope reading
$\hat{\omega}_{c}$ from equation 7 and the measured camera gyroscope reading
$\omega_{c}$ using the mean squared error:

$\mathcal{L}_{gyro\\_phone}(\phi)=\|\hat{\omega}_{c}-\omega_{c}\|^{2}_{2}$
(12)

###  II-J Implementation

Implementation is done using Jax [32], Equinox [33] and Optax. We combine
three different optimizers using Optax’s multi_transform. First, we use an
Adam optimizer with weight decay set to $1e-5$, exponential learning rate
decay starting at $1e-4$ and ending at $1e-5$, and $\beta=(0.9,0.8)$ for the
implicit function. After 10,000 iteration steps, we begin optimizing the
sensor calibrations with an Adam optimizer with a learning rate $1e-5$.
Finally, we optimize the time offsets between the IMU and video data with a
third Adam optimizer with a learning rate $1e-4$, $\beta_{1}=0.85$. We found
that annealing the weights of the sensor losses was crucial for stability. An
optimization for a $\sim$10 second video converges after 20000 optimization
steps and runs at 50 iterations/s on an NVIDIA A6000 GPU. Since our implicit
function $f_{\phi}(t)$ (equation 3) can be evaluated at arbitrary timepoints,
we were able to train our model using the full resolution sampling rate of all
modalities. A training step computes losses over 500 random samples of each
data source.

###  II-K Evaluation

To evaluate the agreement between the camera and sensor system, we calculated
the residuals between our fused model’s prediction and measured values
averaged over the entire recording. This included:

  1. 1.

Keypoint Residual: Mean Euclidean distance between the 3D keypoints detected
from MeTRAbs-ACAE $p_{c}$ and 3D keypoints predicted by our model and rotated
into the camera frame $R_{nc}^{-1}\hat{p}_{n}$.

  2. 2.

Reprojection Residual: Mean distance in pixels between the 2D keypoints
detected from MeTRAbs-ACAE $x$ and the model predictions rotated into the
camera frame and projected through the smartphone camera model
$\Pi(R_{nc}^{-1}\hat{p}_{n}$).

  3. 3.

Phone gyroscope Residual: Mean difference between the phone gyroscope
predicted by our model $\hat{\omega}_{c}$ and gyroscope measured by the
smartphone $\omega_{c}$. Calculated using 12.

  4. 4.

Sensor Gyroscope Residual: Mean difference between the sensor gyroscope
predicted by our model and raw gyroscope readings from the sensors, averaged
across all sensors. Calculated using 11.

  5. 5.

Attitude Residual: Mean difference between attitude predicted by our model
$\hat{R}_{nb_{i}}$ and raw sensor attitude readings $R_{nb_{i}}$ averaged
across all sensors. Calculated using 10.

Next, we synchronized each PBL recording to the MMMC system. As the MMMC and
PBL systems were collected at different sampling rates, we evaluated our model
at the MMMC timepoints for comparison. As sensors were placed above and below
the knee, the knee angle was our primary measure of accuracy between our
fusion solution and multi-camera system. The metrics comparing the fusion and
multi-camera system are:

  1. 1.

Pearson’s correlation coefficient measures the covariation of the flexion
angles calculated from the fusion solution and MMMC system.

  2. 2.

Mean absolute error (MAE) measures the absolute angular error of the flexion
angle between the fusion solution and MMMC system.

  3. 3.

Mean absolute error mean adjusted (MAE-MA) removes the mean value of each
kinematic trace before comparing the comparing the absolute angular error.

Summary statistics are presented as: Median (IQR). To test for significant
differences in these metrics, we used a Wilcoxon signed-rank test, defining
significance as $p<0.05$.

_Occlusion Experiments._ To estimate how well our method tracked kinematics
from the sensors when video is partially occluded, we artificially occluded
the sensorized leg for the middle half of each recording. We allowed our
optimization to proceed as above and compared the same metrics.

Figure 3: Example lower limb kinematics during walking. (A) With no
occlusions, both video and fusion track closely with multi-camera estimates.
In this example, the fusion solution tracks the knee better during heel
strike. (B) The fusion solution continues to track knee kinematics during
occlusion while the video does not.

##  III RESULTS

###  III-A Qualitative Results

Since this system does not require an explicit calibration trial or static
pose, individuals with gait or balance impairments were able to complete this
experiment with relative ease. A recording session consisting of 9+ walking
trials was typically completed in $<$15 minutes, with a majority of the setup
time attributed to EMG placement (data not used in this study). In some
trials, a physical therapist was closely following a participant or even
holding a gait belt, however, the person and keypoint detection done using
[23] and [7] was generally robust to this. While the loss functions were able
to track sensor readings quite closely, the relative weighting/annealing
schedules of the individual loss terms proved to be very important in stable
convergence and accurate results.

TABLE I: Fit residuals (Median (IQR)) across population and fit type |  | Residual
---|---|---
Population | Fit | Keypoint (cm) | Reprojection (pixels) | Phone Gyroscope (deg/s) | Gyroscope (deg/s) | Attitude (deg)
Control | Video | 0.7 (0.08) | 5.35 (1.06) | 1.76 (1.54) | - | -
Fusion | 0.76 (0.08) | 5.52 (0.97) | 1.73 (1.51) | 23.59 (6.17) | 3.68 (1.01)
Stroke | Video | 0.68 (0.09) | 4.72 (0.72) | 2.18 (1.09) | - | -
Fusion | 0.75 (0.07) | 4.93 (0.72) | 2.15 (1.11) | 21.8 (13.6) | 4.68 (1.36)
Prosthetic | Video | 0.69 (0.08) | 4.74 (0.57) | 2.34 (0.54) | - | -
Fusion | 0.78 (0.13) | 4.97 (0.64) | 2.3 (0.54) | 23.03 (6.38) | 4.44 (1.25)

###  III-B Fusion Consistency

Following fusion, residuals errors were quite low (Table I). Our model was
capable of closely following the detected keypoints with an average of $<1$ cm
3D keypoint error across all groups and $\sim 5$ pixels of average
reprojection error. The keypoint residuals were generally lowest in the video-
only fit and greater when the wearable sensors were added to the optimization.
This is not unexpected as the sensors may compensate for incorrectly detected
keypoints. The phone gyroscope estimate tracked within $\sim 2$ deg/s for all
conditions. Figure 1 shows an example of the model disentangling the phone
rotation, supervised with the measured phone gyroscope, and the pelvis
rotation, which, in this case, accurately recovers the global translation
within 6 cm. Sensor gyroscope residuals averaged 21-23 degrees per second. We
note that while this number is quite high, the relative weighting of the
gyroscope error in the optimization was lower than that of other terms.
Attitude residuals were lowest in the control population $\sim 3$ degrees and
higher in the prosthetic and stroke groups.

###  III-C Multi-Camera Comparison

Across all 60 videos recorded of the 38 subjects, video estimates produced a
MAE-MA of 2.42 (1.2) degrees and 3.91 (1.55) degrees for hip and knee,
respectively. These estimates were significantly decreased to 2.16 (0.86)
degrees and 2.9 (1.27) degrees for the hip and knee, respectively. Pearson’s
correlation coefficients were high for video, 0.98 (0.02) for hip and 0.97
(0.06) for knee, with significant increases with the addition of sensors: 0.99
(0.01) for hip and 0.98 (0.02) for knee. This indicates a high temporal
consistency between our estimates and ground truth. A summary of these results
can be found in Table II and example traces in Figure 3. Sensors significantly
decreased knee MAE-MA in all groups. Control populations had the lowest MAE,
MAE-MA and highest correlations.

###  III-D Artificial Video Occlusion

Our results, depicted in Figures 3 and 4, demonstrate the impact of occlusion
on MAE-MA. Specifically, occlusions increased MAE in the video-only approach,
while the integration of sensors mitigated this effect. Fusion fits during
occlusion were not significantly different from video fits. Additionally,
Pearson correlations significantly dropped in the video-only scenario during
occlusion (0.51 (0.41)), but were restored with the use of sensors (0.97
(0.05)).

TABLE II: Comparison of fusion solutions with ground truth (Median (IQR)). Bold for significance vs video. |  |  $\downarrow$MAE (deg) |  $\downarrow$MAE-MA (deg) |  $\uparrow$Pearson’s R
---|---|---|---|---
Population | Fit | Hip | Knee | Hip | Knee | Hip | Knee
Control | Video | 3.49 (3.28) | 4.3 (2.51) | 1.99 (0.7) | 3.35 (1.48) | 0.99 (0.01) | 0.98 (0.02)
Fusion | 4.02 (3.31) | 3.44 (1.73) | 1.94 (0.66) | 2.53 (0.76) | 0.99 (0.0) | 0.99 (0.01)
Stroke | Video | 3.96 (2.35) | 4.83 (2.28) | 2.52 (1.32) | 4.17 (1.57) | 0.98 (0.03) | 0.94 (0.15)
Fusion | 3.96 (1.52) | 4.48 (2.48) | 2.4 (0.96) | 3.06 (0.91) | 0.98 (0.02) | 0.97 (0.08)
Prosthetic | Video | 4.77 (3.04) | 6.08 (3.3) | 3.19 (1.32) | 4.48 (3.4) | 0.97 (0.01) | 0.94 (0.05)
Fusion | 4.64 (2.51) | 5.71 (4.2) | 2.77 (1.52) | 3.65 (1.92) | 0.98 (0.03) | 0.97 (0.05)
Figure 4: Distribution of knee mean adjusted MAE (MAE-MA). Addition of sensors
generally decreases MAE. Artificial occlusion of the sensorized leg increased
MAE in video-only fits, yet fusion fits significantly decreased the error due
to this. Fusion fits during occlusion were not significantly different from
video fits.

##  IV DISCUSSION

While computer vision and wearable sensors offer much promise for
rehabilitation and clinical monitoring, they remain primarily a research tool
that has not translated to standard clinic or home monitoring. Ease of setup,
use, and troubleshooting is fundamental to the adoption of technologies in the
clinic. Clinicians do not have time to utilize a technology unless it is
immediately intuitive and “foolproof” [19]. The same holds for patients using
technology at home for therapy [34]. The objective of this work was to develop
a method for accurate kinematic tracking from a handheld camera/sensor system
with no calibration requirements. We did this by jointly optimizing an
implicit function that fits the pose and camera trajectory as well as the
sensor-body calibration rotations. This approach provides accurate
measurements of lower body kinematics across control, lower limb prosthesis
user, and neurological populations. The accuracy of the knee is within 6
degrees for video alone when comparing to MMMC, with errors less than 4
degrees when removing the bias, which is already quite good. The inclusion of
wearable sensors reduces the measurement error by another degree. Our dataset
contained a large amount of walking data, so we chose to focus on knee
kinematics in this study, however, this method should generalize to other
joints as well.

The results in Table I demonstrate that our algorithm is able to solve for a
solution to Equation 4 that rotates raw IMU data into the body frame. We found
that the phone gyroscope was useful for disentangling the camera vs root
rotation (Figure 1). We were surprised at how high the sensor gyroscope
residuals were, particularly as preliminary experiments with SMPL models
showed lower gyro residuals. We experimented with different
weighting/annealing schedules of the gyroscope loss; however, our model did
not allow for closer fitting without disturbing the kinematic trajectory. An
explanation for this could be that our model did not include varus/valugs
rotation at the knee joint, which could be incorporated in the future.
Alternatively, numerical errors from differentiating through the forward
kinematics could have limited the accuracy (Equations 6,7). All other
residuals we found to be quite low, $<1$ cm and 5 pixels for keypoints, $\sim
2$ deg/s of camera gyroscope error, and $\sim 4$ degrees of attitude error.

Table II shows that our fusion method produces accurate kinematics for our
joint of interest when compared to a multi-camera system [22, 6]. We were
pleased with how small the joint angle error were with a video-only approach
with 3.35, 4.17, and 4.48 degrees of knee MAE-MA for control, stroke and
prosthetic populations, respectively. This indicates that for many
applications, video alone could be sufficient. Despite the high bar set by
video, the addition of the wearable sensors improved kinematic tracking for
all populations by 0.8 - 1.1 degrees. Our synthetic occlusion experiments also
show an additional benefit of sensors, where the fusion allows tracking
through visual occlusions of a leg, and that the sensors recovered the knee
kinematics during this time (Figures 3 and 4).

Previous work has fused sensor and monocular video data; however, using
stationary or calibrated cameras or calibrated IMUs. Work using a moving
camera [35] obtains 12 degrees of angular error across the whole body. This is
not quite a fair comparison to the results presented here as their method
evaluates the mean orientation error of SMPL body segments, which, as they
overparameterize pose space, are likely to not be consistent between
modalities. Our method for evaluating the knee angle does not have this
overparameterization problem so it is not surprising that our method finds
closer agreement. A strength of our method is it allows IMUs to be placed
arbitrarily on arbitrary body segments. Recent work [36] allows for an
arbitrary IMU placement as well and averages 13 degrees of error. OpenCap uses
two calibrated stationary cameras [17] and finds joint orientation errors on
par with our approach (4.1 degrees). OpenSense [18] uses an explicit
calibration step and finds MAE on par with our approach (3-6 degrees) and
lower correlation coefficients (0.60-0.87). Another method, [16] uses three
multi-layer perceptrons to fuse IMU accelerometer, gyroscope and stationary
video for full-body kinematics. While this method achieves results near 4
degrees of angular error, it still requires IMUs to be loosely calibrated to
the body prior the to recording and does not utilize a moving camera. In
summary, we find that our method combines the perks of all these methods:
monocular, moving camera with arbitrarily placed, uncalibrated IMUs to find
accurate, biomechanically grounded kinematics.

While we are excited about the ease of use of this method, it is not without
limitations. First, the inclusion of a loss on the camera gyroscope may make
this method less accessible, however, most smartphones do have a gyroscope and
vision approaches can substitute for raw camera gyroscope values [37]. We
acknowledge that our method only tracks orientation and translation relative
to the smartphone camera. In cases where the experimenter remained stationary
(Figure 1), we saw a mean root translation error of 6 cm. We anticipate future
work using the gait kinematics to inform the trajectory through the world can
provide complementary information to track the person in world coordinates
even in camera-translating conditions [37].

Future work should apply this method to a larger clinical population. In
particular, most lower limb prosthesis users had their device covered with
their pant leg, which likely artificially increased the accuracy in that
population. This work should also be validated on other joint kinematics as
well, including the arm. Another future direction is fitting multiple
trajectories simultaneously and sharing parameters like body scale and sensor
orientation across the trajectories. We use this bilevel optimization approach
[38] in our end-to-end optimization from the MMMC data, and find it is both
accurate and reduces the time to fit multiple trajectories.

##  V CONCLUSIONS

Sensor-body and camera calibration are time-consuming steps that limit the use
of biomechanics and IMUs in routine clinical and at-home settings. We present
a method to fuse handheld smartphone video and uncalibrated IMUs that produces
accurate kinematics while retaining the high-resolution information in IMU
readings. We also find that using an implicit function to reconstruct
biomechanics from monocular video does remarkably well using modern 3D
keypoint detectors and differentiable biomechanical models.

## ACKNOWLEDGMENT

## References

  * [1] Gyrd Thrane, Katharina Stibrant Sunnerhagen and Margit Alt Murphy  “Upper limb kinematics during the first year after stroke: the stroke arm longitudinal study at the University of Gothenburg (SALGOT)”  In _Journal of NeuroEngineering and Rehabilitation_ 17.1, 2020, pp. 76  DOI: 10.1186/s12984-020-00705-2
  * [2] R. Cotton et al.  “Improved Trajectory Reconstruction for Markerless Pose Estimation”  arXiv, 2023  arXiv: http://arxiv.org/abs/2303.02413
  * [3] “Kinematics”  In _Biomechanics and Motor Control of Human Movement_ John Wiley & Sons, Ltd, 2009, pp. 45–81  DOI: https://doi.org/10.1002/9780470549148.ch3
  * [4] Robert M. Kanko, Elise Laende, W. Selbie and Kevin J. Deluzio  “Inter-session repeatability of markerless motion capture gait kinematics” Pages: 2020.06.23.155358 Section: New Results  bioRxiv, 2020  DOI: 10.1101/2020.06.23.155358
  * [5] R. Cotton et al.  “Optimizing Trajectories and Inverse Kinematics for Biomechanical Analysis of Markerless Motion Capture Data”  In _2023 International Conference on Rehabilitation Robotics (ICORR)_ Singapore, Singapore: IEEE, 2023, pp. 1–6  DOI: 10.1109/ICORR58425.2023.10304683
  * [6] R. Cotton  “Differentiable Biomechanics Unlocks Opportunities for Markerless Motion Capture”  arXiv, 2024  arXiv: http://arxiv.org/abs/2402.17192
  * [7] István Sárándi, Alexander Hermans and Bastian Leibe  “Learning 3D Human Pose Estimation from Dozens of Datasets using a Geometry-Aware Autoencoder to Bridge Between Skeleton Formats”  arXiv, 2022  arXiv: http://arxiv.org/abs/2212.14474
  * [8] Matthew Loper et al.  “SMPL: a skinned multi-person linear model”  In _ACM Transactions on Graphics_ 34.6, 2015, pp. 1–16  DOI: 10.1145/2816795.2818013
  * [9] Angjoo Kanazawa, Michael J. Black, David W. Jacobs and Jitendra Malik  “End-to-End Recovery of Human Shape and Pose”, 2018, pp. 7122–7131  URL: https://openaccess.thecvf.com/content_cvpr_2018/html/Kanazawa_End-to-End_Recovery_of_CVPR_2018_paper.html
  * [10] Marilyn Keller et al.  “From Skin to Skeleton: Towards Biomechanically Accurate 3D Digital Humans”  In _ACM Transactions on Graphics_ 42.6, 2023, pp. 253:1–253:12  DOI: 10.1145/3618381
  * [11] Ge Wu et al.  “ISB recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion—part I: ankle, hip, and spine”  In _Journal of Biomechanics_ 35.4, 2002, pp. 543–548  DOI: 10.1016/S0021-9290(01)00222-6
  * [12] Yann Desmarais, Denis Mottet, Pierre Slangen and Philippe Montesinos  “A review of 3D human pose estimation algorithms for markerless motion capture”  In _Computer Vision and Image Understanding_ 212, 2021, pp. 103275  DOI: 10.1016/j.cviu.2021.103275
  * [13] Tong Li and Haoyong Yu  “Visual–Inertial Fusion-Based Human Pose Estimation: A Review”  In _IEEE Transactions on Instrumentation and Measurement_ 72, 2023, pp. 1–16  DOI: 10.1109/TIM.2023.3286000
  * [14] R. Cotton and John Rogers  “Wearable Monitoring of Joint Angle and Muscle Activity”  In _IEEE … International Conference on Rehabilitation Robotics: [proceedings]_ 2019, 2019, pp. 258–263  DOI: 10.1109/ICORR.2019.8779538
  * [15] Robert Mahony, Tarek Hamel and Jean-Michel Pflimlin  “Nonlinear Complementary Filters on the Special Orthogonal Group”  In _IEEE Transactions on Automatic Control_ 53.5, 2008, pp. 1203–1218  DOI: 10.1109/TAC.2008.923738
  * [16] Soyong Shin, Zhixiong Li and Eni Halilaj  “Markerless Motion Tracking with Noisy Video and IMU Data”  In _IEEE Transactions on Biomedical Engineering_ , 2023, pp. 1–11  DOI: 10.1109/TBME.2023.3275775
  * [17] Scott D. Uhlrich et al.  “OpenCap: Human movement dynamics from smartphone videos” Publisher: Public Library of Science  In _PLOS Computational Biology_ 19.10, 2023, pp. e1011462  DOI: 10.1371/journal.pcbi.1011462
  * [18] Mazen Al Borno et al.  “OpenSense: An open-source toolbox for inertial-measurement-unit-based measurement of lower extremity kinematics over long durations”  In _Journal of NeuroEngineering and Rehabilitation_ 19.1, 2022, pp. 22  DOI: 10.1186/s12984-022-01001-x
  * [19] Jessie Mitchell, Camila Shirota and Kelly Clanchy  “Factors that influence the adoption of rehabilitation technologies: a multi-disciplinary qualitative exploration”  In _Journal of NeuroEngineering and Rehabilitation_ 20.1, 2023, pp. 80  DOI: 10.1186/s12984-023-01194-9
  * [20] Anthony Cimorelli, Ankit Patel, Tasos Karakostas and R. Cotton  “Validation of portable in-clinic video-based gait analysis for prosthesis users” Publisher: Nature Publishing Group  In _Scientific Reports_ 14.1, 2024, pp. 3840  DOI: 10.1038/s41598-024-53217-7
  * [21] Anthony Cimorelli, Ankit Patel, Tasos Karakostas and R. Cotton  “Portable in-clinic video-based gait analysis: validation study on prosthetic users”, 2022  DOI: 10.1101/2022.11.10.22282089
  * [22] R. Cotton et al.  “Markerless Motion Capture and Biomechanical Analysis Pipeline”  arXiv, 2023  arXiv: http://arxiv.org/abs/2303.10654
  * [23] R. Cotton  “PosePipe: Open-Source Human Pose Estimation Pipeline for Clinical Research”  arXiv, 2022  arXiv: http://arxiv.org/abs/2203.08792
  * [24] Nicolai Wojke, Alex Bewley and Dietrich Paulus  “Simple online and realtime tracking with a deep association metric”  In _2017 IEEE International Conference on Image Processing (ICIP)_ Beijing: IEEE, 2017, pp. 3645–3649  DOI: 10.1109/ICIP.2017.8296962
  * [25] Saeed Ghorbani et al.  “MoVi: A large multi-purpose human motion and video dataset” Publisher: Public Library of Science  In _PLOS ONE_ 16.6, 2021, pp. e0253157  DOI: 10.1371/journal.pone.0253157
  * [26] Kevin M. Lynch and Frank C. Park  “Modern robotics: mechanics, planning, and control” OCLC: ocn983881868  Cambridge, UK: Cambridge University Press, 2017
  * [27] Emanuel Todorov, Tom Erez and Yuval Tassa  “MuJoCo: A physics engine for model-based control” ISSN: 2153-0866  In _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2012, pp. 5026–5033  DOI: 10.1109/IROS.2012.6386109
  * [28] Vittorio Caggiano et al.  “MyoSuite – A Contact-Rich Simulation Suite for Musculoskeletal Motor Control”, 2022  DOI: 10.48550/arXiv.2205.13600
  * [29] Firas Al-Hafez, Guoping Zhao, Jan Peters and Davide Tateo  “LocoMuJoCo: A Comprehensive Imitation Learning Benchmark for Locomotion”  arXiv, 2023  DOI: 10.48550/arXiv.2311.02496
  * [30] Samuel R. Hamner, Ajay Seth and Scott L. Delp  “Muscle contributions to propulsion and support during running”  In _Journal of Biomechanics_ 43.14, 2010, pp. 2709–2716  DOI: 10.1016/j.jbiomech.2010.06.025
  * [31] Ken Shoemake  “Animating rotation with quaternion curves”  In _Proceedings of the 12th annual conference on Computer graphics and interactive techniques_ , SIGGRAPH ’85  New York, NY, USA: Association for Computing Machinery, 1985, pp. 245–254  DOI: 10.1145/325334.325242
  * [32] James Bradbury  “JAX: composable transformations of Python+NumPy programs” original-date: 2018-10-25T21:25:02Z  Google, 2024  URL: https://github.com/google/jax
  * [33] Patrick Kidger and Cristian Garcia  “Equinox: neural networks in JAX via callable PyTrees and filtered transformations”  arXiv, 2021  arXiv: http://arxiv.org/abs/2111.00254
  * [34] Yu Chen et al.  “Home-based technologies for stroke rehabilitation: A systematic review”  In _International Journal of Medical Informatics_ 123, 2019, pp. 11–22  DOI: 10.1016/j.ijmedinf.2018.12.001
  * [35] Timo Marcard et al.  “Recovering Accurate 3D Human Pose in the Wild Using IMUs and a Moving Camera” Series Title: Lecture Notes in Computer Science  In _Computer Vision – ECCV 2018_ 11214 Cham: Springer International Publishing, 2018, pp. 614–631  DOI: 10.1007/978-3-030-01249-6˙37
  * [36] Tom Van Wouwe et al.  “Diffusion Inertial Poser: Human Motion Reconstruction from Arbitrary Sparse IMU Configurations”  arXiv, 2023  arXiv: http://arxiv.org/abs/2308.16682
  * [37] Soyong Shin, Juyong Kim, Eni Halilaj and Michael J. Black  “WHAM: Reconstructing World-grounded Humans with Accurate 3D Motion”  arXiv, 2023  arXiv: http://arxiv.org/abs/2312.07531
  * [38] Keenon Werling et al.  “AddBiomechanics: Automating model scaling, inverse kinematics, and inverse dynamics from human motion data through sequential optimization”, 2023


2.5cm2.5cm2cm2cm ††affiliationtext: German Institute for Economic Research
(DIW Berlin), Mohrenstr. 58, D-10117 Berlin, Germany<EMAIL_ADDRESS>††affiliationtext: Technische Universität Berlin, D-10623 Berlin, Germany,
<EMAIL_ADDRESS> ††affiliationtext: German Institute for Economic
Research (DIW Berlin), Mohrenstr. 58, D-10117 Berlin, Germany<EMAIL_ADDRESS>††affiliationtext: German Institute for Economic Research (DIW Berlin),
Mohrenstr. 58, D-10117 Berlin, Germany<EMAIL_ADDRESS>
# An open tool for creating battery-electric vehicle time series from
empirical data - emobpy

Carlos Gaete-Morales  Hendrik Kramer  Wolf-Peter Schill  Alexander Zerrahn

###### Abstract

There is substantial research interest in how future fleets of battery-
electric vehicles will interact with the power sector. To this end, various
types of energy models depend on meaningful input parameters, in particular
time series of vehicle mobility, driving electricity consumption, grid
availability, or grid electricity demand. As the availability of such data is
highly limited, we introduce the open-source tool emobpy. Based on mobility
statistics, physical properties of vehicles, and other customizable
assumptions, it derives time series data that can readily be used in a wide
range of model applications. For an illustration, we create and characterize
$200$ battery-electric vehicle profiles for Germany. Depending on the hour of
the day, a fleet of one million vehicles has a median grid availability
between $5$ and $7$ gigawatts, as vehicles are parking most of the time. Four
exemplary grid electricity demand time series illustrate the smoothing effect
of balanced charging strategies.

##  1 Introduction

We introduce emobpy. It is an open-source python-based tool that creates
profiles of battery-electric vehicles (BEV), based on empirical mobility
statistics and customizable assumptions. We additionally provide a first
application of the tool and create vehicle profiles based on representative
German mobility data. An emobpy profile consists of four time series: (i)
vehicle mobility containing the vehicle’s location and distance travelled,
(ii) driving electricity consumption, specifying how much electricity is taken
from the battery for driving; (iii) BEV grid availability, providing
information whether and with which power rating a BEV is connected to the
electricity grid at a certain point in time; and (iv) BEV grid electricity
demand, specifying the actual charging electricity drawn from the grid, based
on different charging strategies.

Such profiles are core input data for a wide range of model applications in
energy, environmental, and economic studies on BEV. Technology developments as
well as energy and climate policy measures drive the deployment of BEV in many
countries [eia_2019]. Growing BEV fleets can have substantial impacts on the
power sector. They increase the electric load, but may also provide temporal
flexibility for integrating variable renewable energy sources and contribute
to decarbonizing transportation [ipcc_2018]. Many model-based analyses
investigate potential power sector interactions of future BEV fleets
[Daina_2017a, mwasilu_2014, richardson_2013, Taljegard_2019] and thus depend
on a meaningful representation of electric vehicles’ mobility patterns.

Yet such data are often not publicly available. In general, empirical data are
scarce because BEV fleets are still small in most countries. And if respective
time series are available, they are often specific to the conditions in which
the data was collected and subject to data protection provisions. Past
approaches make either stylized coarse assumptions [Kempton_2005], derive data
from mobility statistics, but lack documentation, transparency or
reproducibility [Fischer_2019, Kiviluoma_2011, Muratori_2018, robinson_2013,
schaeuble_2017, schill_2015, Heinz_2018], or are idiosyncratic with respect to
geographic characteristics or assumed driver behavior [Chen_2018,
robinson_2013, Wolinetz_2018].

Figure 1: Inputs and outputs for emobpy and sequence of generating four types
of time series. The boxes on the left-hand side show customizable input
assumptions, the boxes on the right-hand side indicate the four types of time
series.

Following [Daina_2017a], we argue that new models are needed to derive
relevant time series in a transparent and flexible way. As a first step in
this direction, the tools Vencopy [Wulff2020] and RAMP-mobility [LOMBARDI2019]
recently emerged. To further fill this gap, we developed emobpy. Our tool
takes empirical mobility statistics, physical properties of vehicles, and
customizable assumptions as inputs and delivers BEV profiles as output. Figure
1 gives a stylized account. We first discuss the outputs, then the required
inputs.

Four output time series constitute one BEV profile. These profiles have a
customizable length and resolution. A handy format for many applications is
all hours of one year. But also other formats are possible by discretion of
the researcher. Likewise, the researcher can choose how many profiles she
wants to create.

The time series of vehicle mobility (i) contains the location of the vehicle
at each time step and the time steps during which the vehicle is driving with
information of the distance traveled. The driving electricity consumption time
series (ii) provides information on how much electricity the vehicle consumes
for driving in each time step. The time series of (iii) grid availability
provides information whether a vehicle is connected to the electricity grid in
a time step and if so, with what power rating for charging or discharging. The
time series of grid electricity demand (iv) provides information on how much
electricity a vehicle demands from the electricity grid in a time step. Time
series (i), (ii) and (iii) are core inputs for models that endogenously
determine the timing of charging (and, potentially, discharging to the grid);
the time series (iv) are core inputs for models that do not endogenously
determine the grid interactions of BEV, but use exogenous input data for this.

The required input data for the time series of vehicle mobility (i) are the
relative frequencies of different driver types, e.g., commuters, of the number
of trips per day, of the destination, distance and duration of trips, and of
the departure hours. Such information can often be derived from national
mobility statistics. If required or desired, a researcher can also make up own
assumptions or resort to the pre-set values from German mobility statistics.
emobpy makes sure that the resulting time series are feasible and consistent.
To this end, a minimum and maximum number of hours at specific locations can
specified, and it is assumed that the last trip of a day heads home. With a
Monte Carlo approach, emobpy ensures variability across profiles.

Based on the vehicle mobility time series, the driving electricity consumption
(ii) time series is derived. This requires further input data, such as
information on nominal motor power, curb weight, drag coefficient, and
dimensions, which the tool includes for several current BEV models. Ambient
temperature is also a significant parameter that affects the consumption of
BEV [Brown_2018, Fischer_2019]. For that reason, emobpy is endowed with a
database of hourly temperature for European countries with a registry of the
last 17 years. Additionally, the vehicle cabin insulation characteristics are
required; this data is not widely available and thus assumed independently of
the BEV models database. Driving cycles are also important input parameters
that are used to simulate every individual trip. The model includes two
driving cycles, Worldwide Harmonized Light Vehicles Test Cycle (WLTC) and
Environmental Protection Agency (EPA). This input data is already provided
within the tool, and the user can select a particular BEV model, country
weather and driving cycle. Additionally, emobpy also allows providing user-
defined custom data.

The required input data for the grid availability time series (iii) is the
driving electricity consumption time series (ii). Further, data or assumptions
on the power rating of charging stations at different generic locations as
well as their availability probabilities are needed. Variability across
profiles is, again, introduced through a Monte Carlo approach, while emobpy
makes sure that the time series (iii) is consistent within each profile .

The required input data for the grid electricity demand time series (iv)
includes the created time series on driving electricity consumption (ii) and
grid availability (iii). Additionally, users can choose a charging strategy,
such as immediate full charging or night-time charging, or make customary
assumptions.

##  2 Results

###  2.1 Application to Germany: parameterization and setup

For a first application of emobpy, we draw on the comprehensive German
mobility survey Mobilität in Deutschland [MID_2017, Mobility in Germany,]. The
survey features mobility data relating to different types of households,
vehicles, individuals, and trips. In this application, we make three general
assumptions: first, we assume that individuals with access to a vehicle carry
out all their trips with the same vehicle; second, we assume that future BEV
drivers have similar mobility patterns as current conventional drivers covered
by the underlying mobility statistics; and third, for simplicity and
tractability, we assume that there are only four BEV models: Hyundai Kona,
Renault Zoe, Tesla Model 3 and Volkswagen ID.3. These models had the largest
market shares in Germany by the time of writing. Again, all of these pre-set
assumptions can easily be modified in emobpy.

We generate $50$ profiles for each BEV model, i.e., $200$ BEV profiles
overall, each consisting of four time series. We focus on two types of
drivers: commuters ($62$% of all drivers) and non-commuters ($38$% of all
drivers). For commuters, we further differentiate between full-time and part-
time employees, with a split of $78$ to $22$% [OCDE_2018]. We exclude
commuting students, apprentices, and trainees, who represent only a small
share of all commuters in the initial dataset. The amount of trips per day
varies between $0$ \- $5$ with different probabilities for weekdays and
weekend days (Table 2.1).

Table 1: Probability distributions (given in %) for the amount of trips per
day by days of the week


where $r,s$ are of the same kind, and that $\lambda(qr)=(\lambda q)r$ for any
scalar $\lambda$ and quantities $q,r$. Although this is technically not part
of the definition of a space of quantities, Raposo also assumes that if $q$
and $r$ are non-zero quantities then $qr$ is a non-zero quantity.

Raposo’s theory of quantity spaces (supplemented with a ’no zero divisors’
condition) is complete and free from anomalies. Compared to the theory
presented in this article, it contains some redundant elements, since there
are more primitive notions. However, the two theories have been shown to be
completely equivalent [26]. This lends credence to both theories, since they
were developed independently.

The table below contains a simplified comparison of some aspects of six of the
approaches to quantity calculus reviewed above.

_Aspects_ | _Authors_
---|---
|  Drobot (1953) |  Quade (1961) |  Whitney (1968) |  Carlson (1979) |  Jonsson (2014) |  Raposo (2016)
Addition of quantities |  Derived |  Primitive |  Derived |  None |  Derived |  Primitive
Scalar product |  |  |  |  |  |
of quantities |  Derived |  Primitive |  Derived |  Derived |  Primitive |  Primitive
Product of quantities |  Primitive |  Primitive |  Primitive |  Derived |  Primitive |  Primitive
Product of dimensions |  Derived |  Derived |  Derived |  (Primitive) |  Derived |  Primitive
Exponent of quantities |  Primitive |  Derived |  Primitive |  Derived |  Derived |  Derived
Exponent of dimensions |  Derived |  Derived |  Derived |  (Primitive) |  Derived |  Derived
Ring of exponents |  $\mathbb{R}$ |  $\mathbb{Z}$ |  $\mathbb{Q}$ or $\mathbb{R}$ |  $\mathbb{Q}$ |  $\mathbb{Z}$ |  $\mathbb{Z}$

Note: For Carlson, we consider products and exponents of pre-units instead of
dimensions.

###  B.2. Rules of quantity calculus

In [2], de Boer lists some fundamental rules for calculation with quantities,
drawn from the literature on quantity calculus. All these rules can be derived
from the theory of quantity spaces. The list below includes references to de
Boer’s rules and relevant definitions or results from this article.

  1. (1)

(A2.1; Definition 2.1, Corollary 3.5). A quantity can be multiplied by a
quantity as in a monoid, and the non-zero quantities form an abelian group
under multiplication.

  2. (2)

(A2.2; Definition 2.1). A quantity $q$ can be multiplied by a number
$\lambda$, and we have the identities $1\cdot q=q$ and
$\alpha\cdot\mathopen{}\mathclose{{}\left(\beta\cdot
q}\right)=\mathopen{}\mathclose{{}\left(\alpha\beta}\right)\cdot q$.

  3. (3)

(A2.2; Definition 2.1). Multiplication of quantities by numbers and by
quantities are related by
$\lambda\cdot\mathopen{}\mathclose{{}\left(pq}\right)=\mathopen{}\mathclose{{}\left(\lambda\cdot
p}\right)q=p\mathopen{}\mathclose{{}\left(\lambda\cdot q}\right)$.

  4. (4)

(A3.1, A6.1; Definition 2.6, Proposition 2.7). A system of quantities can be
partitioned into equivalence classes of quantities of the same kind. (de Boer
makes a distinction between such equivalence classes and dimensions, but in
the present theory these two notions coincide.)

  5. (5)

(A3.2; Proposition 3.9). Quantities of the same kind can be added and form an
abelian group under addition.

  6. (6)

(A3.3; Proposition 3.9). If $p$ and $q$ are of the same kind then
$\lambda\cdot\mathopen{}\mathclose{{}\left(p+q}\right)=\lambda\cdot
p+\lambda\cdot q$ and $\mathopen{}\mathclose{{}\left(\alpha+\beta}\right)\cdot
q=\alpha\cdot q+\beta\cdot q$.

  7. (7)

(A3.4; Remark 3.17.) All ”dimensionless” quantities are of the same kind.

  8. (8)

(A4.1, A6.2; Proposition 3.20). Kinds of quantities can be multiplied, forming
a finitely generated free abelian group under multiplication.

  9. (9)

(A4.2, A6.2; Proposition 2.13, Definition 2.14). If $q$ is a quantity of kind
$K$ and $q^{\prime}$ a quantity of kind $K^{\prime}$ then $qq^{\prime}$ is a
quantity of kind $KK^{\prime}$.

  10. (10)

(A5.1; Proposition 3.8). For every kind of quantities and every quantity $q$
of this kind there is a quantity $u$ of the same kind such that $q=\mu\cdot
u$, where $\mu$ is a uniquely determined number. Such a quantity $u$ is called
a unit.

  11. (11)

(A5.1, A5.2; Proposition 3.10). It is possible to select exactly one unit from
each kind of quantities in such a way that that if $u$ is the unit of kind $K$
and $u^{\prime}$ is the unit of kind $K^{\prime}$ then $uu^{\prime}$ is the
unit of kind $KK^{\prime}$. A set of units satisfying this condition is said
to be coherent.

We can derive some more fundamental rules not considered by de Boer in [2].

  1. (i)

(Corollary 2.11). $\lambda\cdot q$ is a quantity of the same kind as $q$.

  2. (ii)

(Proposition 3.11). If $q$ is a quantity and $r,s$ are quantities of the same
kind then $q\mathopen{}\mathclose{{}\left(r+s}\right)=qr+qs$.

  3. (iii)

(Corollary 3.5.) If $q$ and $r$ are non-zero quantities then $qr$ is non-zero.

###  B.3. On applications of the theory of quantity spaces

At the heart of theoretical metrology is the relationship between measures and
what is measured – quantities. Quantity space theory can clarify this
relationship by elucidating the nature of quantities and systems of
quantities.

In the _International Vocabulary of Metrology_ from 2012 (VIM3) [32], one
sense of ”quantity” (1.1) is a generic one, corresponding to _Grössenart_ or
kind of quantity, while a ”quantity value” (1.19) represents a particular
_Grösse_ or a quantity as defined here. Unit quantities are called
”(measurement) units” (1.9) in VIM3. A set of ”base quantities” (1.4) is in
effect a set $\mathsf{E}$ of selected kinds of quantities, or equivalently
dimensions, which is a basis for some $Q/{\sim}$; the elements of a
corresponding basis
$E=\mathopen{}\mathclose{{}\left\\{e_{1},\ldots,e_{n}}\right\\}$ for $Q$ are
”base units” (1.10). Further, in VIM3 a ”derived unit” (1.11) is some $u\in
Q$, other than a base unit, with an expansion
$u=\mu\cdot\prod_{i=1}^{n}e_{i}^{k_{i}}$, where $\mu\neq 0$, while a ”coherent
derived unit” (1.12) is a derived unit $v$ with an expansion
$v=1\cdot\prod_{i=1}^{n}e_{i}^{k_{i}}$. Following Fourier [9], VIM3 defines a
”quantity dimension” (1.7) of a quantity $q\in Q$ in terms of an integer tuple
$\mathopen{}\mathclose{{}\left(k_{1},\ldots,k_{n}}\right)$ describing how $q$
is expressed as
$\mu_{E}\mathopen{}\mathclose{{}\left(q}\right)\cdot\prod_{i=1}^{n}e_{i}^{k_{i}}$,
where $\mathopen{}\mathclose{{}\left\\{e_{1},\ldots,e{}_{n}}\right\\}$ is a
quantity-space basis for $Q$.

It would be of interest to fully analyze VIM3 (or the upcoming VIM4) in the
light of the theory of quantity spaces. The distinction between concrete and
abstract quantities [15, p. 8–12] should be taken into account in this
connection,

The relationship between measures and the quantities that they represent is
fundamental also in dimensional analysis, which is based on a principle of
covariance: a relation between scalars representing a relation between
quantities relative to a system of unit quantities must continue to hold when
that system is changed in a legitimate way, although individual scalars may
change as unit quantities change. Reference [16] presents an approach to
dimensional analysis explicitly based on this principle and expressed in terms
of quantities, dimensions and quantity functions rather than scalars, units
and scalar functions.

A measure of a quantity is usually assumed to be a real number, but in a
quantity space over $K$ a measure is an element of $K$, where $K$ is any
field. For example, a measure can be a complex number. This makes the present
theory of quantity spaces well suited for applications to problems of quantum
physics.