Iterative Learning Control With Incomplete Information: A Survey

I. Introduction

WITH the improvement of modern life, people pay more and more attention to health and personal safety. Almost every country has been calling for green travel and the shared bicycles are very popular in all big cities around the world, like ofo bicycle and mobike [1]. So the development of intelligent bicycles must be an important trend in the future [2].

Traditional bicycles are widely used in our daily life as an individual transportation tool, which has a merit to relieve the traffic pressure. Although bicycle travel is very environmental friendly, there are still some limitations. Since traditional bicycle is completely driven by cyclist, if the road is in a poor condition with various geometric complexity, the torque output by the cyclist may decrease as well as the speed of bicycle. Besides, it will be very laborious for people to ride on an uphill road.

Once lacking of the power of human input, the bicycle will not maintain running. As far as this major drawback is concerned, man-powered bicycles are not as good as electric bicycles.

Electric bicycle equipped with motor can solve this problem. The energy expenditure of riding an electric bicycle was 31% lower than the conventional bicycle on the uphill [3]. It has greatly extended the scope of application for traditional bicycle by providing assisted power with a motor [4]−[6]. Nowadays, electric bicycle has become one of the main vehicles in the world. Compared with the other vehicles powered by gasoline, electric bicycle can reduce the energy consumption, the air pollution and noise [7]. As we can see from Li et al. [8], the emissions of pollutants by electric bicycles per person per kilometer are several times smaller than those by motorcycles and cars, more or less equivalent to those by buses, and higher than those by bicycles.

For electric bicycles, cyclist can control the output of motor by adjusting the power controller to adapt different road conditions or individual preferences. However, electric bicycle also has its limitations. Since the cyclist can control the motor freely, the speed of electric bicycles may be too fast to be safe. It not only affects the service life of the motor, but also increases pollution.

Moreover, electric bicycles are heavy, when the motor runs out of electricity, it will be very hard for people to ride with his own human power.

In order to protect the environment as well as to achieve the assisted riding, the electric bicycle equipped with torque sensor has been produced. The sensor detects the pedal torque and converts it into an electrical signal, which is the input for the motor controller [9]. However, the torque sensor may be easily damaged since it contacts with the cyclist’s feet directly and the damage costs are not negligible [10].

Besides, the emerging development of connected and automated vehicles imposes a significant challenge on current vehicle control systems [11].

To solve the problems mentioned above, the best way is to design a torque estimation method to measure the pedaling torque in real time and provide assistance [12], [13].

Niki et al. [14] proposed a torque separation method based on high pass filter, they used high pass filter (HPF) to separate the human torque input from running friction. In 2015, Sankaranarayanan and Ravichandran [15] designed a torque sensorless controller for bicycles, but they did not take the disturbance force during cycling into account. The Fourier filtering was also applied to decouple the load torque and pedaling torque [16].

In 2017, Fukushimal and Fujimoto [17] presented a method for pedaling torque estimation using the recursive least square algorithm with multiple forgetting factors and they considered the condition for travelling on upward as well. In 2018, Rallo et al. [18] measured the wheel speed with a magnetic wheel encoder and showed that the pedaling cadence is correlated with the speed oscillation.

In addition, Bertucci et al. [19] analyzed the crank torque in road cycling on level and uphill, and the results show that when the pedaling cadences are the same, the crank torque profile differences are minimal between level ground and uphill road. In this way, we can see that the main factor affecting pedaling torque is cadences. The inertial measurement unit is also used to detect the angular velocity [20], [21].

Nevertheless, none of them focused on the speed control. The purpose of electric assistance is to keep the speed within a certain range in all circumstances. And it is worthwhile to note that although the inertial measurement unit was used, none of the above papers considered the coordinate conversion of angular velocity.

Similarly to the inertial navigation, in order to obtain the accurate positions, angular velocity based on the inertial coordinate system needs to be converted to another coordinate, which is consistent with the movement of the bicycle. For instance, Nilsson et al. [22] presented an open-source wireless foot-mounted inertial navigation module with an intuitive and significantly simplified dead reckoning interface, which can display the data after coordinate conversion.

Disturbance observer is often used to eliminate the effects of interference on system control, considering that there are many kinds of disturbances during cycling [23].

