2.1. Robot configuration

Figure 1 shows the robot structure. The WLM is a mobile robot with 10 degrees of freedom, comprising legs with attached wheels. The robot’s body is designed to accommodate additional equipment. To extend the robot’s workspace, a manipulator with two degrees of freedom, featuring two revolute joints, is mounted on the robot’s platform. The platform is designed and built to navigate diverse terrains and rough surfaces, enabling it to perceive its environment and transport objects. The maximum allowable payload for the manipulator is 5.5 kg. The manipulator’s workspace for the end-effector in the X-Z plane is 0.45 m². Each link is 13.5 cm long and weighs 0.23 kg. The permissible rotation angle for the first link ranges from 0° to 180°, while for the second link, it ranges from −45° to 225°.

Figure 1.

Schematic of the WLM parts.

The mechatronic parts include ten 150-watt Dynamixel motors. These consist of the hip joints, four XM540-W150-R servo motors with RS485 communication protocol, four XL430-W250-T motors with TTL communication protocol for the wheels, and two motors for the manipulator’s links. To determine the WLM’s posture, body position, and manipulator end-effector position, the robot uses MPU-6050 sensors for the main body and end-effector. This sensor is a 6-axis IMU combining an accelerometer and gyroscope. A relay module along with an SMPS2 Dynamixel power supply module is used to switch the robot’s power circuit on and off, which is connected to a 12V power supply that powers the robot’s motors. A buck converter module is employed to generate a stable DC current. The U2D2 board is a converter that establishes communication between the PC and the servo motors, enabling motor commands to be sent and feedback to be received. An interface board is provided to connect the motors to the U2D2 converter using the communication protocols. An Arduino Nano board with the I2C communication protocol is utilized to receive data from the MPU-6050 module and transmit it to the PC.

Table I.

Actual specifications of manipulator and motors.

Table I shows the actual parameters of the manipulator. To maintain the robot’s stability and balance and to prevent it from tipping over, the manipulator is allowed to move only in the X-Z plane. Also, an arm with lightweight links and motors with suitable power has been used for applications in the WLM.

2.2. Robot slippage

As illustrated in Figure 2, the deviation of the WLM from its planned trajectory is mainly caused by two factors. The first is a low slip coefficient at the wheel-ground interface, which leads to skidding. The second is rapid acceleration, wheel spinning, or sudden stops, which result in slipping. Both effects reduce traction at the contact points, causing energy loss and deviation from the intended path. Since the slip and skid vectors act in opposite directions and are difficult to separate, they are measured together as wheel slip and incorporated into the dynamic equations. Furthermore, longitudinal slip and unequal motion of the front wheels induce body rotation, leading to lateral slip and additional path deviation.

Figure 2.

(a) Wheel slipping opposing the WLM’s motion, (b) Wheel skidding in the same direction as the WLM’s motion.

Each robot wheel’s slippage is obtained as follows [Reference Ding, Gao, Deng, Guo and Liu17]:

(1) \begin{equation} \boldsymbol{s}_{\boldsymbol{p}}=\frac{\boldsymbol{r}\boldsymbol{\omega} -\boldsymbol{v}_{\boldsymbol{w}}}{\boldsymbol{r}\boldsymbol{\omega} }\left(\boldsymbol{r}\boldsymbol{\omega} \geq \boldsymbol{v}_{\boldsymbol{w}},0\leq \boldsymbol{s}_{\boldsymbol{p}}\leq 1\right) \end{equation}

Where $\boldsymbol{r}$ is the wheel radius, $\boldsymbol{\omega}$ represents its angular velocity, and $\boldsymbol{v}_{\boldsymbol{w}}$ is the WLM’s wheel forward linear velocity. The skidding rate of each wheel is obtained as follows [Reference Ding, Gao, Deng, Guo and Liu17]:

(2) \begin{equation} \boldsymbol{s}_{\boldsymbol{d}}=\frac{\boldsymbol{v}_{\boldsymbol{w}}-\boldsymbol{r}\boldsymbol{\omega} }{\boldsymbol{v}_{\boldsymbol{w}}}\left(\boldsymbol{r}\boldsymbol{\omega} \leq \boldsymbol{v}_{\boldsymbol{w}},0\leq \boldsymbol{s}_{\boldsymbol{d}}\leq 1\right) \end{equation}

It can be concluded that the direction of the slipping rate and the skidding rate is opposite to each other. Therefore, the combined effects of these two factors are referred to as wheel slip:

(3) \begin{equation} \boldsymbol{s}_{\boldsymbol{p}}=-\boldsymbol{s}_{\boldsymbol{d}}=\frac{\boldsymbol{r}\boldsymbol{\omega} -\boldsymbol{v}_{\boldsymbol{w}}}{\boldsymbol{v}_{\boldsymbol{w}}}\left(\boldsymbol{r}\boldsymbol{\omega} \leq \boldsymbol{v}_{\boldsymbol{w}},-1\leq \boldsymbol{s}_{\boldsymbol{p}}\leq 0\right) \end{equation}

Wheel slip can be obtained from the deviation of the WLM’s actual movement and the WLM’s ideal motion:

(4) \begin{equation} \boldsymbol{s}_{\boldsymbol{l}}=\sum _{i=1}^{n}(| \boldsymbol{r}\boldsymbol{\omega} _{\boldsymbol{i}}-{\boldsymbol{v}_{\boldsymbol{w}}}_{\boldsymbol{i}}| {\unicode[Arial]{x0394}} t) \end{equation}

Where $n$ is the number of measurements and ${\unicode[Arial]{x0394}} t$ is the sampling time. Also, the lateral slip $\boldsymbol{v}_{\boldsymbol{y}}$ is obtained using the slip angle $\hat{\beta }$ from Eq. (5):

(5) \begin{equation} \boldsymbol{v}_{\boldsymbol{y}}=\boldsymbol{v}_{\boldsymbol{x}}\tan(\hat{\beta }) \end{equation}

In the above equation, the slip angle represents the angle between the robot’s velocity $\boldsymbol{v}_{\boldsymbol{x}}$ and its orientation. In measuring lateral slip $\boldsymbol{s}_{\boldsymbol{s}}$ , the amount of deviation from the ideal path at each moment is measured. On a straight path, the lateral velocity $\boldsymbol{v}_{\boldsymbol{d}}$ is equal to zero.

(6) \begin{equation} \boldsymbol{s}_{\boldsymbol{s}}=\sum _{i=1}^{n}(| \boldsymbol{v}_{\boldsymbol{d}}-\boldsymbol{v}_{\boldsymbol{y}}| {\unicode[Arial]{x0394}} t) \end{equation}

3.1. Kinematics of the WLM

By representing the vectors of the links’ position and differentiating them, each respective link’s endpoint velocities, in relation to the WLM’s center of mass, are calculated. This process then allows us to derive the velocities of the WLM.

(7) \begin{equation} \boldsymbol{v}=\boldsymbol{v}_{\boldsymbol{G}}+\left(\dot{\Psi }\hat{\boldsymbol{j}}+\dot{\theta }\hat{\boldsymbol{k}}\right)\times \boldsymbol{r}_{\frac{\boldsymbol{link}}{\boldsymbol{G}}}+\dot{\boldsymbol{r}}_{\frac{\boldsymbol{link}}{\boldsymbol{G}}} \end{equation}

Here, $\boldsymbol{v}_{\boldsymbol{G}}$ is the body’s velocity including the effect of wheel slip terms, while $\boldsymbol{r}_{\frac{\boldsymbol{link}}{\boldsymbol{G}}}$ denotes the link’s endpoint vector to the structure’s center of mass. Additionally, $\dot{\boldsymbol{r}}_{\frac{\boldsymbol{link}}{\boldsymbol{G}}}$ indicates the relative velocity of the link’s endpoint concerning the platform. Due to the base of the manipulator undergoing rotational motion around the Z and Y axes of the body, the velocity’s second term significantly influences the overall velocity dynamics.

Next, the velocity of the body is expressed along the respective coordinate axes as follows [Reference Korayem, Nikseresht and Asadi3]:

(8) \begin{equation} (\boldsymbol{v}_{\boldsymbol{G}})_{x}=\dot{x},(\boldsymbol{v}_{\boldsymbol{G}})_{y}=\dot{y},(\boldsymbol{v}_{\boldsymbol{G}})_{z}=\dot{z}\end{equation}

In the given equation, $\dot{x}$ represents the body’s linear longitudinal velocity, $\dot{y}$ denotes the linear lateral velocity, and $\dot{z}$ corresponds to the linear vertical velocity. Then we write the manipulator’s links position vector.

The manipulator’s first link position vector:

(9) \begin{equation} \boldsymbol{r}_{\boldsymbol{FG}}=(l_{1}\cos (\theta _{1}+\Psi ))\hat{\boldsymbol{i}}+(l_{1}\sin (\theta _{1}+\Psi ))\hat{\boldsymbol{k}} \end{equation}

The manipulator’s second link position vector [Reference Asadi, Nikseresht and Korayem2]:

(10) \begin{equation} \boldsymbol{r}_{\boldsymbol{EG}}=(l_{1}\cos (\theta _{1}+\Psi )+l_{2}\cos (\theta _{1}+\theta _{2}+\Psi ))\hat{i}+(l_{1}\sin (\theta _{1}+\Psi )+l_{2}\sin (\theta _{1}+\theta _{2}+\Psi ))\hat{\boldsymbol{k}} \end{equation}

In the given equations, $l_{1}$ and $l_{2}$ represent the manipulator’s first and second link lengths, in that order. $\theta _{1}$ is the first link angle of rotation, and $\theta _{2}$ is the second link angle of rotation. Also, the position of each manipulator’s link can be determined using position vectors and the platform’s position.

The manipulator’s first link position:

(11) \begin{equation} x_{F}=(x+l_{1}\cos (\theta _{1}+\Psi ))\hat{\boldsymbol{i}}+(y)\hat{\boldsymbol{j}}+(z+l_{1}\sin (\theta _{1}+\Psi ))\hat{\boldsymbol{k}} \end{equation}

The manipulator’s second link position:

(12) \begin{equation}x_{E}=(x+l_{1}\cos (\theta _{1}+\Psi )+l_{2}\cos (\theta _{1}+\theta _{2}+\Psi ))\hat{\boldsymbol{i}}+(y)\hat{\boldsymbol{j}}+(z+l_{1}\sin (\theta _{1}+\Psi )+l_{2}\sin (\theta _{1}+\theta _{2}+\Psi ))\hat{\boldsymbol{k}}\end{equation}

In the given equations, $x$ represents the body’s linear longitudinal position, $x$ denotes the linear lateral position, and $z$ corresponds to the linear vertical position. By deriving from Eqs. (9) and (10), we can define the velocity vectors of each link of the manipulators concerning the platform.

The velocity vector of the manipulator’s first link:

(13) \begin{equation} \dot{\boldsymbol{r}}_{\boldsymbol{FG}}=({-}l_{1}(\dot{\theta }_{1}+\dot{\Psi })\sin (\theta _{1}+\Psi ))\hat{\boldsymbol{i}}+(l_{1}(\dot{\theta }_{1}+\dot{\Psi })\cos (\theta _{1}+\Psi ))\hat{\boldsymbol{k}} \end{equation}

The velocity vector of the manipulator’s second link:

(14) \begin{eqnarray}\dot{\boldsymbol{r}}_{\boldsymbol{EG}}&=&({-}l_{1}(\dot{\theta }_{1}+\dot{\Psi })\sin (\theta _{1}+\Psi )-l_{2}(\dot{\theta }_{1}+\dot{\theta }_{2}+\dot{\Psi })\sin (\theta _{1}+\theta _{2}+\Psi ))\hat{\boldsymbol{i}}+(l_{1}(\dot{\theta }_{1}+\dot{\Psi })\cos (\theta _{1}+\Psi )\nonumber\\&& +\ l_{2}(\dot{\theta }_{1}+\dot{\theta }_{2}+\dot{\Psi })\cos (\theta _{1}+\theta _{2}+\Psi ))\hat{\boldsymbol{k}}\end{eqnarray}

By writing Eq. (7) for each link of the manipulator, the manipulator’s links’ velocities are defined as follows.

The first link’s velocity:

(15) \begin{eqnarray}\boldsymbol{v}_{\boldsymbol{F}} &=& \left( {\dot x + {l_1}\dot \Psi \sin \left( {{\theta _1} + \Psi } \right) - {l_1}({{\dot \theta }_1} + \dot \Psi } \right)\sin \left( {{\theta _1} + \Psi } \right))\hat {\boldsymbol{i}} + \left( {{l_1}\dot \theta \cos \left( {{\theta _1} + \Psi } \right) + \dot y} \right)\hat {\boldsymbol{j}}\nonumber\\[5pt]&& + \left( {\dot z - {l_1}\dot \Psi \dot \theta \cos \left( {{\theta _1} + \Psi } \right) + {l_1}({{\dot \theta }_1} + \dot \Psi )\cos \left( {{\theta _1} + \Psi } \right)} \right)\hat {\boldsymbol{k}}\end{eqnarray}

The second link’s velocity:

(16) \begin{eqnarray}\boldsymbol{v}_{\boldsymbol{E}}&=&(\dot{x}+\dot{\Psi}(l_1 \sin (\theta_1+\Psi)+l_2 \sin (\theta_1+\theta_2+\Psi))-l_1 (\dot{\theta}_1+ \dot{\Psi})\sin (\theta_1+\Psi)-l_2(\dot{\theta}+\dot{\theta}_2 + \dot{\Psi}) \sin\nonumber\\[5pt]&& (\theta_1+\ \theta_2+\Psi))\hat{\boldsymbol{i}}+(\dot{y}+\dot{\theta}(l_1 \cos(\theta_1+\Psi)+l_2 \cos(\theta_1+\theta_2+\Psi))\hat{\boldsymbol{j}}+(\dot{z}-\dot{\Psi}(l_1 \cos(\theta_1+\Psi)\nonumber\\[5pt]&&+\ l_2 \cos(\theta_1+\theta_2+\Psi)+l_1 (\dot{\theta}+\dot{\Psi})\cos(\theta_1+\Psi)+l_2 (\dot{\theta}+\dot{\theta}_2+\dot{\Psi})\cos(\theta_1+\theta_2+\Psi))\hat{\boldsymbol{k}},\end{eqnarray}

Next, we utilize the above relationships to derive the hierarchical dynamic equation of the WLM in slip conditions.

3.2. Dynamics of the WLM

In the following, to derive the hierarchical equations of the WLM’s dynamics, we use the Gibbs–Appell formulation. While the recursive Gibbs–Appell formulation is powerful and has been applied in our prior works [Reference Korayem, Shafei, Doosthoseini, Absalan and Kadkhodaei12, Reference Korayem and Shafei13], in this study, we intentionally used the non-recursive form. This choice was made to simplify the integration of manipulator and base dynamics within a single hierarchical framework. Also, the manipulator’s pseudo-coordinates, $\theta _{1}$ and $\theta _{2}$ , are incorporated into the mechanism’s primary pseudo-coordinates, enabling the system’s accelerations to be expressed in terms of these variables. A Gibbs function specific to the manipulator is added to the robot’s Gibbs functions, resulting in a hierarchical structure for deriving the WLM’s dynamic equations. By applying the Gibbs formulation to the manipulator’s pseudo-coordinates and differentiating them individually, the WLM’s dynamic equations are achieved. The system’s Gibbs functions are calculated independently for each of the two sections of the WLM.

The Gibbs function of the main body [Reference Korayem, Nikseresht and Asadi3]:

(17) \begin{eqnarray}S_{b}&=&\frac{I_{o}}{2}\left(\ddot{\theta }_{rf}^{2}+\ddot{\theta }_{lf}^{2}+\ddot{\theta }_{rr}^{2}+\ddot{\theta }_{lr}^{2}\right)+I_{ly}\left(\ddot{\alpha }^{2}+\ddot{\beta }^{2}\right)+I_{lx}\ddot{\theta }^{2}\left({\kern-1.5pt}\cos ^{2}\alpha +\cos ^{2}\beta \right)+\frac{m_{l}}{2}\left(a_{rf}^{2}+a_{lf}^{2}+a_{rr}^{2}+a_{lr}^{2}\right)\nonumber\\&&+\ \frac{1}{2}\left(I_{Gz}\ddot{\theta }^{2}+I_{Gy}\ddot{\Psi }^{2}+M_{G}\boldsymbol{a}_{\boldsymbol{G}}^{\textbf{2}}\right)\end{eqnarray}

In these equations, $I_{o}$ represents the moment of inertia of the wheel’s axis, while $I_{ly}$ and $I_{lx}$ represent the legs’ moments of inertia with respect to the Y-axis and X-axis, respectively. The parameter $m_{l}$ denotes the mass of the robot’s legs, and $a_{lf}$ , $a_{rf}$ , $a_{lr}$ , and $a_{rr}$ indicate the WLM’s legs accelerations. Additionally, $I_{Gy}$ and $I_{Gz}$ represent the body’s moments of inertia around the Y-axis and Z-axis, respectively, while $M_{G}$ signifies the total mass of the WLM’s body. Finally, $\boldsymbol{a}_{\boldsymbol{G}}$ denotes the WLM’s center of mass acceleration.

The Gibbs function of the links:

(18) \begin{equation} S_{l}=\frac{I_{1}}{2}\left(\ddot{\theta }_{1}+\ddot{\Psi }\right)^{2}+\frac{I_{2}}{2}\left(\ddot{\theta }_{1}+\ddot{\theta }_{2}+\ddot{\Psi }\right)^{2}+\frac{m_{1}}{2}\left(\boldsymbol{a}_{\boldsymbol{F}}\right)^{2}+\frac{m_{2}}{2}\left(\boldsymbol{a}_{\boldsymbol{E}}\right)^{2} \end{equation}

where $I_{1}$ and $I_{2}$ are the first and second manipulator’s links moment of inertia, $m_{1}$ and $m_{2}$ are the first and second manipulator’s links mass and $\boldsymbol{a}_{\boldsymbol{F}}$ and $\boldsymbol{a}_{\boldsymbol{E}}$ are the manipulator’s links accelerations.

The dynamic model of a servo motor is given by:

(19) \begin{equation}\sum\nolimits_{i=1}^{n}\ I_{m}\ \ddot{\theta}_{i}^{2} + \frac{1}{2}\sum\nolimits_{i=1}^{n}\ B_{m}\ \dot{\theta}_{i}^{2} = \tau_{m}\end{equation}

Where, $n$ is the number of motors, $I_{m}$ represents the motor’s moment of inertia, $\ddot{\theta }_{i}$ and $\dot{\theta }_{i}$ denote the angular acceleration and angular velocity of the motors, respectively, $B_{m}$ is the viscous friction coefficient, and $\tau _{m}$ is the torque produced by the motor. The second term of the Eq. (19) is expressed as a Rayleigh dissipation function so that it can be used in the Gibbs–Appell equations [Reference Korayem, Shafei, Doosthoseini, Absalan and Kadkhodaei12]. There are six XL430 motors and four XM540 motors, and Eq. (19) is split into two separate functions accordingly.

The Gibbs function of the motors:

(20) \begin{equation}S_{m}=I_{{m_{1}}}\left(\ddot{\alpha }^{2}+\ddot{\beta }^{2}\right)+\frac{I_{{m_{2}}}}{2}\left(\ddot{\theta }_{rf}^{2}+\ddot{\theta }_{lf}^{2}+\ddot{\theta }_{rr}^{2}+\ddot{\theta }_{lr}^{2}+\ddot{\theta }_{1}^{2}+\ddot{\theta }_{2}^{2}\right)\end{equation}

The viscous damping force of the motors is incorporated into the Rayleigh dissipation function as:

(21) \begin{equation}D_{m}=B_{{m_{1}}}\left(\dot{\alpha }^{2}+\dot{\beta }^{2}\right)+\frac{B_{{m_{2}}}}{2}\left({\dot{\theta }_{rf}}^{2}+{\dot{\theta }_{lf}}^{2}+{\dot{\theta }_{rr}}^{2}+{\dot{\theta }_{lr}}^{2}+{\dot{\theta }_{1}}^{2}+{\dot{\theta }_{2}}^{2}\right)\end{equation}

Then the WLM’s Gibbs function is gained by combining Eq. (17) with Eqs. (18) and (20):

(22) \begin{eqnarray}S_{\textit{total}}&=&\frac{I_{o}}{2}\left(\ddot{\theta }_{rf}^{2}+\ddot{\theta }_{lf}^{2}+\ddot{\theta }_{rr}^{2}+\ddot{\theta }_{lr}^{2}\right)+I_{ly}\left(\ddot{\alpha }^{2}+\ddot{\beta }^{2}\right)+I_{lx}\ddot{\theta }^{2}\left({\kern-1.5pt}\cos ^{2}\alpha +\cos ^{2}\beta \right)+\frac{m_{l}}{2}\left(a_{rf}^{2}+a_{lf}^{2}+a_{rr}^{2}+a_{lr}^{2}\right)\nonumber\\&&+\ \frac{1}{2}\left(I_{Gz}\ddot{\theta }^{2}+I_{Gy}\ddot{\Psi }^{2}+M_{G}\boldsymbol{a}_{\boldsymbol{G}}^{\textbf{2}}\right)+\frac{I_{1}}{2}\left(\ddot{\theta }_{1}+\ddot{\Psi }\right)^{2}+\frac{I_{2}}{2}\left(\ddot{\theta }_{1}+\ddot{\theta }_{2}+\ddot{\Psi }\right)^{2}+\frac{m_{1}}{2}\left(\boldsymbol{a}_{\boldsymbol{F}}\right)^{2}+\frac{m_{2}}{2}\left(\boldsymbol{a}_{\boldsymbol{E}}\right)^{2}\nonumber\\&& +\ I_{{m_{1}}}\left(\ddot{\alpha }^{2}+\ddot{\beta }^{2}\right)+\frac{I_{{m_{2}}}}{2}\left(\ddot{\theta }_{rf}^{2}+\ddot{\theta }_{lf}^{2}+\ddot{\theta }_{rr}^{2}+\ddot{\theta }_{lr}^{2}+\ddot{\theta }_{1}^{2}+\ddot{\theta }_{2}^{2}\right)\end{eqnarray}

Following this, the WLM robot body’s potential energy function is determined using the equation below [Reference Korayem, Nikseresht and Asadi3]:

(23) \begin{equation} V_{b}=gl\left(m_{l}\left({\kern-1.5pt}\sin \left(\alpha \right)+\sin \left(\beta \right)\right)+\frac{M_{G}}{2}\left({\kern-1.5pt}\sin \left(\alpha +\Psi \right)+\sin \left(\beta +\Psi \right)\right)\right) \end{equation}

Here, $g$ denotes gravitational acceleration. The potential function of the manipulator is determined using the following equation.

(24) \begin{eqnarray}V_{m}&=&g\left(m_{1}\left(\frac{1}{2}l_{1}\sin \left(\theta _{1}\right)+\frac{1}{2}l_{x}\sin \left(\Psi \right)+l\sin \left(\alpha \right)\right)\right.\nonumber\\&&+\left. m_{2}\left(l_{1}\sin \left(\theta _{1}\right)+\frac{1}{2}l_{2}\sin \left(\theta _{1}+\theta _{2}\right)+\frac{1}{2}l_{x}\sin \left(\Psi \right)+l\sin \left(\alpha \right)\right)\right)\end{eqnarray}

The function of potential energy of the WLM is then obtained by combining the function of the body and the manipulator’s potential energy. The function is achieved by integrating Eqs. (23) and (24) as shown below:

(25) \begin{eqnarray}V_{total}&=&gl\bigg(m_l ({\kern-1.5pt}\sin(\alpha)+\sin(\beta) )+\frac{M_G}{2}({\kern-1.5pt}\sin(\alpha+\Psi)+\sin(\beta+\Psi))\bigg)\nonumber\\&&+\ g(m_1 \bigg(\frac{1}{2} l_1 \sin(\theta_1 )+\frac{1}{2} l_x \sin(\Psi)+l\sin(\alpha) \bigg)\nonumber\\&&+\ m_2 \bigg(l_1 \sin(\theta_1 )+ \frac{1}{2} l_2 \sin(\theta_1+\theta_2 )+ \frac{1}{2} l_x \sin(\Psi)+l\sin(\alpha) \bigg)\end{eqnarray}

We define the friction function of the wheels considering the Coulomb friction as follows. This function considers the dynamic treatment of the wheels, legs, and manipulator, and it is essential in modeling the system’s motion due to friction [Reference Korayem, Nikseresht and Asadi3]:

(26) \begin{eqnarray}D_w&=& \mu_{dx} (F_{nf} (\dot{x}_{srf}\, sgn(\dot{x}_{srf} )+\dot{x}_{slf}\, sgn(\dot{x}_{slf} ) )+F_{nr} (\dot{x}_{srr}\, sgn(\dot{x}_{srr} )+\dot{x}_{slr}\, sgn(\dot{x}_{slr} ) )\nonumber\\&&+\ 2\mu_{dy} (F_{nf} \dot{y}_f\, sgn(\dot{y}_f )+F_{nr} \dot{y}_r\, sgn(\dot{y}_r ))\end{eqnarray}

In the above equation, $\mu _{dx}$ is the longitudinal coefficient and $\mu _{dy}$ is the lateral coefficient of the WLM’s wheels friction. Additionally, $F_{nf}$ is the front wheels’ normal forces and $F_{nr}$ is the rear wheels’ normal forces. By combining the Eqs. (21) and (26) the total dissipation function of the WLM can be obtained as follows:

(27) \begin{eqnarray}D_{total}&=& \mu_{dx} (F_{nf} (\dot{x}_{srf}\, sgn(\dot{x}_{srf} )+\dot{x}_{slf}\, sgn(\dot{x}_{slf} ) )+F_{nr} (\dot{x}_{srr}\, sgn(\dot{x}_{srr} )+\dot{x}_{slr}\, sgn(\dot{x}_{slr} ) )+2\mu_{dy} (F_{nf} \dot{y}_f\, sgn(\dot{y}_f )\nonumber\\&&+\ F_{nr} \dot{y}_r\, sgn(\dot{y}_r )) +B_(m_1 ) (\dot{\alpha}^2+ \dot{\beta}^2)+\frac{B_{m_2}}{2}\left(\dot{\theta}_{rf}^2+\dot{\theta}_{lf}^2+\dot{\theta}_{rr}^2+\dot{\theta}_{lr}^2+\dot{\theta}_1^2+\dot{\theta}_2^2\right) \end{eqnarray}

The Gibbs equation is stated as a system of equations that incorporates both the dynamic and kinematic behaviors of the mechanism under the influence of slip conditions. This formulation considers the forces and torques generated by wheel slippage, as well as the changes in the acceleration and velocity of the WLM’s components, including the wheels and manipulators, as they interact with the surrounding environment. By factoring in slip, the Gibbs equation offers a clearer depiction of the robot’s movement, enabling a thorough analysis of its performance across various surfaces and terrains [Reference Korayem, Nikseresht and Asadi3]:

(28) \begin{equation} \frac{\partial S_{\textit{total}}}{\partial \ddot{q}}+\frac{\partial D_{\textit{total}}}{\partial \dot{q}}+\frac{\partial V_{\textit{total}}}{\partial q}=\tau \end{equation}

Here, $q$ denotes the system’s pseudo-coordinates, and the generalized torques are represented by $\tau$ . The identification of 14 pseudo-coordinates as independent coordinates has been made, and we focus on the manipulator’s pseudo-coordinates, $\theta _{1}$ and $\theta _{2}$ . The formulation of the remaining 12 pseudo-coordinates can be found in our previous work [Reference Korayem, Nikseresht and Asadi3]. Finally, we derive the torque equations for the manipulator’s links:

Torque of first link’s motor:

(29) \begin{eqnarray}\tau _{\theta 1}&=& I_{1}\left(\ddot{\theta }_{1}+\ddot{\Psi }\right)+I_{2}\left(\ddot{\theta }_{1}+\ddot{\theta }_{2}+\ddot{\Psi }\right)+g\left(\frac{m_{1}l_{1}\cos \left(\theta _{1}\right)}{2}+m_{2}\left(l_{1}\cos \left(\theta _{1}\right)+\frac{l_{2}\cos \left(\theta _{1}+\theta _{2}\right)}{2}\right)\right)\nonumber\\&&+\ I_{{m_{2}}}\ddot{\theta }_{1}+B_{{m_{2}}}\dot{\theta }_{1}\end{eqnarray}

Torque of second link’s motor:

(30) \begin{equation}\tau _{\theta 2}=I_{2}\left(\ddot{\theta }_{1}+\ddot{\theta }_{2}+\ddot{\Psi }\right)+\frac{m_{2}l_{2}g\cos \left(\theta _{1}+\theta _{2}\right)}{2}+I_{{m_{2}}}\ddot{\theta }_{2}+B_{{m_{2}}}\dot{\theta }_{2}\end{equation}

4.1. Slip measurement

To quantify the amount of slippage experienced by the robot, an IMU (MPU-6050) sensor mounted on the robot body was employed [Reference Korayem, Nikseresht and Asadi3]. The sensor provides an orientation accuracy of approximately ±0.5° and an acceleration accuracy of ±0.02 m/s², which were considered when evaluating measurement uncertainty. Slip was estimated by comparing the actual displacement of the robot, obtained via IMU-based integration, with the ideal planned trajectory. Prior to experimentation, the MPU-6050 was calibrated to minimize bias and enhance measurement precision. Gyroscope calibration was performed by keeping the sensor stationary on a flat surface, recording the mean offsets of angular velocity along all three axes, and compensating for these offsets in subsequent measurements to reduce drift. Accelerometer calibration involved positioning the sensor in six known orientations, aligning each axis with gravity in both positive and negative directions, comparing the measured accelerations to the expected ±1 g, and calculating scale factors and biases for correction. This procedure reduced systematic errors and ensured reliable measurements of angular velocity and acceleration during testing. Despite calibration, uncertainty remained in slip estimation due to accelerometer noise (±30–50 mg), gyroscope noise (∼0.005 °/s/ $\surd$ Hz), residual bias drift, sensor misalignment, and limited sampling rate. Variations in surface conditions (dry vs. wet ceramic) further increased measurement error. Overall, slip measurement uncertainty was estimated to be 5–10% of the computed slip ratio. Friction coefficients were derived experimentally from tests conducted on dry and wet (soapy) ceramic surfaces. By integrating IMU-derived acceleration and orientation data, the actual displacement of the robot body was estimated for each time step, and slippage was computed by comparing this displacement with the ideal trajectory. The longitudinal friction coefficient ( $\mu _{dx}$ ) between the robot’s rubber wheels and the ceramic tile was measured as 1.2 on dry, dusty ceramic and 0.65 on wet, soapy ceramic, based on experimental observations of wheel speed and contact conditions [Reference Nikolić, Banić, Stamenković, Simonović, Miltenović and Pavlović31]. The ratio of longitudinal to lateral friction coefficients ( $\mu _{dx}$ / $\mu _{dy}$ ) was assumed to be approximately 1.1 [Reference Xia, Liu and Lu32], resulting in lateral coefficients ( $\mu _{dy}$ ) of 1.09 and 0.59 for dry and wet surfaces, respectively. Wheel slip and skid were modeled using a combination of Coulomb and viscous friction laws, with longitudinal and lateral slip coefficients defined from empirical data and assumed constant in simulations, consistent with prior studies [Reference Nikolić, Banić, Stamenković, Simonović, Miltenović and Pavlović31, Reference Xia, Liu and Lu32]. The Dynamixel XM540-W150-R and XL430-W250-T actuators introduce additional uncertainty due to mechanical tolerances, sensor inaccuracies, communication latency, and environmental factors. Position accuracy is typically within ±0.1° –0.2°, while torque output may fluctuate by ±5%. Communication delay is approximately 1–2 ms. These errors can be mitigated through periodic calibration, maintaining a stable power supply, optimizing communication parameters, and controlling operating conditions such as temperature and electromagnetic interference, thereby improving overall motor reliability during robotic experiments.

4.2. WLM’s first test

The robot’s first test on dry, dusty ceramic demonstrates that the WLM accelerates, rotating its rear legs from 30° to 150°, while the front legs remain fixed at 30° (Figure 6). The WLM’s links begin moving from a 0° angle at the front of the platform without acceleration. The front wheels rotate with acceleration, whereas the rear wheels remain stationary. Overall, the legs exhibit uniform motion, with no additional acceleration.

Figure 6.

Robot and manipulator movement in the first test.

Figure 7.

WLM’s end-effector position during the first test of the robot.

The robot body’s displacement strongly influences the WLM’s end-effector trajectory during the first test, reflecting both the body’s motion and the arm’s own movement (Figure 7). The results also indicate that wheel slip has minimal effect on changes in the end-effector’s height along the z-axis.

Wheel slip causes minor differences between the experimental and simulation speed profiles of the WLM’s end-effector during the first test, although their overall behavior remains similar (Figure 8).

Figure 8.

Velocity of end-effector during the first test.

The manipulator’s end-effector movement is coupled with the robot body’s motion during the first test, resulting in a combined trajectory as the robot moves through space (Figure 9).

Figure 9.

3D position of WLM in the first test.

Figures 7–9 show the manipulator’s end-effector position and velocity during the first test. While the overall trends in simulation and experiment agree, deviations of 8–12% were observed due to wheel slip. Slip had minimal effect on the end-effector’s vertical motion but significantly influenced its horizontal trajectory. Figure 10 presents the arm joint torques. The first joint requires more torque than the second, as it must support the torque of the second link.

Figure 10.

Torque of WLM joints in the first test.

4.3. WLM’s second test

During the second test on wet, soapy ceramic with uneven terrain, the rear legs rotate from 150° to 30° while the front legs move from 30° to 150°. The WLM’s links at the front start from 0° without acceleration, with the front wheels rotating under acceleration and the rear wheels remaining stationary. The legs move at a constant velocity, and the WLM follows a circular path, with the left front wheel speed approaching zero (Figure 11).

Figure 11.

Robot and manipulator circular movement on wet and soapy ceramic in the second test.

Wheel slip and skid during accelerated motion cause some deviation along the x-axis in the WLM’s end-effector position during the second test, although the overall behavior remains similar (Figure 12).

Figure 12.

WLM’s end-effector position during the second test of the robot.

Figure 13 shows the WLM’s end-effector velocity during the second test. The robot body’s speed and acceleration strongly influence the arm’s velocity profile, causing the end-effector to follow both its own motion and the body’s movement.

Figure 13.

Velocity of end-effector in the second test.

Figure 14 illustrates the WLM’s three-dimensional position during the second test. It shows the end-effector’s motion relative to the robot body, demonstrating that the end-effector undergoes a combined movement as the robot moves through space.

Figure 14.

3D Circular position of WLM in the second test.

In this test (Figures 12–14), the WLM followed a circular path. End-effector deviations of 10–14% between simulation and experiment were observed, primarily due to increased slip and skid during accelerated motion. Experimental velocity profiles confirmed the strong coupling between the manipulator and the platform’s body dynamics. Figure 15 shows the arm joint torques during the second test. The torques are similar to those in the first test, despite different body motion, indicating that the manipulator’s dynamics are largely independent of the robot body’s movement.

4.4. Systematic error analysis

To further validate the Gibbs–Appell–based hierarchical model, a quantitative comparison was performed among the theoretical predictions, Webots simulation results, and experimental measurements. The evaluation focused on the end-effector position, velocity, and manipulator joint torque under both dry and wet (soapy ceramic) surface conditions. Error metrics were computed as mean percentage differences between the root mean square (RMS) magnitudes of corresponding data over one complete motion cycle of the robot.

A comprehensive analysis was conducted to identify and quantify the sources of deviation between theoretical models, simulations, and experimental results of the robot in Table III. Overall, the theoretical and Webots simulation results showed close agreement, while larger deviations appeared in physical experiments – particularly on wet or low-traction surfaces. Torque analysis revealed that the first joint consistently required higher torque than the second due to load transmission, with initial deviations of approximately 6.5% and 20%, respectively, converging after 8–14 s. End-effector position errors in experiments ranged between 8 and 12%, primarily attributed to wheel slip, frictional variability, and sensor uncertainty. Slip ratios of 5–10% were observed on dry ( $\mu _{dx}$ ≈ 1.2) and wet or soapy ( $\mu _{dx}$ ≈ 0.65) surfaces, leading to trajectory deviations of about 8–10.5%. Additional discrepancies stemmed from actuator nonidealities – such as torque ripple, backlash, saturation, and communication delay – as well as unmodeled compliance and elastic deformation in legs and wheels, which introduced phase shifts in torque and velocity. Sensor limitations, including IMU drift (e.g., MPU-6050), ±0.02 m/s² noise, and encoder quantization, contributed roughly 5–10% uncertainty in position estimation. Model simplifications, including assumptions of rigid links, constant friction coefficients, linear damping, and neglect of tire deformation and thermal effects, accounted for residual errors of 3–5%. Collectively, these factors explain most of the observed discrepancies between simulation and experiment. Future improvements should focus on closed-loop compensation, online slip and friction estimation, compliance modeling, and incorporation of additional feedback sources such as wheel encoders and visual odometry to enhance model fidelity and robustness.

4.5. Gibbs–Appell vs. Lagrangian formulations

To quantitatively validate the performance of the proposed Gibbs–Appell formulation, a comparative study was conducted against an equivalent Lagrangian-based model using the same setup and data from the two experimental tests presented in this paper. Test 1 involved linear motion on a dry ceramic surface, while Test 2 examined a circular trajectory on a wet ceramic surface coated with a water and soap mixture. In both tests, the models were evaluated using identical input trajectories, link parameters, and friction coefficients as reported in the manuscript. The comparison metrics included the root-mean-square (RMS) error between theoretical predictions, simulation results, and experimental measurements, as well as the symbolic assembly and numerical execution times. For the dry surface, the friction coefficients were set as $\mu _{dx}$ = 1.20 and $\mu _{dy}$ = 1.09, while for the wet surface, they were defined as $\mu _{dx}$ = 0.65 and $\mu _{dy}$ = 0.59. To evaluate the sensitivity of both formulations, a ±20% variation was applied to $\mu _{dx}$ , with $\mu _{dx}$ / $\mu _{dy}=1.1$ .

Table III.

Comparison of average absolute errors.

Figure 15.

Torque of WLM joints in the second test.

Table IV.

Quantitative comparison between the Gibbs–Appell and Lagrangian formulations.

Table IV summarizes the quantitative comparison between the Gibbs–Appell and Lagrangian formulations under varying friction coefficients for both physical tests of the WLM. The Gibbs–Appell formulation consistently achieves lower RMS errors in simulations and experiments. For Test 1 (dry), end-effector RMS errors are 3–4% (simulation) and 8–10% (experimental), while for Test 2 (wet) they are 4–5% and 10–12%, respectively. The symbolic assembly time for Gibbs–Appell ranges from 3.48 to 4.58 s and execution time from 0.39 to 0.49 s, both markedly lower than the Lagrangian model (6.8–8.0 s and 0.55–0.66 s). When μ_dx varies by ±20%, Gibbs–Appell accuracy changes less than ±5%, while the Lagrangian model deviates up to ±15%, demonstrating the superior robustness of Gibbs–Appell dynamics to friction uncertainty. Comparative and sensitivity analyses confirm its higher accuracy and computational efficiency. Symbolic assembly is nearly twice as fast, and prediction accuracy improves, validating earlier qualitative findings. Symbolic differentiation and numerical evaluation were executed in MATLAB software; timing was measured using MATLAB’s tic–toc.

The hierarchical Gibbs–Appell framework provides accurate dynamic prediction and supports various control and planning tasks. By decoupling platform and manipulator dynamics, it enables efficient inverse dynamics computation and real-time feedforward control. Model-based slip estimation enhances adaptive gain tuning and online compensation, improving robustness under changing contact conditions. Its computational efficiency supports real-time trajectory optimization with slip and contact constraints, enabling rapid re-planning in uncertain environments. The hierarchical structure also facilitates residual generation for fault diagnosis, including slip detection, actuator degradation, and sensor faults.