Simple Degree-of-Freedom Modeling of the Random Fluctuation Arising in Human–Bicycle Balance

Abstract

In this study, we propose a new simple degree-of-freedom fluctuation model that accurately reproduces the probability density functions (PDFs) of human–bicycle balance motions as simply as possible. First, we measure the time series of the roll angular displacement and velocity of human–bicycle balance motions and construct their PDFs. Next, using these PDFs as training data, we identify the model parameters by means of particle swarm optimization; in particular, we minimize the Kolmogorov–Smirnov distance between the human PDFs from the participants and the PDFs simulated by our model. The resulting PDF fitnesses were over $98.7\%$ for all participants, indicating that our simulated PDFs were in close agreement with human PDFs. Furthermore, the Kolmogorov–Smirnov statistical hypothesis testing was applied to the resulting human–bicycle fluctuation model, showing that the measured time responses were much better supported by our model than the Gaussian distribution.

1. Introduction

Bicycles provide a useful means of short-distance transportation, and their utilization is expected to contribute to building a healthy and environmentally friendly society [1]. However, at least on Japanese roads, bicycle transportation is not always necessarily safe due to collisions with automobiles. The Japan National Police Agency reported that over 83% of bicycle accidents in Japan in the last ten years have involved automobiles [2]. To avoid such accidents, autonomous vehicle technology will play an important role; if it can predict bicycle motions, the resulting self-driving cars may reduce such accidents. For this purpose, accurate simulation models of bicycle motions are required, and they should be provided as simply as possible for the potential use of electronic control units in self-driving cars.

Bicycle motions with human riders in traffic seem to be broadly classified into two types: voluntary and involuntary. The former comprises purposeful motions such as right and left turning at a street intersection. Google [3] reported that their sensors can detect a cyclists’ hand signals as an indication of an intention to make a turn or shift over. The latter comprises unconscious motions such as human fluctuated balance motions, which have already been found universally in human quiet standing [4,5,6], human stick balancing [7,8], human visuomotor tracking [9,10,11], and so on.

In this study, we propose a simple stochastic model that allows us to simulate the latter type of bicycle motion, i.e., involuntary fluctuated human–bicycle balance motions. To this end, we have conducted an experiment in which each human participant rides a bicycle on bicycle-trainer rollers, allowing it to move without rolling or yawing constraints. During this experiment, we measured the bicycle’s rolling motion using a three-dimensional motion sensor attached to it. The measured time-series of the rolling motion are characterized by their probability density functions (PDFs). Next, we designed our proposed human–bicycle model as a simple pendulum mechanism controlled by our human controller model, which was successfully used in our previous study [11] to simulate random human fluctuations during a visuomotor tracking task. We then identified the model parameters based on the measured PDFs as training data, using particle swarm optimization (PSO) to minimize the Kolmogorov–Smirnov (KS) distance between the measured PDFs and those simulated by our proposed human–bicycle model. The results show that our proposed model successfully reproduces the measured PDFs with fitnesses of over $98.7\%$. Furthermore, we conducted a statistical hypothesis test called the KS test [12,13] on our results to check their stochastic reliability, showing that the measured time series were much better supported by our model than the Gaussian distribution.

Our stochastic human-modeling approach sharply contrasts with other studies in the fields of autonomous or unmanned bicycle-control systems [14,15,16,17,18,19,20,21,22] because their models have been deterministic and not designed to have randomly fluctuating terms. Our approach also differs from Google’s study on voluntary bicycle motions, as mentioned above [3]. Although there has been another study on the stochastic modeling of bicycle fluctuated motion [23], it addressed large-scale bicycle running paths, unlike our study, which deal with small fluctuations.

The rest of the paper is structured as follows: Section 2 describes the experimental test of human–bicycle balance motions. In Section 3, our proposed stochastic model of these motions is presented. In Section 4 and Section 5, we describe the method of parameter identification and the identification results are presented with the KS testing results. Section 6 concludes our study.

2. Human–Bicycle Balance Experiment

