Hybrid nonlinear model predictive and linear quadratic balance control for a two-wheeled self-balancing wheelchair

A two-wheeled self-balancing wheelchair is a highly nonlinear and underactuated wheeled inverted-pendulum system, which faces severe balance control challenges under large tilt angles, load variations, and external disturbances. A hybrid nonlinear model predictive control–linear quadratic regulator (NMPC–LQR) balance control scheme is presented that combines nonlinear optimal recovery and efficient local stabilization while reducing the computational burden. A reduced-order nonlinear dynamic model of the wheelchair–occupant system is derived via the Lagrange method. The proposed controller confines LQR to a locally valid region verified by linearization-error analysis, whereas NMPC is employed for large-deviation recovery. A linearization-error analysis is first used to determine the switching-threshold search range, and a normalized five-index criterion identifies 0.100 rad as the optimal nominal threshold. Mass-sensitivity tests further provide an empirical relation for estimating the optimal threshold within the tested mass range. Comparative simulations under initial tilt angles of 35, 40, and 45° show faster convergence than standalone LQR and standalone NMPC controllers. Benchmark comparisons with SMC demonstrate a larger recoverable initial-angle range. Robustness tests under parameter variations, pulse and step disturbances, and flat-ground friction changes, together with prototype experiments, validate the stability, disturbance rejection, and engineering applicability of the proposed control strategy.

1 Introduction

Two-wheeled self-balancing wheelchairs are becoming increasingly popular due to their compact structure, greater maneuverability in narrow spaces, and better adaptability to complex indoor and outdoor environments (Zhang et al., 2024). They represent a promising solution for improving the mobility of individuals with lower-limb impairment during daily living. However, such wheelchairs must not only provide mobility but also ensure stability and safety under various operating conditions to be truly useful and acceptable (Onozuka et al., 2020). Consequently, stability and safety have become the two key concerns for both researchers and users. The inherent static instability of self-balancing wheelchairs necessitates high gain feedback of fast state estimation, which inevitably amplifies sensor noise and model errors. Kalman filtering has been widely adopted to mitigate this sensitivity, as demonstrated by Huynh et al. (2017) in an optimal control framework, by Zhang et al. (2019) within a PID structure, and by Ma et al. (2024) using a fuzzy Kalman approach. However, these methods share three critical weaknesses: (i) reliance on known model and noise statistics, which do not hold under load variations; (ii) heuristic parameter tuning without online adaptation; and (iii) a decoupled design of estimation and control, leading to degraded performance under large-angle deviations or external disturbances. Thus, developing an adaptive or robust co-design of state estimation and feedback control remains an open challenge.

Model-based controllers have been explored to overcome the poor adaptability existing in conventional methods. Among them, the LQR controller has been widely used in self-balancing systems due to its simple structure and satisfactory local stabilization performance. For instance, Feng et al. (2023) applied LQR to a two-wheel-legged robot and verified its efficient and stable regulation capability near the equilibrium point. Owing to its linear nature, LQR performs well near the equilibrium point but deteriorates rapidly under large-angle deviations. Recent studies have attempted to improve the robustness of LQR by introducing online self-tuning or hierarchical reconfiguration mechanisms. For example, Saleem et al. (2021) proposed several hierarchical self-tuning LQR procedures in which the weighting factors of the quadratic performance index are adjusted online to improve disturbance rejection in under-actuated mechatronic systems. Even with friction feedforward compensation (Zhang et al., 2025), its effective operating region remains fundamentally constrained. Compared with regular self-balancing vehicles, self-balancing wheelchairs are characterized by larger-scale and higher inertia, making safety an overriding concern in control design. Under sudden external disturbances, the system must be stabilized rapidly to avoid overturning while maintaining acceptable recovery performance even at relatively large inclination angles. This inevitably leads to strongly nonlinear dynamic behavior. Analogous control difficulties have been observed in two-wheeled inverted-pendulum systems subject to uncertain rolling resistance and angle limits (Liang et al., 2021). Consequently, the sole use of LQR control is insufficient to achieve satisfactory performance for self-balancing wheelchairs. Hence, the design of advanced control schemes capable of dealing with severe system nonlinearity still poses a considerable challenge.

Predictive, adaptive, and intelligent control approaches have attracted increasing attention to cope with nonlinear behavior during large-angle recovery (Zhang and Mohamad Nor, 2025). For instance, Cui et al. (2015) and Li and Zhang (2010) employed adaptive backstepping and motion/force control strategies, respectively, for wheeled inverted-pendulum systems. Shen and Hong (2022) proposed an MPC strategy for a two-wheeled inverted-pendulum system and demonstrated its capability to handle nonlinear effects through online control. However, its effectiveness under strongly nonlinear large-angle recovery conditions remains insufficiently validated. Liu et al. (2024) introduced NMPC into the control of a self-balancing wheelchair to further improve control performance under highly nonlinear conditions, while Li et al. (2021) demonstrated the applicability of predictive control to nonholonomic wheel-legged robots. These results indicate that predictive control methods are effective in handling nonlinear system dynamics. Nevertheless, it is important to emphasize that model predictive controllers still rely on online optimization at each sampling instant, thereby imposing a substantial computational burden. In addition, the control effectiveness in the near-equilibrium region is not always superior to that of traditional control approaches. An NMPC–LQR hybrid controller is proposed to improve large-angle recovery while reducing the computational burden. When the self-balancing system experiences a relatively large tilt angle, the NMPC controller is activated to provide the control input. As the tilt angle decreases and enters a predefined threshold region, the control law is switched to the LQR controller. The proposed strategy therefore combines the nonlinear recovery capability of NMPC under large-tilt-angle conditions with the computational efficiency and regulation performance of LQR in the near-equilibrium region. Furthermore, a switching-threshold evaluation method is also designed to resolve the problem of switching-angle selection. Subsequent simulations and experiments verify the effectiveness and practical applicability of the proposed hybrid controller for engineering implementation.

The main contributions are summarized as follows:

1. A reduced-order wheelchair–occupant dynamic model is established for longitudinal balance recovery, and the validity range of the modeling assumptions is explicitly discussed.

2. A hybrid NMPC–LQR balance controller is constructed, where NMPC is used for large-deviation nonlinear recovery and LQR is restricted to the local region supported by linearization-validity analysis.

3. A switching-threshold selection procedure is developed by combining linearization-error analysis and a normalized five-index evaluation criterion. In addition, a mass-sensitivity study is conducted to obtain an empirical threshold-estimation relation within the tested equivalent-mass range.

