Integration of Trailer and Tractor Mobile Robot Considering Towing Dynamics

Ünal Dana; Levent Çetin

doi:10.14313/jamris-2026-003

Full Article

1.

Introduction

The advancements in warehouse robotics over the past decade have been significant [1], driving increased attention towards controlling and evaluating various mobile robot designs for warehouse environments [2, 3]. This research focuses primarily on key areas such as energy efficiency, path planning [4, 5], path tracking [6, 7], dynamic behavior analysis [4,8–10] and mechanical design considerations [2,11]. These topics form the foundation for mobile robotics, specifically for autonomous towing of trailers exploiting in transferring products between warehouse and manufacturing facilities or vice versa.

Previous studies in mobile robotics have examined trailer systems [12–14], which share mechanical similarities with non-holonomic two-wheel mobile robots. A significant portion of the literature focuses on trailer dynamics and control for path tracking, often employing Lyapunov-based approaches to address tracking control for non-holonomic wheeled trailer’s [15, 16]. While these methods offer useful insights, they elaborately calculate required torques for actuated wheels by implementing differential drive models concerning only trailers and combined dynamics of tractor and trailer. Nevertheless, in context of mobilization and autonomy of already deployed trailers, these strategies become invalid because of the fact that only towing of unactuated trailers is feasible in their controls. To tow for motion control of the trailer with non-holonomic dynamics, a tractor robot should have fully controlled actuation for planar motion.

Omniwheel robots, known for their maneuverability and holonomic kinematics, have been a significant focus of research. However, their complex kinematic and dynamic analyses present additional challenges [8, 11, 17]. Control strategies for omniwheel robots have been successfully tested, employing a PID control method to adjust wheel accelerations [17], implementing adaptive controllers [9]. This approach calculates required motor torques using inverse dynamic analysis and adjusts the robot’s velocity based on position errors. The desired velocity and planned path are provided as input by the user, with the system controlling the robot accordingly. Although this method is effective, the article does not specify the exact path tracking or planning algorithms tested alongside the control system.

Path planning is a critical aspect that enables autonomous movement in mobile robots and requires an accompanying path tracking algorithm [4, 7]. During path planning, a mobile robot typically relies on map data, which may be preexisting or acquired in real-time using sensors. These data are essential for computing an optimal route from the initial position to the destination through path planning algorithms. The resulting path can be represented as a set of discrete poses or as a mathematical function, both of which require a tracking algorithm to follow the planned route [5, 6]. Path tracking algorithms, such as the Pure Pursuit method, determine the next destination based on a lookahead distance. This approach continuously adjusts the robot’s trajectory by selecting the closest point along the path that lies within the lookahead distance, ensuring smooth and efficient navigation.

This study addresses the challenge of towing a trailer using an omniwheel mobile robot, offering a unique approach to this common problem in mobile robotics.

While previous research has examined the dynamics of omniwheel robots [8, 11, 17–19] and trailer’s [4,7,10,13,19–22], this study re-examinesthem dueto changes in the problem’s requirements. It derives the dynamics of towing and the forces involved, providing new insights into their interaction.

This study introduces a method using cubic Bézier functions with pure pursuit tracking. This approach reduces the need for excessive maneuvers and corrects rotational errors, making it suitable for all trailer operations. The control strategy employs dynamic analysis with Lyapunov stability methods. Though commonly used for non-holonomic robots, this method is adapted here for towing applications, with experiments conducted to assess its effectiveness. This research contributes significantly by enhancing warehouse robotics efficiency in towing operations. Additionally, the findings benefit any system utilizing towed trailer’s, offering valuable insights and methodologies for improved performance.

2.

Dynamical Analysis

The proposed design of the towed trailer mobile robot consists of two parts: the non-holonomic trailer part and the holonomic omniwheel tractor part. The two parts are connected to each other using a revolute joint, this joint allows the rotations on Z direction and transmits only the applied forces without deflection but not the torque. To achieve the desired pose of the trailer, the required towing forces on the revolute joint must be calculated. Then, the corresponding motor torques or motor speeds of the towing robot must be determined based on the calculated external forces of the trailer.

The following assumptions and problem-based constraints were considered throughout the whole research.

- The trailer has 3 degrees of freedom.
- Slip between the wheels and the ground is neglected.
- The center of mass of the towing robot is at the geometrical center of its body, equidistant from all three wheels, as shown in Figure 1a.
- The center of mass of the trailer is located between the middle of the two wheels, as shown in Figure 1b.

2.1.

Towing Dynamics

To model the dynamics of the trailer’s motion the geometrical and mechanical relationship between the trailer’s wheels and body frame should be known. All of the necessary parameters are presented in Figure 2.

The trailer’s linear velocity V_t is determined as the average of the left and right wheel velocities V_L and V_R, while the angular velocity ${\dot{ϕ}}_{t}$ {\dot \phi _t} depends on the difference between these velocities and the robot’s wheelbase L.

These relationships are essential for determining how the robot’s wheels must move to achieve a desired trajectory. Specifically, the relationship between the linear and angular velocities is given by Equation (1).

1

\begin{matrix} V_{t} = \frac{V_{L} + V_{R}}{2}, & {\dot{ϕ}}_{t} = \frac{V_{R} - V_{L}}{L} \end{matrix}

\,\,\matrix{ {{V_t} = {{{V_L} + {V_R}} \over 2},\,} & {{{\dot \phi }_t} = {{{V_R} - {V_L}} \over L}} \cr }

Furthermore, the reverse calculation allows to determine the required wheel velocities for a given linear and angular velocity. By rearranging the equations in Equation (1) the following expressions for the left and right wheel velocities are obtained as shown in Equation (2).

2

\begin{matrix} V_{L} = \frac{2 V_{t} - {\dot{ϕ}}_{t} L}{2}, & V_{R} = \frac{2 V_{t} + {\dot{ϕ}}_{t} L}{2} \end{matrix}

\,\matrix{ {{V_L} = {{2{V_t} - {{\dot \phi }_t}L} \over 2},\,} & {{V_R} = {{2{V_t} + {{\dot \phi }_t}L} \over 2}} \cr }

The velocities of the left and right wheels are related to the angular velocities ${\dot{θ}}_{L}$ {{\dot \theta }_L} and ${\dot{θ}}_{R}$ {{\dot \theta }_R} through the wheel radius r_t, as shown in Equation (3).

3

\begin{matrix} V_{L} = {\dot{θ}}_{L} r_{t}, & V_{R} = {\dot{θ}}_{R} r_{t} \end{matrix}

\matrix{ {{V_L} = {{\dot \theta }_L}{r_t},\,} & {{V_R} = {{\dot \theta }_R}{r_t}} \cr }

To express the angular velocities in terms of the trailer’s linear velocity V_t and angular velocity ${\dot{ϕ}}_{t}$ {\dot \phi _t}, the matrix equation in Equation (4) is derived by using Equation (2) and Equation (3).

4

[\begin{matrix} {\dot{θ}}_{L} \\ {\dot{θ}}_{R} \end{matrix}] = \frac{1}{2 r_{t}} [\begin{matrix} 2 & - L \\ 2 & L \end{matrix}] [\begin{matrix} V_{t} \\ {\dot{ϕ}}_{t} \end{matrix}]

\left[ {\matrix{ {{{\dot \theta }_L}} \cr {{{\dot \theta }_R}} \cr } } \right] = {1 \over {2{r_t}}}\left[ {\matrix{ 2 & { - L} \cr 2 & L \cr } } \right]\left[ {\matrix{ {{V_t}} \cr {{{\dot \phi }_t}} \cr } } \right]

This equation allows to determine the angular velocities of the left and right wheels based on the trailer’s linear and angular velocities. The relationships for the angular accelerations of the wheels in terms of the linear and angular accelerations of the trailer is expressed in Equation (5) by taking the derivative of Equation (4).

5

