Modeling a powered wheelchair with slipping and gravitational disturbances on inclined and non-inclined surfaces

Abstract

Wheelchairs are broadly accepted and are widespread because of the assistance they provide to people with limited mobility. The design of a good controller generally involves the formulation of a comprehensive wheelchair model. Most dynamic models in the literature presume non-inclined planer surfaces within-doors, and therefore fail to take the combined effects of both gravitational forces and rolling friction on the usable-traction into consideration. Wheel-slip situations are also commonly neglected. This paper contributes to wheelchair modeling by proposing and formulating a dynamic model that considers the effects of rolling friction and gravitational potential on the wheelchair’s road-load force, on both inclined and non-inclined surfaces. The dynamic model is derived through the Euler Lagrange procedure, and wheel slip is determined by an approach that reduces the convectional number of slip-detection encoders. In the closed-loop model, the input-output feedback controller is proposed for tracking the user inputs by torque compensation. The optimality of the resulting minimum-phase closed-loop system is then ensured through the performance index of the non-linear continuous-time generalized predictive control with good simulation results.

Keywords

Modeling control feedback linearization friction slipping wheelchair

1. Introduction

Although most wheelchair users are able to comfortably move a joystick and make a fine movement correction when driving, others are only able to click on switches. A number of potential users, on the other hand, are incapable of driving and controlling a powered wheelchair with such interfacing devices, and can only rely upon caretakers to access the environment. The usage and potential users of electric-powered wheelchairs are determined to a large extent by the functionality of their embedded controllers. For controller design therefore, it is necessary that the wheelchair model is comprehensive enough to reflect real situations. Modeling of differential drive wheelchairs and other wheeled-mobile systems has previously been studied from the ideal perspective, where only the kinematic properties are taken into consideration, without regard to the dynamic attributes.^1–5 Such models fail to account for the effects of the system’s mass, inertia and acceleration, and do not consider the contributions of both conservative and non-conservative forces on the system’s motion. However, outdoor wheelchair usage generally demands driving on paths with diverse ground surface characteristics. The slippery and hilly configurations encountered in such situations could complicate the controllability of a wheelchair, and lead to a severe accident with undesirable injuries to the user, if not taken into consideration during design.^6,7 Furthermore, the rolling friction between the wheels and the road surface could be used to realize an optimized transfer of torque to the wheels during acceleration, if properly regarded. This may improve the controllability and lessen the slipping situations during motion. These considerations are therefore very important in the formulation of a wheelchair model.

A few researchers have considered the effects of such dynamic conditions. Fierro and Lewis⁸ and Oubbati et al.⁹ for instance, presented the dynamic models of wheeled-robotic systems with constrained motion on non-inclined planer surfaces. Although the model took into consideration the contributions of mass and inertia, it neither regarded the effects of rolling friction, nor considered the influence of gravitational forces experienced by the system during normal use. Frictional effects have been taken into consideration in various studies,^10–12 while slipping situations in the dynamic model have also been accounted for.^13–15 Although such dynamic models are numerous, they are commonly based on various structures of wheeled-mobile systems, and may therefore not reflect the behavior of the convectional differential drive wheelchair with front castor wheels.

There are relatively few dynamic models formulated specifically to describe the motion of a wheelchair in the literature, and a few models that only take into consideration the influence of mass and inertia have been reported.^16–18 These presentations restrict the motion of a wheelchair to horizontal work-planes and therefore disregard the influence of gravitational potential. Emam et al. presented a wheelchair dynamic model for non-normal driving conditions involving wheel-slip situations.¹⁹ Although the model takes the rolling friction and wheel-slip effects into account, it also constrains the wheelchair’s motion to non-inclined planer surfaces.

A dynamic model that takes the rolling friction as well as the up-hill and down-hill gravitational forces into consideration has been reported.^20–22 However, the authors did not consider the estimation of slipping parameters. The most recent study (to the best of the authors’ knowledge) in this respect has been presented by Chénier et al.^23,24 The study proposed a good dynamic model for manual wheelchair propulsion, usable by a new kind of motorized roller ergometer to simulate the behavior of a wheelchair on both straight and curvilinear level ground paths. Two observations could be made with regards to this model. First, like most dynamic models, it is based on horizontal level grounds. The influence of gravitational potential is therefore assumed to be constant. Second, the model is formulated based on the no-slip condition, which means the rolling friction model does not take slip velocity into account. These are acceptable modeling assumptions for the proposed application of the model. However, the effects of gravitational potential and slipping situations could be significant during the actual outdoor usage. Such conditions therefore need to be taken into consideration during modeling. We could not find a wheelchair model that takes all of the aforementioned dynamic situations into account. This paper therefore contributes by considering the combined effects of extreme dynamic situations accessible to a wheelchair during both indoor and outdoor usage. This involves estimating the slipping parameters of a convectional differential drive wheelchair, by taking the wheelchair’s rolling friction and the varying gravitational potential on both inclined and non-inclined surfaces into consideration during modeling.

Differential drive wheelchair structures with two passive front castor wheels have restricted mobility with a reduced number of actuators, and are therefore both non-holonomic and non-linear. This reduces the feasibility of computing their optimal controllers.^25,26 Some classical controller synthesis methods that have previously been employed for non-holonomic and non-linear systems include feedback linearization and Lyapunov theory.²⁷ The latter has previously been considered in the control of simple non-holonomic models, that did not take into account the effects of rolling friction as well as the varying gravitational potential of inclined grounds.^8,28 This is because it lacks a systematic procedure for constructing the control function, except for passive systems.²⁹ The former method has been employed in this study. Although Bloch and McClamroch²⁵ and Campion et al.²⁶ demonstrated that it was not possible to stabilize non-holonomic systems to a single equilibrium by smooth feedback, Sarkar et al.³⁰ showed that the input-output feedback method could still be used to control such systems. The input-output feedback linearization procedure linearizes the system dynamics between the input and the output, and is considered to be an important control strategy for minimum-phase systems with stable internal dynamics.^31–35 From the conceptual knowledge of wheelchair driving, one can reach every desired location by specifying the linear speed and direction of the vehicle. The proposed close-loop model is therefore considered to track the linear speed and angular orientation from the joystick.

The rest of this document is organized as follows. Dynamic modeling with gravitational forces is presented first by taking the ideal non-holonomic constraints of the wheelchair into account. This is followed by the estimation and incorporation of slipping parameters into the dynamic model. Controller design is then presented, and the simulation results to validate the open-loop and the closed-loop dynamic model are discussed towards the end, together with some concluding remarks.

2. Dynamic model with gravitational forces

2.1. Description of the wheelchair platform

A differential drive wheelchair of the type shown in Figure 1 is considered for modeling. Point O located at distance b from the rear wheels along the Y-axis is the mid-point of rear axle, and is also presumed to be the origin of the body fixed frame ${X Y Z}$ , while point C is the center of mass. The generalized coordinates of the center of mass, stated with respect to the inertial coordinate frame ${x y z}$ are given by $q \in R^{n \times 1} = {[x_{g}, y_{g}, z_{g}, ϕ]}^{T}$ , while the Cartesian components of distance l between C and O are $\bar{x}, \bar{y}$ and $\bar{z}$ , (not shown in Figure 1). Furthermore, the angular positions of the right and left motorized wheels are $γ_{R}$ and $γ_{L}$ respectively. This kind of wheelchair structure has two active rear wheels, usually of radius r, and two passive front castor wheels of radius $r_{C}$ . These are driven by applying the torques $τ_{R}$ and $τ_{L}$ on the right and left rear wheel respectively. It may be necessary to explain that $θ$ is the instantaneous deviation of the Z-axis from the z-axis, whereas the angle between the x-axis and the line of intersection of the moving the XY plane and the stationary xy plane is $ψ$ . Lastly, $ϕ$ is the precession angle about the Z-axis in a counter-clockwise direction as visible in a body fixed frame. The Euler Lagrange’s approach that yields Equation (1) is proposed in the derivation of equations of motion:

\frac{d}{dt} (\frac{\partial L}{\partial {\overset{\cdot}{q}}_{i}}) - (\frac{\partial L}{\partial q_{i}}) - \sum_{j = 1}^{m} A_{ij} (q)^{T} λ_{j} = Q_{i}

(1)

where

$L$ is the Lagrange function formulated in Section 2.3;

${q_{i} |}_{i = 1, 2, \dots, n}$ are the generalized coordinates of the system, with n being the number of generalized coordinates required to describe the time evolution of the dynamical system explicitly;

$A (q)^{T}$ is the transpose of matrix $A (q)$ associated with system constraints, while m is the number of non-holonomic constraints;

$λ$ is a vector of the Lagrange multipliers, while $Q_{i}$ are the external forces aiding or resisting the motion of the system.

Figure 1.

A differential drive wheelchair model.

2.2. System constraints

Considering identical d.c. motors for the right and left rear wheels, the wheelchair’s linear and angular velocities, V and $\overset{\cdot}{ϕ}$ respectively, can be expressed in the body fixed coordinate system, according to Equation (2) in terms of ${\overset{\cdot}{γ}}_{R}$ and ${\overset{\cdot}{γ}}_{L}$ :

\begin{matrix} V = r / 2 ({\overset{\cdot}{γ}}_{R} + {\overset{\cdot}{γ}}_{L}) \\ \overset{\cdot}{ϕ} = r / 2 b ({\overset{\cdot}{γ}}_{R} - {\overset{\cdot}{γ}}_{L}) \end{matrix}

(2)

In order to facilitate the computation of system constraints, the ideal non-slipping condition is initially considered. The body fixed linear velocity vector $[{\overset{\cdot}{X}}_{g} {\overset{\cdot}{Y}}_{g} {\overset{\cdot}{Z}}_{g}]$ , of point C is expressed according to Equation (3), as a linear transformation of the inertial velocity vector. The transformation matrix $B$ of frame ${X, Y, Z}$ relative to ${x, y, z}$ is presented in Equation (43) in terms of Euler angles. The elements of $B$ in the rest of this document will henceforth be referred to as ${B_{XX} |}_{X = 1, 2, 3}$ :

{[{\overset{\cdot}{X}}_{g} {\overset{\cdot}{Y}}_{g} {\overset{\cdot}{Z}}_{g}]}^{T} = [B] {[{\overset{\cdot}{x}}_{g} {\overset{\cdot}{y}}_{g} {\overset{\cdot}{z}}_{g}]}^{T}