Therefore, in order to achieve a safe and convenient riding, it is very important to design an improved highly accurate speed control method for electric bicycles.

This paper proposes an improved speed control method for the electric assisted bicycle. A second-order low-pass filter is designed in disturbance observer to estimate and compensate disturbances during cycling. The quaternion method is used to convert the coordinate and users can set a motor auxiliary ratio which is variable with the real-time speed of bicycle. The result finally comes out that the speed can be self-adaptively stabilized in a set safety range despite of the disturbances during cycling.

The rest of the paper is organized as follows. Section II describes the dynamic model of the electric assisted bicycle system. In Section III, the quaternion method is used to re-coordinate the angular velocity which is measured by inertial measurement units (IMU). Besides, we simplify the calculation of pedaling torque. Section IV presents the design of low-pass filter (LPF) in disturbance observer to estimate and compensate the disturbances during riding. Section V is the introduction of motor assisted control. The stability of the proposed method is analyzed in Section VI. Simulation and results comparison are presented in Section VII. Finally, some conclusions are drawn and summarized in Section VIII.

Table I lists the acronyms and symbols explanations used in this paper.

Table  I.  List of Parameter

Parameter	Description
${E_t}$	Translational kinetic energy
${E_r}$	Rotational kinetic energy
$M$	Mass of an electric assisted bicycle and a cyclist
$V$	Velocity of bicycle
$J$	Moment of inertia of wheel
$\omega $	Angular velocity of wheels
$T$	Total kinetic energy of bicycle
$r$	Radius of rear wheel
${P_h}$	Gravitational potential energy of a bicycle
$h$	Vertical height of cycling
$l$	Length of a wheel movement
$\varphi$	Slope angle of uphill
$\theta$	Rotation angle of wheel
${F_{\rm{man}} }$	Pedaling force of cyclist
${F_{\rm{motor}} }$	Assisted power of motor
${T_{\rm{fric}} }$	Disturbances during cycling
$B$	Friction coefficient
${T_{\rm{man}} }$	Pedaling torque
${J_m}$	Moment of inertia of crank shaft
${\omega _0}$	Pedaling angular velocity
${\omega _1}$	Angular velocity of the large gear
${\omega _2}$	Angular velocity of pinion
${\omega _3}$	Angular velocity of rear wheel
$D$	Quaternion from inertial coordinate to geographic coordinate
$N$	Quaternion from geographic coordinate to carrier coordinate
$C$	Quaternion from inertial coordinate to carrier coordinate
$\Omega $	Absolute angular velocity vectors of carrier coordinate
${\Omega _g}$	Absolute angular velocity vectors of geographic coordinate
${\omega _X}$, ${\omega _Y}$, ${\omega _Z}$	Absolute angular velocity vectors on pedal carrier coordinate
$i$, $j$, $k$	Three axial unit vectors
${\omega _{Xg}}$, ${\omega _{Yg}}$, ${\omega _{Zg}}$	Absolute angular velocity vectors on geographic coordinate
$ \lambda $	Longitude
$\alpha $	Latitude
${\omega _t}$	Absolute angular velocity in carrier coordinate system
${\omega _{gt}}$	Angular velocity of carrier coordinate relative to the geographic coordinate
$d$	Length of crank shaft
$m$	Mass of crank shaft
${G_1}$	Model of bicycle
$\tau $	Disturbance during cycling
$\hat \tau $	Estimated disturbance
${G_{\tau {\omega _{\rm{3}}}}}$	Transfer function from $\tau $ to ${\omega _3}$
${G_{T{\omega _{\rm{3}}}}}$	Transfer function from $T$ to ${\omega _3}$
${Q_1}$	Second-order low-pass filter (LPF)
${i_1}$, ${i_2}$, ${i_3}$, ${i_4}$	Branch current
${V_1}$, ${V_2}$, ${V_3}$	Branch voltage
${A_u}$	Pass-band gain
$Q$	Quality factor
${\omega _c}$	Cut-off angular frequency of LPF
${f_c}$	Cut-off frequency of LPF
$c$	Auxiliary coefficient
$k$	Transmission ratio between gears
${K_t}$	Motor constant
$I$	Current of motor winding
$U$	Average voltage of motor
${K_e}$	Motor constant
$R$	Winding resistance

 | Show Table