[\begin{matrix} {\ddot{θ}}_{L} \\ {\ddot{θ}}_{R} \end{matrix}] = \frac{1}{2 r_{t}} [\begin{matrix} 2 & - L \\ 2 & L \end{matrix}] [\begin{matrix} {\dot{V}}_{t} \\ {\ddot{ϕ}}_{t} \end{matrix}]

\left[ {\matrix{ {{{\ddot \theta }_L}} \cr {{{\ddot \theta }_R}} \cr } } \right] = {1 \over {2{r_t}}}\left[ {\matrix{ 2 & { - L} \cr 2 & L \cr } } \right]\left[ {\matrix{ {{{\dot V}_t}} \cr {{{\ddot \phi }_t}} \cr } } \right]

These relationships helps to understand how changes in the trailer’s velocity affect the angular accelerations of the wheels and by taking inverse of this equation we can understand the inverse relationship shown on Equation (6).

6

[\begin{array}{l} {\dot{V}}_{t} \\ {\ddot{ϕ}}_{t} \end{array}] = [\begin{matrix} \frac{r_{t}}{2} & \frac{r_{t}}{2} \\ - \frac{1}{L} & \frac{1}{L} \end{matrix}] [\begin{array}{l} {\ddot{θ}}_{L} \\ {\ddot{θ}}_{R} \end{array}]

\left[ {\matrix{ {{{\dot V}_t}} \hfill \cr {{{\ddot \phi }_t}} \hfill \cr } } \right] = \left[ {\matrix{ {{{{r_t}} \over 2}} & {{{{r_t}} \over 2}} \cr { - {1 \over L}} & {{1 \over L}} \cr } } \right]\left[ {\matrix{ {{{\ddot \theta }_L}} \hfill \cr {{{\ddot \theta }_R}} \hfill \cr } } \right]

The forces exerted on the wheels are related to the angular accelerations of the wheels through their moments of inertia. Specifically, the forces F_L and F_R are related to the angular accelerations ${\ddot{θ}}_{L}$ {\ddot \theta _L} and ${\ddot{θ}}_{R}$ {\ddot \theta _R} by the system of equations shown in Equation (7).

7

\begin{array}{l} F_{R} = \frac{I_{R} {\ddot{θ}}_{R}}{r_{t}}, & F_{L} = \frac{I_{L} {\ddot{θ}}_{L}}{r_{t}} \end{array}

\matrix{ {{F_R} = {{{I_R}{{\ddot \theta }_R}} \over {{r_t}}},} \hfill & {{F_L} = {{{I_L}{{\ddot \theta }_L}} \over {{r_t}}}} \hfill \cr }

Equation (8) represents the balance of forces and torques acting on the trailer, where ${\vec{F}}_{e} \times {\vec{L}}_{e}$ {\vec F_e} \times {\vec L_e} represents the external torque, and the terms involving ${\vec{F}}_{L}$ {\vec F_L} and ${\vec{F}}_{R}$ {\vec F_R} describe the torque contributions from the wheels. This balance is fundamental for controlling the robot’s rotational motion.

8

{\vec{F}}_{e} \times {\vec{L}}_{e} + {\vec{F}}_{L} \times {\vec{L}}_{t} + {\vec{F}}_{R} \times {\vec{L}}_{t} = I_{t} {\ddot{ϕ}}_{t}

{\vec F_e} \times {\vec L_e} + {\vec F_L} \times {\vec L_t} + {\vec F_R} \times {\vec L_t} = {I_t}{\ddot \phi _t}

Analysis of the forces along the y-axis, which are influenced by the wheel forces and the trailer’s acceleration is given by Equation (9).

9

F_{e y} - F_{L} - F_{R} = m_{t} {\dot{V}}_{t}

{F_{ey}} - {F_L} - {F_R} = {m_t}{\dot V_t}

This equation shows that the net force along the y-axis is the sum of the forces from the left and right wheels, minus the external force F_ey, which is proportional to the robot’s change in velocity. The force F_ex acting along the x-axis can be derived from the moments of inertia of the robot and the geometry of the system. It is expressed in Equation (10).

10

F_{e x} = \frac{F_{L} L}{2 d_{t}} - \frac{F_{R} L}{2 d_{t}} - \frac{I_{t} {\ddot{ϕ}}_{t}}{d_{t}}

{F_{ex}} = {{{F_L}L} \over {2{d_t}}} - {{{F_R}L} \over {2{d_t}}} - {{{I_t}{{\ddot \phi }_t}} \over {{d_t}}}

Further, the force F_ex can also be described in terms of the angular accelerations of the left and right wheels, as shown in Equation (11).

11

F_{e x} = \frac{I_{L} {\ddot{θ}}_{L} L}{2 d_{t} r_{t}} - \frac{I_{R} {\ddot{θ}}_{R} L}{2 d_{t} r_{t}} + \frac{I_{t} {\ddot{θ}}_{L} r_{t} - I_{t} {\ddot{θ}}_{R} r_{t}}{d_{t} L}

{F_{ex}} = {{{I_L}{{\ddot \theta }_L}L} \over {2{d_t}{r_t}}} - {{{I_R}{{\ddot \theta }_R}L} \over {2{d_t}{r_t}}} + {{{I_t}{{\ddot \theta }_L}{r_t} - {I_t}{{\ddot \theta }_R}{r_t}} \over {{d_t}L}}

This shows how the wheel inertia and angular accelerations influence the external force acting on the trailer. The forces along the y-axis is as shown in Equation (12).

12

F_{e y} - F_{L} - F_{R} = m_{t} {\dot{V}}_{t}

{F_{ey}} - {F_L} - {F_R} = {m_t}{\dot V_t}

This equation shows that the net force along the y-axis depends on the forces from both wheels and the external force. The total force along the y-axis can also be expressed as in Equation (13).

13

F_{e y} = F_{L} + F_{R} + m_{t} {\dot{V}}_{t}

{F_{ey}} = {F_L} + {F_R} + {m_t}{\dot V_t}

This shows that the total force along the y-axis is the sum of the forces from the left and right wheels, plus the contribution from the trailer’s change in velocity. The force F_ey can be expressed in terms ofthe angular accelerations of the left and right wheels, as shown in Equation (14).

14

F_{e y} = \frac{I_{L} {\ddot{θ}}_{L}}{r_{t}} + \frac{I_{R} {\ddot{θ}}_{R}}{r_{t}} + m_{t} \frac{{\dot{V}}_{R} + {\dot{V}}_{L}}{2}

{F_{ey}} = {{{I_L}{{\ddot \theta }_L}} \over {{r_t}}} + {{{I_R}{{\ddot \theta }_R}} \over {{r_t}}} + {m_t}{{{{\dot V}_R} + {{\dot V}_L}} \over 2}

The force F_ey can be written as in Equation(15).

15

F_{e y} = \frac{I_{L} {\ddot{θ}}_{L}}{r_{t}} + m_{t} \frac{{\ddot{θ}}_{L} r_{t}}{2} + \frac{I_{R} {\ddot{θ}}_{R}}{r_{t}} + m_{t} \frac{{\ddot{θ}}_{R} r_{t}}{2}

{F_{ey}} = {{{I_L}{{\ddot \theta }_L}} \over {{r_t}}} + {m_t}{{{{\ddot \theta }_L}{r_t}} \over 2} + {{{I_R}{{\ddot \theta }_R}} \over {{r_t}}} + {m_t}{{{{\ddot \theta }_R}{r_t}} \over 2}

This equation gives a more detailed relationship between the forces and the angular accelerations of the left and right wheels. The final equation of model is a matrix representation of the external forces along the x- and y-axes in terms of the angular accelerations of the left and right wheels. This matrix form is useful for numerical computation and simulation of the trailer’s dynamics. The matrix equation is given in Equation (16).

16