(3)

From Equation (2), it is easy to notice that $V = r {\overset{\cdot}{γ}}_{R} - b \overset{\cdot}{ϕ} = r {\overset{\cdot}{γ}}_{L} + b \overset{\cdot}{ϕ}$ . Noting that V in Equation (2) $= {\overset{\cdot}{X}}_{g}$ in Equation (3), the pure rolling constraints of the right and left wheels that express how the longitudinal motion of the wheelchair’s center of mass is restricted by the longitudinal velocity of the wheels may be stated as follows:

\begin{matrix} {\overset{\cdot}{X}}_{g} = r {\overset{\cdot}{γ}}_{R} - b \overset{\cdot}{ϕ} \\ {\overset{\cdot}{X}}_{g} = r {\overset{\cdot}{γ}}_{L} + b \overset{\cdot}{ϕ} \end{matrix}

(4)

Moreover, the non-holonomic restrictions constraining the wheelchair’s motion in the direction perpendicular to the driving axis of the wheels, and on the moving XY plane, (ground surface), could be presented as Equation (5):

\begin{matrix} {\overset{\cdot}{x}}_{g} B_{21} + {\overset{\cdot}{y}}_{g} B_{22} + {\overset{\cdot}{z}}_{g} B_{23} - \overset{\cdot}{ϕ} \bar{x} = 0 \\ {\overset{\cdot}{x}}_{g} B_{31} + {\overset{\cdot}{y}}_{g} B_{32} + {\overset{\cdot}{z}}_{g} B_{33} = 0 \end{matrix}

(5)

Being non-integrable and independent velocity constraints, the equations in (5) may be expressed in terms of Equation (6) with $A (q)$ being a full rank matrix of size $m \times n$ :

A (q) \overset{\cdot}{q} = 0

(6)

2.3. Kinetic and potential energy

Considering the wheelchair to be a single non-elastic rigid body, the kinetic energy, T, of the wheelchair with reference to the center of mass can be computed. Equation (7) observes the wheelchair symmetry, and therefore takes into account the assumption that the wheelchair’s center of mass is likely to lie along the longitudinal axis. This implies that $\bar{y}$ could be considered equal to zero without loss of generality:

\begin{matrix} T = \frac{1}{2} M ({\overset{\cdot}{x}}_{g}^{2} + {\overset{\cdot}{y}}_{g}^{2} + {\overset{\cdot}{z}}_{g}^{2}) + \\ \frac{1}{2} (I_{XX} ω_{X}^{2} + I_{YY} ω_{Y}^{2} + I_{ZZ} ω_{Z}^{2} - 2 I_{XZ} ω_{X} ω_{Z}) + \\ M [\bar{x} (ν_{0 Y} ω_{Z} - ν_{0 Z} ω_{Y}) + \bar{z} (ν_{0 X} ω_{Y} - ν_{0 Y} ω_{X})] \end{matrix}

(7)

where

M is the total mass of the wheelchair including all its components;

$I_{XX}$ , $I_{YY}$ and $I_{ZZ}$ are the moments of inertia of the wheelchair about the X, Y and Z axis respectively through point O, while $I_{XZ}$ is the product inertia about the X and Z axis through the same point;

$ω_{X}$ , $ω_{Y}$ and $ω_{Z}$ are the components of angular velocity $ω$ of the wheelchair given by Equation (8) along the $X, Y$ and Z axis respectively:

\begin{matrix} ω_{X} = \overset{\cdot}{ψ} B_{13} + \overset{\cdot}{θ} \cos ϕ \\ ω_{Y} = \overset{\cdot}{ψ} B_{23} - \overset{\cdot}{θ} \sin ϕ \\ ω_{Z} = \overset{\cdot}{ψ} B_{33} + \overset{\cdot}{ϕ} \end{matrix}

(8)

$ν_{0 X}, ν_{0 Y}$ and $ν_{0 Z}$ are components of inertial velocity of point C given by ${[{\overset{\cdot}{X}}_{g} {\overset{\cdot}{Y}}_{g} {\overset{\cdot}{Z}}_{g}]}^{T}$ along the instantaneous directions of $X, Y$ and Z axis respectively.

Thus with potential energy $U = Mg z_{g}$ , the Lagrange function $L$ can be computed as Equation (9).

\begin{matrix} L = T - U \\ = \frac{1}{2} M ({\overset{\cdot}{x}}_{g}^{2} + {\overset{\cdot}{y}}_{g}^{2} + {\overset{\cdot}{z}}_{g}^{2}) \\ + \frac{1}{2} (I_{XX} ω_{X}^{2} + I_{YY} ω_{Y}^{2} + I_{ZZ} ω_{Z}^{2} - 2 I_{XZ} ω_{X} ω_{Z}) \\ + M [\bar{x} (ν_{0 Y} ω_{Z} - ν_{0 Z} ω_{Y}) + \bar{z} (ν_{0 X} ω_{Y} - ν_{0 Y} ω_{X})] - Mg z_{g} \end{matrix}

(9)

2.4. Dynamic model development

The inclusion of the potential energy of gravitational forces in the Lagrangian function has the logical effect that the resulting equations of motion naturally take into account the gravitational forces subjected to the wheelchair relative to its position in the configuration space. The dynamic equations of motion are computed in accordance with the Lagrangian formalism. This simplifies the modeling process of a rigid body’s dynamics into a straightforward procedure that generates matrices $\bar{M} (q) \in R^{n \times n}$ , $\bar{C} (q, \overset{\cdot}{q}) \in R^{n \times n}$ and $G (q) \in R^{n \times 1}$ (see Appendix 1).

Because the kinetic energy equation (9) is a quadratic function of the generalized velocity vector $\overset{\cdot}{q}$ , $\bar{M} (q)$ could be referred to as the mass matrix, that relates T to $\overset{\cdot}{q}$ according to Equation (45). The symmetric mass matrix $\bar{M} (q)$ depends on the configuration of the dynamic system, and since T is always positive and bounded, matrix $\bar{M} (q)$ is also positive definite and bounded. With these properties, $\bar{M} (q)$ can be considered non-singular with positive eigenvalues. The product of $\overset{\cdot}{q}$ and matrix $\bar{C} (q, \overset{\cdot}{q}) \in R^{n \times n}$ in Equation (51) is the vector of centrifugal and Coriolis forces. It includes all the inertial forces resulting from centrifugal and Coriolis accelerations. The matrices $\bar{M} (q)$ and $\bar{C} (q, \overset{\cdot}{q})$ are related by the skew-symmetric matrix $\overset{\cdot}{\bar{M}} (q) - 2 \bar{C} (q, \overset{\cdot}{q})$ . The $G (q) \in R^{n}$ also in Equation (51) is the vector of gravitational forces, and is often present in a dynamical system from the mechanical point of view. Since it is continuous and depends only on the generalized positions $q$ , $G (q)$ is bounded for each bounded $q$ . Further details regarding these matrices may be found in the literature.^36–38 The ideal dynamic model that does not take the slipping effects into consideration, can therefore be computed as expressed in Equation (10):

\begin{matrix} \bar{M} (q) \overset{\cdot\cdot}{q} + \bar{C} (q, \overset{\cdot}{q}) \overset{\cdot}{q} + G (q) \\ = E (q) τ - F (\overset{\cdot}{q}) + A^{T} (q) λ \end{matrix}

(10)

The matrix $E (q)$ presented in Equation (44) is also an input transformation matrix that transforms the motor torques $τ$ , also in Equation (44), into the applied force components in the generalized coordinates. $F (\overset{\cdot}{q})$ is the vector of frictional forces described in the next section.

The computation of the first order kinematic model in Equation (11) without slip, requires the full rank transformation matrix $S (q) \in R^{n \times (n - m)}$ , which transforms the velocity vector $η \in R^{(n - m) \times 1}$ , to generalized velocity vector $\overset{\cdot}{q}$ . Both $S (q)$ and $η$ are presented in Equation (54). The matrix $S (q)$ is formed by a set of smooth linearly independent vector fields that span the null space of $A (q)$ , such that the product $S^{T} (q) A^{T} (q) = 0$ . It is evident from Equation (6) that $\overset{\cdot}{q}$ is in the null space of $A (q)$ , and it accordingly follows that $\overset{\cdot}{q} \in$ span ${S (q)}$ . However, depending on the physical interpretation given to $η$ , different choices of $S (q)$ may be possible:

\overset{\cdot}{q} = S (q) η (t)

(11)

3. Slipping parameters and frictional force

3.1. Slipping parameters

In this study, a wheelchair is considered slipping if there exists a difference between the computed or theoretical wheel circumferential velocity $r \overset{\cdot}{γ}$ , and the absolute or actual velocity V, at any given time instant. The slip ratio $s_{r}$ is generally computed by Equation (12), where r is the radius, and $\overset{\cdot}{γ}$ is the average rotational velocity of the driving wheel:

s_{r} = {\begin{matrix} \frac{r \overset{\cdot}{γ} - V}{r \overset{\cdot}{γ}} & driving : (r \overset{\cdot}{γ} > V) \\ \frac{r \overset{\cdot}{γ} - V}{V} & braking : (V > r \overset{\cdot}{γ}) \end{matrix}

(12)

Apart from the road surface texture and mechanical load of the wheelchair, the frictional force between the tire and the road surface is also dependent upon the slip ratio $s_{r}$ , because $τ_{R}$ and $τ_{L}$ are translated into wheelchair motion through the road-tire friction. If the slipping condition is allowed beyond a certain threshold, it lowers the traction force and significantly reduces the controllability. Regulation of wheel torques is therefore necessary to maintain traction within the acceptable limit. Slip detection could be a major challenge bearing in mind the slow steering speed of a wheelchair. However, considering that the front castor wheels are relatively far from the center of mass compared to the hind wheels, it could be assumed as a result of force effect that the castor wheels experience no slip. The castor wheels are therefore considered to reflect the real or actual velocity of the wheelchair, while the motorized wheels reflect the theoretical velocity.

Chénier et al. originally proposed an open-loop observer method for estimating the orientational directions of each of the two castor wheels, based on the kinematics of the rear wheels, and without using encoders.³⁹ However, this method is founded on the assumption that none of the wheels will slip. For this reason, a previously published method that requires no additional information regarding the environment or acceleration of the wheelchair is elaborated to take the effects of wheel rotations on linear velocity when driving on inclined paths into consideration.¹⁹ In this method, encoders are used in the determination of slipping velocity by measuring and comparing the actual and expected velocities of the wheelchair. This procedure demands two encoders on each of the front wheels, one for wheelchair orientation measurement, and the other for measuring wheel rotational angle. Similarly, two odometers are required on the driving wheels for absolute velocity and orientation measurement.

3.2. Determination of real velocity

As indicated in the previous section, besides the wheelchair’s geometry, the determination of the real center of mass velocity is based on the utilization of the orientational and rotational velocities of the castor wheels. Considering Figure 2, it is possible to derive Equations (13) to (15), which validate the possibility of representing the direction of the right castor wheel in terms of the left, thus reducing the number of front wheels encoders required for slip detection:

{\overset{\cdot}{ϕ}}_{o} = \frac{v_{Y}}{f + \bar{x}}

(13)

β = \frac{\sin α_{L}}{\cos α_{L} + \frac{b}{\bar{x} + f} \sin α_{L}}

(14)

α_{R} = \frac{\sin α_{L}}{\cos α_{L} + \frac{2 b}{\bar{x} + f} \sin α_{L}}

(15)

where ${\overset{\cdot}{ϕ}}_{o}$ is the rotational velocity of the wheelchair on the XY plane due to $v_{Y}$ and the radial distance $(f + \bar{x})$ along the X-axis. The linear velocities ${\vec{v}}_{R}$ and ${\vec{v}}_{L}$ , of the right and left castor wheels respectively, are considered to originate from the two components of rotational velocities in Equation (16), with each component being computed with respect to the point and axis of rotation. One component ${\vec{v}}_{R_{B}}$ and ${\vec{v}}_{L_{B}}$ , for the right and left castor wheel respectively, represents the effects of rotations about the center of the wheel (point B in Figure 3). The other component ${\vec{v}}_{R_{A}}$ and ${\vec{v}}_{L_{A}}$ represents the effect of rotation about point A in Figure 3, where the castor wheel is attached to the wheelchair. ${\vec{v}}_{R}$ and ${\vec{v}}_{L}$ can therefore be computed according to Equation (16):

\begin{matrix} {\vec{v}}_{R} = {\vec{v}}_{R_{B}} + {\vec{v}}_{R_{A}} \\ {\vec{v}}_{L} = {\vec{v}}_{L_{B}} + {\vec{v}}_{L_{A}} \end{matrix}

(16)

where ${\vec{v}}_{R_{B}}$ , ${\vec{v}}_{R_{A}}$ , ${\vec{v}}_{L_{B}}$ and ${\vec{v}}_{L_{A}}$ are given in Equations (55) to (58). The resultant velocity of the front wheels as conceived at the midpoint between the right and left castor wheels is then computed as Equation (17):

\begin{matrix} \vec{v} = \frac{{\vec{v}}_{R} + {\vec{v}}_{L}}{2} + [\begin{matrix} 0 \\ 0 \\ 0 \\ {\overset{\cdot}{ϕ}}_{o} \end{matrix}] = [\begin{matrix} v_{x} \\ v_{y} \\ v_{z} \\ {\overset{\cdot}{ϕ}}_{o} \end{matrix}] \end{matrix}

(17)

Figure 2.

Geometrical representation of the wheelchair, with the parameters that have been utilized in deriving the velocity of the center of mass from the castor wheels’ velocities.

Figure 3.

(a) Bird view and (b) side view schematic representation of a castor wheel.

Translating $\vec{v}$ to point O (see Figure 1), the new velocity of point O denoted as ${\vec{v}}_{O}$ can be considered equal to $v_{X}$ because of the non-holonomic constraint that restricts the lateral speed of the rear wheels $v_{Y}$ , to zero (see Figure 2). ${\vec{v}}_{O}$ can then be expressed in terms of $\vec{v}$ by Equation (18). The slipping parameters in Equation (19) can be computed from the actual and the relative center of mass velocities:

{\vec{v}}_{O} = \vec{v} {[B_{11} B_{12} B_{13} 1]}^{T}

(18)

ϵ = [\begin{matrix} {\overset{\cdot}{x}}_{f_{o}} \\ {\overset{\cdot}{y}}_{f_{o}} \\ {\overset{\cdot}{z}}_{f_{o}} \\ {\overset{\cdot}{ϕ}}_{o} \end{matrix}] - [\begin{matrix} {\overset{\cdot}{x}}_{o} \\ {\overset{\cdot}{y}}_{o} \\ {\overset{\cdot}{z}}_{o} \\ \overset{\cdot}{ϕ} \end{matrix}]

(19)

where ${\overset{\cdot}{x}}_{f_{o}}, {\overset{\cdot}{y}}_{f_{o}}, {\overset{\cdot}{y}}_{f_{o}}$ and ${\overset{\cdot}{ϕ}}_{o}$ are the components of ${\vec{v}}_{O}$ , while $[{\overset{\cdot}{x}}_{o} {\overset{\cdot}{y}}_{o} {\overset{\cdot}{y}}_{o} \overset{\cdot}{ϕ}]^{T}$ is the velocity of wheelchair at point O translated from $[{\overset{\cdot}{x}}_{g} {\overset{\cdot}{y}}_{g} {\overset{\cdot}{y}}_{g} \overset{\cdot}{ϕ}]^{T}$ of point C. By super-positioning the slipping velocities $ε$ and the generalized velocities without slip $S (q) η$ , presented in Equation (11), the kinematic model with slipping conditions taken into account can be expressed. Upon differentiation of the configuration velocity vector in Equation (20), the acceleration vector in Equation (21) can be expressed as follows:

\overset{\cdot}{q} = S (q) η (t) + ϵ

(20)

\overset{\cdot\cdot}{q} = \overset{\cdot}{S} (q) η + S (q) \overset{\cdot}{η} + \overset{\cdot}{ϵ}

(21)

3.3. Frictional and resistive force modeling

Modeling the rolling friction involves determining the relationship between the driving velocity and the normal force at the area of contact between the wheels and the ground. Because the rolling friction is quite non-linear and dependent upon many parameters, a reduced formulation that depends only upon the new generalized velocity vector, $\overset{\cdot}{q}$ , is considered. The reduced model in Equation (22) is a combination of viscous and Coulomb friction, with a dimensionless coefficient of rolling friction $μ$ , describing the ratio of frictional force responsible for friction coupling between the wheels and the ground:

F (ϵ) = μ_{vis} \overset{\cdot}{q} + μ_{cmb} N s gn (\overset{\cdot}{q})

(22)

where $\overset{\cdot}{q}$ is the new vector of generalized velocities including the slipping parameters as described in Equation (20), N is the normal force at point O, while $μ_{vis}$ and $μ_{cmb}$ denote the coefficient of viscous friction and coefficient of Coulomb friction respectively. Nevertheless, considering the low driving speed involved in wheelchairs, both $μ_{vis}$ and $\overset{\cdot}{q}$ may be considered very small, resulting in much greater Coulomb friction as compared to viscous friction. The viscous friction, $μ_{vis} \overset{\cdot}{q}$ , can therefore be neglected to simplify the model without loss of generality. The frictional force expressed in Equation (22) is not smooth and therefore not differentiable at $\overset{\cdot}{q} = 0$ . Since a continuous and time differentiable friction model is required for simulation, an approximation in Equation (24) based on Equation (23) may be used, with k being a constant that determines the approximation accuracy:

lim_{k \to \infty} \frac{2}{π} \arctan (k \overset{\cdot}{q}) = sgn (\overset{\cdot}{q}) |_{k > > 1}

(23)

s gn (\overset{\cdot}{q}) \approx \frac{2}{π} \arctan (k \overset{\cdot}{q})

(24)

The rolling friction model $(F)$ of the wheelchair is therefore expressed as follows:

F = \frac{2}{π} N μ_{max} [\begin{matrix} \arctan (k {\overset{\cdot}{q}}_{x}) \\ \arctan (k {\overset{\cdot}{q}}_{y}) \\ \arctan (k {\overset{\cdot}{q}}_{z}) \\ \arctan (k {\overset{\cdot}{q}}_{ϕ}) \end{matrix}]

(25)

The normal force N could be determined by solving the Lagrange multiplier related to the vertical velocity constraint at point O in the second component of Equation (5) according to Equation (26):

\begin{matrix} N = λ (2) \\ = - {[A {\bar{M}}^{- 1} A^{T}]}^{- 1} [\overset{\cdot}{A} \overset{\cdot}{q} + A M^{- 1} (E τ - \bar{C} - G)] \end{matrix}

(26)

where $λ (2)$ denotes the second element of the matrix. The general equation of motion (10) can then be transformed into a more appropriate form for control by including Equations (20) and (21). Since matrix $S (q)$ spans the null space of $A (q)$ , multiplying the result by $S^{T} (q)$ eliminates the constraint forces, and produces a dynamic model in Equation (27) which takes the rolling friction and gravitational effects on the wheelchair into account:

\begin{matrix} S^{T} \bar{MS} \overset{\cdot}{η} + S^{T} \bar{M} \overset{\cdot}{S} η + S^{T} \bar{C} S η + S^{T} \bar{M} \overset{\cdot}{ϵ} + S^{T} \bar{C} ϵ \\ + S^{T} G + S^{T} F = S^{T} E τ \end{matrix}

(27)

The dynamic model can be simplified to Equation (28) with $G_{n}$ being full rank and non-singular matrix:

\overset{\cdot}{η} = [F_{n}] + [G_{n}] τ

(28)

where:

\begin{matrix} [F_{n}] = - {(S^{T} \bar{MS})}^{- 1} (S^{T} \bar{M} \overset{\cdot}{S} η + S^{T} \bar{C} S η + S^{T} \bar{M} \overset{\cdot}{ϵ}) \\ - {(S^{T} \bar{MS})}^{- 1} (S^{T} \bar{C} ϵ + S^{T} G + S^{T} F) \\ [G_{n}] = {(S^{T} \bar{MS})}^{- 1} (S^{T} E) \end{matrix}

4. The closed-loop design

In this section, consideration is given to the dynamic model in Equation (28), in an effort to track and stabilize the linear velocity and angular position specified by the user. Expressed in the state space form, system (29) does not satisfy the involutivity condition and is therefore not full-state feedback linearizable.^27,40,30 However, it could be input-output feedback linearizable with appropriate choice of output equation. In this paper, the velocity control problem is translated into torque regulation. This involves taking the dynamic properties into consideration to facilitate the coupling and linearization between the inputs and outputs. The problem is solved by computing the torques to guarantee the tracking of reference inputs on flat, inclined and slippery surfaces, and by optimizing the closed-loop gains through the performance index of non-linear continuous-time generalized predictive control (NCGPC):

\begin{matrix} \overset{\cdot}{x} & = f_{1} (x) + g_{1} (x) u \\ y & = h (x) \end{matrix}

(29)

where:

\begin{matrix} x = [v \overset{\cdot}{ϕ}] h (x) = [v ϕ] \\ f_{1} (x) = [\begin{matrix} - \bar{x} \cos θ \overset{\cdot}{ψ} \overset{\cdot}{ϕ} - g \sin ϕ \sin θ \\ - \frac{Mg \bar{x} \sin θ \cos ϕ}{(I_{ZZ} - M {\bar{x}}^{2})} \end{matrix}] \\ g_{1} (x) = [\begin{matrix} \frac{1}{Mr} & \frac{1}{Mr} \\ \frac{b}{r (I_{ZZ} - M {\bar{x}}^{2})} & \frac{- b}{r (I_{ZZ} - M {\bar{x}}^{2})} \end{matrix}] u = [\begin{matrix} τ_{R} \\ τ_{L} \end{matrix}] \end{matrix}

4.1. Input-output feedback linearization

4.1.1. Navigation task

The design approach of a closed-loop model is generally dependent on the complexity of the control task to be executed. The attention in this study is directed towards wheelchair steering control. The steering task has always been the obligation of a driver, who has to specify the steering signals continuously through the available user interface, in order to direct the wheelchair according to path geometry to a desired destination. A linearized system is necessary in this respect, because it aids the user to make proportional judgments regarding the inputs required to produce a given desired output. The linearization process is accomplished in this study by the input-output feedback control method and with proper consideration of the output variables. To ensure the practicality of the steering task, the output variables of a convectional wheelchair joystick are considered. Because the dynamic model in Equation (28) has two inputs, it is possible to choose any two output variables. The problem of speed and orientation control is therefore considered with an intention of tracking the specified linear velocity and the corresponding angular position. This is accomplished by ensuring that the errors $e = {[\begin{matrix} e_{1}^{1} & e_{1}^{2} \end{matrix}]}^{T}$ specified in Equation (30) results in $lim_{t \to \infty} [e] = 0$ :

\begin{matrix} e_{1}^{1} = v - v_{r} \\ e_{1}^{2} = ϕ - ϕ_{r} \end{matrix}

(30)

4.1.2. System relative degree ( $ρ$ )

From Equation (31), it can be observed that the first derivative of $e$ lacks the control signal $u$ in the second element, which explains why Equation (29) could not be input-state linearized by exact feedback linearization. However, $u$ appears by delaying the first element of $\overset{\cdot}{e}$ while performing the second differentiation on the second element:

\begin{matrix} {\overset{\cdot}{e}}_{1}^{1} & = - \bar{x} \cos θ \overset{\cdot}{ψ} \overset{\cdot}{ϕ} - g \sin ϕ \sin θ - {\overset{\cdot}{v}}_{r} + \frac{1}{Mr} (τ_{R} + τ_{L}) \\ {\overset{\cdot}{e}}_{1}^{2} & = e_{2}^{2} \\ {\overset{\cdot}{e}}_{2}^{2} & = - \frac{Mg \bar{x} \sin θ \cos ϕ}{(I_{ZZ} - M {\bar{x}}^{2})} - {\overset{\cdot\cdot}{ϕ}}_{r} + \frac{b (τ_{R} - τ_{L})}{r (I_{ZZ} - M {\bar{x}}^{2})} \end{matrix}

(31)

The relative degree of $e^{1}$ is therefore one, while $e^{2}$ has a relative degree of two. Since the sum of the components of vector relative degree is greater than the state dimension of the system, state extension is performed to enable computation of the outcome for every state, and a linearized Equation (32) of the system is presented as follows:

\begin{matrix} [\begin{matrix} {\overset{\cdot}{e}}_{1}^{1} \\ {\overset{\cdot}{e}}_{2}^{2} \end{matrix}] = [\begin{matrix} - \bar{x} \cos θ \overset{\cdot}{ψ} \overset{\cdot}{ϕ} - g \sin ϕ \sin θ - {\overset{\cdot}{v}}_{r} \\ - \frac{Mg \bar{x} \sin θ \cos ϕ}{(I_{ZZ} - M {\bar{x}}^{2})} - {\overset{\cdot\cdot}{ϕ}}_{r} \end{matrix}] \\ + [\begin{matrix} \frac{1}{Mr} & \frac{1}{Mr} \\ \frac{b}{r (I_{ZZ} - M {\bar{x}}^{2})} & \frac{- b}{r (I_{ZZ} - M {\bar{x}}^{2})} \end{matrix}] u \end{matrix}

(32)

4.1.3. The control law

Considering $ϑ_{1}$ and $ϑ_{2}$ as the synthetic control inputs, the state feedback law that compensates for input-output non-linearities can be computed, because the two rank decoupling matrix is invertible and non-singular unless $r = 0$ or $I_{ZZ} - M {\bar{x}}^{2} = 0$ . However, since r cannot be equal to zero it is possible to ensure through proper wheelchair design that $I_{ZZ} - M {\bar{x}}^{2} \neq 0$ throughout. In this study therefore, $I_{ZZ}$ , M and $\bar{x}$ are carefully chosen to ensure the decoupling matrix is always invertible. Since system (29) is square, with an equal number of inputs and outputs as well as the non-singular decoupling matrix, it is possible to compute the input $u$ that drives Equation (29) to the desired response. Equation (33) can therefore be expressed as follows:

\begin{matrix} u = {[\begin{matrix} \frac{1}{Mr} & \frac{1}{Mr} \\ \frac{b}{(I_{ZZ} - M {\bar{x}}^{2}) r} & \frac{- b}{(I_{ZZ} - M {\bar{x}}^{2}) r} \end{matrix}]}^{- 1} \times \\ [(\begin{matrix} {\overset{\cdot}{e}}_{1}^{1} \\ {\overset{\cdot}{e}}_{2}^{2} \end{matrix}) - (\begin{matrix} - \bar{x} \cos θ \overset{\cdot}{ψ} \overset{\cdot}{ϕ} - g \sin ϕ \sin θ - {\overset{\cdot}{v}}_{r} \\ - \frac{Mg \bar{x} \sin θ \cos ϕ}{(I_{ZZ} - M {\bar{x}}^{2})} - {\overset{\cdot\cdot}{ϕ}}_{r} \end{matrix})] \end{matrix}

(33)

where ${\overset{\cdot}{e}}_{i} = - \sum_{j = 0}^{ρ_{i} - 1} K_{ij} L_{f}^{j} (h_{i}) |_{i = 1, 2 .}$ with $ρ$ being the relative degree. To guarantee the stability of the closed-loop system, positive values of the tracking coefficients $K_{11}, K_{21}$ and $K_{22}$ are considered in the control law in Equation (34):

(\begin{matrix} {\overset{\cdot}{e}}_{1}^{1} \\ {\overset{\cdot}{e}}_{2}^{2} \end{matrix}) = [\begin{matrix} - K_{11} e_{1}^{1} \\ - K_{21} e_{1}^{2} - K_{22} e_{2}^{2} \end{matrix}]

(34)

Redefining $h (x)$ as ${[e_{1}^{1} e_{2}^{2}]}^{T}$ the internal dynamics of the system can be analyzed. The error $e_{1}^{1} = 0$ implies that ${\overset{\cdot}{e}}_{1}^{1} = 0$ , however $e_{2}^{2} = 0$ means that ${\overset{\cdot}{e}}_{2}^{2} = - K_{21} e_{1}^{2}$ . The error $e_{2}^{2}$ is therefore governed by ${\overset{\cdot}{e}}_{2}^{2} = - K_{21} e_{1}^{2}$ . Nonetheless, since $K_{21}$ is positive, the internal dynamics are stable and the system is well-behaved.

4.2. Non-linear continuous-time generalized predictive control

In this section, the NCGPC is utilized in an attempt to validate and provide an optimal design of the input-output controller gains through its performance index. The output $y = h (x)$ of system (29) with $x \in X \subset R^{n}, y \in Y \subset R^{m} | m = n$ is considered in the receding horizon performance index in Equation (35), where $T_{1}$ is the minimum prediction time and $T_{2}$ is the maximum prediction time with $T_{1} \leq \bar{τ} \leq T_{2}$ :

J = \frac{1}{2} \int_{T_{1}}^{T_{2}} {[y (t + \bar{τ}) - y_{r} (t + \bar{τ})]}^{T} [y (t + \bar{τ}) - y_{r} (t + \bar{τ})] d \bar{τ}

(35)

$y (t + \bar{τ})$ is computed by Taylor series expansion to the order equal to the corresponding relative degree of its elements as given by Equation (36) below:

y_{i} (t + \bar{τ}) = \sum_{k = 0}^{ρ_{i}} \frac{{\bar{τ}}^{k}}{k!} y_{i}^{(k)} (t) + {H (\bar{τ}) |}_{k > ρ_{i}}

(36)

In Equation (36) t is the current time, $t + \bar{τ}$ is the duration of prediction, $ρ_{i}$ is the relative degree of the $i^{th}$ output while ${H (\bar{τ}) |}_{k > ρ_{i}}$ are higher order terms of the expansion. Defining $e (t + \bar{τ}) = y (t + \bar{τ}) - y_{r} (t + \bar{τ})$ and neglecting the higher order terms of Equation (36), the Taylor expansion yields:

e (t + \bar{τ}) = Γ (\bar{τ}) \tilde{E} (t)

(37)

where $Γ (\bar{τ})$ is an $(n \times \sum_{i} (r_{i} + 1))$ matrix whose elements are $(n \times n)$ diagonal matrices and $\tilde{E} (t)$ is a vector of errors expressed in Equation (38).

\begin{matrix} Γ (\bar{τ}) = [\begin{matrix} I & \bar{τ} & {\bar{τ}}^{2} / 2! & \dots & {\bar{τ}}^{p} / p! \end{matrix}] \\ \tilde{E} (t) = [\begin{matrix} ϕ (t) - ϕ_{r} (t) \\ v (t) - v_{r} (t) \\ \overset{\cdot}{ϕ} (t) - {\overset{\cdot}{ϕ}}_{r} (t) \\ \overset{\cdot}{v} (t) - {\overset{\cdot}{v}}_{r} (t) \\ \overset{\cdot\cdot}{ϕ} (t) - {\overset{\cdot\cdot}{ϕ}}_{r} (t) \end{matrix}] \\ = [\begin{matrix} e_{1}^{2} \\ e_{1}^{1} \\ e_{2}^{2} \\ f_{11} - {\overset{\cdot}{v}}_{r} (t) + \frac{τ_{L}}{Mr} + \frac{τ_{R}}{Mr} \\ f_{12} - {\overset{\cdot\cdot}{ϕ}}_{r} (t) - \frac{b τ_{L}}{r (I_{ZZ} - M {\bar{x}}^{2})} + \frac{b τ_{R}}{r (I_{ZZ} - M {\bar{x}}^{2})} \end{matrix}] \end{matrix}

(38)

Equation (35) is then simplified to:

J = \frac{1}{2} {\tilde{E}}^{T} (t) Λ (T_{1}, T_{2}) \tilde{E} (t)

(39)

with:

Λ (T_{1}, T_{2}) = \int_{T_{1}}^{T_{2}} Γ^{T} (\bar{τ}) Γ (\bar{τ}) d \bar{τ}

Minimizing the cost expressed in Equation (39) with respect to the control parameter $u$ results in the following control law:

u = D Λ_{ss} Λ_{s} (\begin{matrix} - e_{1}^{2} \\ - e_{1}^{1} \\ - e_{2}^{2} \\ - f_{11} + {\overset{\cdot}{v}}_{r} \\ - f_{12} + {\overset{\cdot\cdot}{ϕ}}_{r} \end{matrix})

(40)

where $D$ in Equation (32) is the decoupling matrix, and $Λ_{ss} Λ_{s} (T_{1}, T_{2}, ρ) = K$ , see Equations (52) and (53). In this study, $T_{1} = 0$ and $T_{2} = T$ is opted for, for reasons of simplicity. However, this result applies as well to cases with $T_{1} \neq 0$ . The gain matrix $K (T, ρ)$ may thus be expressed:

K (T, r) = (\begin{matrix} 0 & \frac{3}{2 T} & 0 & 1 & 0 \\ \frac{10}{3 T^{2}} & 0 & \frac{5}{2 T} & 0 & 1 \end{matrix})

and the new control law in Equation (41) is obtained as follows:

\begin{matrix} u = {[\begin{matrix} \frac{1}{Mr} & \frac{1}{Mr} \\ \frac{b}{r (I_{ZZ} - M {\bar{x}}^{2})} & \frac{- b}{r (I_{ZZ} - M {\bar{x}}^{2})} \end{matrix}]}^{- 1} \\ \times [\begin{matrix} - f_{11} - \frac{3 e_{1}^{1}}{2 T} + {\overset{\cdot}{v}}_{r} \\ - f_{12} - \frac{10 e_{1}^{2}}{3 T^{2}} - \frac{5 e_{2}^{2}}{2 T} + {\overset{\cdot\cdot}{ϕ}}_{r} \end{matrix}] \end{matrix}

(41)

When Equation (34) is substituted into Equation (33), the generalized predictive control law in Equation (41) resembles the input-output control law in Equation (33) if not for the difference in the optimized gain matrices $K (T, ρ)$ . This characterizes the closed-loop stability of the wheelchair.

4.2.1. Closed-loop stability of the system

Although closed-loop stability is one of the NCGPC’s constraints, the corrected control law procedure that guarantees the stability of the closed-loop system can be applied.⁴¹ However, since $ρ_{i} ≱ 4 | \forall i$ , the corrected control law procedure does not apply, but we demonstrate in this section that the asymptotic stability of the closed-loop system in Equation (41) is guaranteed. Representing the coefficients of $e_{1}^{1}, e_{2}^{1}$ and $e_{2}^{2}$ as $K_{11}, K_{21}$ and $K_{22}$ respectively, and with the invertability of the decoupling matrix in the control law in Equation (41), the system in Equation (29) can be expressed as follows:

[\begin{matrix} - K_{11} e_{1}^{1} + {\overset{\cdot}{v}}_{r} - \overset{\cdot}{v} \\ - K_{21} e_{1}^{2} - K_{22} e_{2}^{2} + {\overset{\cdot\cdot}{ϕ}}_{r} - \overset{\cdot\cdot}{ϕ} \end{matrix}] = 0

or simply as:

\begin{matrix} [\begin{matrix} - K_{11} e_{1}^{1} - {\overset{\cdot}{e}}_{1}^{1} \\ - K_{21} e_{1}^{2} - K_{22} {\overset{\cdot}{e}}_{1}^{2} - {\overset{\cdot\cdot}{e}}_{1}^{2} \end{matrix}] = [\begin{matrix} 3 e_{1}^{1} / 2 T + {\overset{\cdot}{e}}_{1}^{1} \\ 10 e_{1}^{2} / 3 T^{2} + 5 {\overset{\cdot}{e}}_{1}^{2} / 2 T + {\overset{\cdot\cdot}{e}}_{1}^{2} \end{matrix}] \\ = 0 \end{matrix}

(42)

The closed-loop stability of the system can be established by letting the reference signals $v_{r} = 0$ and $ϕ_{r} = 0$ since they do not affect the stability of the linear system in Equation (42). The poles of Equation (42) then become $- 1.5 / T$ and $- 1.25 \pm 1.3307 j / T$ respectively, ensuring that $lim_{t \to \infty} {e_{1}^{1}, e_{1}^{2}} = 0$ . System (29) is therefore asymptotically stable given the control law in Equation (41). Poles of the closed-loop system depend linearly on the prediction time T, thus given a small prediction time, the controlled output will rapidly track the reference signal at the expense of large control torques. Control torques can therefore be regulated by appropriate selection of prediction time T.

5. Simulation and results

The simulation results to validate the proposed open-loop and closed-loop dynamic models are elaborated in this section. Various center of mass trajectories and velocities of the open-loop model are presented first, followed by the closed-loop model results. The dynamic model parameters and the default controller gains used in the entire simulation are presented in Table 1.

Table 1.

Dynamic model and controller parameters used in simulation.

Kinematics	$\begin{matrix} \bar{x} = 0.2 & \bar{y} = 0 \\ b = 0.3 m & f = 0.35 m \\ r_{R} = 0.15 m & r_{L} = 0.15 m \end{matrix}$
Dynamics	$\begin{matrix} M = 120 kg & g = 9.81 m / s^{2} \\ k = 100 & μ_{C} = 0.0143 \end{matrix}$
Default	$\begin{matrix} K_{11} = 3 & K_{22} = 13.33 \\ K_{21} = 5 \end{matrix}$
Surface	$μ_{x} = μ_{y} = μ_{z} = μ_{ϕ} = 0.3$
Prediction horizon	$T = 0.5$ s

5.1. The open-loop model

The inputs considered for the open-loop dynamic model are the wheel torques, $τ_{R}$ and $τ_{L}$ . Given that the two identical rear wheels are independently driven by identical dc motors, straight line trajectories are expected on flat non-slippery surfaces with equal $τ_{R}$ and $τ_{L}$ , while circular trajectories should be envisaged with non-equal torques. The simulation results that follow demonstrate the 20 second behavior of the open-loop model with various torque inputs.

5.1.1. A comparison: The model with and without rolling friction

Figure 4 presents the trajectories of point C, generated by $τ_{R} = τ_{L}$ on a flat surface. The simulation results of the model without rolling friction are shown in sub-plots A to C, while the equivalent with rolling friction are in sub-plots D to F. In sub-plots A and B, and correspondingly D and E, it is notable that longer distances are obtained with larger torques, and in sub-plots C and F, that negative torques generate backward trajectories. This is consistent with the expectation. According to Hoffman et al.⁴² and Chua et al.,⁴³ the coefficients of rolling friction of a manual wheelchair vary significantly with the kind of floor/road surface. Typical values include $μ_{W} = 0.0042$ , $μ_{L} = 0.0061$ and $μ_{C} = 0.0143$ , for wooden, linoleum and carpet floors respectively. However, higher coefficients would be expected from the heavier powered wheelchairs with smaller wheel diameters on non-normal outdoor surfaces. Considering a carpet floor with $μ_{C} = 0.0143$ in sub-plots D to F, an error of 20.15% in terms of the distance covered is observable in sub-plot A without the rolling friction as compared to sub-plot D. Figure 5 presents the change rates of $x_{g}, y_{g}$ and $z_{g}$ for trajectories A and C for clearer perception of the model.

Figure 4.

Straight line trajectories of a wheelchair’s center of mass generated by torques $τ_{R} = τ_{L}$ on a flat surface from coordinate [0 0 0] in 20 seconds.

Figure 5.

The velocity curves for trajectories (A) and (C) respectively in Figure 4.

5.1.2. A comparison: The model with and without gravitation effect

The straight line trajectories and the corresponding change rates of $x_{g}, y_{g}$ and $z_{g}$ , generated from an inclined surface with $θ = 7 . 5^{\circ}$ and $ψ = 0^{\circ}$ are shown in Figure 6. In both plots, the initial position is $[0 0 0]$ with $ϕ = 90^{\circ}$ . The wheelchair rolls backwards in sub-plot A when $τ_{R} = τ_{L} = 0$ due to gravitational potential, while in sub-plot B, the torques are high enough to enable frontwards motion. Comparing sub-plots B in Figure 6 with D in Figure 4, in terms the distance covered given the same amount of torque and coefficient of rolling friction, one observes that the gravitational potential accounts for about 250% of the distance lost. A matching comparison between Chénier et al.’s recent model and the presented simulation model is not possible in the absence of the user’s propulsive-moment data.²⁴ However, it may be observed that although the aforementioned model is based on experimental analysis, which makes its results more valid, its evaluation only regards manual propulsion, and does not consider the varying gravitational influence of inclined planes. Neglecting the contributions of both rolling friction and gravitation potential could therefore be inadequate in the dynamic model.

Figure 6.

Straight line trajectories and rates of change of $x_{g}, y_{g}$ and $z_{g}$ generated from an inclined surface form an initial wheelchair position of $[0 0]$ .

Two circular trajectories with their corresponding velocities are depicted in Figure 7. In the first case (case A) the wheelchair turns leftwards due to high $τ_{R}$ as compared to $τ_{L}$ while in the second case (case B) the opposite is true. Although not presented due to space limitation, simulation results show that the smaller the torque difference between the right and left wheels, the larger the radius of circular trajectories and vice versa. Besides this, the circular trajectories of the same radius generated by equal but opposite $τ_{R}$ and $τ_{L}$ on flat surfaces have also been noticed. These are consistent with the expected behavior.

Figure 7.

Circular wheelchair trajectories and rates of change of $x_{g}, y_{g}$ and $z_{g}$ generated on a flat surface from initial position [0 0] and initial direction $ϕ = 0^{\circ}$ with $τ_{R}$ and $τ_{L}$ equal to 4 and 3, and 3 and 4 in the first and the second sub-plot respectively.

The situations exemplified in Figure 8 arise on inclined surfaces, whenever the initial orientation is neither directly up nor directly down the slope. With $τ_{R} = τ_{L} = 3$ when the initial orientation is perpendicular to the direction of the slope, the S-shaped trajectory depicted in sub-plot A results. The left wheel being raised compared to the right has higher potential energy, which due to the force effect results in faster rotations. The wheelchair thus turns and a curved trajectory is observed. In the new direction, the right wheel is raised. Having a higher potential energy, it rotates faster than the left resulting in the second turning effect. This repetitive process generates the S-shaped trajectory. In sub-plot B, a larger torque is supplied to the right wheel to generate a counter-clockwise circular trajectory. The high torque on the right wheel overcomes the stronger gravitational effect on the raised wheel, the right wheel therefore rotates faster and a circular trajectory begins to form. However, when the right is raised compared to the left, its gravitational force together with its high torque creates a faster left turning effect and the wheelchair turns earlier. This generates the moving circular trajectory which seems to retard as the forces tend to equilibrium.

Figure 8.

Trajectories and velocities generated when initial wheelchair orientation is neither directly up nor directly down the slope. The simulation has been conducted on a surface inclined by $θ = 15^{\circ}$ and $ψ = 0^{\circ}$ from an initial $[0 0]$ wheelchair position.

The explanation given for sub-plots A in Figures 8 and 6 applies to first and second sub-plots in Figure 9 respectively. However it is important to notice the trajectory flip that resulted in the initial $ψ = 90^{\circ}$ .

Figure 9.

The trajectory observed due to $ψ$ on a surface inclined by $θ = 15^{\circ}$ and $ψ = 90^{\circ}$ from an initial $[0 0]$ wheelchair position.

5.1.3. Simulation with slipping disturbance

By introducing a random slipping velocity, the deviations depicted in Figure 10 on a flat surface can be observed. However, further experimental assessment is still necessary, with the use of encoders in the estimation of slipping velocity according the proposed model. This is critical for wheelchair controllability and users’ safety. Although the simulation results of the open-loop model are largely consistent with the expected behavior, it is important to acknowledge the singularity effects of the time derivatives of Euler angles. The analysis and interpretation of simulation results involving such derivatives could sometimes be very complex. As a result, only trivial cases with $\overset{\cdot}{ψ} = 0$ , $ψ = 90^{\circ}$ or $ψ = 0^{\circ}$ have been considered in the analysis. The model however takes the derivatives into consideration.

Figure 10.

Resulting deviation from the intended trajectory on a flat surface with slight slip introduced into the model.

5.2. The closed-loop model

In the closed-loop model simulations, the reference linear velocity $v_{r}$ and angular orientation $ϕ_{r}$ are considered as the inputs. The controller then regulates the torques depending on the navigational situation in order to track the reference signals. Two simulation results are presented to this effect, and a consideration is made in each case to the slippery condition. In Figure 11, by considering a ramp of unit slope as the reference angular orientation with $v_{r} = 1.5$ at $θ = 15^{\circ}$ and $ψ = 90^{\circ}$ , the circular wheelchair trajectory shown in the first sub-plot is generated. During the simulation, a random slip is introduced at time $t = 20$ s for the rest of the simulation time, as depicted in the time series curve for wheel torques, and the controller automatically adjusts the torques to track the specified user inputs. The time series curves for linear velocity and angular orientation error illustrate how well the controller tracks the desired inputs in spite of the slipping and gravitational effects. Similar results are depicted in Figure 12, where the controller tracks a sinusoidal input of unit amplitude, to produce the desired sine wave trajectory. Although a tracking error is evident in the orientation, it is a steady state with a smaller magnitude compared to the desired output. Therefore, the simulation results of the closed-loop model can also be considered consistent with the expectation.

Figure 11.

Circular wheelchair trajectory generated by considering a ramp input for reference angular orientation and $v_{r} = 1.5$ at $θ = 0^{\circ}$ and $ψ = 90^{\circ}$ . As depicted in the time series curve for wheel torques, random slip introduced at time $t = 20$ s for the rest of simulation time does not affect the ability of controller to automatically adjust the torques in order to track the specified user inputs.

Figure 12.

Sinusoidal wheelchair trajectory generated by considering a sine wave input for reference angular orientation and $v_{r} = 1.5$ at $θ = 0^{\circ}$ and $ψ = 90^{\circ}$ . Similarly, the random slip introduced at time $t = 20$ s does not affect the ability of controller to regulate wheel torques.

5.3. Comparison with other differential drive models

Table 2 presents a comparison of the presented model with other previous models, based on model comprehensiveness and the employed modeling approach. The comparison is limited to differential drive structures of wheeled-mobile systems with two front castor wheels.

Table 2.

A comparison of the presented wheelchair model with other wheelchair models.

Study	Dynamic	Configuration surface		Frictional effects	Slip detection
Study	Dynamic	Inclined	Non-inclined	Frictional effects	Slip detection
Johnson and Aylor⁴⁴	✓		✓	✓
Tashiro and Murakami⁴⁵	✓		✓
Wang et al.⁴⁶	✓	✓	✓
Onyango et al.²²	✓	✓	✓
Onyango et al.²⁰	✓	✓	✓	✓
Emam et al.¹⁹	✓		✓		✓
Chénier et al.²³	✓		✓	✓
Chénier et al.²⁴	✓		✓	✓
The presented model	✓	✓	✓	✓	✓

6. Conclusions

This study intended to develop a comprehensive dynamic model that takes the effects of gravitational forces on inclined and non-inclined surfaces as well as the contributions of rolling friction into consideration. This also involved estimating the slipping parameters and formulating the slipping velocities. A new dynamic model is therefore presented. With good simulation results, it is believed that the dynamic model gives a better representation of a real wheelchair, usable under non-normal indoor and outdoor driving conditions. The introduced method of wheel slip detection, with a reduced number of slip detection encoders, is also considered a simple and cost effective solution. The dynamic model is therefore presumed to be comprehensive enough, for testing and validating the wheelchair controllers intended for safety and steering-ease improvements.

7. Recommendations for future work

It is assumed in this study that through the available user interface, the wheelchair user is able to communicate clearly the navigational commands required to make the wheelchair move, stop or turn. Delays in such commands could sometimes be very dangerous to the user. It may be important therefore to consider incorporating assistive systems into the wheelchair, to assist users with steering difficulties. This may involve modeling the driver’s steering behavior, and integrating the steering model in the closed-loop control, to act as a real-time co-driver that predicts and provides on-line solutions to the driver’s possible errors. The authors are currently working on these assistance recommendations.

Footnotes

Appendix 1

The matrix $B$ below is the orthogonal transformation matrix of direction cosines used in Equation (3) to transform the wheelchair velocities in the inertial frame ${x, y, z}$ to the body fixed frame ${X, Y, Z}$ :

(43)

B = [\begin{matrix} \begin{matrix} \cos ϕ \cos ψ - \\ \sin ϕ \sin ψ \cos θ \end{matrix} & \begin{matrix} \cos ϕ \sin ψ + \\ \sin ϕ \cos ψ \cos θ \end{matrix} & \sin θ \sin ϕ \\ \begin{matrix} - \sin ϕ \cos ψ - \\ \cos ϕ \sin ψ \cos θ \end{matrix} & \begin{matrix} - \sin ϕ \sin ψ + \\ \cos ϕ \cos ψ \cos θ \end{matrix} & \sin θ \cos ϕ \\ \sin θ \sin ψ & - \sin θ \cos ψ & \cos θ \end{matrix}]

The applied force on the driving wheels as a result of the motor torques can be obtained by $E (q) τ$ where:

(44)

E (q) = \frac{1}{r} [\begin{matrix} B_{11} & B_{11} \\ B_{12} & B_{12} \\ B_{13} & B_{13} \\ b & - b \end{matrix}] τ = [\begin{matrix} τ_{R} \\ τ_{L} \end{matrix}]

The dynamic equations of motion are computed in accordance with the Lagrangian formalism according to the following procedure given the Lagrangian function in Equation (9):

(45)

L = T - U = \frac{1}{2} \sum_{i, j}^{n} m_{ij} (q) {\overset{\cdot}{q}}_{i} {\overset{\cdot}{q}}_{j} - U (q)

Applying Equation (1), the equations of motion may be computed in the following manner:

(46)

\frac{\partial L}{\partial {\overset{\cdot}{q}}_{k}} = \sum_{j} m_{kj} {\overset{\cdot}{q}}_{j}

(47)

\begin{matrix} \frac{d}{dt} (\frac{\partial L}{\partial {\overset{\cdot}{q}}_{k}}) = \sum_{j} m_{kj} {\overset{\cdot\cdot}{q}}_{j} + \sum_{j} \frac{d}{dt} (m_{kj} {\overset{\cdot}{q}}_{j}) \\ = \sum_{j} m_{kj} {\overset{\cdot\cdot}{q}}_{j} + \sum_{i, j} \frac{\partial m_{kj}}{\partial q_{i}} {\overset{\cdot}{q}}_{i} {\overset{\cdot}{q}}_{j} \end{matrix}

(48)

\frac{\partial L}{\partial q_{k}} = \frac{1}{2} \sum_{i, j} \frac{\partial m_{ij}}{\partial q_{k}} {\overset{\cdot}{q}}_{i} {\overset{\cdot}{q}}_{j} - \frac{\partial U}{\partial q_{k}}

The Euler-Lagrange equations of motion may then be summarized by putting Equations (46–48) together:

(49)

\sum_{j} m_{kj} {\overset{\cdot\cdot}{q}}_{j} + \sum_{i, j} {\frac{\partial m_{kj}}{\partial q_{i}} - \frac{1}{2} \frac{\partial m_{ij}}{\partial q_{k}}} {\overset{\cdot}{q}}_{i} {\overset{\cdot}{q}}_{j} - \frac{\partial U}{\partial q_{k}} = F

In matrix form, the first, second and third terms of Equation (49) may correspondingly be represented by $\bar{M} (q) \overset{\cdot\cdot}{q}$ , $\bar{C} (q, \overset{\cdot}{q}) \overset{\cdot}{q}$ and $G (q)$ in Equations (50) to (51).

(50)

\bar{M} (q) = [\begin{matrix} M & 0 & 0 & - M \bar{x} B_{21} \\ 0 & M & 0 & - M \bar{x} B_{22} \\ 0 & 0 & M & - M \bar{x} B_{23} \\ - M \bar{x} B_{21} & - M \bar{x} B_{22} & - M \bar{x} B_{23} & I_{ZZ} \end{matrix}]

(51)

\begin{matrix} \bar{C} (q, \overset{\cdot}{q}) = \\ [\begin{matrix} 0 & 0 & 0 & M \bar{x} (\overset{\cdot}{ϕ} B_{11} + 2 \overset{\cdot}{ψ} B_{22} - 2 \overset{\cdot}{θ} \sin ψ \sin θ \cos ϕ) \\ 0 & 0 & 0 & M \bar{x} (\overset{\cdot}{ϕ} B_{12} - 2 \overset{\cdot}{ψ} B_{21} + 2 \overset{\cdot}{θ} \sin θ \cos ϕ \cos ψ) \\ 0 & 0 & 0 & M \bar{x} (\overset{\cdot}{ϕ} B_{13} - 2 \overset{\cdot}{θ} \cos ϕ \cos θ) \\ 0 & 0 & 0 & 0 \end{matrix}], \\ G (q) = [\begin{matrix} 0 \\ 0 \\ Mg \\ 0 \end{matrix}] \end{matrix}

(52)

Λ_{ss} = {[\begin{matrix} - \frac{({T_{1}}^{3} - {T_{2}}^{3})}{3} & - \frac{({T_{1}}^{5} - {T_{2}}^{5})}{20} \\ - \frac{({T_{1}}^{3} - {T_{2}}^{3})}{3} & \frac{({T_{1}}^{5} - {T_{2}}^{5})}{20} \end{matrix}]}^{- 1}

(53)

\begin{matrix} Λ_{s} = \\ [\begin{matrix} - \frac{({T_{1}}^{3} - {T_{2}}^{3})}{6} & - \frac{({T_{1}}^{2} - {T_{2}}^{2})}{6} & - \frac{({T_{1}}^{4} - {T_{2}}^{4})}{8} & - \frac{({T_{1}}^{3} - {T_{2}}^{3})}{3} & - \frac{({T_{1}}^{5} - {T_{2}}^{5})}{20} \\ \frac{({T_{1}}^{3} - {T_{2}}^{3})}{6} & - \frac{({T_{1}}^{2} - {T_{2}}^{2})}{6} & \frac{({T_{1}}^{4} - {T_{2}}^{4})}{8} & - \frac{({T_{1}}^{3} - {T_{2}}^{3})}{3} & \frac{({T_{1}}^{5} - {T_{2}}^{5})}{20} \end{matrix}] \end{matrix}

In order to compute the kinematic model in Equation (11) without slip, the full rank transformation matrix $S (q) \in R^{n \times (n - m)}$ is required to transform the velocity vector $η \in R^{(n - m) \times 1}$ to the generalized velocity vector $\overset{\cdot}{q}$ :

(54)

S (q) = [\begin{matrix} B_{11} & \bar{x} B_{21} \\ B_{12} & \bar{x} B_{22} \\ B_{13} & \bar{x} B_{23} \\ 0 & 1 \end{matrix}] η = [\begin{matrix} v \\ \overset{\cdot}{ϕ} \end{matrix}]

Estimation of linear velocities ${\vec{v}}_{R}$ and ${\vec{v}}_{L}$ for the right and left castor wheels respectively, is considered with respect to the two points and axis of rotations A and B in Figures 3(a) and 3(b) respectively. The resultant linear velocity pair due to rotations about point A is presented by vectors ${\vec{v}}_{R_{A}}$ and ${\vec{v}}_{L_{A}}$ , while the pair due to rotations about point B is presented by vectors ${\vec{v}}_{R_{B}}$ and ${\vec{v}}_{L_{B}}$ for the right and left castor wheel respectively:

(55)

{\vec{v}}_{R_{A}} = [\begin{matrix} {\overset{\cdot}{α}}_{R} e (\sin α_{R} \cos ψ - \cos α_{R} \sin ψ \cos θ) \\ {\overset{\cdot}{α}}_{R} e (\sin α_{R} \sin ψ - \cos α_{R} \cos ψ \cos θ) \\ {\overset{\cdot}{α}}_{R} e (\cos α_{R} \sin θ) \\ 0 \end{matrix}]

(56)

{\vec{v}}_{R_{B}} = [\begin{matrix} r_{c} {\overset{\cdot}{γ}}_{f_{R}} (\cos α_{R} \cos ψ - \sin α_{R} \sin ψ \cos θ) \\ r_{c} {\overset{\cdot}{γ}}_{f_{R}} (\cos α_{R} \sin ψ + \sin α_{R} \cos ψ \cos θ) \\ r_{c} {\overset{\cdot}{γ}}_{f_{R}} (\sin α_{R} \sin θ) \\ 0 \end{matrix}]

(57)

{\vec{v}}_{L_{A}} = [\begin{matrix} {\overset{\cdot}{α}}_{L} e (\sin α_{L} \cos ψ - \cos α_{L} \sin ψ \cos θ) \\ {\overset{\cdot}{α}}_{L} e (\sin α_{L} \sin ψ - \cos α_{L} \cos ψ \cos θ) \\ {\overset{\cdot}{α}}_{L} e (\cos α_{L} \sin θ) \\ 0 \end{matrix}]

(58)

{\vec{v}}_{L_{B}} = [\begin{matrix} r_{c} {\overset{\cdot}{γ}}_{f_{L}} (\cos α_{L} \cos ψ - \sin α_{L} \sin ψ \cos θ) \\ r_{c} {\overset{\cdot}{γ}}_{f_{L}} (\cos α_{L} \sin ψ + \sin α_{L} \cos ψ \cos θ) \\ r_{c} {\overset{\cdot}{γ}}_{f_{L}} (\sin α_{L} \sin θ) \\ 0 \end{matrix}]

Funding

The authors of this paper gratefully appreciate the contribution of South African National Research Fund (NRF), Tshwane University of Technology (TUT) and the French South African Institute of Technology (F’SATI) for providing all the relevant and necessary support for this research.

Author Biographies

Stevine Obura Onyango obtained his bachelor’s degree from Egerton University, Njoro, Kenya in 2007, his master of technology degree in electrical engineering from Tshwane University of Technology (TUT), Pretoria, South Africa, in 2010, and his master of science degree in electrical engineering from Ècole Supèrieure d’Ingènieurs en Èlectronique et Èlectrotechnique (ESIEE Paris), Paris, France, in 2011. He is currently a PhD student in the Electrical Engineering Department at Tshwane University of Technology (TUT) at the French South African Institute of Technology (F’SATI), under the supervision of Professor Yskander Hamam, Professor Karim Djouani and Professor Boubaker Daachi. His current research focuses on behavior modeling and system control with human-in-the-loop.

Yskandar Hamam received his bachelor’s degree from the American University of Beirut (AUB) in 1966. He obtained his MSc in 1970 and PhD in 1972 from the University of Manchester Institute of Science and Technology. He also obtained his “Diplôme d’Habilitation á Diriger des Recherches” (equivalent to DSc) from the Université des Sciences et Technologies de Lille in 1998. He has conducted research activities and lectured in England, Brazil, Lebanon, Belgium and France. He was the head of the Control Department and Dean of Faculty at ESIEE, France. He was an active member in modeling and simulation societies and was the president of EUROSIM. He was the scientific director of the French South African Institute of Technology (FSATI) at TUT in South Africa from 2007 to 2012. He is currently professor at the Department of Electrical Engineering of TUT. He has co-authored 3 books. He has authored/co-authored more than 300 papers in archival journals and conference proceedings as well as book contributions. He is a life senior member of the IEEE.

Karim Djouani is a full professor, scientist and technical group supervisor of pattern recognition, soft computing, networking systems and robotics at University Paris Est Créteil (UPEC), France. Since January 2014, he has been the DST/NRF SARChI Chair in Enabled Environment for Assistive Living (). He is leading the SIRIUS team of the LISSI lab, University Paris Est. Since June 2012, he has been the deputy director in charge of research at the IUT-CV (University Institute of Technology) at University Paris Est-Creteil (UPEC). Since January 2011, he has been full professor at the French South African Institute of Technology (FSATI) at Tshwane University of Technology (TUT), Pretoria, South Africa. From July 2008 to December 2010, he was seconded by the French Ministry of Higher Education to the FSATI at TUT, Pretoria, South Africa. Until July 2008 he was in charge of the management of the national and European projects at the LISSI Lab. His current research work focuses on the development of novel and highly efficient algorithms for reasoning systems with uncertainty as well as optimization for distributed systems, networked control systems, wireless ad-hoc network, wireless and mobile communication, and wireless sensors networks as well as robotics. He has authored/co-authored about 200 articles in archival journals and conference proceedings as well as 18 chapters in edited books and 2 books.

Boubaker Daachi is a full professor of computer science at the University of Paris 8. He received his BS degree in computer science from the University of Sétif, Algeria (1995), his MS degree in robotics from the University of Paris 6, France (1998), and his PhD in robotics from the University of Versailles, France (2000). He was an associate professor of robotics and computer science at the University of Paris Est Créteil from September 2003 to August 2014. He evolved as a researcher at CNRS-AIST JRL, Tsukuba, Japan, from September 2014 to August 2015. His research interests include adaptive control of complex systems, soft computing, and BCI. His application fields are robotics and pervasive and distributed systems. He has published more than 60 papers in scientific journals, books, and conference proceedings. He has been involved in the organizing committees of some national and international events.

References

Tarokh

McDermott

Kinematics modeling and analyses of articulated rovers. IEEE T Robot 2005; 21(4): 539–553.

Zhu

Dong

. Robust tracking control of wheeled mobile robots not satisfying nonholonomic constraints. In: International conference on intelligent systems design and applications, Jinan, China, 16–18 October 2006, pp.643–648.

Morales

Feliu

González

. Kinematic model of a new staircase climbing wheelchair and its experimental validation. Int J Robot Res 2006; 25(9): 825–841.

Tian

Shirinzadeh

Zhang

. Design and forward kinematics of the compliant micro-manipulator with lever mechanisms. Precis Eng 2009; 33(4): 466–475.

Zhang

Liu

Kinematics analysis of a novel parallel manipulator. Mech Mach Theory 2009; 44(9): 1648–1657.

Wretstrand

Petzäll

Ståhl

Safety as perceived by wheelchair-seated passengers in special transportation services. Accident Anal Prev 2004; 36(1): 3–11.

Weiss

Frohlich

Zell

Vibration-based terrain classification using support vector machines. In: IEEE/RSJ international conference on intelligent robots and systems, Beijing, China, 9–15 October 2006, pp.4429–4434.

Fierro

Lewis

FL.

Control of a nonholomic mobile robot: Backstepping kinematics into dynamics. J Robotic Syst 1997; 14(3): 149–163.

Oubbati

Schanz

Levi

. Kinematic and dynamic adaptive control of a nonholonomic mobile robot using a RNN. In: IEEE international symposium on computational intelligence in robotics and automation, Espoo, Finland, 27–30 June 2005, pp.27–33.

10.

Kozlowski

Pazderski

Modeling and control of a 4-wheel skid-steering mobile robot. Int J Appl Math Comp 2004; 14(4): 477–496.

11.

Chen

Huang

. 3D dynamical analysis for a caster wheeled mobile robot moving on the frictional surface. In: IEEE/RSJ international conference on intelligent robots and systems, Beijing, China, 9–15 October 2006, pp.3056–3061.

12.

Stonier

Cho

Choi

. Nonlinear slip dynamics for an omniwheel mobile robot platform. In: IEEE international conference on robotics and automation, Rome, Italy, 10–14 April 2007, pp.2367–2372.

13.

Williams

Carter

Gallina

. Dynamic model with slip for wheeled omnidirectional robots. IEEE T Robot 2002; 18(3): 285–293.

14.

Olson

Shaw

Stépán

Nonlinear dynamics of vehicle traction. Vehicle Syst Dyn 2003; 40(6): 377–399.

15.

Sidek

Sarkar

. Dynamic modeling and control of nonholonomic mobile robot with lateral slip. In: Third international conference on systems, Cancun, Mexico, 13–18 April 2008, pp.35–40.

16.

Katsura

Ohnishi

Human cooperative wheelchair for haptic interaction based on dual compliance control. IEEE T Ind Electron 2004; 51(1): 221–228.

17.

Nguyen

. Advanced robust tracking control of a powered wheelchair system. In: Engineering in medicine and biology society. EMBS 2007. 29th annual international conference of the IEEE, Lyon, France, 22–26 August 2007, pp.4767–4770.

18.

Cruz

CDL

Bastos

Carelli

Adaptive motion control law of a robotic wheelchair. Control Eng Pract 2011; 19(2): 113–125.

19.

Emam

Hamam

Monacelli

. Power wheelchair driver behaviour modelling. In: 7th international multi-conference on systems signals and devices (SSD), Amman, Jordan, 27–30 June 2010, pp.1–7.

20.

Onyango

Hamam

Djouani

. Dynamic control of powered wheelchair with slip on an incline. In: 2nd international conference on adaptive science technology, ICAST, Accra, Ghana, 14–16 Januray 2009, pp.278–283.

21.

Onyango

Hamam

Djouani

. Velocity and orientation control in an electrical wheelchair on an inclined and slippery surface. In: Proceedings of the 8th international conference on informatics in control, automation and robotics (ICINCO), Noordwijkerhout, the Netherlands, 28–31 July 2011, pp.112–119.

22.

Onyango

Hamam

Dabo

. Dynamic control of an electrical wheelchair on an incline. In: Proceedings of IEEE AFRICON, Nairobi, Kenya, 23–25 September 2009, pp.1–6.

23.

Chénier

Bigras

Aissaoui

. A new dynamic model of the manual wheelchair for straight and curvilinear propulsion. In: IEEE international conference on rehabilitation robotics (ICORR), Zurich, Switzerland, 29 June–1 July 2011, pp.1–5.

24.

Chénier

Bigras

Aissaoui

A new dynamic model of the wheelchair propulsion on straight and curvilinear level-ground paths. Comput Method Biomech 2015; 18(10): 1031–1043.

25.

Bloch

McClamroch

. Control of mechanical systems with classical nonholonomic constraints. In: Proceedings of the 28th IEEE conference on decision and control, Tampa, FL, USA, 13–15 December 1989, volume 1, pp.201–205.

26.

Campion

d’Andrea Novel

Bastin

Controllability and state feedback stabilizability of non holonomic mechanical systems. In: Canudas de Wit

(ed) Advanced Robot Control, Lecture Notes in Control and Information Sciences. Volume 162. Berlin Heidelberg: Springer, 1991, pp.106–124.

27.

Isidori

Nonlinear Control Systems II (Communications and Control Engineering). London: Springer, 1999.

28.

Yang

Wang

Adaptive tracking control for uncertain dynamic nonholonomic mobile robots based on visual servoing. J Contr Theor Appl 2012; 10(1): 56–63.

29.

Spong

Vidyasagar

Robot dynamics and control. Hoboken, NJ, USA: John Wiley and Sons, 2008.

30.

Sarkar

Yun

Kumar

Control of mechanical systems with rolling constraints application to dynamic control of mobile robots. Int J Robot Res 1994; 13(1): 55–69.

31.

Lewis

Vandegrift

Flexible robot arm control by a feedback linearization/singular perturbation approach. In: IEEE international conference on robotics and automation, Atlanta, GA, USA, 2–6 May 1993, volume 3, pp.729–736.

32.

d’Andréa Novel

Campion

Bastin

. Control of nonholonomic wheeled mobile robots by state feedback linearization. Int J Robot Res 1995; 14(6): 543–559.

33.

Kim

JH.

Tracking control of a two-wheeled mobile robot using inputoutput linearization. Control Eng Pract 1999; 7(3): 369–373.

34.

Caracciolo

De Luca

Iannitti

. Trajectory tracking control of a four-wheel differentially driven mobile robot. In: IEEE international conference on robotics and automation, Detroit, MI, USA, 10–15 May 1999, volume 4. pp.2632–2638.

35.

Oriolo

De Luca

Vendittelli

WMR control via dynamic feedback linearization: Design, implementation, and experimental validation. IEEE T Contr Syst T 2002; 10(6): 835–852.

36.

Spong

Hutchinson

Vidyasagar

Robot modeling and control. Hoboken, NJ: John Wiley and Sons, 2006.

37.

Taghirad

Parallel Robots: Mechanics and Control. Boca Raton, FL, USA: Taylor and Francis, 2013.

38.

Kelly

Davila

Perez

Control of Robot Manipulators in Joint Space. London: Springer, 2006.

39.

Chénier

Bigras

Aissaoui

An orientation estimator for the wheelchair’s caster wheels. IEEE T Contr Syst T 2011; 19(6): 1317–1326.

40.

Yun

Kumar

Sarkar

. Control of multiple arms with rolling constraints. In: IEEE international conference on robotics and automation, Nice, France, 12–14 May 1992, volume 3, pp.2193–2198.

41.

Dabo

Chafouk

Langlois

. Unconstrained NCGPC with a guaranteed closed-loop stability: Case of nonlinear siso systems with the relative degree greater than four. In: Proceedings of the 48th IEEE conference on decision and control, Shanghai, China, 15–18 December 2009, pp.1980–1985.

42.

Hoffman

Millet

Hoch

. Assessment of wheelchair drag resistance using a coasting deceleration technique. Am J Phys Med Rehab 2003; 82(11): 880–889.

43.

Chua

Fuss

Subic

Rolling friction of a rugby wheelchair. Procedia Eng 2010; 2(2): 3071–3076.

44.

Johnson

Aylor

Dynamic modeling of an electric wheelchair. IEEE T Ind Appl 1985; IA-21(5): 1284–1293.

45.

Tashiro

Murakami

Step passage control of a power-assisted wheelchair for a caregiver. IEEE T Ind Electron 2008; 55(4): 1715–1721.

46.

Wang

Salatin

Grindle

. Real-time model based electrical powered wheelchair control. Med Eng Phys 2009; 31(10): 1244–1254.

Modeling a powered wheelchair with slipping and gravitational disturbances on inclined and non-inclined surfaces

Abstract

Keywords

1. Introduction

2. Dynamic model with gravitational forces

2.1. Description of the wheelchair platform

2.2. System constraints

2.3. Kinetic and potential energy

2.4. Dynamic model development

3. Slipping parameters and frictional force

3.1. Slipping parameters

3.2. Determination of real velocity

3.3. Frictional and resistive force modeling

4. The closed-loop design

4.1. Input-output feedback linearization

4.1.1. Navigation task

4.1.2. System relative degree ( ρ )

4.1.3. The control law

4.2. Non-linear continuous-time generalized predictive control

4.2.1. Closed-loop stability of the system

5. Simulation and results

5.1. The open-loop model

5.1.1. A comparison: The model with and without rolling friction

5.1.2. A comparison: The model with and without gravitation effect

5.1.3. Simulation with slipping disturbance

5.2. The closed-loop model

5.3. Comparison with other differential drive models

6. Conclusions

7. Recommendations for future work

Footnotes

Appendix 1

Funding

Author Biographies

References

4.1.2. System relative degree ( $ρ$ )