1. Introduction
Since the appearance of spherical robots (SR) in 1996, they have received extensive attention from scholars at home and abroad due to their compact structure, flexible movement, and ability to withstand various task loads [Reference Halme, Suomela, Schonberg and Wang1]. The SR is a fully enclosed robot in which the drive mechanism, control system, and battery are all enclosed in the spherical shell. This prevents adverse damage to the external environment, which makes the SR suitable to be applied in humid, dusty, radiation, or toxic environment.
At present, the drive mechanism of the SR is mainly divided into omnidirectional wheel eccentric torque drive [Reference Armour and Vincent2], frictional internal drive [Reference Chen, Chen, Tsai, Yang and Lin3], and eccentric torque drive of heavy pendulum [Reference Yang, Tao, Huo, Neusypin, Zhang, Zhang and Guo4]. The pendulum drive is favored by industry and academia due to its simple implementation and significant driving torque. However, during the pendulum SR movement, the internal suspension is in the under actuation state, which affects its stability and accuracy, and even causes control failure.
In recent years, most studies on the pendulum SR’s dynamic modeling and motion control have mainly focused on the speed control of the SR. For instance, a motion control approach was proposed in ref. [Reference Joshi and Banavar5] for a spherical mobile robot rolling on a plane, in which only one Euler angle was varied in each step, improving computational efficiency. In refs. [Reference Sun, Xiao, Jia and Wang6–Reference Kayacan, Bayraktaroglu and Saeys8], the motion equation was decomposed into horizontal and vertical subsystems. References [Reference Sun, Xiao, Jia and Wang6, Reference Liu, Sun, Jia and Wang7] designed full-state feedback and PI controllers for rolling and forward motions, respectively. In ref. [Reference Kayacan, Bayraktaroglu and Saeys8], a feedback linearized PD-type fuzzy controller was employed for the closed-loop control of robot motion. Roozegar M. et al. [Reference Roozegar and Mahjoob9] utilized the heuristic fuzzy controller and PID controllers to realize the SR’s motion control on an inclined plane. The above approaches did not consider the system’s robustness to parametric uncertainties and the reduction of stationary errors in detail. Wang Y. et al. [Reference Wang, Guan, Hu, Zhang, Wang, Wang, Liu and Li10] constructed a yaw angle prediction-based fuzzy PID controller to adjust the attitude of a SR, which can perform fuzzy reasoning and autonomously modify control parameters based on the robot’s current state information. However, the controller is based on a kinematics design, which cannot effectively control the SR under the uncertain parameters of its model. Xie L. et al. [Reference Xie, Yu and Chen11] studied the dynamics modeling and joint trajectory tracking control of robotic systems with external disturbances and uncertain parameters and proposed a robust fuzzy sliding mode control to achieve desired trajectory tracking of joints. The synovial controller is independent of object model parameters and disturbances, the sliding modes of which can be designed.
In order to solve the problems of inaccurate dynamic modeling, parameter changes, and interference in the SR control that led to poor performance of traditional controllers, Cai Y. et al. [Reference Cai, Zhan and Xi12] proposed a speed control method combining fuzzy logic and SMC, while Kayacan E. et al. [Reference Kayacan, Ramon and Saeys13] utilized the adaptive neuro-fuzzy controller and the learning algorithm of the SMC theory for the SR’s speed control. Liu Y. et al. [Reference Liu, Wang, Guan, Wang, Jin, Hu, Ren, Hao, Zhang and Li14, Reference Liu, Wang, Guan, Hu, Zhang, Jin, Wang, Hao and Li15] designed a hierarchical SMC (HSMC) for speed and orientation control of the SR, improving control efficiency and stability. Nevertheless, it is challenging to fit sensors in the limited space of the sphere to measure variables such as rolling speed in spherical robots. Therefore, applying the HSMC method for the motion control of SRs is challenging due to the unmeasurable state and the uncertain rolling resistance torque generated by the ground. Yue M. et al. [Reference Yue, Liu, An and Sun16] combined HSMC and the state expansion observer to promote the designed closed-loop system’s robustness and adaptability. As we all know, since the state observer can reconstruct the system with the measurable state and obtain the value of the internal state variable, the uncertain rolling resistance can be obtained in real time.
Most of the above studies are in the simulation step, and even if there are physical experiments, they are mostly confined to smooth ground applications. SR control systems are susceptible to unmodeled dynamics, parameter changes, uncertainties, and disturbances. For example, the parameters are susceptible to disturbance when the spherical rolling robot rolls on uneven ground. Active disturbance rejection controller (ADRC) [Reference Han17] is independent of the system model and can estimate and compensate for different disturbances to the system in real time during motion. Tao T. et al. [Reference Tao, Yang, Huo, Zhang, Neusypin, Guo, Yan, Duan and Yu18] and Lin R. et al. [Reference Lin, Liu, Huo, Zhang, Yang and Guo19] proposed an ADRC-based motion controller for SRs moving on uneven ground, which enhanced their anti-interference capability. Liu M. et al. [Reference Liu, Lin, Yang, Nazarova and Huo20] designed a motion controller (PSO-ADRC) that employs particle swarm optimization (PSO) to adjust the parameters
$b_{0}$
,
$\,\beta_{01}$
,
$\beta_{02}$
, and
$\beta_{03}$
in ADRC to enhance the SRs’ anti-interference efficiency in actual motion. However, the ADRC controller has many parameters where their tuning can considerably influence the SR’s motion control efficiency.
There are other concerns about improving the control capability of disturbed parameters during the motion of SRs. For example, Yue M. [Reference Yue and Deng21], Sun H. [Reference Sun, Zheng and Jia22], and Loh R.N.K. [Reference Loh and Karsiti23] designed the optimal control law of a linear quadratic regulator (LQR) with the change of state variables as the objective function to reduce the disturbance caused by the dynamic change of the state effects on SRs. Zadeh F. K. et al. [Reference Zadeh, Moallem, Asiri and Zadeh24] designed an efficient LQR to eliminate the oscillation in a SR’s motion and improve the system’s stability. Zhan Q. et al. [Reference Zhan, Liu and Cai25–Reference Zheng, Zhan, Liu and Cai29] utilized backstepping, the neurodynamic shunt model, the cerebellar model joint control model, and RBF neural network to improve the controller’s disturbance immunity.
Although the above control methods ensure the SR’s stable control to a certain extent, the motion parameters of SRs vary significantly due to the lack of their environmental perception ability and susceptibility to external and internal uncertainties. Therefore, designing the control system based on an accurate model has certain limitations, especially when the system is greatly disturbed, degrading the system’s control effect significantly. Model reference adaptive control (MRAC) is a powerful tool for solving the mentioned problem [Reference Kireçci, Topalbekiroglu and Eker30]. Typically, an MRAC comprises two distinct loops, where the inner one, including the plant and the controller, is connected in parallel with the reference model, and the outer one is an adaptation mechanism that adjusts the controller parameters. Since the reference signal applied to the adjustable controller also acts on the reference model, the system can effectively solve the reduction of the control effect due to the parameter change of the system model so that the system can obtain better control performance. For example, Roozegar M. et al. [Reference Roozegar, Mahjoob and Ayati31] developed an MRAC-based adaptive method for the control of pendulum-driven SRs. Bomfim M. et al. [Reference Bomfim, Monteiro and Lima II32] designed a Lyapunov-based hybrid MRAC for a pneumatic artificial muscle-actuated manipulator.
Inspired by the above research, a Lyapunov-based MRAC control algorithm is proposed. The essential contributions of the current work are (1) According to the actual situation of SR motion, taking the pendulum’s swing angle and the spherical shell’s rolling angle as generalized coordinates, the Euler–Lagrange approach is adopted to establish the SR learning model’s omnidirectional dynamics. (2) The dynamic model is simplified sequentially to design the controller. According to the steering principle of the car during driving, the XK-I SR’s motion control model is decoupled into linear and steering, and the state space expressions of the SR’s linear and steering motions are obtained. (3) The proposed MRAC control algorithm is based on the Lyapunov method. MRAC is utilized to design the SR’s adaptive motion control law, and the controller parameters are adjusted through the Lyapunov approach to overcome the poor control effect induced by the SR’s inaccurate dynamic model and parameter uncertainty. Finally, the paper conducts experimental verification and compares MRAC with LQR, Fuzzy-PID, PSO-ADRC, and HSMC.
The remainder of this article is arranged as follows. Section 2 describes the XK-I SR’s mathematical model, including the uncoupled dynamic model and state space equations. Section 3 designs the adaptive motion controller for the SR. The simulation and experimental results are given in Section 4. In the end, the conclusions of this study and future research perspectives are presented in Section 5.
2. Description of the pendulum SR system
2.1. Dynamic model of the pendulum SR
The XK-I SR in this study is composed of a two-degree-of-freedom pendulum driven by two sets of motors. As described in Fig. 1, a Cartesian coordinate system
$o - x_{1}y_{1}z_{1}$
is established, considering the sphere’s center as the origin
$o$
. The driving torque produced by the motor makes the pendulum swing around the
$ox_{1}$
and
$oy_{1}$
axes by angles
$\alpha$
and
$\beta$
, respectively. In Fig. 1,
$o - xyz$
describes an inertial coordinate system fixed on the moving ground, and
$o - x_{0}y_{0}z_{0}$
describes a moving coordinate system with the center of the sphere as the origin. The
$o - x_{0}y_{0}z_{0}$
coordinate system is only translated relative to the
$o - xyz$
coordinate system, and the
$o - x_{1}y_{1}z_{1}$
coordinate system is only rotated with respect to the
$o - x_{0}y_{0}z_{0}$
coordinate system.
$\theta$
,
$\varphi$
, and
$\psi$
are the rolling angles of the sphere around
$ox$
,
$oy$
, and
$oz$
, respectively. From the actual motion of the XK-I type, only considering the sphere rotation around the lateral axis ox and the longitudinal axis
$oy$
, the rotation matrix between
$o - x_{1}y_{1}z_{1}$
and
$o - x_{0}y_{0}z_{0}$
is
\begin{equation} \boldsymbol{R} = \boldsymbol{R}(x,\theta )\boldsymbol{R}(y,\varphi ) = \begin{bmatrix}{\cos\varphi } &\quad 0 &\quad {-\!\sin\varphi } \\ \\[-7pt]{\sin\theta \sin\varphi } &\quad {\cos\theta } &\quad {\sin\theta \cos\varphi } \\ \\[-7pt]{\cos\theta \sin\varphi } &\quad {-\!\sin\theta } &\quad {\cos\theta \cos\varphi } \\ \end{bmatrix} \end{equation}
XK-I spherical robot coordinate system.
In order to construct the XK-I SR’s kinematics and dynamics models, the following assumptions are made:
-
1. There is no sliding friction when the ball moves on the ground;
-
2. The spherical shell’s density and thickness are uniformly distributed, and the mass is
$m_{1}$
; -
3. The square frame of the XK-I type is like a pendulum, and the mass is
$m_{2}$
; -
4. The pendulum’s center is connected to the sphere’s center with a rigid rod regardless of mass, and the length is
$l$
; -
5. The whole system’s center of mass is located at the spherical shell’s geometric center;
-
6. The sphere’s rotational motion around the vertical
$oz$
axis is ignored; -
7. The swing direction is specified as the positive direction.
Assuming the spherical shell center of mass’s angular velocity is
$\boldsymbol{\omega }_{1} = (\!- \!\dot{\theta }, - \dot{\varphi },0)^{T}$
, then the corresponding velocity is
$\boldsymbol{v}_1= r{\boldsymbol{\omega }_1} ={(\!\!- \!r\dot \theta, - r\dot \varphi,0)^T}$
, where r is sphere radius. Thus, the spherical shell’s Lagrangian function is
\begin{equation}{L_1} ={E_1} -{P_1} \end{equation}
where
$E_1$
describes the sum of the spherical shell’s horizontal displacement and rotational kinetic energies. Since the spherical shell is equivalent to a uniform thin spherical shell, its moment of inertia is
${I_1} = \dfrac{2}{3}{m_1}{r^2}$
. Thereby,
\begin{equation} E_{1} = \frac{1}{2}m_{1}\left \| \boldsymbol{v}_{1} \right \|^{2} + \frac{1}{2}I_{1}\left \|\boldsymbol{\omega }_{1} \right \|^{2} \end{equation}
where
$P_1$
is the gravitational potential energy of the spherical shell, and the rolling plane is a zero-potential energy surface, then:
${P_1} ={m_1}gr$
.
Assuming that the pendulum’s angular velocities in the
$o -{x_0}{y_0}{z_0}$
and
$o -{x_1}{y_1}{z_1}$
coordinate systems are denoted by
$\boldsymbol{\omega }_2$
and
$\boldsymbol{\omega }$
, respectively, then,
\begin{equation}{\boldsymbol{\omega }_2} ={\boldsymbol{\omega }_1} + R \boldsymbol{\omega} ={(\dot \alpha \cos\varphi - \dot \theta,\dot \alpha \sin\theta \sin\varphi + \dot \beta \cos\theta - \dot \varphi,\dot \alpha \cos\theta \sin\varphi - \dot \beta \sin\theta )^T} \end{equation}
where
$\boldsymbol{\omega } ={(\dot \alpha,\;\dot \beta,0)^T}$
. The position vector of the pendulum’s center of mass is described as
$\boldsymbol{r}_2$
and
$\boldsymbol{r}$
in the
$o -{x_0}{y_0}{z_0}$
and
$o -{x_1}{y_1}{z_1}$
coordinate systems, respectively, then
\begin{equation}{\boldsymbol{r}_2} = \boldsymbol{R}\boldsymbol{r} \end{equation}
Since the pendulum can swing around the
${\rm{o}}{x_1}$
and
${\rm{o}}{y_1}$
axes, the position of the pendulum’s center of mass changes in the three-dimensional space, then
\begin{equation} \boldsymbol{r} = x\boldsymbol{i} + y\boldsymbol{j} + z\boldsymbol{k} = \left ({\! - l\sin\beta } \right )\!\boldsymbol{i} + \left ({l\sin\alpha } \right )\boldsymbol{j} - \left ({l\cos\alpha \cos\beta } \right )\!\boldsymbol{k} \end{equation}
where
$ ({i,j,k} )$
describes a unit vector. The linear velocity of the pendulum’s center of mass is
\begin{equation} \boldsymbol{v}_{2} = \boldsymbol{v}_{1} + \boldsymbol{\omega }_{2} \times \boldsymbol{r}_{2} \end{equation}
The Lagrangian function of the pendulum is
${L_2} ={E_2} -{P_2}$
, where
$E_2$
is the sum of the pendulum’s horizontal movement and rotational kinetic energies, which is
\begin{equation} E_{2} = \frac{1}{2}m_{2}\left \|{\boldsymbol{v}_{2}} \right \|^{2} + \frac{1}{2}I_{2}\left \|{\boldsymbol{\omega }_{2}} \right \|^{2} \end{equation}
where
$I_2$
is the pendulum’s moment of inertia,
${I_2} = \dfrac{1}{{12}}{m_2}{l^2} + \dfrac{1}{4}{m_2}{l^2}$
.
$P_2$
describes the sum of the pendulum’s potential energy,
${P_3} ={m_2}g ({r - l\cos \alpha \cos \beta } )$
. Combining Eqs. (1)–(8), the whole XK-1 SR system’s Lagrangian function is attained as
\begin{equation} L = \frac{1}{2}m_{1}\left \| \boldsymbol{v}_{1} \right \|^{2} + \frac{1}{2}I_{1}\left \| \boldsymbol{\omega }_{1} \right \|^{2} - m_{1}gr + \frac{1}{2}m_{2}\left \| \boldsymbol{v}_{2} \right \|^{2} + \frac{1}{2}I_{2}\left \| \boldsymbol{\omega }_{2} \right \|^{2} - m_{2}g\left ( r - l{\cos \alpha }{\cos \beta } \right )\, \end{equation}
From the Lagrangian approach, the pendulum swing angles
$\alpha$
,
$\beta$
and spherical shell rolling angles
$\theta$
,
$\varphi$
are taken as four generalized coordinates, that is,
$\boldsymbol{q} ={(\theta,\alpha,\varphi,\beta )^T}$
. The SR movement is generally driven by the two motors torques,
$\tau_x$
and
$\tau_y$
output by the two sets of motors. Therefore, the input torque of the system is
$\boldsymbol{Q} ={({\tau_x},{\tau_x},{\tau_y},{\tau_y})^T}$
. Now, the motion’s Euler–Lagrangian equation for the simplified model is attained as
\begin{equation} \left \{ \begin{matrix} \begin{matrix}{\dfrac{d}{dt}\left ( \dfrac{\partial L}{\partial \dot{q_{1}}} \right ) - \dfrac{\partial L}{\partial q_{1}} = \tau_{x}} \\ \\[-8pt]{\dfrac{d}{dt}\left ( \dfrac{\partial L}{\partial \dot{q_{2}}} \right ) - \dfrac{\partial L}{\partial q_{2}} = \tau_{x}} \\ \\[-8pt] \end{matrix} \\ \\[-8pt] \begin{matrix}{\dfrac{d}{dt}\left ( \dfrac{\partial L}{\partial \dot{q_{3}}} \right ) - \dfrac{\partial L}{\partial q_{3}} = \tau_{y}} \\ \\[-8pt] {\dfrac{d}{dt}\left ( \dfrac{\partial L}{\partial \dot{q_{4}}} \right ) - \dfrac{\partial L}{\partial q_{4}} = \tau_{y}} \\ \\[-8pt] \end{matrix} \\ \end{matrix} \right. \end{equation}
According to Eq. (10), the SR’s dynamic equation is described with the following matrix form:
\begin{equation} \boldsymbol{M}(\boldsymbol{q}) \ddot{\boldsymbol{q}}+\boldsymbol{V}(\boldsymbol{q}, \dot{\boldsymbol{q}}) \dot{\boldsymbol{q}}=\boldsymbol{u} \end{equation}
\begin{equation} \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} m_{1} & a & b & c \\ \\[-7pt] a & m_{2} & d & e \\ \\[-7pt] b & d & m_{3} & f \\ \\[-7pt] c & e & f & m_{4} \\ \\[-7pt] \end{array}\right) \cdot \begin{bmatrix} \begin{matrix} \ddot{q_{1}} \\ \\[-7pt] \ddot{q_{2}} \\ \\[-7pt] \ddot{q_{3}} \\ \\[-7pt] \ddot{q_{4}} \\ \end{matrix} \\ \end{bmatrix} + \begin{bmatrix} \begin{matrix} V_{11} \\ \\[-7pt] V_{21} \\ \\[-7pt] V_{31} \\ \\[-7pt] V_{41} \\ \end{matrix} \\ \end{bmatrix} = \begin{bmatrix} \begin{matrix} \tau_{x} \\ \\[-7pt] \tau_{x} \\ \\[-7pt] \tau_{y} \\ \\[-7pt] \tau_{y} \\ \\[-7pt] \end{matrix} \\ \end{bmatrix} \end{equation}
where
$\boldsymbol{M} \in \mathbb{R}^{4 \times 4}$
,
$\boldsymbol{V} \in \mathbb{R}^{4 \times 1}$
. The specific elements of matrices
$\boldsymbol{M}$
and
$\boldsymbol{V}$
are in Appendix 1.
2.2. State-space equations of the pendulum SR
The dynamic equation (12) is a coupled nonlinear system. It is challenging to construct the motion control system of type XK-I described by Eq. (12). According to the method of [Reference Yang, Tao, Huo, Neusypin, Zhang, Zhang and Guo4, Reference Kayacan, Bayraktaroglu and Saeys8, Reference Liu, Lin, Yang, Nazarova and Huo20], the XK-I SR’s rolling and steering motion control can be separated; that is, the driving force produced by a group of motors makes the XK-I move in a straight line. If the XK-I type needs steering motion during linear motion, it is only necessary to input control signals to another group of motors. Therefore, the SR’s linear motion is driven by the torque
$\tau_x$
output by one set of motors. In contrast, the torque
$\tau_y$
output by another set of motors derives the steering movement. This decoupling approach converts the nonholonomic system into two decoupled subsystems, simplifies the SR dynamics equations, and improves the maneuverability of the system.
Besides, in the XK-I type’s actual movement process, the pendulum’s swing angle is slight, and the Taylor series is utilized to linearize the sine and cosine functions of
$\alpha$
and
$\beta$
.
When the SR moves in a straight line, only
${u_1} ={\tau_x}$
is the torque provided by the linear motor. Therefore, the pendulum’s swinging angle is only
$\alpha$
, and the spherical shell’s rolling angle is only
$\theta$
. Then, the state variable of the SR’s linear motion is
${\boldsymbol{x}_1} ={ [{{q_1},{{\dot q_1}},{q_2},{{\dot q_2}}} ]^T}$
. Since the pendulum swing angle
$\alpha$
is small during the XK-I spherical robot’s actual movement,
$\sin\alpha$
and
$\cos\alpha$
are described with their corresponding Taylor series. Since for small swinging angles, we have
${\sin\alpha }=\alpha$
and
${\cos\alpha }=1$
, the state-space expression of the SR’s linear motion can be attained as
\begin{equation} \left \{ \begin{matrix}{\dot{\boldsymbol{x}_{1}} = \boldsymbol{A}_{1}\boldsymbol{x}_{1} + \boldsymbol{B}_{1}u_{1}} \\ \\[-7pt] {\boldsymbol{y}_{1} = \boldsymbol{C}_{1}\boldsymbol{x}_{1} + \boldsymbol{D}_{1}u_{1}} \\ \end{matrix} \right. \end{equation}
where
$\boldsymbol{A}_{1} \in \mathbb{R}^{4 \times 4}$
,
$\boldsymbol{B}_{1} \in \mathbb{R}^{4 \times 1}$
,
$\boldsymbol{C}_{1} \in \mathbb{R}^{4 \times 4}$
, and
$\boldsymbol{D}_{1} = \lbrack \textbf{0}_{4 \times 1} \rbrack$
. The specific elements of matrices
$\boldsymbol{A}_1$
,
$\boldsymbol{B}_1$
and
$\boldsymbol{C}_1$
are in Appendix 2.
In the same way, when the spherical robot turns and moves, only
${u_2} ={\tau_y}$
is the torque provided by the motor. Therefore, the pendulum’s swing angle is only
$\beta$
, and the spherical shell’s rolling angle is only
$\varphi$
. Since the pendulum swing angle
$\beta$
is small during the actual movement of the XK-I SR,
$\sin\beta$
and
$\cos\beta$
can be represented with their corresponding Taylor series. The state variable of the SR turning motion is
${\boldsymbol{x}_2} ={ [{{q_3},{{\dot q_3}},{q_4},{{\dot q_4}} } ]^T}$
. Then, the SR turning motion’s state-space equations are
where
$\boldsymbol{A}_{2} \in \mathbb{R}^{4 \times 4}$
,
$\boldsymbol{B}_{2} \in \mathbb{R}^{4 \times 1}$
,
$\boldsymbol{C}_{2} \in \mathbb{R}^{4 \times 4}$
, and
$\boldsymbol{D}_{2} = \lbrack \textbf{0}_{4 \times 1} \rbrack$
. The specific elements of matrices
${\boldsymbol{A}_2}$
,
${\boldsymbol{B}_2}$
, and
${\boldsymbol{C}_2}$
are in Appendix 3.
3. Design of the XK-I adaptive motion controller
The design block diagram of the model reference adaptive controller of the XK-I SR is presented in Fig. 2. The controller mainly includes a reference model, adaptive regulator, and parameter debugger. Therefore, the XK-I SR controller’s reference model is first constructed according to the state space equation and the LQR. Here, we focus on analyzing the XK-I SR’s linear motion control, and the steering motion can be obtained similarly. In the XK-I SR’s motion control, we give the desired motion state
$\boldsymbol{x}_{1}^{\prime}$
, such as the angular velocity and angle of the spherical shell and the pendulum. It is hoped that the SR output is the same as the expected value at each moment. However, the actual process is imperfect, and there will be errors in the control at every moment. When only the error at the current moment is utilized in the controller design, the steady-state error is inevitable, degrading the control system’s performance. The desired system performance can be attained by minimizing the accumulated state errors of the XK-I SR in the entire movement time. Then, the controller is designed to minimize the cumulative error when the SR’s motion state reaches the desired value at a particular moment. Therefore, the quadratic performance cost function of the error state and control signals should be constructed for the state-space equations (13) or (14). Finally, the optimal control signal
$u_1$
or
$u_2$
is determined to minimize the following quadratic performance index
\begin{equation} J=\frac{1}{2} \int_{0}^{\infty }\left (\boldsymbol{e}_{1}^{T} \boldsymbol{Q}_{1} \boldsymbol{e}_{1}+R u_{1}^{2}\right ) d t \rightarrow \min \end{equation}
where
$\boldsymbol{e}_{1} \in \mathbb{R}^{4 \times 1}$
is the error vector,
$\boldsymbol{e}_{1} = \boldsymbol{x}_{1}^{\prime} - \boldsymbol{x}_{1}$
, and
$\boldsymbol{x}_{1}^{\prime}$
describes the given state vector,
$\boldsymbol{Q}_{1} \in \mathbb{R}^{4 \times 4}$
is a positive definite matrix and
$R$
describes a weighting constant. The control law that minimizes the performance index Eq. (15) is
\begin{equation}{u_1} ={R^{ - 1}}{\boldsymbol{B}_1}^T\boldsymbol{P}{e_1} = \boldsymbol{K}{\boldsymbol{e}_1} \end{equation}
where
$P$
is the Riccati equation’s positive definite symmetrical solution, then,
$\boldsymbol{A}_{1}^{\boldsymbol{T}}\boldsymbol{P} + \boldsymbol{P}\boldsymbol{A}_{1} + \boldsymbol{Q}_{1} - \boldsymbol{P}\boldsymbol{B}_{1}R^{-1}{\boldsymbol{B}_{1}^{\boldsymbol{T}}}\boldsymbol{P} = \textbf{0}$
,
$K$
is the control parameter of the optimal controller. Therefore, the reference model’s state-space equations are designed as
Adaptive control system of spherical robot.
\begin{equation} \left \{ \begin{matrix}{{\dot{\boldsymbol{x}}}_{1r} = \boldsymbol{A}_{1}\boldsymbol{x}_{1r} + \boldsymbol{B}_{1}u_{1}} \\ \\[-7pt] {\boldsymbol{y}_{1r} = \boldsymbol{C}_{1}\boldsymbol{x}_{1r}} \\ \end{matrix} \right. \end{equation}
Bringing Eq. (16) into the state equation (17) of the reference model, and assume that the state
$\boldsymbol{x}_{1}$
can be consistent with the state
$\boldsymbol{x}_{1r}$
of the reference model, we get,
\begin{equation} \left \{\begin{array}{c} \dot{\boldsymbol{x}}_{1 r}=\boldsymbol{A}_{1 r} \boldsymbol{x}_{1 r}+\boldsymbol{B}_{1 r} \boldsymbol{x}_{1}^{\prime } \\ \\[-8pt] \boldsymbol{y}_{1 r}=\boldsymbol{C}_{1} \boldsymbol{x}_{1 r} \end{array}\right. \end{equation}
where
${\boldsymbol{A}_{1r}} = ({\boldsymbol{A}_1} -{\boldsymbol{B}_1}\boldsymbol{K} )$
,
${\boldsymbol{B}_{1r}} ={\boldsymbol{B}_1}\boldsymbol{K}$
. In practice, it is difficult for
$\boldsymbol{x}_{1r}$
to be consistent with
$\boldsymbol{x}_{1}$
, then the error between the reference model output
$\boldsymbol{x}_{1r}$
and the state vector
$\boldsymbol{x}_{1}$
is
${ \boldsymbol{e}_{1r}} ={\boldsymbol{x}_{1r}} -{\boldsymbol{x}_1}$
. Combining Eqs. (13) and (18), we can get,
\begin{equation}{\dot{\boldsymbol{e}}}_{1r} ={\dot{\boldsymbol{x}}}_{1r} -{\dot{\boldsymbol{x}}}_{1} = \boldsymbol{A}_{1r}\boldsymbol{e}_{1r} + \left ({\boldsymbol{A}_{1r} - \boldsymbol{A}_{1}} \right )\boldsymbol{x}_{1} + \boldsymbol{B}_{1r}\boldsymbol{x}_{1}^{\prime} - \boldsymbol{B}_{1}u_{1} \end{equation}
The goal is to make
$\boldsymbol{x}_{1r}$
approach
$ \boldsymbol{x}_{1}$
, then
${\boldsymbol{e}_{1r}} = 0$
. Thus,
$u_1$
can be calculated by Eq. (19),
\begin{equation} u_{1} = \left ( \boldsymbol{B}_{1}^{T}\boldsymbol{B}_{1} \right)^{- 1}\boldsymbol{B}_{1}^{T}\left ({\boldsymbol{A}_{1r} - \boldsymbol{A}_{1}} \right )\boldsymbol{x}_{1} + \left ( \boldsymbol{B}_{1}^{T}\boldsymbol{B}_{1} \right)^{- 1}\boldsymbol{B}_{1}^{T}\boldsymbol{B}_{1r}\boldsymbol{x}_{1}^{\prime} = \boldsymbol{\rho }\boldsymbol{x}_{1} + \boldsymbol{\eta }\boldsymbol{x}_{1}^{\prime} \end{equation}
Thus, Eq. (16) can be rewritten as Eq. (20). Since the reference model is an ideal control system,
$\boldsymbol{\rho }$
and
$\boldsymbol{\eta }$
are ideal tuning parameters, and Eq. (20) is a desired control law. Additionally, since the error term is not included in the control law (20), it is challenging to eliminate the error by adjusting the parameters. Therefore, the actual control law of the XK-I robot’s linear motion under complex road conditions is considered as follows:
\begin{equation}{\hat{u}}_{1} = \hat{\boldsymbol{\rho }}\boldsymbol{x}_{1} + \hat{\boldsymbol{\eta }}\boldsymbol{x}_{1}^{\prime} \end{equation}
Multiplying Eq. (20) by matrix
$\boldsymbol{B}_1$
gives
\begin{equation} \boldsymbol{B}_{1}u_{1} = \left ({\boldsymbol{A}_{1r} - \boldsymbol{A}_{1}} \right )\boldsymbol{x}_{1} + \boldsymbol{B}_{1r}\boldsymbol{x}_{1}^{\prime} = \boldsymbol{B}_{1}\!\left ( \boldsymbol{\rho }\boldsymbol{x}_{1} + \boldsymbol{\eta }\boldsymbol{x}_{1}^{\prime} \right ) \end{equation}
Substitute formula (21) into deviation Eq. (19), extracting matrix
$\boldsymbol{B}_{1}$
, and combining formula (20), we can get
\begin{equation}{\dot{\boldsymbol{e}}}_{1r} = \boldsymbol{A}_{1r}\boldsymbol{e}_{1r} + \boldsymbol{B}_{1}\!\left ( \tilde{\boldsymbol{\rho }}\boldsymbol{x}_{1} + \tilde{\boldsymbol{\eta }}\boldsymbol{x}_{1}^{\prime} \right ) \end{equation}
where
$\boldsymbol{\tilde \rho }$
and
$\boldsymbol{\tilde \eta }$
are the differences between the ideal and actual adjustment coefficients; that is,
$\boldsymbol{\tilde \rho } = \boldsymbol{\rho } - \boldsymbol{\hat \rho }$
,
$\boldsymbol{\tilde \eta } = \boldsymbol{\eta } - \boldsymbol{\hat \eta }$
. In order to obtain the actual control law, the following Lyapunov function candidate is defined in terms of the errors
$\boldsymbol{e}_{1r}$
and
$\boldsymbol{\tilde \rho }$
,
$\boldsymbol{\tilde \eta }$
are
\begin{equation} \boldsymbol{V} = \boldsymbol{e}_{1r}^{T}{\boldsymbol{H}\boldsymbol{e}}_{1r} + \left | \boldsymbol{B}_{1} \right |\left ( \frac{1}{2}{\tilde{\boldsymbol{\rho }}}^{\boldsymbol{T}}\boldsymbol{N}_{1}\tilde{\boldsymbol{\rho }} +{\tilde{\boldsymbol{\eta }}}^{\boldsymbol{T}}\boldsymbol{N}_{2}\tilde{\boldsymbol{\eta }} \right ) \end{equation}
where
$\boldsymbol{N}_{1},{\boldsymbol{N}_2}$
are diagonal matrices, where their diagonal elements are positive. The matrix
$H$
satisfies the following relationship:
$\boldsymbol{A}_{1r}^{T}\boldsymbol{H} + \boldsymbol{H}\boldsymbol{A}_{1r} = - \boldsymbol{I}$
. Differentiating Eq. (24) gives
\begin{equation} \dot{\boldsymbol{V}} = \boldsymbol{e}_{1r}^{T}\left ( \boldsymbol{A}_{1r}^{T}\boldsymbol{H} + \boldsymbol{H}\boldsymbol{A}_{1r} \right )\boldsymbol{e}_{1r}+2\boldsymbol{e}_{1r}^{T}\boldsymbol{H}\boldsymbol{B}\left ( \tilde{\boldsymbol{\rho }}\boldsymbol{x}_{1} + \tilde{\boldsymbol{\eta }}\boldsymbol{x}_{1}^{\prime} \right ) + \left | \boldsymbol{B}_{1} \right |\left ({\tilde{\boldsymbol{\rho }}}^{\boldsymbol{T}}\boldsymbol{N}_{1}\dot{\tilde{\boldsymbol{\rho }}} +{\tilde{\boldsymbol{\eta }}}^{\boldsymbol{T}}\boldsymbol{N}_{2}\dot{\tilde{\boldsymbol{\eta }}} \right ) \end{equation}
Since the first term in Eq. (25) is negative, the following equality should be satisfied to ensure
$ \dot{\boldsymbol{V}} \le \textbf{0}$
,
\begin{equation} 2\boldsymbol{e}_{1r}^{T}\boldsymbol{H}\boldsymbol{B}\!\left ({\tilde{\boldsymbol{\rho }}\boldsymbol{x}_{1} + \tilde{\boldsymbol{\eta }}\boldsymbol{x}_{1}^{\prime}} \right ) + \left | \boldsymbol{B}_{1} \right |\left ({{\tilde{\boldsymbol{\rho }}}^{\boldsymbol{T}}\boldsymbol{N}_{1}\dot{\tilde{\boldsymbol{\rho }}} +{\tilde{\boldsymbol{\eta }}}^{\boldsymbol{T}}\boldsymbol{N}_{2}\dot{\tilde{\boldsymbol{\eta }}}} \right ) = 0 \end{equation}
Equation (26) has
$\boldsymbol{\tilde \rho },\boldsymbol{\tilde \eta }$
variables. If
$\boldsymbol{\tilde \rho }$
matches with
$\boldsymbol{\tilde \eta }$
, then Eq. (26) has a solution. Therefore, assuming
${\boldsymbol{N}_1} = {\boldsymbol{N}_2} = \boldsymbol{N}$
,
$- \boldsymbol{\rho }= \boldsymbol{\eta }$
, and
$ - \boldsymbol{\hat \rho } = \boldsymbol{\hat \eta }$
we can get
$ - \boldsymbol{\tilde \rho } = \boldsymbol{\tilde \eta }$
, simplifies Eq. (26) as
\begin{equation} - 2\boldsymbol{e}_{1r}^{T}\boldsymbol{H}\boldsymbol{B}\tilde{\boldsymbol{\rho }}\boldsymbol{e}_{1} + 2\left | \boldsymbol{B}_{1} \right |{\tilde{\boldsymbol{\rho }}}^{\boldsymbol{T}}\boldsymbol{N}\dot{\tilde{\boldsymbol{\rho }}} = 0 \end{equation}
When
$ - \dot{\hat{\boldsymbol{\rho }}} = \dot{\tilde{\boldsymbol{\rho }}}$
, the control law and coefficient adjustment rule are obtained as
\begin{equation}{\hat{u}}_{1} = \boldsymbol{\hat \rho }\boldsymbol{e}_{1} \end{equation}
\begin{equation} \dot{\hat{\boldsymbol{\rho }}} = - \text{sgn}\boldsymbol{B}_{1} \cdot \boldsymbol{N}^{- 1}\boldsymbol{e}_{1r}^{T}{\boldsymbol{H}\boldsymbol{e}}_{1} \end{equation}
4. Experimental analysis and discussion
In order to evaluate the performance of the designed adaptive motion controller, simulation and real experiments of the XK-I SR are performed. The simulation experiments in this study are divided into two groups: (1) The ideal situation without disturbance and noise, in which the designed controller speed is verified and (2) random interference noise is added to the feedback loop to evaluate the anti-interference capability of the constructed controller. In order to evaluate the efficiency of the MRAC, the simulation experiments are compared with LQR, Fuzzy-PID, PSO-ADRC, and HSMC.
4.1. Simulation results
Table I shows the parameter settings in the simulation process based on the XK-I SR’s physical prototype (as presented in Fig. 3). During the simulation, the linear and steering motions are carried out independently; that is, the motor groups M1 and M2 control the linear and steering motions, respectively. Besides, the rotational speed of the XK-I SR’s linear and steering motions is set to 3 rad/s, and random signals are selected as the interference source, where the interference amplitude is less than 1 rad/s.
Robot parameters.
XK-I spherical robot prototype.
In the simulations, when the LQR controller is utilized to control the motion effect, the cumulative amount of the state change of the system should be as small as possible during the change of the state variable of the robot from the initial value to the expected value. In this way, the influence of the disturbance caused by the dynamic change in the robot’s state can be reduced. Therefore, Eq. (15) is adopted as the system’s objective function.
When designing the Fuzzy-PID controller, the input error is quantified in the range (−6,6). Further, the error is fuzzified to determine the fuzzy subset as {NL, NM, NS, ZE, PS, PM, PL}. Besides, the triangular membership function is selected. The specific fuzzy rules are presented in Tables II–IV.
Fuzzy rules for Kp parameters in Fuzzy-PID.
Fuzzy rules for Ki parameters in Fuzzy-PID.
Fuzzy rules for Kd parameters in Fuzzy-PID.
The PSO-ADRC controller adopts the method provided in ref. [Reference Liu, Lin, Yang, Nazarova and Huo20]. In the particle swarm algorithm, both the population size and the number of iterations are set to 50 so that the overall number of iterations of the ADRC controller reaches 2500, and the random generation of parameters ranges from 0.0001 to 100000 values. The controller tuning parameters mainly select the velocity factor r and the system integration time h in the tracking differentiator (TD) in ADRC, and the compensation coefficient
$b_{0}$
and the observer parameters
$b_{01}$
,
$b_{02}$
,
$b_{03}$
in the extended state observer (ESO) can be considered a coupling of two second-order subsystems. The HSMC controller selects the relationship between the spherical shell and the pendulum in the state-space equations as the target-controlled object and calculates the difference between the output velocity and acceleration of the spherical shell and the ideal signal.
During simulation, the ideal model’s output speed is selected as the reference speed signal (3 rad/s). The linear motion simulation results are presented in Figs. 4 and 5. As shown in Figs. 4 and 5, the presented MRAC strategy makes the output velocity approach the reference signal quickly under ideal and interference conditions. In terms of the overall speed, in the ideal case, we have MRAC>PSO-ADRC>Fuzzy-PID>LQR>HSMC; in the case of interference, we have MRAC>PSO-ADRC>Fuzzy-PID>HSMC>LQR. In terms of error volatility, Fuzzy-PID, LQR, and HSMC have noticeable fluctuations with a specific magnitude under ideal conditions; in the case of interference, in the error change at 0.5 s,
$\;{e_{MRAC}} \lt \;{e_{PSO - ADRC}} \lt \;{e_{Fuzzy - PID}} \lt \;{e_{HSMC\;}} \lt{e_{LQR\;}}$
. Fuzzy-PID, LQR, and HSMC have evident error fluctuations in the next few seconds.
Linear motion simulation of the SR: (a) no interference environment; (b) interference environment.
Errors of each controller and reference signal in linear motion of the SR: (a) no interference environment; (b) interference environment.
Similarly, the steering motion simulation results can be obtained, as shown in Figs. 6 and 7. Under ideal conditions, in terms of fast performance, we have MRAC>PSO-ADRC>Fuzzy-PID>LQR>SMC. After adding interference, the five controllers all have significant fluctuations in the initial period at 0.5 s, while the MRAC controller fluctuation is close to the impulse signal. It can be seen that the reference model can be quickly tracked even under the action of disturbance, which can efficiently counter the disturbance and restore the system’s stability quickly.
Steering motion simulation of SR: (a) interference-free environment; (b) interference environment.
Errors between the controllers and the reference signal in the spherical robot’s linear motion: (a) no interference environment; (b) interference environment.
4.2. Linear motion experiment
Further, real experiments of linear motion control are performed on smooth and uneven ground. The motor velocity’s expected values of the SR are set as 50 rps, 100 rps, and 150 rps, respectively. The experimental results are presented in Figs. 8 and 9, and the analysis data are presented in Table V. Figures 10 and 11 show the process of using the MRAC controller to control the XK-I linear motion under smooth and uneven ground conditions.
The spherical robot’s linear motion control under the smooth ground condition.
The spherical robot’s linear motion control under the uneven ground condition.
Process diagram of the XK-I SR moving on smooth ground.
Process diagram of the XK-I SR moving on uneven ground.
Control error and stability comparison table of four controllers under different conditions.
Figure 8 presents the SR test results on smooth ground. Due to the error of modeling and model decoupling, the system parameters may change during the SR’s actual motion. Figure 8 shows that Fuzzy-PID, MRAC, and HSMC can control the movement speed to reach the desired value within 1s, while the transient response efficiency of the LQR and PSO-ADRC controllers is low. Among the Fuzzy-PID, MRAC, and HSMC controllers, the Fuzzy-PID controller leads to a considerable overshoot in the transient response. The transient response fluctuations of the Fuzzy-PID controller are more significant than the MRAC. The HSMC controller always has evident static error problems and is slower than the MRAC and Fuzzy-PID controllers. Besides, Table V indicates that the movement speed variance provided by the MRAC is 0.196, 0.195, and 0.296, respectively, which is the smallest variance compared to the other four controllers. These features demonstrate the good control effect of the MRAC controller.
Figure 9 presents the SR test results on the bumpy ground. Under strong interference, all five controllers exhibit significant fluctuations. In terms of speed, MRAC, Fuzzy-PID, and HSMC controllers all provide good results. Besides, the designed MRAC’s anti-interference capability is at a medium level, which is stronger than the LQR controller and the same as the HSMC controller but weaker than the PSO-ADRC and Fuzzy-PID controllers.
In summary, in the simulation and real experiments, the MRAC controller can quickly approach the control motor’s expected value of the SR and quickly restore stability in the presence of disturbances. The simulation and smooth ground test results demonstrate that the MRAC still has good control ability under motion disturbances.
The data in Table V are divided into two cases: smooth ground and uneven ground. Under smooth ground, the mean value of MRAC is closest to the target value. The error is within 2 rps. The control effect of other controllers is quite different from the set value. The control signal error follows MRAC<Fuzzy-PID<HSMC<PSO-ADRC<LQR. In terms of variance, the volatility of MRAC is also the smallest. The fluctuation of the control signal follows MRAC<Fuzzy-PID<HSMC<PSO-ADRC<LQR. The overall reflection of the fast performance of MRAC is the best.
Under uneven ground, there is external interference. As presented in Table V, both the mean and variance of MRAC are in the middle of several controllers. Figure 8 intuitively shows that the fast performance of MRAC is very prominent, while its anti-interference performance is not optimal. Besides, the noise immunity performances of Fuzzy-PID, HSMC, and PSO-ADRC are almost the same.
The ideal and interference environments in Table V reflect the controller’s fast and anti-interference performance, respectively. The data reflect that the designed MRAC has excellent fast performance, while the anti-interference is insufficient. In the real environment, various disturbances, such as the changes in the external environment, will inevitably affect the SRs’ motion control performance. In the future, combined with the MRAC, a motion control strategy for SRs that can effectively suppress external disturbances, such as a disturbance observer-based control method [Reference Chen, Yang, Guo and Li33], will be studied. In this method, the difference between the actual object and the nominal model induced by external disturbances and model parameter variations is equivalent to the control input, and equivalent compensation is incorporated into the controller to eliminate the disturbances.
4.3. Rolling of the XK-I SR over curvilinear trajectories
In order to verify the practical application of the presented control approach, an experiment was performed in which the MRAC controlled the XK-I SR to move along a curve on smooth ground. The trajectory diagram of moving along the curve is presented in Fig. 12, and the actual scene diagrams of the XK-I SR moving along the curve are shown in Figs. 13 and 14. Figure 12 shows a particular deviation when the MRAC controls the XK-I SR to move along the curve, indicating that it cannot accurately track the curve movement. During the MRAC design, controlling the angular velocity and angle of the XK-I SR’s pendulum swing and spherical shell rolling is considered while ignoring its position and attitude control, deviating from the control of the XK-I along the curve. In the future, the expected value in the MRAC can be set to the trajectory coordinate value moving along the curve to reconstruct the control law and realize the XK-I SR’s path following control.
Experimental results of rolling over curvilinear trajectories.
The actual scene of the XK-I SR moving along a circular trajectory.
The actual scene of the XK-I SR moving along the trajectory.
5. Conclusions
The current work presented a Lyapunov-based MRAC control algorithm. The SR’s omnidirectional dynamic model is established via the Euler–Lagrange method by taking the pendulum’s swing angle and the spherical shell’s rolling angle as the generalized coordinates. According to the steering principle of the car during driving, the XK-I SR’s dynamic model is decoupled into linear and steering models. Besides, MRAC is utilized to design the SR’s adaptive motion control law, and a Lyapunov-based adaptation mechanism for the controller parameters is designed to overcome the poor control effect induced by the SR’s inaccurate dynamic model and parameter uncertainty. Finally, to evaluate the effect of the designed motion controller, the simulation and real experiments are compared with LQR, Fuzzy-PID, PSO-ADRC, and HSMC. Experiments indicate that the designed controller can be successfully applied to the SR control system. In the case of interference, it can quickly restore stability. Especially when the system parameters are uncertain, the SR can still be well controlled to achieve the expected effect after decoupling the SR. In the future, the SR’s control system should be further studied, mainly to solve the anti-interference ability problem and trajectory tracking of the MRAC controller.
Author contributions
All authors proposed the research. Jianwen Huo conceived and designed the study. Rui Lin conducted data gathering and performed statistical analyses. Jianwen Huo and Rui Lin wrote the article.
Financial support
This research was funded by the Research Fund of Southwest University of Science and Technology (No. 22zx7109).
Conflicts of interest
The authors declare that they have no conflicts of interest.
Ethical approval
None.
Acknowledgement
The authors would like to thank the financial support from the Research Fund of Southwest University of Science and Technology.
Appendix 1
$m_{1} = ({{m\_1} +{m\_2}} ){r^2} +{I\_1} +{I\_2} + 2{m\_2}{l^2}\sin \varphi \cos \varphi \cos \alpha \sin \beta \cos \beta +{m\_2}{l^2} (\!\sin{\alpha ^2} $
$+ \sin{\varphi ^2}\sin{\beta ^2} + \cos{\varphi ^2}\cos{\alpha ^2}\cos{\beta ^2} );$
$m_{2} ={I_2} +{m_2}{l^2} (\!{\sin{\alpha ^2} + \cos{\alpha ^2}\cos{\beta ^2}} );$
$m_{3} = ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2} (\!\cos{\varphi ^2}\sin{\beta ^2} + \sin{\theta ^2}\sin{\alpha ^2} - 2\sin \varphi \cos \varphi \cos \alpha \sin \beta \cos \beta $
$+ \cos{\alpha ^2}\sin{\varphi ^2}\cos{\beta ^2} + \cos{\theta ^2}\sin{\varphi ^2}\sin{\beta ^2} + 2\sin \theta \cos \theta \sin \alpha \cos \alpha \cos \beta \cos \varphi$
$+ \cos{\theta ^2}\cos{\alpha ^2}\cos{\varphi ^2}\cos{\beta ^2} + 2\sin \theta \cos \theta \sin \varphi \sin \alpha \sin \beta + 2\cos{\theta ^2}\sin \varphi \cos \varphi \cos \alpha \sin \beta \cos \beta );$
$m_{4} ={I_2} +{m_2}{l^2} (\!{\sin{\beta ^2} + \cos{\alpha ^2}\cos{\beta ^2}} );$
$ a = -{I_2}\cos \varphi +{m_2}lr\sin \varphi \sin \alpha +{m_2}{l^2} ({ -\!\sin \varphi \cos \alpha \sin \beta \cos \beta -\!\cos \varphi \cos{\alpha ^2}\cos{\beta ^2} -\!\cos \varphi \sin{\alpha ^2}} );$
$b ={m_2}{l^2} (\!\sin \theta \sin{\varphi ^2}\sin \beta \cos \beta \cos \alpha -\!\sin \theta \sin \varphi \cos \varphi \sin{\beta ^2} -\!\cos \theta \sin \varphi \sin \alpha \cos \alpha \cos \beta $
$+ \cos \theta \cos \varphi \sin \alpha \sin \beta -\!\sin \theta \cos{\varphi ^2}\cos \alpha \sin \beta \cos \beta + \sin \theta \sin \varphi \cos \varphi \cos{\alpha ^2}\cos{\beta ^2} );$
$ c ={m_2}lr (\!{\sin \varphi \sin \beta + \cos \varphi \cos \alpha \cos \beta } ) +{m_2}{l^2} ({ + \sin \varphi \sin \alpha \cos \alpha \cos \beta -\!\cos \varphi \sin \beta \sin \alpha } );$
$d = -{I_2}\sin \theta \sin \varphi +{m_2}lr ({ -\!\cos \theta \cos \alpha \cos \beta -\!\sin \theta \cos \varphi \sin \alpha } ) +{m_2}{l^2} (\!-\!\cos \theta \sin \alpha \sin \beta $
$-\!\sin \theta \sin \varphi \sin{\alpha ^2} + \sin \theta \cos \varphi \cos \alpha \sin \beta \cos \beta -\!\sin \theta \sin \varphi \cos{\alpha ^2}\cos{\beta ^2} );$
$ e ={m_2}{l^2} (\!{\sin \alpha \sin \beta } );$
$f= -{I_2}\cos \theta +{m_2}lr ({ -\!\sin \theta \cos \varphi \sin \beta -\!\sin \theta \sin \varphi \cos \alpha \cos \beta } ) +{m_2}{l^2} (\!-\!\cos \theta \sin{\beta ^2} $
$-\!\cos \theta \cos{\alpha ^2}\cos{\beta ^2} -\!\sin \theta \sin \varphi \sin \alpha \sin \beta -\!\sin \theta \sin \varphi \sin \alpha \cos \alpha \cos \beta .);$
$V_{11} = m_{2}rl ({-{\dot \alpha } ({{\sin \varphi }{\sin \alpha }} )^{\prime} + \dot{\beta } ({{\sin \varphi }{\sin \beta } +{\cos \varphi }{\cos \alpha }{\cos \beta }} )^{\prime}} )$
$+ m_{2}l^{2} ( \dot{\alpha } ({-{\sin \varphi }{\cos \alpha }{\sin \beta }{\cos \beta } -{\cos \varphi }{\cos \alpha }^{2}{\cos \beta }^{2} -{\cos \varphi }{\sin \alpha }^{2}} )^{\prime} $
$+ \dot{\beta } ({-{\cos \varphi }{\sin \alpha }{\sin \beta } +{\sin \varphi }{\sin \alpha }{\cos \alpha }{\cos \beta }} )^{\prime} + \dot{\theta } ( 2{\sin \varphi }{\cos \varphi }{\cos \alpha }{\sin \beta }{\cos \beta } $
$+{\cos \varphi }^{2}{\cos \alpha }^{2}{\cos \beta }^{2} +{\sin \varphi }^{2}{\sin \beta }^{2} +{\sin \alpha }^{2} )^{\prime} + \dot{\varphi } (\!-{\sin \theta }{\sin \varphi }{\cos \varphi }{\sin \beta }^{2} $
$+{\sin \theta }{\sin \varphi }^{2}{\cos \alpha }{\sin \beta }{\cos \beta }{{+ \cos }\theta }{\cos \varphi }{\sin \alpha }{\sin \beta } -{\cos \theta }{\sin \varphi }{\sin \alpha }{\cos \alpha }{\cos \beta }$
$-{\sin \theta }{\cos \varphi }^{2}{\cos \alpha }{\sin \beta }{\cos \beta } +{\sin \theta }{\sin \varphi }{\cos \varphi }{\cos \alpha }^{2}{\cos \beta }^{2} )^{\prime}$
$+ I_{2} ({\dot{\alpha }\dot{\varphi }{\sin \varphi } + \dot{\beta }\dot{\varphi }{\sin \theta } - \dot{\alpha }\dot{\varphi }{\cos \theta }{\sin \varphi }} .)+ m_{2} ( Z_{1} \cdot (\!- \dot{\varphi }l{\sin \theta }{\sin \varphi }{\sin \beta } + \dot{\varphi }l{\cos \theta }{\sin \alpha } $
$- \dot{\varphi }l{\sin \theta }{\cos \varphi }{\cos \alpha }{\cos \beta } ) + Z_{2} \cdot (\!- \dot{\alpha }l{\sin \theta }{\cos \alpha }{\cos \beta }+ \dot{\alpha }l{\cos \theta }{\cos \varphi }{\sin \alpha } $
$+ \dot{\beta }l{\cos \theta }{\cos \varphi }{\sin \beta } - \dot{\beta }l{\cos \theta }{\sin \varphi }{\cos \alpha }{\cos \beta }+ \dot{\theta }l{\sin \theta }{\sin \varphi }{\sin \beta }$
$- \dot{\theta }l{\cos \theta }{\sin \alpha } + \dot{\theta }l{\sin \theta }{\cos \varphi }{\cos \alpha }{\cos \beta } )+ Z_{3} \cdot (\!- \dot{\alpha }l{\sin \theta }{\cos \varphi }{\sin \alpha } $
$- \dot{\alpha }l{\cos \theta }{\cos \alpha }{\cos \beta }+ \dot{\theta }l{\cos \theta }{\sin \varphi }{\sin \beta } + \dot{\theta }l{\sin \theta }{\sin \alpha } + \dot{\theta }l{\cos \theta }{\cos \varphi }{\cos \alpha }{\cos \beta }$
$- \dot{\beta }l{\sin \theta }{\cos \varphi }{\sin \beta } + \dot{\beta }l{\sin \varphi }{\sin \theta }{\cos \alpha }{\cos \beta } );$
${V_{21}} ={I_2} ({\dot \theta \dot \varphi \sin \varphi - \dot \theta \dot \varphi \cos \theta \sin \varphi -{{\dot \varphi }^2}\sin \theta \cos \varphi } )$
$ +{m_2}rl ({\dot \theta{{ (\!{\sin \varphi \sin \alpha } )}^{\rm{\prime}}} + \dot \varphi{{ ({ -\!\cos \theta \cos \alpha \cos \beta -\!\sin \theta \cos \varphi \sin \alpha } )}^{\rm{\prime}}}} )$
$+{m_2}{l^2} ( \dot \alpha{{ (\!{\sin{\alpha ^2} + \cos{\alpha ^2}\cos{\beta ^2}} )}^{\rm{\prime}}} + \dot \beta{{ (\!{\sin \alpha \sin \beta } )}^{\rm{\prime}}} $
$ + \dot \theta{{ ({ -\!\sin \varphi \cos \alpha \sin \beta \cos \beta -\!\cos \varphi \sin{\alpha ^2} -\!\cos \varphi \cos{\alpha ^2}\cos{\beta ^2}} )}^{\rm{\prime}}}$
$+ \dot \varphi (\!-\!\cos \theta \sin \alpha \sin \beta -\!\sin \theta \sin \varphi \sin{\alpha ^2} + \sin \theta \cos \varphi \cos \alpha \sin \beta \cos \beta $
$-\!\sin \theta \sin \varphi \cos{\alpha ^2}\cos \beta )^{\rm{\prime}}$
$+ m_{2} ( Z_{1} \cdot ({- \dot{\alpha }l{\sin \varphi }{\cos \alpha } + \dot{\beta }l{\cos \varphi }{\sin \alpha }{\cos \beta } + \dot{\varphi }l{\sin \theta }{\cos \alpha } - \dot{\varphi }l{\cos \theta }{\cos \varphi }{\sin \alpha }{\cos \beta }} ) $
$+ Z_{2} \cdot (\!- \dot{\alpha }l{\cos \theta }{\sin \alpha }{\cos \beta } + \dot{\alpha }l{\sin \theta }{\cos \varphi }{\cos \alpha } + \dot{\beta }l{\sin \theta }{\sin \varphi }{\sin \alpha }{\cos \beta } $
$- \dot{\theta }l{\sin \theta }{\cos \alpha } + \dot{\theta }l{\cos \theta }{\cos \varphi }{\sin \alpha }{\cos \beta } )$
$+ Z_{3} \cdot ( \dot{\alpha }l{\cos \theta }{\cos \varphi }{\cos \alpha } + \dot{\alpha }l{\sin \theta }{\sin \alpha }{\cos \beta } $
$ - \dot{\theta }l{\cos \theta }{\cos \alpha } - \dot{\theta }l{\sin \theta }{\cos \varphi }{\sin \alpha }{\cos \beta } + \dot{\beta }l{\sin \varphi }{\cos \theta }{\sin \alpha }{\cos \beta }$
$- \dot{\varphi }l{\sin \varphi }{\sin \alpha }{\cos \beta } )- m_{2}gl{\sin \alpha }{\cos \beta } ;$
${V_{31}} ={I_2} ({ - \dot \alpha \dot \theta \cos \theta \sin \varphi - 2\dot \alpha \dot \varphi \sin \theta \cos \varphi + \dot \beta \dot \theta \sin \theta + \dot \theta \dot \alpha \sin \varphi } )$
$+{m_2}rl ( \dot \alpha{{ ({ -\!\cos \theta \cos \alpha \cos \beta -\!\sin \theta \cos \varphi \sin \alpha } )}^{\rm{\prime}}} $
$+ \dot \beta{{ ({ -\!\sin \theta \cos \varphi \sin \beta + \sin \theta \sin \varphi \cos \alpha \cos \beta } )}^{\rm{\prime}}} )$
$+{m_2}{l^2} ( \dot \alpha (\!\sin \theta \cos \varphi \cos \alpha \sin \beta \cos \beta -\!\cos \theta \sin \alpha \sin \beta -\!\sin \theta \sin \varphi \sin{\alpha ^2} $
$ -\!\sin \theta \sin \varphi \cos{\alpha ^2}\cos{\beta ^2}$
$+ \dot \beta{ ({ -\!\cos \theta \sin{\beta ^2} -\!\cos \theta \cos{\alpha ^2}\cos{\beta ^2} -\!\sin \theta \sin \varphi \sin \alpha \sin \beta -\!\sin \theta \cos \varphi \sin \alpha \cos \alpha \cos \beta } )^{\rm{\prime}}}$
$+ \dot \varphi (\!{\cos{\theta ^2}\sin{\varphi ^2}\sin{\beta ^2} + 2\sin \theta \cos \theta \sin \varphi \sin \alpha \sin \beta + 2\cos{\theta ^2}\sin \varphi \cos \varphi \cos \alpha \sin \beta \cos \beta } $
$ + \sin{\theta ^2}\sin{\alpha ^2} + 2\sin \theta \cos \theta \cos \varphi \sin \alpha \cos \alpha \cos \beta + \cos{\theta ^2}\cos{\varphi ^2}\cos{\alpha ^2}\cos{\beta ^2} + \cos{\varphi ^2}\sin{\beta ^2}$
$- 2\sin \varphi \cos \varphi \cos \alpha \sin \beta \cos \beta + \sin{\varphi ^2}\cos{\alpha ^2}\cos{\beta ^2} )^{\rm{\prime}}$
$+ \dot \theta (\!-\!\sin \theta \sin \varphi \cos \varphi \sin{\beta ^2} + \cos \theta \cos \varphi \sin \alpha \sin \beta -\!\sin \theta \cos{\varphi ^2}\cos \alpha \sin \beta \cos \beta $
$+ \sin \theta \sin{\varphi ^2}\cos \alpha \sin \beta \cos \beta -\!\cos \theta \sin \varphi \sin \alpha \cos \alpha \cos \beta + \sin \theta \sin \varphi \cos \varphi \cos{\alpha ^2}\cos{\beta ^2} )^{\rm{\prime}})$
$+{m_2} ({Z_1} \cdot (\!- \dot \alpha l\cos \varphi \sin \alpha - \dot \beta l\cos \varphi \sin \beta + \dot \beta l\sin \varphi \cos \alpha \cos \beta + \dot \varphi l\cos \theta \cos \varphi \sin \beta $
$- \dot \varphi l\cos \theta \sin \varphi \cos \alpha \cos \beta ) +{Z_2}$
$\cdot (\!- \dot \alpha l\cos \theta \sin \alpha \cos \beta - \dot \alpha l\sin \theta \sin \alpha \sin \varphi - \dot \beta l\sin \theta \sin \varphi \sin \beta - \dot \beta l\sin \theta \cos \varphi \cos \alpha \cos \beta $
$- \dot \theta l\cos \theta \cos \varphi \sin \beta + \dot \theta l\cos \theta \sin \varphi \cos \alpha \cos \beta ) +{Z_2}$
$\cdot (\!- \dot \alpha l\cos \theta \sin \varphi \sin \alpha + \dot \theta l\sin \theta \cos \varphi \sin \beta - \dot \theta l\sin \theta \sin \varphi \cos \alpha \cos \beta - \dot \beta l\cos \theta \sin \varphi \sin \beta $
$- \dot \beta l\cos \theta \cos \varphi \cos \alpha \cos \beta + \dot \varphi l\sin \varphi \sin \beta + \dot \varphi l\cos \varphi \cos \alpha \cos \beta ) );$
${V_{41}} ={I_2}\dot \theta \dot \varphi \sin \theta +{m_2}rl ({\dot \theta{{ (\!{\sin \varphi \sin \beta + \cos \varphi \cos \alpha \cos \beta } )}^{\rm{\prime}}} + \dot \varphi{{ ({ -\!\sin \theta \cos \varphi \sin \beta + \sin \theta \sin \varphi \cos \alpha \cos \beta } )}^{\rm{\prime}}}} )$
$+{m_2}{l^2} ( \dot \alpha{{ (\!{\sin \alpha \sin \beta } )}^{\rm{\prime}}} + \dot \beta{{ (\!{\cos{\alpha ^2}\cos{\beta ^2} + \sin{\beta ^2}} )}^{\rm{\prime}}} $
$+ \dot \theta{{ (\!{\sin \varphi \sin \alpha \cos \alpha \cos \beta -\!\cos \varphi \sin \alpha \sin \beta } )}^{\rm{\prime}}}$
$+ \dot \varphi{{ ({ -\!\cos \theta \sin{\beta ^2} -\!\sin \theta \sin \varphi \sin \alpha \sin \beta -\!\sin \theta \cos \varphi \sin \alpha \cos \alpha \cos \beta -\!\cos \theta \cos{\alpha ^2}\cos{\beta ^2}} )}^{\rm{\prime}}} )$
$+{m_2} ({Z_1} \cdot ({ - \dot \beta l\sin \varphi \cos \beta + \dot \beta l\cos \varphi \cos \alpha \sin \beta + \dot \varphi l\cos \theta \sin \varphi \cos \beta - \dot \varphi l\cos \theta \cos \varphi \cos \alpha \sin \beta } ) $
$+{Z_2} \cdot (\!- \dot \alpha l\cos \theta \cos \alpha \sin \beta + \dot \beta l\sin \theta \cos \varphi \cos \beta + \dot \beta l\sin \theta \sin \varphi \cos \alpha \sin \beta $
$- \dot \theta l\cos \theta \sin \varphi \cos \beta + \dot \theta l\cos \theta \cos \varphi \cos \alpha \sin \beta ) +{Z_3}$
$\cdot ( \dot \alpha l\sin \theta \cos \alpha \sin \beta + \dot \theta l\sin \theta \sin \varphi \cos \beta - \dot \theta l\sin \theta \cos \varphi \cos \alpha \sin \beta + \dot \beta l\cos \theta \cos \varphi \cos \beta $
$+ \dot \beta l\cos \theta \sin \varphi \cos \alpha \sin \beta - \dot \varphi l\cos \varphi \cos \beta - \dot \varphi l\sin \varphi \cos \alpha \sin \beta ) -{m_2}gl\cos \alpha \sin \beta ;$
${Z_1} = - \dot \alpha l\sin \varphi \sin \alpha - \dot \beta l\sin \varphi \sin \beta - \dot \beta l\cos \varphi \cos \alpha \cos \beta + \dot \varphi l\cos \theta \sin \varphi \sin \beta + \dot \varphi l\sin \theta \sin \alpha$
$+ \dot \varphi l\cos \theta \cos \varphi \cos \alpha \cos \beta - r\dot \theta ;$
${Z_2} = \dot \alpha l\cos \theta \cos \alpha \cos \beta + \dot \alpha l\sin \theta \cos \varphi \sin \alpha + \dot \beta l\sin \theta \cos \varphi \sin \beta - \dot \beta l\sin \theta \sin \varphi \cos \alpha \cos \beta$
$- \dot \theta l\cos \theta \sin \varphi \sin \beta - \dot \theta l\sin \theta \sin \alpha - \dot \theta l\cos \theta \cos \varphi \cos \alpha \cos \beta - r\dot \varphi ;$
${Z_3} = \dot \alpha l\cos \theta \cos \varphi \sin \alpha - \dot \alpha l\sin \theta \cos \alpha \cos \beta + \dot \theta l\sin \theta \sin \varphi \sin \beta - \dot \theta l\cos \theta \sin \alpha$
$+ \dot \theta l\sin \theta \cos \varphi \cos \alpha \cos \beta + \dot \beta l\cos \theta \cos \varphi \sin \beta - \dot \beta l\sin \varphi \cos \theta \cos \alpha \cos \beta$
$- \dot \varphi l\cos \varphi \sin \beta + \dot \varphi l\sin \varphi \cos \alpha \cos \beta ;$
Appendix 2
\begin{equation*} A_{1}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 0 & 1 & 0 & 0 \\ \\[-8pt] 0 & 0 & F_{1} H_{1}^{-1} & 0 \\ \\[-8pt] 0 & 0 & 0 & 1 \\ \\[-8pt] 0 & 0 & F_{2} H_{1}^{-1} & 0 \end {array}\right];\quad B_{1}=\left[\begin{array}{c} 0 \\ \\[-8pt] F_{3} H_{1}^{-1} \\ \\[-8pt] 0 \\ \\[-8pt] F_{4} H_{1}^{-1} \end {array}\right[;\quad C_{1}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 0 & 1 & 0 & 0 \\ \\[-8pt] 0 & 0 & 0 & 0 \\ \\[-8pt] 0 & 0 & 0 & 0 \\ \\[-8pt] 0 & 0 & 0 & 0 \end {array}\right]; \end{equation*}
${F_1} = ({{{ ({{I_2} +{m_2}{l^2}} )}^3} -{{ ({{I_2} +{m_2}{l^2}} )}^2} ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ) + ({{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2}} ){m_2}gl;$
${F_2} = ({\rm{\raise .5ex\hbox{$\scriptstyle 1$}\kern -.1em/ \kern -.15em\lower .25ex\hbox{$\scriptstyle 4$} }}{m_2}{\rm{}}{{\rm{\raise .5ex\hbox{$\scriptstyle 1$}\kern -.1em/ \kern -.15em\lower .25ex\hbox{$\scriptstyle 2$} }}^2}{l^2}{r^2}(({m_1} +{m_2}){r^2} +{I_1} +{I_2}) - (({m_1} +{m_2}){r^2} +{I_1} +{I_2})(({m_1} +{m_2}){r^2}$
$+{I_1} +{I_2} +{m_2}{l^2})({I_2} +{m_2}{l^2}) + (({m_1} +{m_2}){r^2} +{I_1} +{I_2}){({I_2} +{m_2}{l^2})^2}){m_2}gl ;$
${H_1} = ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2}} ){ ({{I_2} +{m_2}{l^2}} )^3} - ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2}} ) ({{m_1} +{m_2}} ){r^2}$
$ +{I_1} +{I_2} +{m_2}{l^2} ){ ({{I_2} + {m_2}{l^2}} )^2}$
$+ ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2}} ) ({{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2}$
$+{ ({{I_2} +{m_2}{l^2}} )^4} ({ ({{m_1} +{m_2}} ){r^2}+{I_1} +{I_2} +{m_2}{l^2}} ) -{ ({{I_2} +{m_2}{l^2}} )^4} - 2{ ({{I_2} +{m_2}{l^2}} )^2}{m_2}^2{l^2}{r^2}$
$ + ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ) ({{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2} -{m_2}^4{l^4}{r^4};$
${F_3} ={ ({{I_2} +{m_2}{l^2}} )^2} ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ) -{ ({{I_2} +{m_2}{l^2}} )^2} - ({{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2};$
${F_4} = - ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2}} ) ({{I_2} +{m_2}{l^2}} ){m_2}lr +{ ({{I_2} +{m_2}{l^2}} )^2}{m_2}lr +{m_2}^3{l^3}{r^3}$
$-{ ({{I_2} +{m_2}{l^2}} )^2} ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ) +{ ({{I_2} +{m_2}{l^2}} )^3} + ({{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2};$
Appendix 3
\begin{equation*} \begin {aligned} \boldsymbol{A}_{2} & =\left [\begin {array}{c@{\quad}c@{\quad}c@{\quad}c} 0 & 1 & 0 & 0 \\ \\[-8pt] 0 & 0 & K_{1} H_{2}^{-1} & 0 \\ \\[-8pt] 0 & 0 & 0 & 1 \\ \\[-8pt] 0 & 0 & K_{2} H_{2}^{-1} & 0 \end {array}\right ] ; \quad \boldsymbol{B}_{2}=\left [\begin {array}{c} 0 \\ \\[-8pt] K_{3} H_{2}^{-1} \\ \\[-8pt] 0 \\ \\[-8pt] K_{4} H_{2}^{-1} \end {array}\right ] ; \quad \boldsymbol{C}_{2}=\left [\begin {array}{c@{\quad}c@{\quad}c@{\quad}c} 0 & 1 & 0 & 0 \\ \\[-8pt] 0 & 0 & 0 & 0 \\ \\[-8pt] 0 & 0 & 0 & 0 \\ \\[-8pt] 0 & 0 & 0 & 0 \end {array}\right ] ; \\ \end {aligned} \end{equation*}
${K_1} = ({{{ ({{I_2} +{m_2}{l^2}} )}^3} - ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ){{ ({{I_2} +{m_2}{l^2}} )}^2} +{m_2}^2{l^2}{r^2} ({{I_2} +{m_2}{l^2}} )} ){m_2}gl;$
${K_2} = ({{({({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} )}^2} - ({({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2} + ( ({{m_1} +{m_2}} ){r^2}$
$+{I_1} +{I_2} +{m_2}{l^2} )^2 ({{I_2} +{m_2}{l^2}} ) ){m_2};$
${K_3} = ({ - ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ){ ({{I_2} +{m_2}{l^2}} )^2} +{ ({{I_2} +{m_2}{l^2}} )^3} + ({{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2}$
$- ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ) ({{I_2} +{m_2}{l^2}} ){m_2}lr +{ ({{I_2} +{m_2}{l^2}} )^2}{m_2}lr +{m_2}^3{l^3}{r^3};$
${K_4} = - ({{I_2} +{m_2}{l^2}} ){ ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} )^2} + ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2}$
$+ ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ){ ({{I_2} +{m_2}{l^2}} )^2};$
${H_2} = ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ){ ({{I_2} +{m_2}{l^2}} )^3} -{ ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} )^2}{ ({{I_2} +{m_2}{l^2}} )^2}$
$+ ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ) ({{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2} + ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ){ ({{I_2} +{m_2}{l^2}} )^3}$
$-{ ({{I_2} +{m_2}{l^2}} )^4} - 2{ ({{I_2} +{m_2}{l^2}} )^2}{m_2}^2{l^2}{r^2} + ({ ({{m_1} +{m_2}} ){r^2} +{I_1} +{I_2} +{m_2}{l^2}} ) ({{I_2} +{m_2}{l^2}} ){m_2}^2{l^2}{r^2}$
$-{m_2}^4{l^4}{r^4};$