DownLoad: CSV

II. System Modeling

Assuming that an electric assisted bicycle is running on an uphill road, the longitudinal dynamic model of the electric bicycle is shown in Fig. 1 and the forces during riding are indicated. We consider a rear-wheel-assisted electric bicycle, the motor is installed on the rear wheel, so as to collect the angular velocity and provide assisted power directly.

Figure  1.  Electric assisted bicycle model

DownLoad: Full-Size Img PowerPoint

The kinetic energy of a moving rigid body consists of translational kinetic energy and rotational kinetic energy.

The expression of translational kinetic energy is

where $ M $ is the mass of the electric assisted bicycle and a cyclist, $ V $ is the velocity of the bicycle.

Considering there are two wheels on the bicycle and the expression of rotational kinetic energy is

where $ J $ is the moment of inertia of the wheel and $ \omega $ is the angular velocity of the wheel.

Therefore, the total kinetic energy of the electric assisted bicycle can be expressed as

where $ T $ is the total kinetic energy of the bicycle.

Considering $V = \omega \times r$, where $ r $ is the radius of rear wheel, here we suppose that the front wheel has the same radius as the rear wheel, T can be written as

Due to the uphill road, the gravitational potential energy of the bicycle $ {P_h} $ can be expressed as

where $ h $ is the vertical height of cycling.

Fig. 2 shows the uphill road riding, $ l $ is the length of the wheel movement, $ \varphi $ is the slope angle of the uphill, and $ \theta $ is the rotation angle of the wheel.

Figure  2.  Uphill road riding

DownLoad: Full-Size Img PowerPoint

With the equation of arc length calculation, we have

The vertical height of cycling can be rewritten as

The potential energy of the bicycle can be rewritten as

Based on the kinetic equation and the Lagrange function $ L = T - {P_h} $, we get the Lagrange function

The Lagrange function satisfies $ \dfrac{d}{{dt}}\left(\dfrac{{\partial L}}{{\partial \theta }}\right) - \dfrac{{\partial L}}{{\partial \theta }} = Q $. Considering there are two major powers during cycling, which are the pedaling force of cyclist ${F_{\rm{man}}}$ and the assisted power of motor ${F_{\rm{motor}}}$, we have

At the same time, considering the disturbances during cycling ${T_{\rm{fric}}}$ and the impedance force between tires and ground, $ B $ is friction coefficient, the final Lagrange equation can be written as follows:

Based on $ \ddot \theta = \dot \omega $ and $ \dot \theta = \omega $, the equation (11) can be rewritten as

Then the transfer function can be expressed as

Assuming $ {J_{\rm{1}}} =M{r^2}{\rm{ + 2}}J $, the block diagram of a transfer function of the bicycle is shown in Fig. 3.

Figure  3.  Block diagram of a transfer function

DownLoad: Full-Size Img PowerPoint

III. Coordinate Conversion

To get the pedaling torque ${T_{\rm{man}}}$, we can simplify the pedaling process into a circular motion and use moment formula $ M = {J_m}\dfrac{{d\omega }}{{dt}} $ to replace it, where $ {J_m} $ is the moment of inertia and $ \omega $ is the angular velocity of the pedal rotation.

Inertial measurement unit (IMU) is used to collect the pedal angular velocity, it is mainly composed of a three axes of gyroscope and an three directional accelerometer, which can measure the data of the three axes of the $ OXYZ $ coordinate system.

In order to measure the angular velocity of the pedal and the rear wheel, one IMU is installed directly on the right pedal of the bicycle, and the other IMU is firmly connected with the rear wheel to detect the angular velocity of the electric assisted bicycle.

Firstly, the coordinate must be mentioned. In this paper, there are three kinds of coordinate systems. The first one is the earth centered inertial coordinate system ($ i $ system). The origin of coordinate is located at the earth’s center, $ {O_i}{Z_i} $ is along the polar axis of the earth, $ {O_i}{X_i} $ and $ {O_i}{Y_i} $ are on the equatorial plane of the earth.