2.1. Experimental Setup and Procedure

Figure 1 shows a photograph of our experimental device, a participant, and an experimenter. The experimental device consists of four units: a bicycle (BE-ELL03, Panasonic, Osaka, Japan), a set of bicycle-trainer rollers (E-MOTION, Elite, Italy), a three-dimensional motion sensor (CSM-MG100, Tokyo Aircraft Instrument, Tokyo, Japan), and a computer.

The experimental participants were eight healthy males in their early twenties. They were first instructed on the operation of the experimental device, the number of trials, and the duration of each trial. The experiment was performed according to the principles of the Declaration of Helsinki, and informed consent was obtained.

In each trial, the participant rode a bicycle on trainer rollers for significantly more than 180 s at a speed of about 15 km/h. Several practice trials were performed prior to measurement. The air pressure of the tires was set to 300 kPa.

Note that, in our setup, the bicycle is constrained on the stationary trainer. This may cause our measured motion to differ from free motions on the road. Such potential discrepancy, however, did permit a three-minute safety trial without external disturbances and was minimized by using a roller-type trainer known to have fewer constraints than other types of trainers.

2.2. Experimental Data

During each trial, the time series of the measurement vector,

$$\mathit{x}\left(t\right)={(x_1\left(t\right),x_2\left(t\right))}^T:={(\theta \left(t\right),\dot{\theta }\left(t\right))}^T$$

(1)

(hereafter, ${(·)}^T$ denotes a transpose), was obtained by the motion sensor. Here, $\theta$ [rad] is the roll angle from the vertical line to the bicycle’s vertical axis, and $\dot{\theta }:=d\theta /dt$ [rad/s] is the corresponding angular velocity. Figure 2 schematically shows the definition of the roll angle in the front view of the bicycle. The bicycle’s vertical axis was nominally determined based on the direction of the bicycle frame. The motion sensor was initially calibrated to output $\theta =0$ when the bicycle’s vertical axis was parallel to the vertical line. Therefore, $\theta =0$ does not indicate the upright equilibrium of the unmanned bicycle. Nevertheless, since the resulting calibration setting was commonly maintained for all participants and trials, the obtained datasets can be compared with each other.

Throughout the experiment, $\mathit{x}\left(t\right)$ were stored in the computer in the following form:

$$\{{\mathit{x}}_{hum}^{(s,n)}\left(t_0\right),\cdots ,{\mathit{x}}_{hum}^{(s,n)}\left(t_i\right),\cdots ,{\mathit{x}}_{hum}^{(s,n)}\left(t_{I-1}\right)\},i=0,\cdots ,I-1,s=1,\cdots ,S,n=1,\cdots ,N,$$

(2)

where $t_i:=i\Delta t$ [s] is a discrete time with a sampling period of $\Delta t$ [s], I is the length of the time series, s and S are an index and the number of participants, respectively, and n and N are an index and the number of trials, respectively.

For this study, we chose $\Delta t={10}^{-2}$ s and $I=\mathrm{18,001}$ to obtain the physical data length $(I-1)\Delta t=180$ s. The number of participants was $S=8$, and the number of trials undertaken by each participant was $N=5$.

Figure 3 shows the measured time series for $(s,n)=(1,1)$, i.e., for the first participant’s first trial. The result clearly exhibits fluctuations specific to human balancing motions in which large-amplitude spikes intermittently arise among the moderate-amplitude fluctuation process, which has already been recognized as temporal intermittency in the field of nonlinear physics [6,7,8]. Thus, we obtained the time series for all $s=1,\cdots ,S$ and $n=1,\cdots ,N$.

As shown in the upper graph of Figure 3, the vibrational center of $\theta \left(t\right)$ is shifted from the origin, $\theta =0$. This is mainly because the calibrated $\theta =0$ does not indicate the vertical equilibrium of the bicycle, as mentioned above. We statistically evaluate these shifts by identifying them as the mean roll angles given by