[\begin{array}{l} F_{e x} \\ F_{e y} \end{array}] = [\begin{matrix} \frac{I_{L} L}{2 d_{t} r_{t}} + \frac{I_{t} r_{t}}{d_{t} L} & - \frac{I_{R} L}{2 d_{t} r_{t}} - \frac{I_{t} r_{t}}{d_{t} L} \\ \frac{I_{L}}{r_{t}} + m_{t} \frac{r_{t}}{2} & \frac{I_{R}}{r_{t}} + m_{t} \frac{r_{t}}{2} \end{matrix}] [\begin{array}{l} {\ddot{θ}}_{L} \\ {\ddot{θ}}_{R} \end{array}]

\left[ {\matrix{ {{F_{ex}}} \hfill \cr {{F_{ey}}} \hfill \cr } } \right] = \left[ {\matrix{ {{{{I_L}L} \over {2{d_t}{r_t}}} + {{{I_t}{r_t}} \over {{d_t}L}}} & { - {{{I_R}L} \over {2{d_t}{r_t}}} - {{{I_t}{r_t}} \over {{d_t}L}}} \cr {{{{I_L}} \over {{r_t}}} + {m_t}{{{r_t}} \over 2}} & {{{{I_R}} \over {{r_t}}} + {m_t}{{{r_t}} \over 2}} \cr } } \right]\left[ {\matrix{ {{{\ddot \theta }_L}} \hfill \cr {{{\ddot \theta }_R}} \hfill \cr } } \right]

2.2.

Omniwheel Robot Dynamics

The designed omniwheel mobile robot, used for towing the trailer, calculates the external forces affecting the trailer to determine the needed equivalent motor angular velocities. In this context, the geometrical relationships between the motor vectors resulting from the mechanical connections to the mobile robot body play a crucial role. These relationships are fundamental to calculating the necessary velocities. All required geometrical and mechanical parameters are specified in the conceptual design in Figure 3.

The velocity vector of the omniwheel robot’s center of mass, denoted by ṗ_o, is given by the following equation.

17

\begin{array}{l} {\dot{p}}_{o} = {[\begin{array}{l} V_{x o} & V_{y o} & {\dot{ϕ}}_{o} \end{array}]}^{T}, & {\dot{q}}_{o} = {[\begin{array}{l} {\dot{θ}}_{1} & {\dot{θ}}_{2} & {\dot{θ}}_{3} \end{array}]}^{T} \end{array}

\matrix{ {{{\dot p}_o} = {{\left[ {\matrix{ {{V_{xo}}} \hfill & {{V_{yo}}} \hfill & {{{\dot \phi }_o}} \hfill \cr } } \right]}^T},} \hfill & {{{\dot q}_o} = {{\left[ {\matrix{ {{{\dot \theta }_1}} \hfill & {{{\dot \theta }_2}} \hfill & {{{\dot \theta }_3}} \hfill \cr } } \right]}^T}} \hfill \cr }

Equation (17) defines the velocity vector ṗ_o ofthe omniwheel robot’s center of mass, where V_xo and V_yo represent the linear velocities along the X_o and Y_o axes, and ${\dot{ϕ}}_{o}$ {\dot \phi _o} is the angular velocity of the omniwheel robot.

To calculate the linear velocity along the X_o direction V_xo, we use the following equation.

18

V_{x o} = \frac{2 V_{2} - V_{1} - V_{3}}{3}

{V_{xo}} = {{2{V_2} - {V_1} - {V_3}} \over 3}

Equation (18) combines the contributions of the individual motor velocities V₁, V₂, and V₃, where the motors are arranged in a triangular configuration. The equation shows how the omniwheel robot’s X_o direction velocity is influenced by these motor velocities.

Similarly, the linear velocity along the Y_o direction V_yo is computed as:

19

V_{y o} = \frac{\sqrt{3} V_{3} - \sqrt{3} V_{1}}{3}

{V_{yo}} = {{\sqrt 3 {V_3} - \sqrt 3 {V_1}} \over 3}

Equation (19) shows that the Y_o direction velocity V_yo depends on the motor velocities V₁ and V₃, with a factor of $\sqrt{3}$ \sqrt 3 from the geometry of the robot’s wheels.

20

{\dot{ϕ}}_{o} = \frac{V_{1} + V_{2} + V_{3}}{3 d_{o}}

{\dot \phi _o} = {{{V_1} + {V_2} + {V_3}} \over {3{d_o}}}

Equation (20) provides the angular velocity at the robot’s center of mass based on the sum of the motor velocities and the distance d_o from the center of mass to the wheels.

The individual motor velocities V₁, V₂, and V₃ are related to the motor angular velocities ${\dot{θ}}_{1}, {\dot{θ}}_{2},$ {\dot \theta _1},{\dot \theta _2}, and ${\dot{θ}}_{3}$ {\dot \theta _3} by the following relations.

21

\begin{array}{l} V_{1} = {\dot{θ}}_{1} r_{o}, & V_{2} = {\dot{θ}}_{2} r_{o}, & V_{3} = {\dot{θ}}_{3} r_{o} \end{array}

\matrix{ {{V_1} = {{\dot \theta }_1}{r_o},} \hfill & {{V_2} = {{\dot \theta }_2}{r_o},} \hfill & {{V_3} = {{\dot \theta }_3}{r_o}} \hfill \cr }

In equation (21), r_o represents the radius of the omniwheels, and each omniwheel’s velocity is proportional to its angular velocity.

Rearranging the system of equations for linear and angular velocities to write the matrix equation that transforms the motor angular velocities into the linear velocities of the omniwheel robot’s center of mass.

22

[\begin{matrix} V_{x o} \\ V_{y o} \\ {\dot{ϕ}}_{o} \end{matrix}] = r_{o} [\begin{matrix} - \frac{1}{3} & \frac{2}{3} & - \frac{1}{3} \\ - \frac{\sqrt{3}}{3} & 0 & \frac{\sqrt{3}}{3} \\ \frac{1}{3 d_{o}} & \frac{1}{3 d_{o}} & \frac{1}{3 d_{o}} \end{matrix}] [\begin{matrix} {\dot{θ}}_{1} \\ {\dot{θ}}_{2} \\ {\dot{θ}}_{3} \end{matrix}]

\left[ {\matrix{ {{V_{xo}}} \cr {{V_{yo}}} \cr {{{\dot \phi }_o}} \cr } } \right] = {r_o}\left[ {\matrix{ { - {1 \over 3}} & {{2 \over 3}} & { - {1 \over 3}} \cr { - {{\sqrt 3 } \over 3}} & 0 & {{{\sqrt 3 } \over 3}} \cr {{1 \over {3{d_o}}}} & {{1 \over {3{d_o}}}} & {{1 \over {3{d_o}}}} \cr } } \right]\left[ {\matrix{ {{{\dot \theta }_1}} \cr {{{\dot \theta }_2}} \cr {{{\dot \theta }_3}} \cr } } \right]

Equation (22) expresses the relationship between the motor angular velocities and the linear velocities of the robot’s center of mass in matrix form. This matrix allows for the direct calculation of the robot’s velocities from the motor angular velocities.

The inverse of this relationship is found to calculate the motor angular velocities ${\dot{θ}}_{1}, {\dot{θ}}_{2},$ {\dot \theta _1},{\dot \theta _2}, and ${\dot{θ}}_{3}$ {\dot \theta _3} from the linear velocities V_xo, V_yo, and ${\dot{ϕ}}_{o}$ {\dot \phi _o}.

23

[\begin{array}{l} {\dot{θ}}_{1} \\ {\dot{θ}}_{2} \\ {\dot{θ}}_{3} \end{array}] = \frac{1}{r_{o}} [\begin{matrix} - \frac{1}{2} & - \frac{\sqrt{3}}{2} & d_{o} \\ 1 & 0 & d_{o} \\ - \frac{1}{2} & \frac{\sqrt{3}}{2} & d_{o} \end{matrix}] [\begin{array}{l} V_{x o} \\ V_{y o} \\ {\dot{ϕ}}_{o} \end{array}]