The second one is the geographic coordinate system ($ t $ system), where the coordinate origin is also the barycenter of the carrier. $ {O_t}{Z_t} $ axis points to the zenith along the local reference ellipsoid’s normal line, $ {O_t}{X_t} $ axis is along the local latitude line to the east and the $ {O_t}{Y_t} $ axis points to the north along the local meridian.

The third is the carrier coordinate system ($ b $ system), which is firmly connected with a carrier. In this paper, the inertial measurement unit is loaded on the pedal, so the pedal is the carrier. The origin of this coordinate is located at the barycenter of pedal. The $ {O_b}{X_b} $ axis points to the direction of bicycle heading on the plane of the carrier, the $ {O_b}{Y_b} $ axis is perpendicular to the $ {O_b}{X_b} $ axis and points to the right side. The $ {O_b}{Z_b} $ axis is perpendicular to the plane that consists of the $ {O_b}{X_b} $ axis and $ {O_b}{Y_b} $ axis.

$ {\omega _0} $ is the pedaling angular velocity, $ {\omega _1} $ is the angular velocity of the large gear and $ {\omega _2} $ is the angular velocity of the pinion. The angular velocity $ \omega $ that we mentioned in Section II is the angular velocity $ {\omega _3} $ of the rear wheel, which can be detected by the IMU. However, IMU measures the projection of the absolute angular velocity vector on the axis of the carrier coordinate system. To be exact, they are not equal to the previous $ {\omega _0} $ and $ {\omega _3} $.

In order to obtain the accurate angular velocity and acceleration to control the speed of the electric bicycle, it is necessary to get the real-time data of the carrier coordinate system (pedal) relative to the geographic coordinate system. In this paper, the quaternion method is used to convert the data measured in the carrier coordinate system to the geographic coordinate system [24].

Matveev et al. [25] proposed a coordinate transformation method in 2009. Assume that the quaternion of the transformation from the inertial coordinate system to the geographic coordinate system is $ D $, the quaternion of the transformation from the geographic coordinate system to the carrier coordinate system is $ N $, and $ C $ is the quaternion that only needs to be transformed once from the inertial coordinate system to the carrier coordinate system. It can be composed of the two intermediate coordinate transformations. Based on the principle of quaternion, we get

where “$ \circ $” is the product symbol of the quaternions.

The purpose of the paper is to get the data of the coordinate system relative to the geographic coordinate system, so we need to calculate the quaternion $ N $.

Left multiplying both sides of (14) by the conjugate quaternion $ \bar D $, we can get

Considering $ \bar D \circ D = 1 $ and reversing the equation, we get

Take the derivative of both sides of (16)

Using the dynamic equation of quaternion $ 2\dot N = N \circ \Omega $ and $ 2\dot {\bar N} = - \Omega \circ \bar N $, $ \dot C $ and $ \dot {\bar D} $ are written as follows

Substituting (18) and (19) into (17), we get

Substituting (16) into (20), the final expression can be written as

where $ \Omega $ and $ {\Omega _g} $ are the absolute angular velocity vectors of the carrier coordinate system and geographic coordinate system. The three components of $ \Omega $ can be measured directly by the three axes gyroscopes in IMU.

Assuming $ {\omega _X} $, $ {\omega _Y} $ and $ {\omega _Z} $ are the absolute angular velocity vectors on the three axes of the pedal carrier coordinate system, $ \Omega $ can be expressed as

where $ i $, $ j $ and $ k $ are the three axial unit vectors, and the quaternion matrix for the corresponding $ \Omega $ can be described as

Assuming that $ {\omega _{Xg}} $, $ {\omega _{Yg}} $ and $ {\omega _{Zg}} $ are the absolute angular velocity vectors on the three axes of the geographic coordinate system, $ {\Omega _g} $ can be expressed as

where $ {\omega _{ie}} $ is the rotational angular velocity of the earth and the value is $7.292 \times {10^{ - 5}}\;{\rm{rad/s}}$, $ \lambda $ is the longitude and $ \alpha $ is the latitude.

Assuming $ N = {\lambda _0} + i{\lambda _1} + j{\lambda _2} + k{\lambda _3} $ and substituting $ \Omega $ and $ {\Omega _g} $ into (21), we can get