$$E^{\left(s\right)}\left[\theta \right]:=\frac{1}{N\times I}{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{i=0}^{I-1}}{x}_1^{(s,n)}\left(t_i\right),s=1,\cdots ,S,$$

(3)

which is the temporal average of the sth participant’s $\theta \left(t\right)$, further averaged over all of his trials. The resulting values are listed in

Table 1. The mean roll angle of the sth participant.

s	1	2	3	4
$E^{\left(s\right)}\left[\theta \right]$	$1.71\times {10}^{-2}$	$1.55\times {10}^{-2}$	$1.71\times {10}^{-2}$	$1.69\times {10}^{-2}$
$S{D}^{\left(s\right)}\left[\theta \right]$	$1.46\times {10}^{-2}$	$1.33\times {10}^{-2}$	$1.63\times {10}^{-2}$	$1.07\times {10}^{-2}$
$\mathit{s}$	5	6	7	8
$E^{\left(s\right)}\left[\theta \right]$	$1.59\times {10}^{-2}$	$1.65\times {10}^{-2}$	$1.83\times {10}^{-2}$	$2.31\times {10}^{-2}$
$S{D}^{\left(s\right)}\left[\theta \right]$	$1.40\times {10}^{-2}$	$1.37\times {10}^{-2}$	$1.30\times {10}^{-2}$	$1.42\times {10}^{-2}$

with standard deviations given by

$$S{D}^{\left(s\right)}\left[\theta \right]=\sqrt{V^{\left(s\right)}\left[\theta \right]},V^{\left(s\right)}\left[\theta \right]:=\frac{1}{N\times I+1}{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{i=0}^{I-1}}{\left(x_1^{(s,n)}\left(t_i\right)-E^{\left(s\right)}\left[\theta \right]\right)}^2,s=1,\cdots ,S.$$

(4)

Viewing the average over all participants’ values, $E^{\left(s\right)}\left[\theta \right]$ takes about $1.75\times {10}^{-2}$, indicating that the actual equilibrium angle of the bicycle was approximately at $\theta \approx 1.75\times {10}^{-2}$ rad (or ${1.01}^{\circ }$) in our setup. In addition, each $E^{\left(s\right)}\left[\theta \right]$ value is slightly different due to the respective riding forms of the participants.

2.3. Construction of Measured PDFs

First, we obtain $P_{hum}^{(s,n)}\left(x_1,x_2\right)$, the joint PDF with respect to the components of the time series in (2) for the sth participant’s nth trial, by normalizing the two-dimensional histogram of ${\{{\mathit{x}}^{(s,n)}\left(t_i\right)\}}_{i=0}^{I-1}$ with bin width $(\overline{x_k}-\underline{x_k})/N_{\mathrm{bin}}$ ($k=1,2$). Here, $N_{\mathrm{bin}}$ is the number of histogram bins and $\overline{x_k}$ and $\underline{x_k}$ are the upper and lower limits of $x_k$, respectively.

Next, the resulting $P_{hum}^{(s,n)}\left(x_1,x_2\right)$ is averaged over all trials $n=1,\cdots ,N$ to obtain the sth participant’s joint PDF as

$$P_{hum}^{\left(s\right)}\left(x_1,x_2\right)=\frac{1}{N}{\displaystyle \sum _{n=1}^N}{P}_{hum}^{(s,n)}\left(x_1,x_2\right).$$

(5)

We call (5) the measured PDF of the sth participant.

Figure 4 shows the measured PDF of all participants ($s=1,\cdots ,8$). In this study, we commonly set $N_{\mathrm{bin}}=40$, $(\underline{x_1},\overline{x_1})=(-0.04,0.08)$, and $(\underline{x_2},\overline{x_2})=(-0.3,0.3)$ for all PDFs, as was done in Figure 4.

3. Fluctuation Model of the Human–Bicycle Balance

In this section, we propose a new fluctuation model that allows us to reproduce the measured PDFs obtained above as simply as possible.

3.1. A Human–Bicycle Fluctuation Model