\left[ {\matrix{ {{{\dot \theta }_1}} \hfill \cr {{{\dot \theta }_2}} \hfill \cr {{{\dot \theta }_3}} \hfill \cr } } \right] = {1 \over {{r_o}}}\left[ {\matrix{ { - {1 \over 2}} & { - {{\sqrt 3 } \over 2}} & {{d_o}} \cr 1 & 0 & {{d_o}} \cr { - {1 \over 2}} & {{{\sqrt 3 } \over 2}} & {{d_o}} \cr } } \right]\left[ {\matrix{ {{V_{xo}}} \hfill \cr {{V_{yo}}} \hfill \cr {{{\dot \phi }_o}} \hfill \cr } } \right]

Equation (23) provides the inverse matrix that calculates the motor angular velocities from the omniwheel robot’s linear velocities and angular velocity.

The total torque T_zo acting on the omniwheel robot due to the forces exerted by the motors can be expressed as: 24 $T_{z o} = F_{1} d_{o} + F_{2} d_{o} + F_{3} d_{o} + I_{o} {\dot{w}}_{o}$ {T_{zo}} = {F_1}{d_o} + {F_2}{d_o} + {F_3}{d_o} + {I_o}{\dot w_o}

Equation (24) describes the total torque acting on the robot. It includes the contributions from the individual motor forces F₁, F₂, and F₃, as well as the rotational inertia term I_oẇ_o.

25

F_{x o} = \frac{2 F_{2} - F_{3} - F_{1}}{3} + m_{o} a_{x}

{F_{xo}} = {{2{F_2} - {F_3} - {F_1}} \over 3} + {m_o}{a_x}

26

F_{y o} = \frac{\sqrt{3} F_{3} - \sqrt{3} F_{1}}{3} + m_{o} a_{y}

{F_{yo}} = {{\sqrt 3 {F_3} - \sqrt 3 {F_1}} \over 3} + {m_o}{a_y}

Equations (25) and (26) calculate the forces along the X_o and Y_o axes.

The forces exerted by each motor can be calculated as: 27 $\begin{array}{l} F_{1} = \frac{I_{w o} {\ddot{θ}}_{1}}{r_{o}}, & F_{2} = \frac{I_{w o} {\ddot{θ}}_{2}}{r_{o}}, & F_{3} = \frac{I_{w o} {\ddot{θ}}_{3}}{r_{o}} \end{array}$ \matrix{ {{F_1} = {{{I_{wo}}{{\ddot \theta }_1}} \over {{r_o}}},} \hfill & {{F_2} = {{{I_{wo}}{{\ddot \theta }_2}} \over {{r_o}}},} \hfill & {{F_3} = {{{I_{wo}}{{\ddot \theta }_3}} \over {{r_o}}}} \hfill \cr }

In equation (27), the forces generated by the individual motors are related to the motor angular accelerations ${\ddot{θ}}_{1}, {\ddot{θ}}_{2},$ {\ddot \theta _1},{\ddot \theta _2}, and ${\ddot{θ}}_{3}$ {\ddot \theta _3} and depend on the moment of inertia I_wo and the radius r_o.

To express the forces F_xo and F_yo in terms of the motor angular accelerations, we use the following equations.

28

F_{x o} = A_{11} {\ddot{θ}}_{1} + A_{12} {\ddot{θ}}_{2} + A_{13} {\ddot{θ}}_{3}

{F_{xo}} = {A_{11}}{{\ddot \theta }_1} + {A_{12}}{{\ddot \theta }_2} + {A_{13}}{{\ddot \theta }_3}

29

F_{y o} = A_{21} {\ddot{θ}}_{1} + A_{22} {\ddot{θ}}_{2} + A_{23} {\ddot{θ}}_{3}

{F_{yo}} = {A_{21}}{{\ddot \theta }_1} + {A_{22}}{{\ddot \theta }_2} + {A_{23}}{{\ddot \theta }_3}

30

T_{z o} = A_{31} {\ddot{θ}}_{1} + A_{32} {\ddot{θ}}_{2} + A_{33} {\ddot{θ}}_{3}

{T_{zo}} = {A_{31}}{{\ddot \theta }_1} + {A_{32}}{{\ddot \theta }_2} + {A_{33}}{{\ddot \theta }_3}

Equations (28) (29) and (30) are expressed with coefficients to shorten the equations.

These relationships can be written in matrix form: 31 $[\begin{array}{l} F_{x o} \\ F_{y o} \\ T_{z o} \end{array}] = [\begin{array}{l} A_{11} & A_{12} & A_{13} \\ A_{21} & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{array}] [\begin{array}{l} {\ddot{θ}}_{1} \\ {\ddot{θ}}_{2} \\ {\ddot{θ}}_{3} \end{array}]$ \left[ {\matrix{ {{F_{xo}}} \hfill \cr {{F_{yo}}} \hfill \cr {{T_{zo}}} \hfill \cr } } \right] = \left[ {\matrix{ {{A_{11}}} \hfill & {{A_{12}}} \hfill & {{A_{13}}} \hfill \cr {{A_{21}}} \hfill & {{A_{22}}} \hfill & {{A_{23}}} \hfill \cr {{A_{31}}} \hfill & {{A_{32}}} \hfill & {{A_{33}}} \hfill \cr } } \right]\left[ {\matrix{ {{{\ddot \theta }_1}} \hfill \cr {{{\ddot \theta }_2}} \hfill \cr {{{\ddot \theta }_3}} \hfill \cr } } \right]

Equation (31) expresses the forces and torques in terms of the motor angular accelerations in matrix form.

The coefficients A_ij are defined as follows: 32 $A_{11} = A_{13} = - \frac{I_{w o}}{3 r_{o}} - m_{o} \frac{r_{o}}{3}$ {A_{11}} = {A_{13}} = - {{{I_{wo}}} \over {3{r_o}}} - {m_o}{{{r_o}} \over 3} 33 $A_{12} = \frac{2 I_{w o}}{3 r_{o}} + m_{o} \frac{2 r_{o}}{3}$ {A_{12}} = {{2{I_{wo}}} \over {3{r_o}}} + {m_o}{{2{r_o}} \over 3} 34 $A_{21} = - A_{23} = - \frac{\sqrt{3} I_{w o}}{3 r_{o}} - m_{o} \frac{\sqrt{3} r_{o}}{3}$ {A_{21}} = - {A_{23}} = - {{\sqrt 3 {I_{wo}}} \over {3{r_o}}} - {m_o}{{\sqrt 3 {r_o}} \over 3} 35 $A_{22} =0$ {A_{{\rm{22}}}}{\rm{ = 0}} 36 $A_{31} = A_{32} = A_{33} = \frac{I_{w o} d_{o}}{r_{o}} + I_{o} \frac{r_{o}}{3 d_{o}}$ {A_{31}} = {A_{32}} = {A_{33}} = {{{I_{wo}}{d_o}} \over {{r_o}}} + {I_o}{{{r_o}} \over {3{d_o}}}

Equations (32) through (36) define the coefficients A_ij. The angular accelerations of the motors can be expressed as a function of the forces and torques.

37

[\begin{array}{l} {\ddot{θ}}_{1} \\ {\ddot{θ}}_{2} \\ {\ddot{θ}}_{3} \end{array}] = [\begin{array}{l} B_{11} & B_{12} & B_{13} \\ B_{21} & B_{22} & B_{23} \\ B_{31} & B_{32} & B_{33} \end{array}] [\begin{array}{l} F_{x o} \\ F_{y o} \\ T_{z o} \end{array}]

\left[ {\matrix{ {{{\ddot \theta }_1}} \hfill \cr {{{\ddot \theta }_2}} \hfill \cr {{{\ddot \theta }_3}} \hfill \cr } } \right] = \left[ {\matrix{ {{B_{11}}} \hfill & {{B_{12}}} \hfill & {{B_{13}}} \hfill \cr {{B_{21}}} \hfill & {{B_{22}}} \hfill & {{B_{23}}} \hfill \cr {{B_{31}}} \hfill & {{B_{32}}} \hfill & {{B_{33}}} \hfill \cr } } \right]\left[ {\matrix{ {{F_{xo}}} \hfill \cr {{F_{yo}}} \hfill \cr {{T_{zo}}} \hfill \cr } } \right]