Then, (28) can be simplified to

By solving the first order equation group (29), we can get the values of $ {\lambda _0} $,$ {\lambda _1} $,$ {\lambda _2} $,$ {\lambda _3} $, and the expression $ N = {\lambda _0} + i{\lambda _1} + j{\lambda _2} + k{\lambda _3} $.

$ N $ is the quaternion that transforms from the geographic coordinate system to the carrier coordinate system, while our purpose is to convert the coordinate reversely. Therefore, after measuring the absolute angular velocity $ {\omega _t} $ in the carrier coordinate system, the angular velocity $ {\omega _{gt}} $ of the carrier coordinate system relative to the geographic coordinate system can be obtained by

For the convenience of calculation, we simplify the human pedal model in this paper and the pedaling torque formula can be described as

where $ {\omega _0} $ is the angular velocity of pedal, $ {J_m} $ is the moment of inertia of the crank shaft. $ d $ and $ m $ are the length and mass of crank shaft respectively.

The chain transmission model of electric bicycle is shown in Fig. 4. With the principle of chain transmission, the rotational angular velocity of the rear wheel $ {\omega _3} $ can be calculated by the pedaling angular velocity $ {\omega _0} $. However, $ {\omega _3} $ is not proportional to $ {\omega _0} $ because of the disturbance factors such as inertia and friction during cycling. When the pedal stops rotating, $ {\omega _3} $ will still exist. That is why in this paper, an inertial measurement unit is loaded on the rear wheel to detect the angular velocity $ {\omega _3} $ in real time.

Figure  4.  Chain transmission model of the electric bicycle

DownLoad: Full-Size Img PowerPoint

IV. Observer Design

As we mentioned in Section II, the electric bicycle will be disturbed by many kinds of frictions during riding. This paper mainly considers three kinds of frictions, namely the friction between tires and road, the wind resistance and the impedance force due to the gravitational potential energy.

On the basis of kinematic model, Huang et al. [26] used a sensor fusion algorithm as an observer. In this paper, in order to compensate the disturbances of frictions during riding, a disturbance observer is used to estimate the interference, as shown in Fig. 5 [27]. Considering the variability of frictions, we set a special cut-off frequency in LPF so that the majority of disturbances can be observed and compensated.

Figure  5.  Block diagram of the proposed method

DownLoad: Full-Size Img PowerPoint

$ {G_1} $ is the model of the bicycle, we assume that the bicycle system is a first order model defined as $ {G_1}(s) = \dfrac{{\rm{1}}}{{{J_1}s + B}} $, $ {G_{\rm{2}}}(s) = {G^{{\rm{ -1}}}(s)} $, where $ {J_1} = M{r^2} + 2J $.

$ \tau $ is the disturbance during cycling and $\tau = - {T_{\rm{fric}}}$, it can be observed and compensated by disturbance observer described in this paper. In this way, the transfer function from $ \tau $ to the angular velocity of the wheel $ {\omega _3} $ can be written as

$ T $ is the input of the observer. With Mason formula, the transfer function from $ T $ to $ {\omega _3} $ can be expressed as

$ \hat \tau $ is the estimated disturbance, it can be observed to compensate $ \tau $, the equation is as follow:

$ {Q_1} $ is a second-order low-pass filter (LPF), it is used in observer to filter out high-frequency [28], [29]. As shown in (36), when $ {Q_{\rm{1}}}(s) = 1 $ in the low-frequency stage, we get $\hat \tau = \tau$. The disturbance can be estimated in this way. Fig. 6 shows the circuit diagram of LPF.

Figure  6.  The second-order low-pass filter (LPF)

DownLoad: Full-Size Img PowerPoint

Firstly, with Kirchhoff’s law, the following can be obtained

Therefore, the transfer function of LPF can be expressed as

where $ {A_u} $ is the pass-band gain, $ Q $ is quality factor, $ {\omega _c} $ is the cut-off angular frequency of LPF and $ {f_c} $ is the cut-off frequency.