In view of our setup in Figure 2, we model the human–bicycle mechanics by a simple inverted pendulum about the contact point of the bicycle wheel. Based on the time series statistics in Table 1, we also assume that the roll angle $\theta$ about the equilibrium is sufficiently small; the maximal three standard deviation indicates $3\times {max}_{s}S{D}^{\left(s\right)}\left[\theta \right]\approx 4.88\times {10}^{-2}$ rad (or ${2.80}^{\circ }$ deg). Therefore, we model the bicycle’s rolling motion by a linearized inverted pendulum of the form:

$$m{r}^2\ddot{\theta }-mgr(\theta -\widehat{\theta })=T,$$

(6)

where m [kg] and r [m] are the mass and length of the pendulum, respectively, $g:=9.81$ m/s${}^2$ is the gravitational acceleration, T [Nm] is a torque input, and $\widehat{\theta }$ [rad] is the equilibrium angle. The equation of motion (6) is nondimensionalized and represented in the following state-space form:

$${\dot{x}}_1=x_2,{\dot{x}}_2=k(x_1-\widehat{\theta })+u,$$

(7)

where $\mathit{x}={(x_1,x_2)}^T:={(\theta ,\dot{\theta })}^T$, $k:=g/r$, and $u:=T/m{r}^2$.

Next, we specify u to simulate the human fluctuation during human–bicycle motion. As successfully demonstrated in our previous study [11], some human fluctuations can be accurately simulated by the following state-feedback mechanism:

$$u=F_1\{1+{\sigma }_1{\xi }_1\left(t\right)\}(x_1-\widehat{\theta })+F_2{x}_2+{\sigma }_2{\xi }_2\left(t\right),$$

(8)

where ${\xi }_1\left(t\right)$ and ${\xi }_2\left(t\right)$ are independent white Gaussian noises with zero mean and unit variance. The first coefficient $F_1\{1+{\sigma }_1{\xi }_1\left(t\right)\}$ represents a random proportional gain with mean $F_1$ and variance ${\left(F_1{\sigma }_1\right)}^2$. The second coefficient $F_2$ is a deterministic derivative gain. The third term ${\sigma }_2{\xi }_2\left(t\right)$ represents an additive random perturbation.

Finally, we substitute (8) into (7), rewrite it to reduce the dependency between parameters, and propose a human–bicycle fluctuation model of the form

$${\dot{x}}_1=x_2,{\dot{x}}_2=\{p_1+p_2{\xi }_1\left(t\right)\}(x_1-\widehat{\theta })+p_3{x}_2+p_4{\xi }_2\left(t\right),$$

(9)

with the generalized parameter vector

$$\mathit{p}=(p_1,p_2,p_3,p_4):=(g/r+F_1,F_1{\sigma }_1,F_2,{\sigma }_2)\in {\mathbb{R}}^4,$$

(10)

which parameterizes the equivalent properties of human–bicycle motion. Here, $p_1$ and $p_2$ are a mean and a standard deviation of the randomly fluctuating stiffness of the human–bicycle motion, respectively, $p_3$ is a deterministic viscous damping, and $p_4$ is an additive random fluctuation strength.

3.2. Calculation of Simulated PDFs

Using given $\widehat{\theta }$ and $\mathit{p}$, we obtain $N^{\prime }$ samples of the stationary numerical solution of (9) as

$${\{{\mathit{x}}_A^{\left(n^{\prime }\right)}(t_i;\widehat{\theta },\mathit{p})\}}_{i=0}^{I-1},n^{\prime }=1,\cdots ,N^{\prime },$$

(11)

by means of a fourth-order Runge–Kutta–Gill method with time step $\Delta t={10}^{-2}$ s, the same as the experimental sampling period. To generate these samples, $N^{\prime }$ different sequences of normal pseudo-random numbers [24],

$${\{{\nu }_i^{\left(n^{\prime }\right)}\}}_{i=0}^{I-1},n^{\prime }=1,\cdots ,N^{\prime },$$

(12)

are used to simulate the independent white Gaussian noises ${\xi }_1\left(t\right)$ and ${\xi }_2\left(t\right)$ by