Equation (37) is the inverse of Equation (31). This inverse equation calculates ${\ddot{θ}}_{1}, {\ddot{θ}}_{2},$ {\ddot \theta _1},{\ddot \theta _2}, and ${\ddot{θ}}_{3}$ {\ddot \theta _3} based on the equivalent forces and torques on the robot. The coefficients Bij are defined as follows: 38 $B_{11} = B_{31} = - 2 B_{21} = - \frac{1}{2} (I_{w o} r_{o} + m_{o} r_{o})$ {B_{11}} = {B_{31}} = - 2{B_{21}} = - {1 \over 2}\left( {{I_{wo}}{r_o} + {m_o}{r_o}} \right) 39 $B_{22} =0$ {B_{{\rm{22}}}}{\rm{ = 0}} 40 $B_{12} = - B_{32} = \frac{- 3}{2 (\sqrt{3} I_{w o} r_{o} + \sqrt{3} m_{o} r_{o})}$ {B_{12}} = - {B_{32}} = {{ - 3} \over {2\left( {\sqrt 3 {I_{wo}}{r_o} + \sqrt 3 {m_o}{r_o}} \right)}} 41 $B_{13} = B_{12} = B_{11} = \frac{d_{o} r_{o}}{3 I_{w o} d_{o}^{2} + I_{o} r_{o}^{2}}$ {B_{13}} = {B_{12}} = {B_{11}} = {{{d_o}{r_o}} \over {3{I_{wo}}d_o^2 + {I_o}r_o^2}}

Equations (38) through (41) define the coefficients of B_ij.

2.3.

Mobile Robot Integration

The integration of two systems, the trailer and the omniwheel mobile robot, can be modeled using the dynamic analysis of both systems. When these systems are integrated as shown in the Figure 4, the external forces acting on the trailer, denoted F_ex and F_ey, are equal to the tractor’s resulting tracking forces F_xo and F_yo, as expressed in Equation (28) and Equation (29).

By multiplying the Jacobian from the dynamic models of both the trailer and the omniwheel mobile robot, we can establish the relationship between the wheel accelerations of the omniwheel mobile robot and the trailer wheel velocities. The resulting equation is provided below, which models this relationship.

42

[\begin{array}{l} {\ddot{θ}}_{1} \\ {\ddot{θ}}_{2} \\ {\ddot{θ}}_{3} \end{array}] = K_{o} K_{t} [\begin{array}{l} {\ddot{θ}}_{L} \\ {\ddot{θ}}_{R} \end{array}]

\left[ {\matrix{ {{{\ddot \theta }_1}} \hfill \cr {{{\ddot \theta }_2}} \hfill \cr {{{\ddot \theta }_3}} \hfill \cr } } \right] = {K_o}{K_t}\left[ {\matrix{ {{{\ddot \theta }_L}} \hfill \cr {{{\ddot \theta }_R}} \hfill \cr } } \right]

Equation (42) defines how the wheel accelerations of the omniwheel mobile robot are related to the trailer wheel velocities through the product of the Jacobian matrices K_o and K_t.

The matrix K_t, representing the trailer’s dynamic model, is given by: 43 $K_{t} = [\begin{matrix} \frac{I_{w t} L}{2 d_{t} r_{t}} + \frac{I_{t} r_{t}}{d_{t} L} & - \frac{I_{w t} L}{2 d_{t} r_{t}} - \frac{I_{t} r_{t}}{d_{t} L} \\ \frac{I_{w t}}{r_{t}} + m_{t} \frac{r_{t}}{L} & \frac{I_{w t}}{r_{t}} + m_{t} \frac{r_{t}}{L} \\ 0 & 0 \end{matrix}]$ {K_t} = \left[ {\matrix{ {{{{I_{wt}}L} \over {2{d_t}{r_t}}} + {{{I_t}{r_t}} \over {{d_t}L}}} & { - {{{I_{wt}}L} \over {2{d_t}{r_t}}} - {{{I_t}{r_t}} \over {{d_t}L}}} \cr {{{{I_{wt}}} \over {{r_t}}} + {m_t}{{{r_t}} \over L}} & {{{{I_{wt}}} \over {{r_t}}} + {m_t}{{{r_t}} \over L}} \cr 0 & 0 \cr } } \right]

Equation (43) provides the matrix K_t, which takes into account the physical properties of the trailer, including moments of inertia I_wt and I_t, the wheel radius r_t, and the distance between the wheels. This matrix provides the relationship between the trailer’s wheel accelerations and its dynamics.

The matrix K_o, which models the dynamic properties of the omniwheel mobile robot, is defined as: 44 $K_{o} = [\begin{array}{l} B_{11} & B_{12} & B_{13} \\ B_{21} & B_{22} & B_{23} \\ B_{31} & B_{32} & B_{33} \end{array}]$ {K_o} = \left[ {\matrix{ {{B_{11}}} \hfill & {{B_{12}}} \hfill & {{B_{13}}} \hfill \cr {{B_{21}}} \hfill & {{B_{22}}} \hfill & {{B_{23}}} \hfill \cr {{B_{31}}} \hfill & {{B_{32}}} \hfill & {{B_{33}}} \hfill \cr } } \right]

Equation (44) shows the matrix K_o, where the elements B_ij represent coefficients that account for the robot’s geometry and dynamics, including the effects of forces and torques on the wheel accelerations. This matrix is essential in transforming the trailer’s wheel velocities into the mobile robot’s motor accelerations.

3.

Control Method

As it is shown in dynamical analysis. The control task for motion control of trailer realized by applying the towing forces which are generated via the omni-wheel tractor robot. Therefore, control requirements are given as: path planning between way points concerning non-holonomic trailer kinematics, derivation of required forces using its differential kinematics based control strategy and actuator level control of omniwheel robot motion for controlled towing.

3.1.

Path Planning

Path planning in trailer-tractor robots are essential in order to achieve autonomous movement [4, 5, 20]. This autonomous movement can be planned in different ways, they can aim choosing the most efficient, feasible or wanted set of movements for this reason the path planning of mobile robots is a complicated and important task. Path planning of mobile robots can be done using both mathematical functions and algorithms depending on the path planning strategy. A path planning strategy for a mobile robot is often determined by considering a planned scenario or a set of movements, the planned scenario includes the work environment and all of the constraints of the mobile robot. An example for this is if a mobile robot is in a closed and blocked environment, the movement of the mobile robot is unconstrained and full autonomy of the mobile robot is desired from the path planner, the most common path planning strategies uses local and global maps [23,24] to determine the most efficient movement sets by help of many different algorithms [25–28]. These path planning strategies will aim to be the most beneficent for that certain hypothetical scenarios. In this research a path planning method is required in order to test the path tracking success of the trailer tractor robot.

Due to this need, a hypothetical scenario and all of the defined constraints of the trailer tractor robot are identified and based on this the designed path planning is simulated and tested in MATLAB Simulink.

In an ideal environment without obstacles, the simplest method to achieve the desired pose is to move along the shortest linear path. However, this straight path is not suitable for the given problem because it would result in rotations around the center of mass of the trailer. These rotations around the center of mass of the trailer is harder to apply because of the non-holonomic structure of the trailer. To solve this, a cubic Bézier curve (with n = 3) is employed to generate a path that avoids these rotations. The general form of the Bézier curve [29] is given by Equation (45).

45

\begin{array}{l} f (t) = \sum_{i = 0}^{n} b_{i} B_{i}^{n} (t), & 0 \leq t \leq 1 \end{array}

\matrix{ {f(t) = \sum\limits_{i = 0}^n {{b_i}} B_i^n(t),} \hfill & {0 \le t \le 1} \hfill \cr }

This equation expresses the Bézier curve as a weighted sum of basis functions $B_{i}^{n} (t)$ B_i^n(t), where b_i are the control points, and t is the parameter that varies between 0 and 1. This sum produces a curve that interpolates between the control points.