The filtering effect on a certain frequency section can be achieved by selecting the different cut-off angular frequencies. In the low-frequency stage, $ {Q_{\rm{1}}}(s) = 1 $, $ {G_{T{w_{\rm{3}}}}} = {G_{\rm{2}}}^{{\rm{ - 1}}} $ and $ {G_{\tau {\omega _{\rm{3}}}}} = 0 $, the effect of disturbance can be compensated by $ {G_{\rm{2}}} $.

Similarly, in the high-frequency stage, $ {Q_{\rm{1}}}(s) = 0 $, the disturbance which is higher than $ {f_c} $ is effectively suppressed.

In general, people often pedal 65−80 times in a minute when they ride bicycles in a normal frequency. The corresponding frequency is from 1.083 Hz to 1.33 Hz. But the frequency can not be cut off immediately in the low-pass filter, because of the transition-band. Therefore, the cut-off frequency of LPF is set to be 4 Hz in this paper.

Taking the characteristics of the circuit into account, the resistance value cannot be too large, and the $ Q $ value is preferably less than 3. We set $ Q = 0.7 $, then the value of $ {R_1} $, $ {R_2} $, $ {C_1} $, $ {C_2} $ can be calculated.

The expression of LPF can be described as follows:

Fig. 7 shows the amplitude-frequency characteristic of designed LPF. We can see clearly that the curve shows a declining trend after 4 Hz. As indicated in the label, when it is at the cut-off frequency, the gain value is −2.93 dB, which is in accordance with the design of the low-pass filter.

Figure  7.  The amplitude-frequency characteristic of LPF

DownLoad: Full-Size Img PowerPoint

In order to test and verify the accuracy of disturbance observer $ \tau $, we consider the disturbance as sine wave, the expression is $ \tau = {\rm{2}}0\sin t $.

The simulation results are shown in Fig. 8. The dotted line represents the simulated disturbance $ \tau $ and the solid line is the estimated disturbance $ \hat \tau $ through observation. As we can see, there are small phase and amplitude differences between the two waveforms.

Figure  8.  Disturbance and observation

DownLoad: Full-Size Img PowerPoint

Fig. 9 captures and amplifies the waveforms from 0.6 s to 2 s. The two dashed lines indicate the amplitude values of the waveform vertexes respectively. It is clear that the difference in amplitude is 0.4, which proves that the observer can estimate the majority of disturbance.

Figure  9.  Difference between disturbance and observation

DownLoad: Full-Size Img PowerPoint

V. Motor-Assisted Control

To achieve the purpose of the assisted riding, a permanent magnet brushless DC motor is used in this paper and installed in the rear wheel.

When the riding speed decreases and riding becomes more laborious, the motor can provide assistance based on the real-time angular velocity of wheel. As shown in Fig. 5, $ {\omega _{\rm{3}}} $ is fed back to the motor-assist module, and the auxiliary equation can be expressed as

where $ c $ is the auxiliary coefficient determined in real-time by the wheel angular velocity.

Japanese industry proposed the PAP method, which has been adopted by regulation. The PAP method can be formulated by the linear piecewise equation [30]. Considering the motor assistance, we rewrite the auxiliary coefficient, and the specific relationship can be expressed as

where $ {\omega _{\min }} $ and $ {\omega _{\max }} $ are the pre-defined angular velocity limits. Once the wheel angular velocity $ {\omega _{\rm{3}}} $ is lower than the minimum limit $ {\omega _{\min }} $, the motor will provide assistance and the auxiliary coefficient $ c = 1 $. The motor and cyclist will provide power together with a proportion of 1:1.

When $ {\omega _{\rm{3}}} $ exceeds the maximum limit $ {\omega _{\max }} $, the motor will be reversed and $ c = {\rm{ - }}1 $. In this way, the electric bicycle will not be overspeed and become too dangerous.

As the permanent magnet brushless DC motor is used in this paper, it is equipped with position sensors. According to the detected difference of the rotor position, each phase coil will be electrified in turns, making the direction change of magnetic field, which is produced by stators. Therefore, we can use logic circuits to control the conduction sequence, so as to achieve motor inversions.

In the same way, if $ {\omega _{\rm{3}}} $ is in a range of $ {\omega _{\min }} \leq \omega \leq {\omega _{\max }} $, the auxiliary coefficient will vary between 0 and 1, so as to provide assistance intelligently.