$${\xi }_l\left(t_i\right)\approx {\nu }_i^{\left(n^{\prime }\right)}{(\Delta t)}^{-1/2},l=1,2,$$

(13)

where ${(\Delta t)}^{-1/2}$ is the numerical factor required for integrating stochastic differential equations [25].

From the simulated time series in (11), we construct the $n^{\prime }$th sample’s joint PDF $P_{sim}^{\left(n^{\prime }\right)}(x_1,x_2;\widehat{\theta },\mathit{p})$ using the same procedure and conditions as applied for deriving the measured PDF in Section 2.3. We also take the average $P_{sim}^{\left(n^{\prime }\right)}(x_1,x_2;\widehat{\theta },\mathit{p})$ over all samples by

$$\begin{array}{c}\hfill P_{sim}(x_1,x_2;\widehat{\theta },\mathit{p})=\frac{1}{N^{\prime }}\sum _{n^{\prime }=1}^{N^{\prime }}{P}_{sim}^{\left(n^{\prime }\right)}(x_1,x_2;\widehat{\theta },\mathit{p}).\end{array}$$

(14)

We call (14) the simulated PDFs for $\widehat{\theta }$ and $\mathit{p}$, which are to be compared with the measured $P_{hum}^{\left(s\right)}\left(x_1,x_2\right)$.

4. Method of Parameter Identification

In this section, we formulate the identification problem of an unknown parameter vector $\mathit{p}$ that allows the simulated $P_{sim}(x_1,x_2;\widehat{\theta },\mathit{p})$ to reproduce the measured $P_{hum}^{\left(s\right)}\left(x_1,x_2\right)$.

4.1. Parameter Identification Problem

We solve the optimization problem

$$\underset{\mathit{p}}{\mathrm{Minimize}}E\left(\mathit{p}\right)$$

(15)

with the cost function

$$\begin{array}{c}E\left(\mathit{p}\right):=\underset{(x_1,x_2,i)}{max}\left\{|F_{sim}^i(x_1,x_2;\widehat{\theta },\mathit{p})-F_{hum}^{i\left(s\right)}(x_1,x_2)|\right\}.\end{array}$$

(16)

Equation (16) is known as a two-dimensional Kolmogorov–Smirnov (KS) distance [12] and is used for two-dimensional goodness-of-fit testing between the empirical distribution of data and a hypothetical density law or between two distributions of separate data [12,13]. Here, $F_{sim}^i(x_1,x_2)$ and $F_{hum}^{i\left(s\right)}(x_1,x_2)$ ($i=1,\cdots ,4$) are the cumulative distribution functions (CDFs) with respect to the four quadrants about $(x_1,x_2)$ on the $(x_1,x_2)$-plane, i.e.,

$$\begin{array}{c}\hfill F_{sim}^i(x_1,x_2;\widehat{\theta },\mathit{p}):={\int \int }_{R_i}{P}_{sim}(x_1,x_2;\widehat{\theta },\mathit{p})d{x}_{1}d{x}_2,F_{hum}^{i\left(s\right)}(x_1,x_2;\widehat{\theta },\mathit{p}):={\int \int }_{R_i}{P}_{hum}^{\left(s\right)}(x_1,x_2)d{x}_{1}d{x}_2,\end{array}$$

(17)

with their domains $R_1:=[\underline{x_1},x_1]\times [\underline{x_2},x_2]$, $R_2:=[x_1,\overline{x_1}]\times [\underline{x_2},x_2]$, $R_3:=[\underline{x_1},x_1]\times [x_2,\overline{x_2}]$, and $R_4:=[x_1,\overline{x_1}]\times [x_2,\overline{x_2}]$, respectively. The cost function (16) evaluates the CDFs’ reproduction error and satisfies $E\left(\mathit{p}\right)=0$ if $P_{sim}(x_1,x_2;\widehat{\theta },\mathit{p})=P_{hum}^{\left(s\right)}\left(x_1,x_2\right)$. Hence, it also indicates the PDFs’ reproduction error.