The basis functions B_i(t) for the Bézier curve are defined in Equation (46).

46

\begin{array}{l} B_{i} (t) = (\begin{matrix} n \\ i \end{matrix}) t^{i} {(1 - t)}^{n - i}, & i = 0, \dots, n \end{array}

\matrix{ {{B_i}(t) = \left( \matrix{ n \cr i \cr} \right){t^i}{{(1 - t)}^{n - i}},} \hfill & {i = 0, \ldots ,n} \hfill \cr }

These functions determine the influence of each control point on the final curve. The coefficients are computed using binomial expansions, ensuring that the curve smoothly interpolates between the given control points.

Algorithm 1 Path Generation Procedure

1: procedure PathGenerated (P_f, P_i)

2: P₀ ← [P_i(1); P_i(2)]

3: P₃ ← [P_f(1); P_f(2)]

4: L ← norm(P₃ – P₀)

5: P₁ ← (P₀ + 0.5L)[cos(P_i(3)); sin(P_i(3))]

6: P₂ ← P₃ – 0.5L [cos(P_f(3)); sin(P_f(3))]

7: N ← 50

8: t ← linspace(0,1, N)

9: patd_generated ← zeros(2,N)

10: for i = 1 to N do

11: patft_generated (:, i) ← (1 – t(i))³P₀ + 3(1 – t(i))²t(i)P₁ + 3(1 – t(i))t(i)²P₂ + t(i)³P₃

12: end for

13: end procedure

For cubic Bézier curves, where n = 3, the curve is determined by four control points. The path planning algorithm utilizes two control points P₁ and P₂ along with the initial and desired poses of the trailer. Control point P₁ is positioned towards the front of the trailer’s initial pose, while control point P₂ is placed in the opposite direction, towards the front of the trailer’s desired pose as shown on the Figure 5.

The cubic Bézier curve ensures that the trailer follows a smooth path without excessive rotation around its center of mass, providing an efficient solution for path planning.

This approach avoids the complexities of more advanced algorithms while still achieving the desired motion. The generated path is calculated using the path generation procedure, described in Algorithm 1.

3.2.

Path Tracking

Path tracking is the process that allows autonomous systems to travel along a desired trajectory. The process typically consists of several sets of rules; hence, these are referred to as algorithms.

The path tracking algorithms are responsible for determining the next waypoint on the planned path and calculating the desired velocities to reach the way-point. Many different path tracking algorithms exist, but one of the most suitable algorithms for the trailer tractor system is the pure pursuit algorithm [6].

The pure pursuit algorithm, implemented in MATLAB Simulink, uses a specific distance called the lookahead point distance to decide the next waypoint, as shown in Figure 6. The waypoints are generated from the path planner by dividing the generated path into a number of waypoints and listing them from start to finish. Each waypoint consists of the position and orientation of the vehicle in task space. The next way-point is the farthest point within the radius of the lookahead distance.

The waypoints and the actual pose of the vehicle are used to calculate the steering angle. The algorithm can generate a constant linear velocity, and the angular velocity can be bounded to the desired limits.

The success of trajectory tracking depends on the following factors:

- The number of waypoints in the path and the lookahead distance are important. Within the lookahead distance radius, there must always be at least one waypoint.
- The vehicle’s desired linear velocity and angular velocity should be proportional to the lookahead distance. A lower lookahead distance causes the vehicle to oscillate around the planned path, while a higher lookahead distance causes the vehicle to converge to the planned path over a longer time. The parameters should be tuned with a constant distance between each waypoint for the algorithm to work regardless of the planned path.
- The algorithm always has a slight error due to the lookahead distance; therefore, the tuning of parameters is done by considering oscillations on the planned path.
- The pure pursuit algorithm does not stop the vehicle at the last waypoint; therefore, a stopping algorithm is needed.

Algorithm 2 Pure Pursuit Path Tracking Algorithm

1: procedure PurePursuit(P, P_v, θ_v, L)

2: d_min ← Max

3: P_lh ← None

4: for each P_i in P do

5: d ← Dist(P_v,P_i)

6: if d ≥ L and d < d_min then

7: d_min ← d

8: P_lh ← P_i

9: end if

10: end for

11: dx ← P_lh(1) – P_v(1)

12: dy ← P_lh(2) – P_v(2)

13: α ← atan2(dy, dx) – θ_v

14: $k \leftarrow \frac{2 \cdot \sin (α)}{L}$ k \leftarrow {{2\cdot\sin (\alpha )} \over L}

15: δ ← atan(k)

16: return δ

17: end procedure

3.3.

Dynamic Tracking Control

The problem definition states that one of the objectives of this study is to control the trailer pose p_c and velocity V_t, ${\dot{ϕ}}_{t}$ {\dot \phi _t} in the task space. This means that trailer dynamics need to be controlled with a tracking control method. The dynamic tracking control method should maintain the velocity at the referenced value while tracking continues, by calculating the forces or torques that need to be applied.

The dynamic tracking control equations for the trailer are given by: 47 $F_{x} = - I_{t} {\dot{w}}_{c} - K_{b} w_{c},$ {F_x} = - {I_t}{\dot w_c} - {K_b}{w_c}, 48 $F_{y} = m_{t} {\dot{V}}_{c} + K_{a} V_{c},$ {F_y} = {m_t}{\dot V_c} + {K_a}{V_c},

The dynamic tracking control based on the Lyapunov function [15, 16] works well for the non-holonomic wheeled mobile robot. As shown above Towing a trailer is a non-holonomic motion control problem which results in nonlinear mathematical model. Hence, Lyapunov control design method results in deployable control laws compatible with kinematics, it is adopted in our problem as a complementary part of motion planning [15, 16]. In Figure 8, the simulation setup for dynamic tracking control is shown.

Where ${\dot{V}}_{c}$ {\dot V_c}{\rm{ }} is the trailer’s acceleration, K_a is the tangential velocity control coefficient, V_c is the trailer’s desired tangential velocity, ẇ_c is the trailer’s angular acceleration, w_c is the trailer’s desired angular velocity, and K_b is the angular velocity control coefficient. The Equations (47) and (48) describe the forces that need to be applied to maintain the desired velocity and orientation of the trailer. The control coefficients K_a and K_b are tuned to ensure that the path tracking is performed successfully. The method is based on Lyapunov stability, which guarantees that the system will remain stable while following the desired path.

3.4.

Motor Tuning

The actuator is the selected motor Pololu 37Dx73L, 12 V 131:1 gear motor with an internal 64 CPR two-channel encoder. This motor has a gearbox in front of it and has an encoder at the back, which causes the inaccuracies in the encoder readings to be less effective on the motor torque or position.

This motor is a gray box system model, where system models are called gray box when the theoretical structure is known, but the parameters are unknown. Gray box model systems are theoretical white box models because, with proper data, these system parameters can be identified, and this process is called system identification.

The system identification starts with data acquisition. To collect data from the three DC gear motors, a setup containing an NI myDAQ card, Arduino Mega, and L298N motor driver is made. With this setup, different PWM values are given to the motor, and their encoder data are collected from one channel. With the help of a myDAQ card used as a switch and data collector, the start of the motor is coordinated with the given input voltage. The encoder data is collected with a sampling rate of 50 kHz. All of the collected samples are gathered for 5 seconds. The PWM voltage is given in the first second and cut off at the exact third second. All of the encoder data collected is stored in an Excel file. This collected data must be processed to be used in system identification.

The theoretical structure of this system is most similar to the transfer function of a basic armature-controlled DC motor. The armature-controlled DC motor transfer function has the input as voltage and the output as angular velocity. System identification data collected from encoder channels need to be converted to motor angular velocities, and the input PWM should be converted to the equivalent voltages. These conversions are made in MATLAB using MATLAB scripts. The conversion for equivalent voltage of PWM is simple: multiplying the 12 V with the duty percentage works. The encoder data collected provides the rotation information, with each rising pulse counted, and 16 counts equals one rotation. With the time data collected, we can find the RPM values by dividing the total rotation counted over seconds. Disturbed encoder data is cleaned and filtered with a low-pass filter. The change of motor angular velocity and input voltage over time is given in Figure 9.