4. Comparative simulations, SMC benchmark tests, disturbance tests, friction-coefficient tests, computational-efficiency tests, and prototype experiments are conducted to evaluate the proposed method.

However, three gaps remain:

1. switching thresholds in existing hybrid schemes are often selected empirically, and their searchable range or parameter-dependent estimation is rarely quantified;

2. the validity range of LQR linearization is rarely quantified when large initial tilt angles are considered;

3. limited studies validate wheelchair-oriented hybrid control under load variations, friction mismatch, external disturbances, and prototype experiments.

The remainder of this paper is organized as follows. Section 2 presents the reduced-order dynamic model of the wheelchair–occupant system and discusses the validity range of the modeling assumptions. Section 3 introduces the hybrid NMPC–LQR controller, the linearization-error-based threshold search range, the five-index threshold evaluation method, the mass-dependent threshold estimation, and the stability considerations. Section 4 reports comparative simulations, SMC benchmark tests, pulse- and step-disturbance tests, friction-coefficient tests, computational-efficiency evaluations, and prototype experiments. Section 5 concludes the paper and discusses the limitations and future research directions.

2 Modeling

Since both LQR and NMPC are model-based control methods, the accuracy of the dynamic model is of great importance. The wheelchair system is simplified as an inverted-pendulum system, following the classical modeling framework commonly adopted for nonholonomic two-wheeled inverted-pendulum robots (Kim et al., 2005). The main symbols and physical parameters used in the dynamic model are summarized in Table 1. The three-dimensional wheelchair–occupant configuration and the planar reduced model used in the subsequent derivation are shown in Fig. 1.

Table 1 Main symbols and physical parameters of the wheelchair–occupant system.

Figure 1 Dynamic modeling configuration (left) and equivalent mass distribution (right).

The following assumptions are adopted:

• Item 1. Only the motion of the wheelchair system in the longitudinal plane is considered, while yaw, roll, and left–right asymmetry are neglected.

• Item 2. The two driving wheels are assumed to have identical structural and inertial parameters.

• Item 3. Pure rolling without slipping is assumed between the driving wheels and the ground.

• Item 4. The relative posture between the wheelchair body and the occupant is assumed unchanged and is treated as a rigidly connected whole.

• Item 5. The effects of tire elastic deformation and local ground unevenness are neglected.

The above assumptions are adopted to obtain a tractable reduced-order model for real-time controller design. Therefore, the proposed model mainly describes longitudinal balance recovery and pitch stabilization during straight-line motion. Yaw, roll, lateral motion, tire deformation, wheel slip, and active occupant motion are not explicitly considered. The rigid-body assumption between the wheelchair and the occupant is reasonable when the occupant maintains a relatively fixed seated posture, but it may introduce modeling errors when the occupant posture or mass distribution changes. In addition, friction and external forces are not included in the nominal prediction model but are treated as unmodeled dynamics and external disturbances. Their influence is further evaluated through parameter-variation and disturbance-rejection simulations.

The two driving wheels are connected by a light linkage, and the total mass of the system is represented by an equivalent concentrated mass. This equivalent mass is connected to the midpoint of the two wheels. These uncertainties are further investigated in the robustness tests (Dai et al., 2015).

2.1 Determination of the overall center of mass

The equivalent pendulum-body mass of the upper structure is defined as

$$\begin{array}{}\text{(1)}& m_{\mathrm{p}}=m_{\mathrm{b}}+m_{\mathrm{h}}.\end{array}$$

The total mass of the system consists of the wheelchair body, excluding the two driving wheels, and the occupant mass. The overall center of mass of the wheelchair–occupant system is determined by mass-weighted averaging:

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{rl}{\displaystyle }{x}_{\mathrm{c}}& {\displaystyle }={\displaystyle \frac{m_{\mathrm{b}}{x}_{\mathrm{b}}+m_{\mathrm{h}}{x}_{\mathrm{h}}}{m_{\mathrm{b}}+m_{\mathrm{h}}}},\\ {\displaystyle }{z}_{\mathrm{c}}& {\displaystyle }={\displaystyle \frac{m_{\mathrm{b}}{z}_{\mathrm{b}}+m_{\mathrm{h}}{z}_{\mathrm{h}}}{m_{\mathrm{b}}+m_{\mathrm{h}}}},\end{array}\right.\end{array}$$

where (x_b,z_b) and (x_h,z_h) denote the center-of-mass coordinates of the wheelchair body and the occupant, respectively, and (x_c,z_c) denotes the overall center-of-mass coordinates of the wheelchair–occupant system.

2.2 Unconstrained dynamic modeling

The total kinetic energy of the system consists of the translational and rotational kinetic energy of the equivalent upper body T_p and the two driving wheels T_ω:

$$\begin{array}{}\text{(3)}& {\displaystyle }{T}_{\mathrm{p}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}_{\mathrm{p}}{v}_{\mathrm{p}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_{\mathrm{p}}{\dot{\mathit{\theta }}}^{\mathrm{2}},\text{(4)}& {\displaystyle }{T}_{\mathit{\omega }}& {\displaystyle }=\mathrm{2}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}_{\mathit{\omega }}{\dot{y}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_{\mathit{\omega }}{\dot{\mathit{\alpha }}}^{\mathrm{2}}\right).,\end{array}$$

where v_p is the translational velocity of the equivalent upper-pendulum body, and $\dot{\mathit{\alpha }}$ denotes the common angular velocity of the two driving wheels under straight-line motion, namely ${\dot{\mathit{\alpha }}}_r={\dot{\mathit{\alpha }}}_l=\dot{\mathit{\alpha }}$. The potential energy is determined mainly by the gravitational position of the equivalent upper body:

$$\begin{array}{}\text{(5)}& V=m_{\mathrm{p}}gL\cos \mathit{\theta },\end{array}$$

and the Lagrangian is

$$\begin{array}{}\text{(6)}& \mathcal{L}=T-V.\end{array}$$

Substituting the above expressions into the Lagrange equations yields

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\partial \mathcal{L}}{\partial \dot{q}}}\right)-{\displaystyle \frac{\partial \mathcal{L}}{\partial q}}=\mathcal{Q},\end{array}$$

where $\mathit{q}=[y,\mathit{\theta },{\mathit{\alpha }}_{\mathrm{r}},{\mathit{\alpha }}_{\mathrm{l}}{]}^{\mathrm{T}}$ is the generalized coordinate vector, and 𝒬 denotes the generalized force vector. For convenience in subsequent analysis, the unconstrained nonlinear dynamics can be written in the compact form