4.2. Particle Swarm Optimization (PSO)

We employ PSO [26] to solve (15). Consider a swarm of M candidate solutions,

$$\begin{array}{c}\hfill \{{\mathit{p}}^1,{\mathit{p}}^2,\cdots ,{\mathit{p}}^i,\cdots ,{\mathit{p}}^M\},{\mathit{p}}^i\in {\mathbb{R}}^4,\end{array}$$

(18)

which are called particles. Each component of ${\mathit{p}}^i$ is recursively updated by

$$\begin{array}{c}\left\{\begin{array}{c}{v}_j^i(k+1)=\omega v_j^i\left(k\right)+c_1{\eta }_{1j}^i\left(k\right)\left(p{b}_j^i\left(k\right)-p_j^i\left(k\right)\right)+c_2{\eta }_{2j}^i\left(k\right)\left(g{b}_j\left(k\right)-p_j^i\left(k\right)\right),\hfill \\ p_j^i(k+1)=p_j^i\left(k\right)+v_j^i(k+1),(k=0,1,\cdots ,K),\hfill \end{array}\right.\hfill \end{array}$$

(19)

where $p_j^i\left(k\right)$ denotes the jth component of ${\mathit{p}}_i$ at iteration k; $v_j^i\left(k\right)$ is the corresponding velocity; $\omega$, $c_1$ and $c_2$ are system parameters of PSO; ${\eta }_{1j}^i\left(k\right)$ and ${\eta }_{2j}^i\left(k\right)$ are random numbers independently generated by [24] for each i, j, and k with a uniform distribution in $[0,1]$; and $p{b}_j^i\left(k\right)$ and $g{b}_j\left(k\right)$ are the jth components of the vectors ${\mathit{pb}}_i\left(k\right)$ and $\mathit{gb}\left(k\right)\in {\mathbb{R}}^4$, respectively.

${\mathit{pb}}_i\left(k\right)$ is the position of the particle taking the lowest cost among those at ${\mathit{p}}_i\left(0\right),\cdots ,{\mathit{p}}_i\left(k\right)$; this is called the personal best. $\mathit{gb}\left(k\right)$ is the position of the particle with the lowest cost among all particles for all iterations up to k; this is called the global best. For sufficiently large K, $\mathit{gb}\left(K\right)$ is expected to be close to the optimal solution ${\mathit{p}}^{*}$.

5. Identification Results

5.1. Identification Condition

In our PSO application, the number of particles was set to $M=32$, and the initial particles ${\mathit{p}}^i\left(0\right)$, $i=1,\cdots ,M$ were given by random points uniformly distributed within the four-dimensional hyperrectangle

$${\mathit{D}}_0:=\{(p_1,p_2,p_3,p_4)|-30\le p_1\le 0,0\le p_2\le 3,-10\le p_3\le 0,0\le p_4\le 3\}.$$

(20)

The number of iterations is $K=500$. We set the number of the simulated samples in (14) to $N^{\prime }=5$, the same as the number of experimental trials par participant.

We hereafter denote by ${\mathit{p}}^{\left(s\right)}$ the optimized solution $\mathit{gb}\left(K\right)$ obtained from the sth participant’s data: $P_{hum}^{\left(s\right)}\left(x_1,x_2\right)$ and $\widehat{\theta }={\widehat{\theta }}^{\left(s\right)}:=E^{\left(s\right)}\left[\theta \right]$. We also use the notation $P_{sim}^{\left(s\right)}\left(x_1,x_2\right):=P_{sim}(x_1,x_2;{\widehat{\theta }}^{\left(s\right)},{\mathit{p}}^{\left(s\right)})$ for the simulated PDFs derived from the sth participant’s data.

5.2. Identification Results

Table 2. Identified ${\mathit{p}}^{\left(s\right)}=(p_1,p_2,p_3,p_4)$ and its cost value for the sth participant. “**” denotes the best result and “*”, the worst. The last two columns show the corresponding Gaussian results.

Our Proposed Fitting							Gaussian Fitting
$\mathit{s}$	${\mathit{p}}_{\mathbf{1}}$	${\mathit{p}}_{\mathbf{2}}$	${\mathit{p}}_{\mathbf{3}}$	${\mathit{p}}_{\mathbf{4}}$	$\mathit{E}\left({\mathit{p}}^{\left(\mathit{s}\right)}\right)$	Fitness	$\mathit{E}$	Fitness
1	$-32.6272$	$3.35867$	$-5.84823$	$2.74614$	$9.611\times {10}^{-3}$	$99.04$%	$1.134\times {10}^{-1}$	$88.66$%
2	$-39.6902$	$1.21198$	$-1.72440$	$1.53507$	$7.189\times {10}^{-3}$	$99.28$%	$1.684\times {10}^{-1}$	$83.16$%
3	$-44.5088$	$1.08744$	$-5.97971$	$3.70913$	$7.056\times {10}^{-3}$	$99.29$%	$1.234\times {10}^{-1}$	$87.66$%
4	$-46.6330$	$1.17335$	$-6.31916$	$2.58083$	$7.478\times {10}^{-3}$	$99.25$%	$1.652\times {10}^{-1}$	$83.48$%
5	$-37.5965$	$1.61197$	$-1.77080$	$1.54042$	$1.228\times {10}^{-2}$	$98.77$% *	$1.640\times {10}^{-1}$	$83.60$%
6	$-34.6138$	$0.993204$	$-7.01156$	$2.96107$	$5.689\times {10}^{-3}$	$99.43$% **	$1.351\times {10}^{-1}$	$86.49$%
7	$-38.5872$	$1.21683$	$-7.31069$	$3.07624$	$6.011\times {10}^{-3}$	$99.40$%	$9.134\times {10}^{-2}$	$90.87$%
8	$-32.5560$	$2.55268$	$-2.15329$	$1.63554$	$8.500\times {10}^{-3}$	$99.15$%	$1.225\times {10}^{-1}$	$87.75$%

lists the identified vector components of ${\mathit{p}}^{\left(s\right)}$ by PSO for all participants $s=1,\cdots ,8$ and the corresponding KS cost value $E\left({\mathit{p}}^{\left(s\right)}\right)$. The seventh column indicates the cost value as a PDF fitness value of the form

$$\mathrm{Fitness}:=(1-E)\times 100\%,$$

(21)

which indicates the accuracy of our human–bicycle fluctuation model (9) in terms of the reproducibility of PDFs. The best and worst results are indicated by “**” and “*”, respectively. The second-last and the last columns show the cost value and the corresponding fitness value, respectively, between the measured $P_{hum}^{\left(s\right)}\left(x_1,x_2\right)$ and the mathematical two-dimensional Gaussian PDF $P_{Gauss}^{\left(s\right)}\left(x_1,x_2\right)$ with the same mean vector and covariance matrix as those of the measured $P_{hum}^{\left(s\right)}\left(x_1,x_2\right)$.

The results clearly show that our proposed model (9) successfully achieved over 98.7% PDF fitness, even in the worst case ($s=5$). This implies that it provides much better fitness than conventional Gaussian models.

Next, the left column of Figure 5 shows the difference of our simulated $P_{sim}^{\left(s\right)}(x_1,x_2)$ from the measured $P_{hum}^{\left(s\right)}(x_1,x_2)$, and the right column shows that for the Gaussian PDF $P_{Gauss}^{\left(s\right)}(x_1,x_2)$, for all s. The plot ranges are the same as those of all plots.

As already indicated by the fitness values in Table 2, our simulated PDFs agree well with the measured PDFs, compared to the Gaussian PDFs. As the shapes of Gaussian PDFs are symmetric by definition, the right column’s results imply that our measured PDFs have asymmetric shapes, as indicated by the deep-colored peaks and troughs.

Therefore, our proposed model properly reproduced such asymmetric shapes, which were barely reproduced by the Gaussian distribution.

5.3. KS Test

Up to this point, we have found that our simple model (9) provides high reproducibility of the PDF shapes of the human–bicycle fluctuation. However, this does not directly support the stochastic reliability of our model. Therefore, in this section, we conduct a statistical-hypothesis test called the KS test [12,13] on our results to check their stochastic reliability.

As described in [12], the statistic of the one-sample KS test is provided by

$$Z(n,E)=n^{1/2}E,$$

(22)

where E is the KS distance already given in (16), and n is the number of samples. On the other hand, the two-sample KS test employs the following statistic:

$$Z(n_1,n_2,E)={\left(\frac{n_1{n}_2}{n_1+n_2}\right)}^{1/2}E,$$

(23)

where $n_1$ and $n_2$ are the numbers of two independent samples. Under the null hypothesis (the sample follows a given distribution or the two samples follow the same distribution) and for large n (or $n_1$ and $n_2$), the random variable Z follows the following CDF [12,13]:

$$F\left(z\right):=\mathrm{Prob}(Z>z)=2{\displaystyle \sum _{i=1}^{\infty }}{(-1)}^{i-1}exp(-2i^2{z}^2),$$

(24)

where the lowercase z denotes a value of Z. The value of $F\left(z\right)$ is called the p-value of the observed z. If this p-value is smaller than a pre-defined $\alpha$ called a significance level, the null hypothesis is rejected.

In Figure 6, the solid curve shows $F\left(z\right)$. The small circles plot the p-values between $P_{hum}^{\left(s\right)}$ and $P_{sim}^{\left(s\right)}$ at $z=Z(n_1,n_2,E)$ using $E=E\left({\mathit{p}}^{\left(s\right)}\right)$ listed in Table 2, the measured data length $n_1=I\times N=\mathrm{18,001}\times 5$, and the simulated $n_2=I\times N^{\prime }=n_1$. Under the significance level $\alpha =0.01$, the null hypothesis is rejected for $s=1,5$, and 8; that is to say that three of the eight measurements cannot be said to follow our simulated distributions. On the other hand, the cross marks show the results between $P_{hum}^{\left(s\right)}$ and $P_{Gauss}^{\left(s\right)}$ at $z=Z(n,E)$ using the E values listed in Table 2 and $n=n_1$. In this Gaussian case, the null hypothesis is rejected for all s; i.e., no measurements can be said to follow these Gaussian distributions.

Given the above, we conclude in terms of statistical hypothesis testing that our simple model (9) can simulate the time series of human–bicycle fluctuations much better than the Gaussian distribution.

Note that our proposed model can separately fit the measurements of individual participants; it does not provide a general description of human–bicycle fluctuation. Conversely, it provides a means of mechanical parameterization of individual difference. The obtained parameter vector $\mathit{p}$ is useful for comparing individuals and seeing how they differ, in a mechanical manner.

6. Conclusions

In this study, we have constructed a simple degree-of-freedom human–bicycle fluctuation model that accurately reproduces the PDFs of experimentally measured human–bicycle balance motions.

First, we measured the time series of the roll angular displacement and the velocity of the human–bicycle balance motions and constructed their PDFs. Using these PDFs as the training data, we identified the model parameters by PSO, minimizing the KS distance between the measured PDF from the participant and the simulated PDF from our model. The resulting PDF fitnesses were over $98.7\%$, indicating that the simulated PDFs were in close agreement with the measured ones.

Next, we applied the KS statistical hypothesis test to our results, showing that our model simulated the time series of human–bicycle fluctuation much better than the Gaussian distribution.

The above result leads to the conclusion that our proposed model can provide an accurate single-degree-of-freedom model of human–bicycle fluctuations.

In future work, using our model parameters, we plan to compare various cyclists of different ages and genders who ride different types of bicycles in different environments. We also plan to develop a multi-degree-of-freedom fluctuation model of human–bicycle balance motions, making it possible to simulate fluctuating bicycle running paths based on physically identified human–bicycle parameters.