These data conversions are made for every five input voltage duty percentages, and encoder data are acquired from all three motors. When the conversions are done, the data are used in the MATLAB system identification toolbox. The transfer function is obtained using an iterative curve-fitting method.

The System Identification Toolbox in MATLAB requires the user to know the number of poles and zeros of the transfer function. The required pieces of information are known since the system is a gray box model, the system has two poles and no zeros because it works similarly to basic armature-controlled DC motors.

Figure 10 shows the curve fitting of the transfer function. The created transfer function’s system response is %83.63 accurate according to the processed encoder data. The motor transfer function is known, which means we can use other MATLAB tools to find the optimal PID coefficients. Finding the PID coefficients can be done with another toolbox called PID Tuner. While tuning the PID, the control effort should be saturated. This saturation is necessary to apply these PID values in real-life applications and to see the same results in simulation. The control effort is the voltage input of this system, which means it can only be between 12 V and -12 V. The output can only be between 1050 rad/s and -1050 rad/s. These saturation values come from the specs of the components used. The sample time of the simulation must be the same as the sample time of the used hardware, and since the hardware is digital, a fixed-step solver should be used.

The motor transfer function is controlled using the tuned PID parameters in Figure 11–12. The control system reaches a steady state in 0.5 seconds. Steadystate error and overshoot of the system response are zero.

49

G (s) = \frac{62840}{s^{2} + 47.1 s + 738.4}

G(s) = {{62840} \over {{s^2} + 47.1s + 738.4}}

50

u (s) = P + I \frac{1}{s} + D \frac{N}{1 + N \frac{1}{s}}

u(s) = P + I{1 \over s} + D{N \over {1 + N{1 \over s}}}

This method is applied to all three motors, and different PID parameters are tuned. Table 1 shows the PID parameters tuned for all three motors. Equation (50) describes the effect of PID control parameters on the control effort.

Table 1.

Tuned PID Parameters

	P	I	D	N
Motor 1	0.231	0.358	0.007	15796.2
Motor 2	0.158	0.493	0.011	5750.9
Motor 3	0.166	0.535	0.012	6962.8

4.

Simulation Procedure and Results

The proposed mobile robot is drawn in a 3D CAD program, and the mobile robot is exported to MATLAB Simulink using Simscape. Simscape is a Simulink add-on that can generate mechanical and electromechanical models inside a simulation world, working together with MATLAB. This tool allows this research to be tested in an ideal and known simulation setting. Thus, the allowed degrees of freedom of the world, the friction coefficient, and many other variables can be adjusted. Additionally, simulation results or any variable can be measured, viewed, and logged.

Before simulating the full towing robot, the trailer dynamics, dynamic tracking control, path planning, and path tracking are tested using a simplified Simscape model displayed on Figure 13–14 that includes only the trailer portion of the mechanical system shown on Figure 19. This approach allows for a more manageable analysis by dividing the problem into two parts, facilitating the evaluation of the control methods and algorithms effectiveness [30].

The path generator subsystem shown in Figure 15 uses the pseudocode in Algorithm 1 and creates the path of the cubic Bézier function. The path is divided into waypoints and a lookahead distance is calculated. The waypoints are then used in the path tracking algorithm to calculate the desired velocities. The desired output velocities are then used in the dynamic tracking control subsystem shown in Figure 16. This subsystem uses the control equations and produces the desired external towing forces.

The Simscape Multibody Link Plugin exports the trailer’s geometrical and physical attributes to MATLAB Simulink from SolidWorks assembly and saves attributes in a MATLAB file, Simscape schematic in a Simulink file, and model appearance in STL or STEP format. After the export, some variables can be changed to adjust for the assumptions made for the convenience of the analysis in the first chapter. The Simscape mechanical model is created in an empty world when first exported. The movement plane, main body constraints with the world frame, wheel contact parameters, body frames, and sensors are added to the subsystem of the mechanical Simscape model shown in Figure 17. The generated path is shown in simulation Figures 23–24 with a spline. The frame created on the towing point is used for applying external towing forces, and the frame on the center of mass is used for position and velocity measurements.

The pure pursuit algorithm is insufficient when the mobile robot needs to stop. Therefore, the current and final positions of the trailer are used to create a stopping algorithm. The distance between the center of mass of the trailer and the final waypoint is calculated.

When the trailer is within a radius of 0.01 m of the final position, the trailer stops moving and the simulation ends. Figure 19 illustrates the simulation settings from isometric, front, top, and bottom views.

After the trailer control simulation approves the success of the control methods and algorithms, the omniwheel mobile robot is added to the Simscape mechanical model, and the external towing forces are replaced with the inverse dynamics of the omniwheel mobile robot as shown in Figure 20. The stopping, path tracking, and path planning algorithms remain the same; therefore, the implementation of the mechanical model and inverse dynamic analysis is simpler.

Figure 21 illustrates the same method of dynamic tracking control used in Figure 16, but the overall dynamics of the system has changed, so the coefficients K_a and K_b are adjusted again.

The simulation is shown on Figure 22. The analysis of the simulations are done by investigating the graphs generated from the simulation data. The dynamic tracking control parameters, path tracking parameters, implementation of selected motors to the simulation and PID tuning of the motors are done by analyzing the results. Analysis of the success of path tracking requires the simulation path to be generated in different forms such as square, circle and archimedian curved paths rather than using the path generation algorithm to generate Bézier curved paths.

4.1.

Simulation Results for Tuned Dynamic Tracking Control Parameters

The tuning of the dynamic tracking control parameters (K_a, K_b) are done by analyzing the settling time, rise time, overshoot and steady state error data gained from simulations. The parameter values are kept when a desired result is reached. In this case the control parameters are K_a =200 and K_b = 300.

Figure 23 demonstrates the path taken by trailer when the tuned dynamic tracking control parameters are used. The path taken by the trailer remained correct. The trailer control simulation has an RMS error of 29.7 mm, and the tractor robots simulation has a 3.23 mm RMS error. Although the error seems lower in the tractor robot simulation, the path taken by the trailer is closer to the desired path in the trailer control simulation. The lookahead distance in the tracking algorithm of both simulations was the same; therefore, the difference must be due to the dynamic tracking control coefficients. The reason for the error on way-points 49 and 50 is that in both of the simulations, the trailer starts to shift its orientation irrelevant to the path planning when the lookahead distance radius cannot detect any more waypoints ahead. The path tracking algorithm needs to be adjusted so that the last point can be detected. The solution to the path tracking algorithm is to adjust the lookahead distance according to the stopping radius in the stopping algorithm.

The dynamic tracking control method successfully controls the trailer with the desired velocities in the trailer control simulation shown in Figures 24 and 25. The tractor robot simulation is controlled with the same dynamic tracking control method, but the results have high steady-state errors. This steady-state error can result from the hand-tuned dynamic tracking control parameters. The settling time of the trailer control simulation is faster than the tractor robot simulation. The trailer control simulation’s linear velocity controls are slow and robust with no overshoot, while the tractor robot simulation’s linear velocity control is fast and aggressive and has overshoots. The simulated angular velocities of both simulations have difficulty settling, and the desired angular velocity is always changing. The rapid change of the desired value causes the simulated values to have errors. The tractor robot requires a faster response to rapid changes, or another solution is to reduce the rate of rapid changes. Therefore, better tuning of the dynamic tracking control coefficients is needed to achieve faster settling times.

Figures 24–25 shows the measured and desired velocities of the trailer when simulation is started with the tuned dynamic tracking control parameters.