$$\begin{array}{}\text{(8)}& M\left(q\right)\ddot{q}+C_{\mathrm{d}}(q,\dot{q})\dot{q}+G\left(q\right)={\mathbf{B}}_{\mathrm{d}}\left(q\right)u+{\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}\left(q\right)\mathit{\lambda },\end{array}$$

where u denotes the single-wheel driving torque under symmetric actuation, with ${\mathit{\tau }}_{\mathrm{r}}={\mathit{\tau }}_{\mathrm{l}}=u$. Therefore, the total longitudinal driving force generated by the two wheels is $({\mathit{\tau }}_{\mathrm{r}}+{\mathit{\tau }}_{\mathrm{l}})/r=\mathrm{2}u/r$. A_c(q) is the nonholonomic constraint matrix, B_d(q) is the dynamic input matrix, and λ denotes the Lagrange multipliers associated with the wheel–ground motion constraints.

2.3 Nonholonomic constraint equations and model reduction

Due to the pure rolling assumption of the two driving wheels, the wheel–ground motion constraints can be written as

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{rl}{\displaystyle }y-{\displaystyle \frac{r}{\mathrm{2}}}({\mathit{\alpha }}_{\mathrm{r}}+{\mathit{\alpha }}_{\mathrm{l}})& {\displaystyle }=\mathrm{0},\\ {\displaystyle }{\mathit{\alpha }}_{\mathrm{r}}-{\mathit{\alpha }}_{\mathrm{l}}& {\displaystyle }=\mathrm{0},\end{array}\right.\end{array}$$

where y denotes the longitudinal displacement of the wheelchair body; r is the wheel radius; and α_r and α_l are the angular displacements of the right and left wheels, respectively.

Taking the time derivative of Eq. (9), the velocity-level nonholonomic constraints are obtained as

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{rl}{\displaystyle }\dot{y}-{\displaystyle \frac{r}{\mathrm{2}}}({\dot{\mathit{\alpha }}}_r+{\dot{\mathit{\alpha }}}_l)& {\displaystyle }=\mathrm{0},\\ {\displaystyle }{\dot{\mathit{\alpha }}}_r-{\dot{\mathit{\alpha }}}_l& {\displaystyle }=\mathrm{0}.\end{array}\right.\end{array}$$

Let the full generalized coordinate vector be defined as

$$\begin{array}{}\text{(11)}& \mathit{q}={\left[\begin{array}{cccc}y& \mathit{\theta }& {\mathit{\alpha }}_{\mathrm{r}}& {\mathit{\alpha }}_{\mathrm{l}}\end{array}\right]}^{\mathrm{T}}.\end{array}$$

Then Eq. (10) can be expressed in matrix form as

$$\begin{array}{}\text{(12)}& {\mathbf{A}}_{\mathrm{c}}\left(q\right)\dot{q}=\mathrm{0},\end{array}$$

where the constraint Jacobian is

$$\begin{array}{}\text{(13)}& {\mathbf{A}}_{\mathrm{c}}\left(q\right)=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& -\frac{r}{\mathrm{2}}& -\frac{r}{\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& -\mathrm{1}\end{array}\right].\end{array}$$

Before eliminating the constraint forces, the constrained Lagrange equations can generally be written as

$$\begin{array}{}\text{(14)}& {\mathbf{M}}_{\mathrm{f}}\left(q\right)\ddot{q}+C_{\mathrm{f}}(q,\dot{q})\dot{q}+G_{\mathrm{f}}\left(q\right)={\mathbf{B}}_{\mathrm{f}}\left(q\right)u+{\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}\left(q\right)\mathit{\lambda },\end{array}$$

where M_f(q) is the full mass matrix, $C_{\mathrm{f}}(q,\dot{q})\dot{q}$ denotes the Coriolis and centrifugal terms, G_f(q) is the gravitational term, and B_f(q) is the full input distribution matrix. The input u denotes the single-wheel driving torque under symmetric actuation, namely ${\mathit{\tau }}_{\mathrm{r}}={\mathit{\tau }}_{\mathrm{l}}=u$, and λ is the vector of Lagrange multipliers associated with the rolling constraints.

To avoid explicitly solving the constraint reaction forces in the control model, a null-space projection method is employed. The following derivation gives the projection procedure used to eliminate the Lagrange multipliers and obtain the reduced-order dynamics. Let S(q) be a basis matrix of the null space of A_c(q), satisfying

$$\begin{array}{}\text{(15)}& {\mathbf{A}}_{\mathrm{c}}\left(q\right)\mathbf{S}\left(q\right)=\mathrm{0}.\end{array}$$

For straight-line symmetric motion, the independent velocity vector is selected as

$$\begin{array}{}\text{(16)}& \mathit{v}={\left[\begin{array}{cc}\dot{y}& \dot{\mathit{\theta }}\end{array}\right]}^{\mathrm{T}}.\end{array}$$

The full generalized velocity can then be parameterized as

$$\begin{array}{}\text{(17)}& \dot{q}=\mathbf{S}\left(q\right)v,\end{array}$$

where

$$\begin{array}{}\text{(18)}& \mathbf{S}\left(q\right)=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\\ \frac{\mathrm{1}}{r}& \mathrm{0}\\ \frac{\mathrm{1}}{r}& \mathrm{0}\end{array}\right].\end{array}$$

It can be verified from Eqs. (13) and (18) that

$$\begin{array}{}\text{(19)}& {\mathbf{A}}_c\left(q\right)\mathbf{S}\left(q\right)=\mathrm{0},\end{array}$$

which indicates that all generalized velocities generated by S(q)v automatically satisfy the nonholonomic rolling constraints.

Taking the time derivative of Eq. (17) gives

$$\begin{array}{}\text{(20)}& \ddot{q}=\mathbf{S}\left(q\right)\dot{v}+\dot{S}\left(q\right)v.\end{array}$$

Substituting Eqs. (17) and (20) into Eq. (14) yields

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\mathbf{M}}_{\mathrm{f}}\left(q\right)\left[\mathbf{S}\left(q\right)\dot{v}+\dot{S}\left(q\right)v\right]& +C_{\mathrm{f}}(q,\dot{q})\mathbf{S}\left(q\right)v+G_{\mathrm{f}}\left(q\right)\\ & ={\mathbf{B}}_{\mathrm{f}}\left(q\right)u+{\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}\left(q\right)\mathit{\lambda }.\end{array}\end{array}$$