As shown in Fig. 10, the predefined wheel angular velocity limits are $ {\omega _{\min }}{\rm{ = 12}} $ rad/s and $ {\omega _{\max }} = 16 $ rad/s respectively. The auxiliary coefficient values at different wheel angular velocities are shown in the figure.

Figure  10.  Auxiliary coefficient

DownLoad: Full-Size Img PowerPoint

Because of the gear transmission between motor and the rear wheel, the formula of speed transmission can be expressed as

where ${\omega _{\rm{wheel}}}$ is the angular velocity of the rear wheel, $ k $ is the transmission ratio between gears and ${\omega _{\rm{motor}}}$ is the angular velocity of motor.

In this paper, $ {\omega _{\rm{3}}} $ is ${\omega _{\rm{wheel}}}$, and the value of $ k $ is related to the radii of the gears.

For the permanent magnet brushless DC motor, the electromagnetic torque of the motor can be expressed as

where $ {K_t} $ is the motor constant, and $ I $ is the current of the motor winding.

The motor winding current $ I $ can be calculated as

where $ U $ is the average voltage of motor, $ {K_e} $ is the motor constant, ${\omega _{\rm{motor}}}$ is the angular velocity of motor and $ R $ is the winding resistance.

So, ${T_{\rm{motor}}}$ can be rewritten as

Substituting ${T_{\rm{motor}}} = c{T_{\rm{man}}}$ into (46), we have

Considering (43), the angular velocity of the bicycle wheel can be rewritten as

In this way, the speed of the electric bicycle can be controlled and dynamically stabilized in a certain safe range.

VI. System Analysis

Based on the controller (41) and the coefficient function $ c(\omega_3) $, we have

Combining (11), the differential equation of the system can be written as

The system dynamics can be written as

where $ \omega^{*} $ is the pre-set expected angular velocity, and we have

Consider the Lyapunov function $ V(x) $ as

Its time derivative is given by

where $ A = \left[ c(x+\omega ^*)-c(\omega ^{*}) \right ] $. Since $ c $ is a function of $ \omega_3 $, $ \omega ^* $ is set between $\omega_{\rm{min}}$ and $\omega_{\rm{max}}$. Then $ x\neq 0 $, $ Ax <0 $. We also have $ Ax = 0 $ when $ x = 0 $.

Due to the fact that $J_m > 0,\; \dot{\omega}_0>0,\;B>0$, then $ \dot{V}(x) $ is negative-definte and $ \dot{V}(x) = 0 $ only when $ x = \omega_3-\omega^* = 0 $. Therefore, the system is asymptotically stable at $ \omega^* $.

VII. Simulation

To show the validity of the proposed method, several simulations are implemented. For the convenience of this research, the signals for the input angular velocities of the three directions $ XYZ $ are all assumed to be sinusoid.

The expressions are depicted as $ {\omega _x} = 10\sin (2\pi t) + 10 $, $ {\omega _y} = 12\sin (2\pi t) + 12 $ and ${\omega _\textit{z}} = 14\sin (2\pi t) + 14$, respectively. The waveform is shown in Fig. 11.

Figure  11.  Angular velocity before transformation

DownLoad: Full-Size Img PowerPoint

Besides, the dynamic parameters used in simulation are shown in Table II.

Table  II.  Simulation Parameter Values

Parameter	Value
$M(man + cyclist)$	$100\;{\rm{kg}}$
$J$	$0.15\;{\rm{N \cdot m \cdot {s^2}/rad}}$
$r$	$0.3\;{\rm{m}}$
$B$	$0.2$
${J_m}$	$0.12\;{\rm{kg \cdot {m^2}}}$

 | Show Table

DownLoad: CSV

Fig. 12 shows the angular velocity and three components curve diagram after coordinate transformation. We can see clearly that the value of angular velocity has changed from the original carrier coordinate system to the geographic coordinate system, so as to provide the subsequent calculation of the pedaling torque.

Figure  12.  Angular velocity after transformation

DownLoad: Full-Size Img PowerPoint