The use of cubic Bézier curves for path generation and the implementation of the path tracking algorithm both functions effectively. The resulting simulation path aligns closely with the desired trajectory and the trailer successfully reaches its target pose with minimal discrepancies. Both simulations achieve the samedesired pose by employing identical path planning and tracking algorithms. However, in the trailer control simulation, the trailer reaches the desired pose more quickly compared to the tractor robot’s simulation. This difference arises from variations in the trailer’s average speed; specifically, the slower average speed observed in the tractor robot is attributed to steady-state errors and overshoots within the control mechanism. As shown in Figure 23, at the coordinate [1, 1.6], the simulation of the tractor robot reveals that the trailer deviates from its intended path by 0.05 meters. This deviation persists until the target point is finally reached.

After tuning the dynamic tracking control coefficients, the simulations are repeated, and the desired linear velocity is reached with a 3.96% steady-state error, as shown in Figure 24, while the error occurring in the desired angular velocity is decreased (Figure 25). The tuned dynamic tracking control coefficients are used in other paths, and the trajectory control is analyzed.

4.2.

Simulation Results for Different Paths

The success of the path tracking is further analyzed by generating non-optimal paths for path tracking in Simulink. The tuned dynamic tracking control parameters are used and the effect of lookahead distance is analyzed on a square path. The lookahead distance is initially 0.25 m, for the square path analysis two lookahead distances are examined, 0.75 m and 1.25 m.

Figures 26–28 are the result of simulations when the lookahead distance is 0.75 m.

The square path analysis in Figures 26–31 shows the importance of the lookahead distance. Errors on the track appear when a sharp turn is ahead and a higher lookahead distance is needed to track the path. The desired angular velocities change with the lookahead distance in both of the lookahead distances.

The desired velocities are achieved in Figures 27 and 30, therefore, the path error only happens because the desired angular velocities are incorrect. Figures 29–31 are the result of simulations when the lookahead distance is 1.25 m and the trailer is following a square path.

Circular path analysis uses the same tuned parameters and analyzes the behavior of the tracking algorithm in circular motion. Figure 32 shows the circular path taken by the towed trailer in simulation.

The circular path in Figure 32 is tracked with the same parameters as the Bézier curve path planning simulations. The circular path should have a fixed desired angular velocity, but the small oscillations in Figure 34 are due to the low resolution of the path. The desired linear velocity is reached within an acceptable margin of error.

The circular path and the archimedian curve path are analyzed with the tuned dynamic tracking control parameters (K_a = 200 and K_b = 300.) and the initial lookahead distance (0.25 m).

The Archimedean spiral path in Figure 35 is tracked with an error at the beginning of the path. The curve at the beginning of the path is sharp, so the maximum value of the angular velocity is desired in Figure 37. It is observed that the limit angular velocity of 0.1 rad/s is not enough to take the turn without errors. Even with a larger curved turn, the trailer finds its track and reaches the desired velocity as shown in Figure 36.

4.3.

Simulation Results for Implementation of PID Controlled DC Motors

The three motors selected for the towing robot are integrated in the simulation in this section. In the simulation the estimated transfer function for the three motors are used to simulate the electromechanical effect motors. In order to simulate the electromechanical limits of the motor, stall torques are defined using saturation blocks. PID control of the motors can now be implemented in the simulation with the found PID parameters.

The implemented on Simulink subsystem is shown in Figure 39.

The characteristics of the selected dc motors and the PID control parameters determined with tests are integrated into the simulation to apply a feasible simulation setting in real life application, as shown in Figures 38–39. The integration of the motors and PID control in the simulation caused the trailer to have small oscillations around the desired path shown in Figure 40. The linear and angular velocities of the trailer in Figures 41–42 are closer to the desired inputs when the PID control is applied.

4.4.

Comparison of Results

The performance results of the simulations are presented in Table 2 for each step of different tuning and implemented control system for tractor trailer mobile robot. The performance criteria for the simulation results are; steady state linear velocity error (V_ess,%) to analyze the closeness of the final linear velocity to the desired value the in order to tune the controllers, linear velocity overshoot (%) to analyze the initial maximum linear velocity error of the mobile robot for the determination of success of the tuned controllers, steady state time (T_ss, s) in order to determine the time that the controller takes to be in the %5 range of the desired output, maximum and average angular velocity error $({\dot{ϕ}}_{emax}, {\dot{ϕ}}_{eavg}, rad / s)$ \left( {{{\dot \phi }_{{\rm{emax }}}},{{\dot \phi }_{{\rm{eavg }}}},{\rm{rad}}/{\rm{s}}} \right), execution time (s), tracking error (Max, RMS, m) and tracked and planned path length (m) to decide if the tracking control algorithms and path planning are effective. Each performance criteria has a purpose and effects the tuning of control parameter. All of the results in the Table 2 are simulated using the same path planning and path tracking algorithms to achieve the same pose from the same initial conditions and with equal conditions. The difference between non-tuned, tuned and motor integrated PID tuned results is, non-tuned results refer to the results when there is no PID controller, dynamic tracking parameters K_a is 50 and K_b is 120 and they are not tuned for integrated trailer tractor mobile robot, tuned results are from the Figures 23–25 when dynamic tracking parameters K_a is 200 and K_b is 300 without PID control and motor integration, the motor integrated PID tuned results are from the Figures 38–42 when the motors transfer function and the tuned PID controller is added to the control system, working together with the tuned dynamic tracking.According to these results comments about success of the methods can be made.

Table 2.

Performance Comparison of Non-Tuned, Tuned, and Motor Integrated PID Tuned Controllers

Metric	Non-Tuned	Tuned	Motor Integrated PID Tuned
Steady State Linear Velocity Error (V_ess,%)	13.29	6.26	4.32
Linear Velocity Overshoot (%)	8.68	17.48	17.17
Steady State Time (T_ss, s)	N/A	N/A	61.9
Maximum Angular Velocity Error $({\dot{ϕ}}_{emax}, rad / s)$ \left( {{{\dot \phi }_{{\rm{emax }}}},{\rm{rad}}/{\rm{s}}} \right)	0.104	0.092	0.117
Average Angular Velocity Error $({\dot{ϕ}}_{eavg}, rad / s)$ \left( {{{\dot \phi }_{{\rm{eavg }}}},{\rm{rad}}/{\rm{s}}} \right)	0.010	0.006	0.008
Execution Time (s)	74.7720	68.7573	67.2488
Tracking Error (Max, m)	0.0777	0.0755	0.0756
Tracking Error (RMS, m)	0.0425	0.0408	0.0404
Tracked Path Length (m)	6.620	6.620	6.630
Planned Path Length (m)	6.627	6.627	6.627

In order to compare the success of the simulated control method, results regarding control performance collected from references that uses different methods for path tracking and control of tractortrailer mobile robots [21, 31, 32]. The methods such as adaptive control, backwards motion control, LQR (Linear Quadratic Regulator) algorithm, model predictive control and PID control are used to control the trajectory successfully. These methods are just as effective as the dynamic tracking control and PID control used in this study. To make this claim, the simulation and test results of these methods shown in different figures are compared and analyzed.

5.

Conclusion

The simulation results show the trailer controlled by external forces on the joint successfully tracks the generated path. The simulation results indicate that the designed path planning algorithm for trailer dynamics is working as intended. Path planning can now be developed when the trailer detects obstacles and the path planning is adjusted according to the map data. The pure pursuit tracking algorithm works well with trailer dynamics. The algorithm should have a built-in stopper when the last waypoint is the next target destination stopper should activate. When the trailer starts to move away from the last waypoint, it means that the trailer has passed by, the stopper algorithm should activate and stop the trailer by force.

The trailer’s dynamic tracking control is performed effectively and as for tractor robot, the control system is fast and accurate. Steady state error and settling time of the tractor robot can be improved by the better tuning of the dynamic tracking, PID and path tracking control parameters. Oscillations of the simulated angular velocities can be weakened by further tuning of the dynamic tracking controller coefficients or by increasing the resolution of the path generation. The overall trajectory control of the omniwheel tractor robot performs successfully.

Integration of Trailer and Tractor Mobile Robot Considering Towing Dynamics

Full Article

Paradigm

My account