Premultiplying both sides of Eq. (21) by S^T(q) gives

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{S}^{\mathrm{T}}{\mathbf{M}}_{\mathrm{f}}S\dot{v}& +S^{\mathrm{T}}\left({\mathbf{M}}_{\mathrm{f}}\dot{S}+C_{\mathrm{f}}S\right)v+S^{\mathrm{T}}{G}_{\mathrm{f}}\\ & =S^{\mathrm{T}}{\mathbf{B}}_{\mathrm{f}}u+S^{\mathrm{T}}{\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}\mathit{\lambda }.\end{array}\end{array}$$

According to Eq. (15),

$$\begin{array}{}\text{(23)}& S^{\mathrm{T}}{\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}=({\mathbf{A}}_{\mathrm{c}}\mathbf{S}{)}^{\mathrm{T}}=\mathrm{0}.\end{array}$$

Therefore, the term containing the Lagrange multipliers is eliminated by the null-space projection, and the reduced-order dynamic model can be written as

$$\begin{array}{}\text{(24)}& M_{\mathrm{r}}\left(q\right)\dot{v}+C_{\mathrm{r}}(q,v)v+G_{\mathrm{r}}\left(q\right)=B_{\mathrm{r}}\left(q\right)u,\end{array}$$

where

$$\begin{array}{}\text{(25)}& \left\{\begin{array}{rl}{\displaystyle }{M}_{\mathrm{r}}\left(q\right)& {\displaystyle }=S^{\mathrm{T}}\left(q\right){\mathbf{M}}_{\mathrm{f}}\left(q\right)S\left(q\right),\\ {\displaystyle }{C}_{\mathrm{r}}(q,v)& {\displaystyle }=S^{\mathrm{T}}\left(q\right)\left[{\mathbf{M}}_{\mathrm{f}}\left(q\right)\dot{S}\left(q\right)+C_f(q,\dot{q})S\left(q\right)\right],\\ {\displaystyle }{G}_{\mathrm{r}}\left(q\right)& {\displaystyle }=S^{\mathrm{T}}\left(q\right)G_{\mathrm{f}}\left(q\right),\\ {\displaystyle }{B}_{\mathrm{r}}\left(q\right)& {\displaystyle }=S^{\mathrm{T}}\left(q\right){\mathbf{B}}_{\mathrm{f}}\left(q\right).\end{array}\right.\end{array}$$

Since S(q) in Eq. (18) is constant for the present longitudinal rolling model, $\dot{S}\left(q\right)=\mathrm{0}$.

For the reduced longitudinal balance model, Eq. (24) can be further expressed as

$$\begin{array}{}\text{(26)}& M\left(q\right)\dot{v}+C_{\mathrm{d}}(q,v)v+G\left(q\right)=B_{\mathrm{d}}\left(q\right)u,\end{array}$$

where

$$\begin{array}{}\text{(27)}& v=\left[\begin{array}{c}\dot{y}\\ \dot{\mathit{\theta }}\end{array}\right].\end{array}$$

The reduced matrices are given by

$$\begin{array}{}\text{(28)}& \left\{\begin{array}{rl}{\displaystyle }M\left(q\right)& {\displaystyle }=\left[\begin{array}{cc}{m}_{\mathrm{p}}+\mathrm{2}{m}_{\mathit{\omega }}+\frac{\mathrm{2}{J}_{\mathit{\omega }}}{r^{\mathrm{2}}}& m_{\mathrm{p}}L\cos \mathit{\theta }\\ m_{\mathrm{p}}L\cos \mathit{\theta }& J_{\mathrm{p}}+m_{\mathrm{p}}{L}^{\mathrm{2}}\end{array}\right],\\ {\displaystyle }{C}_{\mathrm{d}}(q,v)& {\displaystyle }=\left[\begin{array}{cc}\mathrm{0}& -m_{\mathrm{p}}L\sin \mathit{\theta }\dot{\mathit{\theta }}\\ \mathrm{0}& \mathrm{0}\end{array}\right],\\ {\displaystyle }G\left(q\right)& {\displaystyle }=\left[\begin{array}{c}\mathrm{0}\\ -m_{\mathrm{p}}gL\sin \mathit{\theta }\end{array}\right],\\ {\displaystyle }{B}_{\mathrm{d}}\left(q\right)& {\displaystyle }=\left[\begin{array}{c}\frac{\mathrm{2}}{r}\\ \end{array}\right].\end{array}\right.\end{array}$$

Here, m_p denotes the equivalent mass of the wheelchair–occupant body, m_ω denotes the mass of a single wheel, J_ω is the wheel moment of inertia, J_p is the pitch moment of inertia of the wheelchair–occupant body, L is the distance from the wheel axle to the equivalent center of mass, and g is the gravitational acceleration.

For numerical implementation in the NMPC prediction model and Simulink simulation, Eq. (26) is solved explicitly as

$$\begin{array}{}\text{(29)}& \dot{v}=M^{-\mathrm{1}}\left(q\right)\left[{\mathbf{B}}_{\mathrm{d}}\left(q\right)u-C_{\mathrm{d}}(q,v)v-G\left(q\right)\right].\end{array}$$

Let

$$\begin{array}{}\text{(30)}& \left\{\begin{array}{rl}{\displaystyle }{M}_{\mathrm{11}}& {\displaystyle }=m_{\mathrm{p}}+\mathrm{2}{m}_{\mathit{\omega }}+{\displaystyle \frac{\mathrm{2}{J}_{\mathit{\omega }}}{r^{\mathrm{2}}}},\\ {\displaystyle }{M}_{\mathrm{12}}& {\displaystyle }=m_{\mathrm{p}}L\cos \mathit{\theta },\\ {\displaystyle }{M}_{\mathrm{22}}& {\displaystyle }=J_{\mathrm{p}}+m_{\mathrm{p}}{L}^{\mathrm{2}},\\ {\displaystyle }{\mathrm{\Delta }}_M& {\displaystyle }=M_{\mathrm{11}}{M}_{\mathrm{22}}-M_{\mathrm{12}}^{\mathrm{2}}.\end{array}\right.\end{array}$$

Then the longitudinal acceleration and pitch angular acceleration are obtained as

$$\begin{array}{}\text{(31)}& \left\{\begin{array}{rl}{\displaystyle }\ddot{y}& {\displaystyle }={\displaystyle \frac{M_{\mathrm{22}}\left(\frac{\mathrm{2}u}{r}+m_{\mathrm{p}}L\sin \mathit{\theta }{\dot{\mathit{\theta }}}^{\mathrm{2}}\right)-M_{\mathrm{12}}\left(m_{\mathrm{p}}gL\sin \mathit{\theta }\right)}{{\mathrm{\Delta }}_M}},\\ {\displaystyle }\ddot{\mathit{\theta }}& {\displaystyle }={\displaystyle \frac{-M_{\mathrm{12}}\left(\frac{\mathrm{2}u}{r}+m_{\mathrm{p}}L\sin \mathit{\theta }{\dot{\mathit{\theta }}}^{\mathrm{2}}\right)+M_{\mathrm{11}}\left(m_{\mathrm{p}}gL\sin \mathit{\theta }\right)}{{\mathrm{\Delta }}_M}}.\end{array}\right.\end{array}$$

For the balance control problem, the state vector is selected as

$$\begin{array}{}\text{(32)}& x={\left[\begin{array}{cccc}y& \mathit{\theta }& \dot{y}& \dot{\mathit{\theta }}\end{array}\right]}^{\mathrm{T}}.\end{array}$$

The nonlinear state-space model used in the controller design is then written as

$$\begin{array}{}\text{(33)}& \dot{x}=f(x,u)=\left[\begin{array}{c}\dot{y}\\ \dot{\mathit{\theta }}\\ \ddot{y}\\ \ddot{\mathit{\theta }}\end{array}\right].\end{array}$$

3 Hybrid NMPC–LQR control method

3.1 Overall construction of the control method

The hybrid controller is designed to address the online optimization problem in self-balancing control. The controller provides the driving torque command to the PD-actuated system. It consists of two components, namely an NMPC controller and an LQR controller. When the tilt angle exceeds the predefined switching threshold θ_sw, the NMPC controller is activated to cope with the nonlinear operating condition. In the controller design, u denotes the single-wheel torque command, and the same command is applied to the left and right wheels under straight-line symmetric motion. Once the tilt angle falls below the switching threshold, the control law is switched to the LQR controller:

$$\begin{array}{}\text{(34)}& u=\left\{\begin{array}{ll}{u}_{\mathrm{NMPC}},& \left|\mathit{\theta }\right|\ge {\mathit{\theta }}_{\mathrm{sw}},\\ u_{\mathrm{LQR}},& \left|\mathit{\theta }\right|<{\mathit{\theta }}_{\mathrm{sw}}.\end{array}\right.\end{array}$$

The overall hybrid control architecture is illustrated in Fig. 2.

Figure 2 Overall hybrid NMPC–LQR control framework for the self-balancing wheelchair.

3.2 LQR controller design

The LQR controller is used for local stabilization when the system approaches the equilibrium point. Based on the linearized state-space model and with reference to Chacko and Abraham (2023), the cost function of the LQR controller is constructed as

$$\begin{array}{}\text{(35)}& J_{\mathrm{LQR}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left[(x-x_{\mathrm{ref}}{)}^{\mathrm{T}}{\mathbf{Q}}_L(x-x_{\mathrm{ref}})+u^{\mathrm{T}}{\mathbf{R}}_{L}u\right]\mathrm{d}t,\end{array}$$

where Q_L and R_L are the weighting matrices for the state error and the control input, respectively. They are selected as

$$\begin{array}{}\text{(36)}& \left\{\begin{array}{rl}{\displaystyle }{\mathbf{Q}}_L& {\displaystyle }=\mathrm{diag}(\mathrm{20},\mathrm{200},\mathrm{20},\mathrm{50}),\\ {\displaystyle }{\mathbf{R}}_L& {\displaystyle }=\mathrm{1}.\end{array}\right.\end{array}$$

By solving the algebraic Riccati equation, the optimal feedback gain matrix is obtained, and the LQR control law becomes

$$\begin{array}{}\text{(37)}& u_{\mathrm{LQR}}=-\mathbf{K}(x-{\mathit{x}}_{\mathrm{ref}}),\end{array}$$

where K is the LQR feedback gain matrix and x_ref denotes the reference state vector.

3.3 NMPC controller design

The NMPC controller is employed to address the highly nonlinear behavior of the system. When the tilt angle of the wheelchair exceeds the switching threshold, the NMPC controller is activated. At each control update, an optimal control sequence over a finite prediction horizon is solved online on the basis of the current system state and the nonlinear dynamic model. Only the first control input is then applied to the system, thereby forming a receding-horizon optimization framework.

With the state at the kth sampling instant, prediction horizon N_p, and sampling time defined, the finite-horizon objective function is written as

$$\begin{array}{}\text{(38)}& \begin{array}{rl}{J}_{\mathrm{NMPC}}& =\sum _{i=\mathrm{0}}^{N_{\mathrm{p}}-\mathrm{1}}\left[(x_{k+i}-x_{\mathrm{ref}}{)}^{\mathrm{T}}{Q}_n(x_{k+i}-x_{\mathrm{ref}})\right.\\ & \left.+u_{k+i}^{\mathrm{T}}{R}_n{u}_{k+i}\right],\end{array}\end{array}$$

where k is the current sampling instant, i is the prediction-step index, N_p is the prediction horizon, and u_k+i is the predicted equivalent torque input at the (k+i)th step. The weighting matrices are selected by prioritizing attitude safety. Larger weights are assigned to the pitch angle and pitch angular velocity, while the input term penalizes excessive wheel torque:

$$\begin{array}{}\text{(39)}& \left\{\begin{array}{rl}{\displaystyle }{Q}_n& {\displaystyle }=\mathrm{diag}(\mathrm{500},\mathrm{1000},\mathrm{50},\mathrm{300}),\\ {\displaystyle }{R}_n& {\displaystyle }=\mathrm{0.1}.\end{array}\right.\end{array}$$

The nonlinear state equations are directly used as prediction-model constraints, and actuator limits are imposed as input constraints:

$$\begin{array}{}\text{(40)}& \left\{\begin{array}{rl}{\displaystyle }{x}_{k+i+\mathrm{1}}& {\displaystyle }=f(x_{k+i},u_{k+i}),\\ {\displaystyle }{u}_{min}& {\displaystyle }\le u_{k+i}\le u_{max}.\end{array}\right.\end{array}$$

Here, u_min and u_max are the lower and upper bounds of the single-wheel torque command.

3.4 Switching-threshold evaluation and parameter sensitivity

The switching threshold determines the boundary between the nonlinear recovery stage and the local stabilization stage. According to the switching law in Eq. (34), NMPC is activated when $\left|\mathit{\theta }\right|\ge {\mathit{\theta }}_{\mathrm{sw}}$, whereas LQR is used when $\left|\mathit{\theta }\right|<{\mathit{\theta }}_{\mathrm{sw}}$. Therefore, if the threshold is too small, the NMPC controller remains active until the system is very close to the equilibrium, which increases the computational burden and may delay the use of efficient local regulation. Conversely, if the threshold is too large, the controller may switch to LQR while the system is still in a region with strong nonlinear effects, which may weaken the recovery capability and reduce the validity of the local linear controller. Hence, the switching threshold must be selected by considering both nonlinear recovery and local stabilization performance.

3.4.1 Determination of the threshold search range

Before evaluating the switching threshold, a linearization-error analysis was conducted to determine a reasonable search interval. The nonlinear model and the linearized model were compared under different local angular ranges. The tested ranges were selected as 0.015, 0.030, 0.050, 0.080, 0.100, and 0.150 rad, which correspond to the candidate switching thresholds used in the following evaluation. An additional range of 0.200 rad was included as an outside reference.

Figure 3 Linearization mismatch under different angular ranges. The shaded region denotes the selected switching-threshold search interval of 0.015–0.150 rad, while 0.200 rad is used as an outside reference.

The results in Fig. 3 show that the pitch angular–acceleration mismatch increases as the angular range expands. Within the interval of 0.015–0.150 rad, the mismatch increases gradually, whereas a more pronounced increase is observed when the range is enlarged to 0.200 rad. Therefore, the interval of 0.015–0.150 rad is adopted as the switching-threshold search range. It should be noted that the linearization-error analysis is used to determine the candidate range, while the final threshold is selected according to the five-index performance evaluation.

3.4.2 Five-index evaluation under nominal parameters

Based on the above search interval, six switching thresholds are selected for simulation tests:

$${\mathit{\theta }}_{\mathrm{sw}}=\mathrm{0.015},\mathrm{0.030},\mathrm{0.050},\mathrm{0.080},\mathrm{0.100},\mathrm{0.150}\mathrm{rad}.$$

Under identical initial conditions, controller parameters, and input constraints, the five performance indices corresponding to each threshold are extracted, as shown in Table 2.

Table 2 Performance indices under different switching thresholds under nominal model parameters.

The min–max normalization method is employed to remove the influence of different physical dimensions and numerical ranges among the evaluation indices:

$$\begin{array}{}\text{(41)}& {\mathit{\eta }}_i^{\ast }={\displaystyle \frac{{\mathit{\eta }}_i-min\left(\mathit{\eta }\right)}{max\left(\mathit{\eta }\right)-min\left(\mathit{\eta }\right)}},\end{array}$$

where η_i denotes the original value of a performance index and ${\mathit{\eta }}_i^{\ast }$ denotes its normalized value.

A five-index composite cost function is then established to comprehensively evaluate the threshold-dependent performance:

$$\begin{array}{}\text{(42)}& J_{\mathrm{5}}={\mathit{\nu }}_{\mathrm{1}}{t}_{\mathrm{s}}^{\ast }+{\mathit{\nu }}_{\mathrm{2}}{e}_{\mathrm{ss}}^{\ast }+{\mathit{\nu }}_{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{s}}^{\ast }+{\mathit{\nu }}_{\mathrm{4}}{d}_{\mathrm{r}}^{\ast }+{\mathit{\nu }}_{\mathrm{5}}{e}_{max}^{\ast },\end{array}$$

where $t_{\mathrm{s}}^{\ast }$, $e_{\mathrm{ss}}^{\ast }$, ${\mathit{\sigma }}_{\mathrm{s}}^{\ast }$, $d_{\mathrm{r}}^{\ast }$, and $e_{max}^{\ast }$ denote the normalized settling time, steady-state error, steady-state standard deviation, drift rate, and maximum steady-state deviation, respectively. The weights are selected as

$$\begin{array}{}\text{(43)}& [{\mathit{\nu }}_{\mathrm{1}},{\mathit{\nu }}_{\mathrm{2}},{\mathit{\nu }}_{\mathrm{3}},{\mathit{\nu }}_{\mathrm{4}},{\mathit{\nu }}_{\mathrm{5}}]=[\mathrm{0.18},\mathrm{0.30},\mathrm{0.20},\mathrm{0.15},\mathrm{0.17}].\end{array}$$

A larger weight is assigned to the steady-state error because accurate balance near the equilibrium is a primary requirement for wheelchair stabilization. A lower value of J₅ indicates better overall control performance.

The calculated five-index costs are summarized in Table 3. The results show that smaller thresholds lead to shorter settling time in some cases, but they may also produce larger steady-state error and fluctuation. As the threshold increases, the steady-state behavior improves within a certain range, while the settling time gradually increases. Among the tested candidates, θ_sw=0.100 rad yields the minimum five-index cost, indicating the best overall trade-off between transient recovery and steady-state regulation under the nominal model parameters.

Table 3 Five-index costs under different switching thresholds under nominal model parameters.

3.4.3 Parameter sensitivity of the selected threshold

To examine whether the selected switching threshold strongly depends on nominal model parameters, a one-factor-at-a-time sensitivity analysis was further conducted. The equivalent upper-body mass m_p and the equivalent pendulum length L were varied independently. For each parameter factor, the threshold was scanned within the local range of 0.08–0.12 rad, and the optimal threshold was selected according to the minimum five-index cost.

For the equivalent mass variation, the mass factor is defined as

$$\begin{array}{}\text{(44)}& {\mathit{\rho }}_{\mathrm{m}}={\displaystyle \frac{m_{\mathrm{p}}}{m_{p\mathrm{0}}}}.\end{array}$$

The optimal threshold decreases as the equivalent mass increases. A linear fitting gives

$$\begin{array}{}\text{(45)}& {\mathit{\theta }}_{\mathrm{sw},\mathrm{opt}}=\mathrm{0.212}-\mathrm{0.110}{\mathit{\rho }}_{\mathrm{m}},\end{array}$$

with R²=0.945. This negative correlation indicates that a larger equivalent upper-body mass tends to require a smaller switching threshold, allowing NMPC to remain active over a wider angular region before switching to LQR. It should also be noted that the optimal thresholds for the 0.8 and 0.9 mass factors are located at the upper boundary of the tested interval, suggesting that the optimal threshold may shift slightly upward when the equivalent mass decreases.

For the equivalent pendulum-length variation, the optimal threshold does not show a strictly monotonic trend. Under several length-variation cases, the minimum J₅ occurs at the lower boundary of the tested interval, namely θ_sw=0.08 rad. This indicates that when the center-of-mass-related parameter deviates from the nominal value, a smaller switching threshold is generally beneficial because NMPC remains active over a wider angular region and provides stronger nonlinear recovery capability.

Table 4 Optimal switching thresholds under representative parameter variations.

Overall, the nominal threshold evaluation identifies θ_sw=0.100 rad as the best fixed switching threshold under the nominal model parameters. The mass-sensitivity analysis further provides an empirical estimate of the optimal threshold within the tested mass range. For ${\mathit{\rho }}_{\mathrm{m}}=m_{\mathrm{p}}/m_{\mathrm{p}\mathrm{0}}\in [\mathrm{0.8},\mathrm{1.2}]$, the fitted relation ${\mathit{\theta }}_{\mathrm{sw},\mathrm{opt}}=\mathrm{0.212}-\mathrm{0.110}{\mathit{\rho }}_{\mathrm{m}}$ indicates that a larger equivalent mass generally requires a smaller switching threshold. Therefore, θ_sw=0.100 rad is adopted as the nominal fixed threshold, while the fitted mass–threshold relation provides a reference for parameter-dependent threshold adjustment within the tested range.

3.5 Stability considerations

The stability consideration of the proposed controller is based on the functional division between NMPC and LQR. LQR is not used over the entire nonlinear operating range but is restricted to the near-equilibrium region verified by the linearization-error analysis in Sect. 3.4. In this region, the nonlinear model is linearized around the equilibrium point, and the LQR gain is designed so that the closed-loop matrix A_l−B_lK is Hurwitz. Thus, the LQR stage provides local asymptotic stability for the linearized system.

For states outside the local region, NMPC is activated to handle large-deviation nonlinear recovery under input constraints. The switching threshold therefore serves as a boundary between nonlinear recovery and local stabilization, preventing premature use of LQR in regions where the linear approximation is less reliable.

The above discussion provides a local stability guarantee for the LQR stage and a practical stability argument for the hybrid framework. Rigorous global Lyapunov proof for the completely switched nonlinear system is beyond the scope of this study. Therefore, the overall closed-loop behavior is validated through large-initial-angle simulations, disturbance tests, parameter-sensitivity tests, friction-coefficient tests, and prototype experiments.

4 Experiments

4.1 Prototype

The desktop-scale two-wheeled self-balancing prototype used for experimental validation was developed based on the WHEELTEC B570 balancing vehicle architecture and employed an STM32F103C8T6 microcontroller together with an MPU6050 sensor for attitude measurement and closed-loop control implementation.

4.2 Simulation

Simulink in MATLAB is used to construct the simulation platform. The model includes the wheelchair body, driving wheels, center-of-mass and inertia parameters, external-disturbance input, and state-output modules. The body is modeled as a rigid body pitching about the wheel axis, and the driving torque is applied through the wheel-actuator module. The outputs include the pitch angle, pitch angular velocity, wheel displacement, and wheel velocity, which are used by the LQR, NMPC, and hybrid controllers. Disturbance input ports are reserved for robustness tests. The corresponding signal–flow relationship is shown in Fig. 4.

Figure 4 Signal–flow diagram of the hybrid NMPC–LQR balance control system.

4.2.1 Numerical discretization of the NMPC model

The continuous-time nonlinear dynamic model is numerically discretized before being incorporated into the finite-horizon optimization problem. The state evolution is approximated with high accuracy by jointly utilizing the state values and dynamic information at both ends of the interval as well as the midpoint state in each sampling interval:

$$\begin{array}{}\text{(46)}& \left\{\begin{array}{rl}{\displaystyle }{x}_{\mathrm{m},k}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(x_k+x_{k+\mathrm{1}}\right)+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{8}}}\left[f\left(x_{k+\mathrm{1}},u_{k+\mathrm{1}}\right)-f\left(x_k,u_k\right)\right],\\ {\displaystyle }{x}_{k+\mathrm{1}}& {\displaystyle }=x_k+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{6}}}\left[f\left(x_k,u_k\right)+\mathrm{4}f\left(x_{\mathrm{m},k},u_{\mathrm{m},k}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+f\left(x_{k+\mathrm{1}},u_{k+\mathrm{1}}\right)\right],\end{array}\right.\end{array}$$

where Δt is the sampling interval; x_m,k is the midpoint state; and u_m,k is the midpoint control input, which can be approximated as $u_{\mathrm{m},k}=(u_k+u_{k+\mathrm{1}})/\mathrm{2}$. The resulting discrete state–transition relationship enables the construction of the discrete prediction model required by NMPC. Based on the above discrete state–transition relationship, the NMPC problem is formulated as a nonlinear programming problem over a finite prediction horizon. The optimization variables include the state sequence and control input sequence over the prediction horizon, while the constraints include consistency between the current state and the measured state, discrete state–transition constraints, and upper and lower bounds on the control input. To guarantee real-time performance and closed-loop robustness, a receding-horizon optimization mechanism is adopted. At each sampling instant, the optimization problem is reformulated using the current measured state as the initial condition. The optimal future control sequence is computed, and only the first control input is implemented.

4.3 Balance verification

4.3.1 Comparative simulation under large initial tilt angles

The recovery motion of the self-balancing wheelchair is simulated under initial tilt angles of 35, 40, and 45°. The standalone LQR controller, the standalone NMPC controller, and the proposed hybrid NMPC–LQR controller are compared under the same system parameters, input constraints, and simulation settings. The representative response curves are shown in Fig. 5, and the quantitative performance indices are summarized in Table 5.

Figure 5 Comparative simulation responses of the LQR, NMPC, and hybrid controllers under initial tilt angles of 35, 40, and 45°.

Table 5 Performance comparison of LQR, NMPC, and hybrid controllers under large initial tilt angles.

4.3.2 Benchmark comparison with sliding-mode control

The SMC controller was tested using the same simulation model and input constraint as the proposed controller. Since sliding-mode control is widely applied to nonlinear systems, it provides a meaningful benchmark for evaluating the proposed hybrid strategy. Recent studies have also used SMC as a robust benchmark to evaluate intelligent LQR-based controllers for inverted-pendulum-type systems (Saleem et al., 2025).

Table 6 Performance of the SMC benchmark under different recoverable initial tilt angles.

Table 7 Recoverability comparison between the proposed hybrid NMPC–LQR controller and the SMC benchmark under different initial tilt angles.

4.3.3 Comparative analysis of control performance

Table 5 shows that the proposed hybrid controller achieves the shortest settling time under all three large-initial-angle conditions. At 35 and 40°, LQR gives slightly smaller steady-state error and fluctuation, reflecting its local regulation advantage. However, the standalone LQR controller has weaker large-deviation recovery capability, whereas the standalone NMPC controller shows longer settling time and larger residual oscillation.

At 45°, the hybrid controller achieves the fastest recovery and the smallest steady-state error, standard deviation, and maximum steady-state deviation. This indicates that NMPC improves large-angle nonlinear recovery, while LQR enhances local stabilization after the state enters the near-equilibrium region.

The SMC benchmark results in Tables 6 and 7 show that SMC fails under 40 and 45°, whereas the proposed hybrid controller remains stable under all tested initial angles. Therefore, the hybrid NMPC–LQR controller provides a better balance between recoverability, local stabilization, and robustness under the same input constraint.

4.3.4 Computational efficiency test

A computational time test was conducted under the same simulation conditions to evaluate the computational efficiency of the proposed hybrid controller. Five repeated runs were performed under the same simulation settings. The average computation time decreases from 48.04 min for standalone NMPC to 35.00 min for the hybrid controller, indicating that the proposed switching strategy reduces the online computational burden. The detailed results are provided in Appendix B.

4.3.5 External disturbance test

The disturbance-rejection capability of the hybrid controller is evaluated using both pulse and step disturbances. The disturbance magnitudes are 50, 65, 70, and 80 N, and the response curves are provided in Appendix A.

The system exhibits a secondary pitch deflection after the pulse is applied, but no sustained divergence is observed. Under 50 and 65 N disturbances, the wheelchair quickly returns to the stable region. Under 70 and 80 N disturbances, the recovery process becomes longer and the residual oscillation increases, but the response remains bounded. These results indicate that the proposed controller retains disturbance-rejection capability under representative pulse impacts.

In addition to the pulse-disturbance test, a step-disturbance test was conducted to evaluate the response of the proposed controller under persistent external force. The step disturbance was applied at t=2 s after the initial recovery process, and the tested disturbance magnitudes were 10, 20, and 30 N. The detailed response curves are provided in Appendix D. The results show that the post-step peak deviation increases with the disturbance magnitude. Specifically, the maximum post-step pitch deviations are 3.62, 5.89, and 8.11° under 10, 20, and 30 N step disturbances, respectively. Nevertheless, the responses remain bounded and gradually return to the near-equilibrium region. This indicates that the proposed hybrid controller can maintain recovery capability under persistent external disturbances.

4.3.6 Flat-ground friction-coefficient robustness test

Flat-ground simulations were further conducted to examine the influence of wheel–ground friction uncertainty. The slope angle was set to zero, and all controller parameters were kept unchanged. Only the static friction coefficient in the spatial contact model was varied. The tested coefficients were 0.03, 0.10, and 0.30, corresponding to low-, medium-, and high-friction contact conditions, respectively.

Table 8 Flat-ground robustness test under different static friction coefficients.

The results in Table 8 show that the wheelchair does not overturn under the tested friction coefficients, indicating bounded recovery under friction-parameter variations. The steady-state error changes slightly from 0.03 to 0.10, while the high-friction case of 0.30 produces larger steady-state fluctuation. This suggests that stronger wheel–ground contact interaction may amplify residual pitch oscillations.

Since the responses do not enter and remain within the prescribed error band during the 10 s simulation, these results are interpreted as bounded recovery rather than exact convergence to the upright target. The detailed response curves are provided in Appendix C, where the final steady-state mean angle is also marked for each friction condition.

4.3.7 Prototype balance-recovery experiments

Prototype balance-recovery experiments were conducted under initial tilt angles of 25, 30, and 35°. The measured responses are shown in Fig. 6. In all cases, the prototype recovers to a small neighborhood around the upright equilibrium, indicating the physical realizability of the proposed controller.

Figure 6 Prototype balance-recovery experiment results under initial tilt angles of 25, 30, and 35°.

When the initial angle increases, the transient rebound and oscillation become more evident. Nevertheless, the response remains bounded and eventually approaches the upright region. These experimental results are consistent with the simulation trends and further support the engineering applicability of the proposed hybrid control strategy.

5 Conclusions

This paper develops a hybrid NMPC–LQR balance control strategy for two-wheeled self-balancing wheelchairs. NMPC is used for large-deviation nonlinear recovery, while LQR is restricted to the local stabilization region verified by linearization-error analysis. Comparative simulations, SMC benchmark tests, disturbance tests, friction-coefficient tests, and prototype experiments show that the proposed controller improves large-angle recovery capability while reducing the computational burden compared with the standalone NMPC.

For switching-threshold selection, the linearization-error analysis is first used to determine the threshold search range. A normalized five-index evaluation criterion is then adopted to compare different threshold candidates. Under the nominal model parameters, θ_sw=0.100 rad gives the lowest comprehensive cost and is selected as the nominal fixed threshold. The mass-sensitivity analysis further provides an empirical estimate of the optimal threshold within the tested mass range. For ${\mathit{\rho }}_{\mathrm{m}}=m_{\mathrm{p}}/m_{p\mathrm{0}}\in [\mathrm{0.8},\mathrm{1.2}]$, the fitted relation ${\mathit{\theta }}_{\mathrm{sw},\mathrm{opt}}=\mathrm{0.212}-\mathrm{0.110}{\mathit{\rho }}_{\mathrm{m}}$ indicates that a larger equivalent upper-body mass generally requires a smaller switching threshold.

Several limitations remain. The model is a planar reduced-order model and does not include yaw, roll, tire deformation, lateral instability, or active human-body motion. The stability discussion is based on local LQR analysis and closed-loop simulations rather than rigorous global Lyapunov proof for the switched nonlinear system. In addition, the prototype experiments are conducted on a desktop-scale platform. Future work will focus on human–wheelchair coupled modeling, multi-parameter threshold fitting, adaptive switching-threshold calculation, explicit friction compensation, robust NMPC with terminal constraints, and full-scale experiments on slopes and uneven terrains.

Appendix A: Pitch-angle responses under external pulse disturbances

Figure A1 Pitch-angle responses under different external pulse disturbances.

Appendix B: Detailed results of the computational time tests

Table B1 Detailed results of the computational time tests.

Appendix C: Pitch-angle responses under different flat-ground friction coefficients

Figure C1 Pitch-angle responses under different flat-ground static friction coefficients. The horizontal dash-dotted line denotes the final steady-state mean angle θ_ss, calculated from the last 10 % of the simulation time.

Appendix D: Pitch-angle responses under step disturbances

Figure D1 Pitch-angle responses under different step disturbances. The shaded region denotes the period after the step disturbance is applied, and the red marker indicates the maximum post-step pitch deviation.