Take a set of data of the initial angular velocity as an example, $ {\omega _x} = 10 $, $ {\omega _y} = 12 $ and ${\omega _\textit{z}} = 14$. They make up an angular velocity matrix [0 10 12 14]^T and the transformed matrix we get is [0 5.836 9.814 17.68]^T. Thus, the quaternion from the geographic coordinate system to the carrier coordinate system can be calculated and expressed as N = [–0.9265 –0.216 –0.08526 –0.3018]^T.

As we have got the value of transformed angular velocity, the pedaling torque can be obtained by the torque formula. The result is shown in Fig. 13.

Figure  13.  Pedaling torque

DownLoad: Full-Size Img PowerPoint

With ${T_{\rm{man}}}$ and the angular velocity of the wheel, the motor module can provide assistance in real time.

Article 58 of China’s road traffic safety law stipulates that the maximum velocity of an electric bicycle riding on a non-motorized lane must not exceed 15 km/h.

In 2013, Zhao et al. [31] studied that the average velocity of a bicycle is 15 km/h on an isolated lanes between motor vehicles and non-motor vehicles with better riding conditions in China.

The electric bicycle safety specification stipulates that the maximum velocity of an electric bicycle shall not exceed 25 km/h, which is issued in China in May 2018. Before this specification, the Chinese national standard stipulated that the maximum velocity of electric bicycles was 20 km/h.

Considering that the electric assisted bicycle discussed in this paper is not a fully electric powered bicycle, we believe that velocity should be controlled according to the standard of the bicycle.

Therefore, the safe velocity of ordinary people riding a bicycle is between 12 km/h and 16 km/h. Assuming the radius of wheel $ r $ is 0.3 m and the formula $v = \omega \times r$, we can get the range of angular velocity is 11.1−14.8 rad/s.

Accordingly, the minimum limited angular velocity for the wheel can be set to 10 rad/s and the maximum limit is 13 rad/s. If the wheel angular velocity is lower than 10 rad/s, the motor and manpower will provide power together by the ratio of 1:1. If the wheel angular velocity is higher than 13 rad/s, motor will be reversed, so as to avoid danger.

The angular velocity variation curves of the rear wheel are shown in Figs. 14 and 15.

Figure  14.  Angular velocity of the wheel in previous methods

DownLoad: Full-Size Img PowerPoint

Figure  15.  Angular velocity of wheel in previous and new proposed methods

DownLoad: Full-Size Img PowerPoint

Three kinds of methods are simulated in our system, two of which are the previous method. A low pass filter is used in method 1. In method 2, low pass filter and band pass filter are used in combination to estimate disturbance and human pedaling torque.

In Fig. 14, the dotted line indicates the result of previous method 1, and the solid line represents the previous method 2. After setting the same safety angular velocity range for both of methods, the rapidity of method 2 is faster but fluctuation is more turbulent.

Fig. 15 shows the contrast curve between the results of this paper and the previous method 1. It can be seen that they can both well control the angular velocity but the rapidity and stability are not the same. Obviously, the method proposed in this paper is much faster and more stable, it takes reaction in the first 5 seconds.

As the cyclist begins to step on, the angular velocity value is the largest at the beginning in Fig. 15. When the value exceeds 13 rad/s, motor will turn to the reverse direction to decrease the value and the curve will show a downward trend. Later, when the value is between 10 rad/s and 13 rad/s, the curve rises up due to the motor’s assistance.

With the control and assistance, we can see clearly that the wheel angular velocity can still be effectively stabilized in a safe range, even though there are disturbances such as frictions and slopes during riding, so the same as the velocity of an electric bicycle.

VIII. Conclusion

The contribution of this paper lies in the presenting of an improved speed control method for the electric assisted bicycle without torque sensor installed.

On the basis of dynamic model system, we combine the disturbance observer, coordinate transformation and motor assistance to realize the speed control.

Based on the simulation results above, the proposed method can well estimate and compensate the majority of disturbances during cycling, and improve the measurement and control accuracy by means of coordinate system conversion. Finally, the bicycle speed can be intelligently stabilized within a safe range with the motor assistance.

Compared with the previous methods, this new proposed method has a higher rapidity and stability than that of the previous methods. Meanwhile, the stability of the method is also analyzed. It further validates that the speed control method can not only assist the cycling under different road conditions, but also improve the accuracy.