Kinematic Modeling of a Combined System of Multiple Mecanum-Wheeled Robots with Velocity Compensation

In industry, combination configurations composed of multiple Mecanum-wheeled mobile robots are adopted to transport large-scale objects. In this paper, a kinematic model with velocity compensation of the combined mobile system is created, aimed to provide a theoretical kinematic basis for accurate motion control. Motion simulations of a single four-Mecanum-wheeled virtual robot prototype on RecurDyn and motion tests of a robot physical prototype are carried out, and the motions of a variety of combined mobile configurations are also simulated. Motion simulation and test results prove that the kinematic models of single- and multiple-robot combination systems are correct, and the inverse kinematic correction model with velocity compensation matrix is feasible. Through simulations or experiments, the velocity compensation coefficients of the robots can be measured and the velocity compensation matrix can be created. This modified inverse kinematic model can effectively reduce the errors of robot motion caused by wheel slippage and improve the motion accuracy of the mobile robot system.


Introduction
In recent years, intelligent and flexible manufacturing has motivated the development of autonomous mobile robots for workpiece and equipment handling and transportation [1,2], and in particular, automated guided vehicle (AGV) technology is widely studied and applied [3][4][5][6][7]. The omni-directional mobile AGV with Mecanum-wheeled robot platform, which has good driving force and easy control performance, can improve the utilization of workshop space and the efficiency of workshop transportation, and is very appropriate for transporting heavy goods in complex industrial environments [8,9]. The AGV that adopts a four-Mecanum-wheeled mobile robot platform with symmetrical structure is the most basic form and most widely used in industry [10][11][12]. In the field of large-scale equipment manufacturing, such as electric multiple unit (EMU) and aircraft, the objects are very heavy and large, and the carrying capacity and size of the four-Mecanum-wheeled AGV cannot meet the transportation requirements. In order to solve the problem, two schemes can be adopted: (1) Using a heavy-duty omnidirectional mobile platform with more Mecanum wheels, such as 8-or 12-wheeled platforms. These robot platforms can also be used in tandem to carry larger loads. For example, German company CLAAS (which has been acquired by MBB Industries AG, and is now Aumann Beelen GmbH) developed an omnidirectional mobile heavy-duty mobile handling The latest research on cooperative transportation mainly focuses on path planning and navigation algorithms in working environments, and control methods in the cooperative transportation process. Research on the Mecanum-wheeled robot platform is mainly focused on single-Mecanum-wheeled robots [8,9,[36][37][38][39], while there is less research on the kinematic of a combination system and motion compensation of multiple-Mecanum-wheeled robot platforms. However, studying the kinematic and characteristics of the combined omnidirectional mobile system is the basis for studying the motion control, path planning, and navigation of the combination system. In this research work, we study the combined mobile system of multiple Mecanum-wheeled robots as a whole, not as a cooperative multi-robot system, and the main research interest is the kinematic of this combined system and velocity compensation for the kinematic model, aiming to provide a theoretical basis for motion control of the combined system. This paper is organized as follows: In Section 2, on the basis of studying the kinematic constraints of a single Mecanum wheel in an omnidirectional mobile system, kinematic models of a four-Mecanumwheeled robot platform with symmetrical structure and a combined robot system composed of multiple Mecanum-wheeled robots are established, and a motion compensation model is proposed. In Section 3, in order to verify the correctness and feasibility of the kinematic and motion compensation models, the typical motion modes of a virtual prototype of a four-Mecanum-wheeled robot are simulated on RecurDyn. In Section 4, translation motion tests are carried out using a four-Mecanum-wheeled robot physical prototype, and the test results with and without motion compensation are compared. The tests in Section 4 verify the correctness of the motion and motion compensation models and the feasibility of the simulation method in Section 3. In Section 5, the motions of a variety of combined mobile system of a multiple Mecanum-wheeled robot are simulated on RecurDyn, and the kinematic and motion compensation models of a multi-robot system are verified.

Kinematic Constraint Model of a Single Mecanum Wheel
The kinematic constraints of the i-th Mecanum wheel of robot system O XYZ − consisting of n Mecanum wheels are shown in Figure 1 [40][41][42].  In Figure 1, the Cartesian coordinate systems of the i-th Mecanum wheel and the roller of the wheel are O wi − X wi Y wi Z wi and P i − g i h i z i ; r w and r r are the radius of the wheel and the roller, respectively; Sensors 2020, 20, 75 4 of 37 (l i , α i ) is used to describe the relative installation orientation of origin O of the body coordinate system and center O i of the wheel; and the angle between the X wi axis and l i is β i , which is defined as the installation attitude angle of the local coordinate system of the wheel. The velocity of the motion center is v o in the current state, and the angle between v o and the X axis is θ o ; . θ is the rotation angular velocity of the system when moving in the plane. The angle between the projection of X wi and h i on the plane is the tilt angle γ i (0 • < γ i < 90 • ) of the roller; and v gi is the velocity of the roller touching the ground on the i-th Mecanum wheel.
It is assumed that the movement between the roller and the ground is pure rolling, the contact point of the roller with the ground does not slip, and the instantaneous velocity is 0. According to the constraints of rolling and sliding, the following formulas can be obtained [20,43].
Because the rollers rotate passively, the velocity of roller v gi is an uncontrollable variable, which is usually not taken into account. By eliminating v gi from Equation (1), we obtain the following [21]: Let . ζ = .
x . y . θ T , then the form of the matrix of Equation (2) is as shown in Equation (3), that is, the inverse kinematic equation of any (i-th) Mecanum wheel.
The motion state of the robot in local coordinate system O − XYZ can be mapped to global coordinate system O I − X I Y I Z I , which is expressed as: where . ζ I = .
x I . y I

Kinematic Model of Four-Mecanum-Wheeled Mobile Robot with Symmetrical Structure
In this section, a four-Mecanum-wheeled robot with symmetrical structure is taken as the research object, as shown in Figure 2, and the coordination motion relationship of each wheel is discussed [44][45][46][47]. Choosing the geometrically symmetrical center as system origin O, the rectangular coordinate system XOY fixed with the mobile platform is established. The coordinate systems of Mecanum wheels X wi O wi Y wi (i = 1, 2, 3, 4) are established, taking the wheel centers as the system origins. According to the geometric characteristics of the robot, for any Mecanum wheel O wi , α i + β i = 0. According to Equation (3), the following formula can be obtained. . ϕ i = −1 r w sin γ i cos γ i sin γ i l sin(β i + γ i ) Sensors 2020, 20, 75 5 of 37 where roller angle γ 1 = γ 3 = −γ = −45 • , γ 3 = γ 4 = γ = 45 • , then the inverse kinematic equation for the four-Mecanum-wheeled robot is [17][18][19]21]: The inverse kinematic Jacobian matrix J is expressed as: In the global coordinate system X I O I Z I , the inverse kinematic equation of the four-Mecanumwheeled robot is Sensors 2020, 20, x FOR PEER REVIEW 5 of 37 where roller angle for the four-Mecanum-wheeled robot is [17][18][19]21]:

Kinematic Model of Multiple Mecanum-Wheeled Mobile Robot System
In Figure 3, in the global coordinate system is X I O I Y I , when a multiple Mecanum-wheeled-robot system composed of m (m = 1, 2, · · · j · · · , k) robots co-transports an object, the poses and relative positions of these robots remain unchanged. In order to better describe the motion of the multi-robot system, a common coordinate system X S O S Y S is established at a specified location. X j O j Y j is the local coordinate system of the j-the robot of this system, which consists of n j (n j = 1, 2, · · · i · · · , k) Mecanum wheels. The relative position of coordinate systems X S O S Y S and X j O j Y j of each robot is determined by the geometric parameters of the multi-robot system. Let . Φ jn j , J j , and . ζ j be the wheel rotation Sensors 2020, 20, 75 6 of 37 matrix, Jacobian matrix, and motion state in the local coordinate system of the j-th robot, respectively. The inverse kinematic equation of the j-th robot is shown in Equation (9): Sensors 2020, 20, x FOR PEER REVIEW 7 of 37  The robot configuration composed of m same four-Mecanum-wheeled robots connected end-toend and side-by-side is a common multi-robot system, as shown in   After transforming the local coordinate system X j O j Y j of the j-th robot to the designated coordinate , and the inverse kinematic equation of the robot is as follows: where The inverse kinematic equation of the j-th robot in global coordinate system In global coordinate system X I O I Y I , the inverse kinematic equation of the multi-robot mobile system composed of m (m = 1, 2, · · · j · · · , k) robots is shown as follows: Sensors 2020, 20, 75 7 of 37 If the multiple-robot system only performs translational motion, its inverse kinematic equation is as follows The robot configuration composed of m same four-Mecanum-wheeled robots connected end-to-end and side-by-side is a common multi-robot system, as shown in Figure 4. The coordinate system X S O S Y S of the robot system is usually established on the structurally symmetric center line.
. Φ j4 , J, and K j are the wheel velocity matrix, Jacobian matrix, and motion state transformation matrix, respectively, relative to the specified reference coordinate system X S O S Y S of the j-th four-Mecanum-wheeled robot. The inverse kinematic equation of the multi-robot configuration is as follows: Sensors 2020, 20, x FOR PEER REVIEW 8 of 37

Velocity Compensation of the Robot System
During the movement of the robot, due to manufacturing errors, wheel slippage, and other factors, there are errors between the actual velocity and the desired velocity, and the relative errors of longitudinal, lateral, and rotational angular speed are also different. Normally, the actual velocity of the robot is less than the set desired velocity. Different robots have different wheel slip rates and velocity errors due to different structural dimensions, manufacturing precision, number of wheels, deformation rate of wheel rollers, and friction coefficient between wheel and ground. For example, the precision of the roller angle of the Mecanum wheel has a great influence on the motion accuracy and velocity of the robot. In the following, taking a four-Mecanum-wheeled robot as an example, the influence of the precision of the roller angle on velocity is analyzed.
In order to analyze the influence of roller angle γ on the velocity of the robot, the partial

Velocity Compensation of the Robot System
During the movement of the robot, due to manufacturing errors, wheel slippage, and other factors, there are errors between the actual velocity and the desired velocity, and the relative errors of longitudinal, lateral, and rotational angular speed are also different. Normally, the actual velocity of the robot is less than the set desired velocity. Different robots have different wheel slip rates and velocity errors due to different structural dimensions, manufacturing precision, number of wheels, Sensors 2020, 20, 75 8 of 37 deformation rate of wheel rollers, and friction coefficient between wheel and ground. For example, the precision of the roller angle of the Mecanum wheel has a great influence on the motion accuracy and velocity of the robot. In the following, taking a four-Mecanum-wheeled robot as an example, the influence of the precision of the roller angle on velocity is analyzed.
In order to analyze the influence of roller angle γ on the velocity of the robot, the partial derivative of Equation (6) with respect to γ is calculated.
If there is an error in roller angle ∆γ and γ = 45 • , the wheel velocity will be compensated according to Equation (17): According to Equation (17), the longitudinal motion velocity of the robot is not related to the roller angle, so it is unnecessary to compensate the wheel velocity, the lateral velocity and rotation velocity need to be compensated; and the rotation velocity compensation rate is related to the structural parameters of the robot. If the angle error of each roller of each Mecanum wheel is different, the situation will be more complicated. Because the velocity error is the result of many factors, the contribution of each factor is difficult to measure accurately. For a specific robot, its geometric size, manufacturing error, wheel deformation, and friction coefficient with the designated ground have been determined. The influence of these factors on the speed error can be determined. The velocity error may have a certain functional relationship with the motion variables, such as velocity of the robot. At different velocities, the relative errors of the robot's longitudinal, transverse, and rotational velocity can be obtained by experiments, so that the velocity compensation coefficients with different velocities can be obtained. For the multi-robot system in Figure 3, the velocity compensation equation matrix of the j-th robot is: Then, the inverse kinematic equation of the j-th robot with velocity compensation is The inverse kinematic equation of this multi-robot system composed of m robots with velocity compensation is Sensors 2020, 20, 75 9 of 37 If the multi-robot system only performs translational motion, the inverse kinematic equation with velocity compensation is For the four-Mecanum-wheeled robot in Figure 2, according to Equations (8) and (19), the inverse kinematic equation with velocity compensation is

Motion Simulation of Four-Mecanum-Wheeled Robot Virtual Prototype and Motion Test of Physical Robot
In order to further verify the correctness of the above theoretical research, the Mecanum-wheeled robot was simulated and analyzed by virtual prototyping technology. In this paper, the virtual prototype of the Mecanum-wheel mobile platform was modelled in SolidWorks, and the simulation was carried out on RecurDyn. Before importing into RecurDyn, the robot model should be simplified. The simplified assembly model built in SolidWorks should be imported into RecurDyn in Parasolid format. In the process of modeling, a single Mecanum wheel is defined as a subsystem in RecurDyn, which not only facilitates the hierarchical management of the model, but also enables the subsystem to be established, modified, imported, and exported independently. That can improve the reuse rate of the model, reduce the chance of model construction errors, and improve the efficiency.
According to the actual movement of the robot platform, the constraint relationship between components is created based on determining the assembly and motion relationships of components. In the robot virtual prototype simulation model, 36 revolute joints are created, as shown in Figure 5a. There are four revolute joints between the hub of the Mecanum wheel and the main body, and eight revolute joints between the rollers and the hub in a single Mecanum wheel subsystem, as shown in Figure 5b. The following are the basic dimensions of the robot model: γ = 45 • , W = 370 mm, H = 350 mm, and R = 152.4 mm. The weight of the robot is 200 kg. By setting the angular velocity motion functions on the revolute joints between the four Mecanum wheels and the main body, the robot platform can move in simulation at desired speed. In order to prevent system anomalies caused by sudden changes in speed during the simulation, the drive function is usually applied to the revolute joints using the system's predefined STEP function. The STEP function is a third-order polynomial function built into the software and is used to define a relatively smooth load curve, which is suitable for smoothly loading drivers for mobile robot platforms. motion functions on the revolute joints between the four Mecanum wheels and the main body, the robot platform can move in simulation at desired speed. In order to prevent system anomalies caused by sudden changes in speed during the simulation, the drive function is usually applied to the revolute joints using the system's predefined STEP function. The STEP function is a third-order polynomial function built into the software and is used to define a relatively smooth load curve, which is suitable for smoothly loading drivers for mobile robot platforms.

Creating Contacts and Setting Parameters
Only by adding contacts between the rollers and the ground and setting reasonable contact parameters can the virtual prototype be simulated in accordance with the predetermined scheme. The RecurDyn contact toolkit provides rich contact models, from which Solid Contact was selected. Usually, the outer layer of the roller is vulcanized with a layer of rubber or polyurethane material. The higher the surface hardness of the roller, the stronger the bearing capacity and the higher the elastic modulus. The surface hardness of rollers with Mecanum wheels for small mobile platforms is usually 75 A Shaw hardness, and 90 A for the rollers of heavy-duty AGV. Moreover, we assume that the workplace has C30 grade concrete ground. The estimated values of the attribute parameters of the roller and the ground are shown in Table 1.

Creating Contacts and Setting Parameters
Only by adding contacts between the rollers and the ground and setting reasonable contact parameters can the virtual prototype be simulated in accordance with the predetermined scheme. The RecurDyn contact toolkit provides rich contact models, from which Solid Contact was selected. Usually, the outer layer of the roller is vulcanized with a layer of rubber or polyurethane material. The higher the surface hardness of the roller, the stronger the bearing capacity and the higher the elastic modulus. The surface hardness of rollers with Mecanum wheels for small mobile platforms is usually 75 A Shaw hardness, and 90 A for the rollers of heavy-duty AGV. Moreover, we assume that the workplace has C30 grade concrete ground. The estimated values of the attribute parameters of the roller and the ground are shown in Table 1. Based on Hertz contact theory, the stiffness coefficient between two objects can be expressed as: where ρ is the comprehensive radius of curvature, 1 ρ = 1 ρ 1 + 1 ρ 2 ; and ρ 1 and ρ 2 are the radii of curvature of the two objects at the contact; and E * is the comprehensive elastic modulus; 1 , E 2 and υ 1 , υ 2 are the elastic modulus and Poisson's ratio of the two materials, respectively.
The material property parameters of the ground and the rollers in Table 1, the radius (ρ 1 = 10 10 mm) of curvature of the horizontal ground, and the radius (ρ 2 = 33.5 mm) of curvature of the roller are substituted into Equation (23), and the stiffness coefficient between the two objects can be calculated. According to the theoretical calculation and feedback of simulation results, the contact parameters are set as follows: spring coefficient is 1200 N/mm, damping coefficient is 4, dynamic friction coefficient is 0.35, stiffness exponent is 2, and the rest are kept as default. On each roller of the Mecanum wheel, solid contact with the ground is set. Each Mecanum wheel is provided with 8 solid contacts on the ground, so 32 solid contacts are established for the mobile platform. The virtual prototype simulation model established through the above steps is shown in Figure 5a.

Motion Simulation and Analysis of a Single Four-Mecanum-Wheeled Robot
The motion simulations of the four-Mecanum-wheeled robot platform were carried out using 11 modes to analyze the motion characteristics. In order to better show the motion state of the robot platform during the simulation process, the trajectory tool in the software was used to display the motion track of the robot center or the specified marker point. The motion situation and trajectories of the robot in the simulation using 11 motion modes are shown in Figure 6. The simulation trajectories of the center of the robot platform in longitudinal, lateral, and oblique motion at 45 • are shown in Figure 6a-c, respectively.  Table 2.  Figure 7 shows the load curves of a wheel and a roller on the wheel during longitudinal movement of the mobile platform. By observing these two load curves, it can be seen that both the wheel and the roller bear periodic loads. The periodic load on the roller is derived from the roller being in contact with the ground every revolution of the wheel, and the contact duration is short and the load changes drastically. When the robot platform is in a stable motion state, the periodicity of  Table 2. Trajectory equations of robot moving in circular, figure-8-shaped, and harmonic curves.

Names of Curve
Trajectory Equations in Figure 6 Circular curve x = r i sin ω i t, y = r i cos ω i t, θ = 0, where r 1 = 2500 mm, r 2 = r 3 = 1350 mm, ω 1 = π/10, ω 2 = ω 3 = π/5 Figure Figure 7 shows the load curves of a wheel and a roller on the wheel during longitudinal movement of the mobile platform. By observing these two load curves, it can be seen that both the wheel and the roller bear periodic loads. The periodic load on the roller is derived from the roller being in contact with the ground every revolution of the wheel, and the contact duration is short and the load changes drastically. When the robot platform is in a stable motion state, the periodicity of the wheel load is mainly affected by the change of the state of contact between the roller and the ground on the body. In addition, the load on the roller is significantly greater than the load on the wheel, since the load on the entire robot is all applied to the roller in contact with the ground. The longitudinal and lateral motions of the robot are simulated with a series of different velocities, and the average velocities of the robot in stable state are obtained, as shown in Table 3. These simulation velocities are less than the corresponding set velocities. From the simulation data, the relative error of longitudinal velocity is very small; the average relative error is about −0.22%, so its correlation with velocity is not obvious. Five sets of simulation results show that the velocity loss of lateral translation is significantly greater than that of longitudinal translation when the robot runs at the same set velocity. Moreover, there is a negative correlation between the error of the robot and the set velocity in the lateral translation motion. The relative error curves of the velocities of the robot in the two motion modes are shown in Figure 8.         Figure 9a shows the displacement curve of the robot in the x-direction in lateral motion. Figure  9c shows the x-direction velocity curve of the robot platform when it moves laterally at a set velocity of 500 mm/s. In the figure, the velocity value rises gradually from zero, which reflects the process of accelerating the robot platform from static state to stable moving state, and then the value fluctuates steadily in a range. In the simulation of lateral movement at a set speed of 500 mm/s, the maximum  Figure 9a shows the displacement curve of the robot in the x-direction in lateral motion. Figure 9c shows the x-direction velocity curve of the robot platform when it moves laterally at a set velocity of 500 mm/s. In the figure, the velocity value rises gradually from zero, which reflects the process of accelerating the robot platform from static state to stable moving state, and then the value fluctuates steadily in a range. In the simulation of lateral movement at a set speed of 500 mm/s, the maximum lateral velocity is 497.5 mm/s, the minimum velocity is 460.4 mm/s, and the average velocity is 483.2 mm/s. Figure 9b shows the longitudinal displacement and velocity curves in the y-direction during lateral motion. When the robot moves laterally in the x-direction, it fluctuates slightly in the y-direction and deviates from the x-axis direction. The velocity of the robot along the y-axis fluctuates around the zero line, as shown in Figure 9d. In 10 seconds of lateral movement, that is, when the robot moves about 5000 mm in the lateral direction, it deviates from the x-axis direction by about 4 mm. So, the longitudinal offset has little effect on the lateral motion of the robot platform.
(2) Simulation of Oblique Motion at 45 • When the robot moves in a straight line with an oblique direction of 45 • , only two Mecanum wheels on one diagonal line rotate; the wheels on the other diagonal line do not rotate, but the bottom rollers on the wheels that do not turn passively rotate. Table 4 shows the simulation velocities and relative errors of the 45 • oblique translation of the robot. According to the simulation data, there is loss of velocity in motion, and the velocity loss in the x-and y-directions is almost equal. The relative error of velocity in the x-direction is smaller than that of the lateral translation set at the same velocity. lateral motion. When the robot moves laterally in the x-direction, it fluctuates slightly in the ydirection and deviates from the x-axis direction. The velocity of the robot along the y-axis fluctuates around the zero line, as shown in Figure 9d. In 10 seconds of lateral movement, that is, when the robot moves about 5000 mm in the lateral direction, it deviates from the x-axis direction by about 4 mm. So, the longitudinal offset has little effect on the lateral motion of the robot platform. (2) Simulation of Oblique Motion at 45° When the robot moves in a straight line with an oblique direction of 45°, only two Mecanum wheels on one diagonal line rotate; the wheels on the other diagonal line do not rotate, but the bottom rollers on the wheels that do not turn passively rotate. Table 4 shows the simulation velocities and relative errors of the 45° oblique translation of the robot. According to the simulation data, there is loss of velocity in motion, and the velocity loss in the x-and y-directions is almost equal. The relative error of velocity in the x-direction is smaller than that of the lateral translation set at the same velocity.  Figure 10a shows that the velocity of the x-and y-axis increases rapidly from 0 to a predetermined value within 0.6 s when the mobile platform moves along the oblique direction of 45°. Then the velocities in the two directions are basically the same, and fluctuate within a certain range.   Figure 10a shows that the velocity of the x-and y-axis increases rapidly from 0 to a predetermined value within 0.6 s when the mobile platform moves along the oblique direction of 45 • . Then the velocities in the two directions are basically the same, and fluctuate within a certain range. In the simulation, the velocity in the x-and y-axis directions is set at 500 mm/s, and the average test velocities are 490.164 mm/s and 490.061 mm/s, respectively. The curves of the two velocities are symmetrical relative to the zero line. The combination of these two directions enables the robot platform to move along the 45 • direction. The trajectory curve of the robot in the oblique 45 • motion is shown in Figure 10b.
In the simulation, the velocity in the x-and y-axis directions is set at 500 mm/s, and the average test velocities are 490.164 mm/s and 490.061 mm/s, respectively. The curves of the two velocities are symmetrical relative to the zero line. The combination of these two directions enables the robot platform to move along the 45° direction. The trajectory curve of the robot in the oblique 45° motion is shown in Figure 10b. (3) Turning on the Spot Table 5 shows the relative errors of the angular velocity of the robot in the turning on the spot motion simulation at different angular velocities, which are positively correlated with increased velocity. The trajectory of the robot center in this simulation in five cycles is shown in Figure 11a. The trajectory has great randomness. In the first cycle, the displacement of the robot center in both x-and y-axis directions is within 1 mm. However, with the increased number of rotation cycles, the displacement tends to increase. After five cycles of rotation, the maximum displacement of the robot has exceeded 4 mm. The displacement curves of the robot center in the x-and y-axis directions in five cycles are shown in Figure 11b. It can be seen that the offset of the center coordinate of the robot increases with the increased number of rotation cycles.  (3) Turning on the Spot Table 5 shows the relative errors of the angular velocity of the robot in the turning on the spot motion simulation at different angular velocities, which are positively correlated with increased velocity. The trajectory of the robot center in this simulation in five cycles is shown in Figure 11a. The trajectory has great randomness. In the first cycle, the displacement of the robot center in both x-and y-axis directions is within 1 mm. However, with the increased number of rotation cycles, the displacement tends to increase. After five cycles of rotation, the maximum displacement of the robot has exceeded 4 mm. The displacement curves of the robot center in the x-and y-axis directions in five cycles are shown in Figure 11b. It can be seen that the offset of the center coordinate of the robot increases with the increased number of rotation cycles.  According to the errors between simulation results and theoretical values of the velocities shown in Tables 3-5, the simulation velocities of the robot are less than the theoretical velocities. In these motion modes, the relative errors of simulation velocity vary slightly with different velocities but have a positive correlation on the whole, especially in lateral translation, oblique 45° translation, and in situ According to the errors between simulation results and theoretical values of the velocities shown in Tables 3-5, the simulation velocities of the robot are less than the theoretical velocities. In these motion modes, the relative errors of simulation velocity vary slightly with different velocities but have a positive correlation on the whole, especially in lateral translation, oblique 45 • translation, and in situ rotation. Comparing these three kinds of motion simulation, the velocity error of lateral motion is the largest, and that of longitudinal motion is the smallest. The relative error of velocity of 45 • oblique motion is smaller than that of lateral motion and more than that of longitudinal motion. The relative error of angular velocity of in situ rotation is equivalent to that of 45 • oblique translation motion.
Possible reasons for the errors between the simulation and theoretical set velocity are as follows: First, the wheel skidding phenomenon occurs in the process of rotation, which is common in wheeled mobile robots; slipping is most serious in the process of lateral movement. Second, the rim of the multiple rollers of the Mecanum wheel is theoretically a circle; however, in the simulation movement, deformation of the rubber layer at different positions of the rollers under load is different, and the effective radius of the wheel changes with rotation, thereby causing periodic changes even if the angular velocity is the same. Third, the roller cuts into the ground under load, which causes the radius of the wheel to be less than the theoretical radius, which is similar to the compression of the elastic roller in actual motion.
It can be seen that there are errors between the robot's velocity in motion simulation and the set velocity, especially in lateral translation motion. In order to ensure precise motion of the robot along the desired trajectory, it is necessary to compensate its velocity.

Velocity Compensation in Motion Simulation of Complex Trajectories
According to the simulation results of simple motion in the previous section, it can be seen that there are velocity losses in longitudinal translation, lateral translation, and rotational motion. When the robot translates along a complex trajectory, such as circular, 8-shaped, and simple harmonic motion curves, its motion can be decomposed into motion in both x-and y-directions. When the robot performs simulated motion along complex trajectories, if a theoretical velocity value is assigned to each Mecanum wheel, the robot will not be able to move to the expected position, so there is an error between the simulated and theoretical trajectory. The trajectories of the robot in the motion simulation along circular, 8-shaped, and harmonic curves are plotted with blue lines in Figure 12. These simulation curves are basically consistent with the theoretical curves, plotted with dashed red lines in Figure 12, but because of the velocity error of the robot, there are deviations between the simulation and theoretical curves. The circular and 8-shaped simulation curves are closed, so the curves deviate inward. The displacement errors of circular curves in Figure 12a in the x-and y-axis directions are 3.48% and 0.14%, respectively. The displacement errors of the 8-shaped trajectory in Figure 12b in the xand y-axis are 2.92% and 0.28%, respectively. Relative errors of displacement of x-and y-axis directions in motion simulation along the circular and 8-shaped curves are shown in Table 6. With the increased motion period, the deviation of the simple harmonic trajectory in the x-direction increases, as shown in Figure 12c. If the robot keeps a certain tilt angle during the translational motion, there is a large displacement error in the trajectory in the lateral direction. Figure 12d,e show the trajectories of the robot in motion simulation shown in Figure 6h,i, maintaining a 60 • angle during translating motion. According to Figure 12d,e, the deviation between the robot's trajectory and the theoretical curve in the oblique direction is the largest. The oblique direction is the x-direction in the local coordinate system of the robot. Sensors 2020, 20, x FOR PEER REVIEW 18 of 37

(c)
Trajectory after velocity compensation X (mm) Figure 12. Trajectories of robot in motion simulation shown in Figure 6: (a) circular trajectory in Figure 6e; (b) 8-shaped trajectory in Figure 6f; (c) simple harmonic motion trajectory in Figure 6g; (d) circular trajectory in Figure 6h; (e) 8-shaped trajectory in Figure 6i; (f) centripetal circular motion in Figure 6j; (g) centripetal motion of 45 • angle in Figure 6l. Compensating the robot's velocity can reduce the trajectory error and make the trajectory closer to the theoretical trajectory. In the translational motion of the robot, only the velocity in the x-and y-axis directions needs to be compensated. According to the simulation data in Table 3, as the velocity changes, the relative error of the robot's velocity changes, but the amplitude is small. When the robot makes a curve motion, the velocities of the robot and the wheels change with the moment of change at all times. So a fixed compensation coefficient can be set for the virtual prototype model of the robot. In this simulation test, the average values of the velocity compensation coefficients of the longitudinal and lateral translation motions in Table 3 are used as the compensation coefficients in the x-and y-axis directions. The motion shown in Figure 6e-i is translation motion; according to Table 3, the velocity compensation matrix is set to C 1 . The centripetal circular motion shown in Figure 6j-m is a combination of translation and rotation; therefore, it is also necessary to compensate for its rotational motion. According to Table 5, the compensation coefficient of rotational angular velocity is set to 1.0195, and the compensation coefficient of lateral velocity is set to 1.04 due to its large simulated lateral velocity. Then, the compensation matrix of velocity is set to C 2 .
The trajectories of the robot in the motion simulation along the circular, 8-shaped, and simple harmonic motion curves after the velocity compensation are plotted with solid black lines in Figure 12. It can be seen from Figure 12a-e that the trajectories after velocity compensation are very consistent with the theoretical curves. The trajectories of the centripetal circular motions in Figure 6j-m are shown in Figure 12f,g, and the simulation trajectories after velocity compensation are better than those without velocity compensation. The simulation results show that this velocity compensation method is effective.

Motion Test System of Mecanum-Wheeled Robot Using Optitrack Optical Motion Capture System
In Section 3, RecurDyn software was used to simulate the robot's motion, and the motion simulations verified the correctness of the kinematic model of the four-Mecanum-wheeled robot and the feasibility of velocity compensation. However, they were carried out to simulate an ideal robot on ideal ground, which is quite different from real situations. The movement of a robot in a real environment is affected by many factors, including machining error of the Mecanum wheels, angle error of rollers, installation error of the four wheels, flatness of the ground, and control performance of the motor. In order to verify the correctness of the motion and motion compensation models, motion tests of the robot physical prototype were carried out. In order to measure the coordinates and trajectories of the robot during movement, the motion test system was constructed using the OptiTrack optical motion capture system, as shown in Figure 13 [48]. Three OptiTrack Prime 13 cameras (high-speed motion capture cameras), shown in Figure 13a, were arranged on each side of the test area. The cameras were connected to the data and power supply by using a Gigabit Ethernet GigE/PoE interface. All cameras were connected to a Gigabit network hub with Ethernet cables. An installed workstation with Motive optical motion capture was connect to the hub with a cable. The Motive software was used for recording, presentation, playback, and remote data services of the position data. The test used a four-Mecanum-wheeled robot prototype, shown in Figure 13b. The distance between the center lines of the front and rear wheels is 400 mm, and of the left and right wheels is 450 mm. The robot is controlled by another computer with a human-computer interaction (HCI) system. A Hand Rigid Bodies Marker Set was fixed on the robot prototype to test its space coordinates in the test space. error of rollers, installation error of the four wheels, flatness of the ground, and control performance of the motor. In order to verify the correctness of the motion and motion compensation models, motion tests of the robot physical prototype were carried out. In order to measure the coordinates and trajectories of the robot during movement, the motion test system was constructed using the OptiTrack optical motion capture system, as shown in Figure 13 [48]. Three OptiTrack Prime 13 cameras (high-speed motion capture cameras), shown in Figure 13a, were arranged on each side of the test area. The cameras were connected to the data and power supply by using a Gigabit Ethernet GigE/PoE interface. All cameras were connected to a Gigabit network hub with Ethernet cables. An installed workstation with Motive optical motion capture was connect to the hub with a cable. The Motive software was used for recording, presentation, playback, and remote data services of the position data. The test used a four-Mecanum-wheeled robot prototype, shown in Figure 13b. The distance between the center lines of the front and rear wheels is 400 mm, and of the left and right wheels is 450 mm. The robot is controlled by another computer with a human-computer interaction (HCI) system. A Hand Rigid Bodies Marker Set was fixed on the robot prototype to test its space coordinates in the test space.

Motion Shots of the Robot in Test
During the testing process, the Mecanum-wheeled robot was controlled to move in these motion modes: longitudinal translation, lateral translation, oblique 45° translation, turning on the spot, circular and 8-shaped curves, and translation motion along a simple harmonic curve. The OptiTrack optical motion capture system captured and recorded the spatial coordinates of the robot. In order to visualize the movement process, the Motion Shot app was used to create a motion picture of a series of coherent images of the robot in the same test environment. The motion picture recorded the position of the robot during movement and showed its movement track. Figure 14 shows six motion

Motion Shots of the Robot in Test
During the testing process, the Mecanum-wheeled robot was controlled to move in these motion modes: longitudinal translation, lateral translation, oblique 45 • translation, turning on the spot, circular and 8-shaped curves, and translation motion along a simple harmonic curve. The OptiTrack optical motion capture system captured and recorded the spatial coordinates of the robot. In order to visualize the movement process, the Motion Shot app was used to create a motion picture of a series of coherent images of the robot in the same test environment. The motion picture recorded the position of the robot during movement and showed its movement track. Figure 14 shows six motion pictures of the robot in motion, including longitudinal translation (Figure 14a (Figure 14g,h), and translation along a simple harmonic curve (Figure 14d). The arrows in Figure 14 indicate the moving direction of the robot.  (Figure 14g,h), and translation along a simple harmonic curve (Figure 14d). The arrows in Figure 14 indicate the moving direction of the robot.

Analysis and Motion Tests of Robot in Four Simple Motions
(1) Longitudinal and Lateral Motion Tests Regardless of the longitudinal or lateral translational motion, there is a certain velocity error of the robot, so the actual velocity is less than the set value. The test velocities and velocity errors are shown in Table 7. When the moving velocity of the robot is set to 100 mm/s, the longitudinal velocity error is 2.18% and the lateral translation velocity error is 12.38%. It can be seen that the lateral movement velocity error is much larger than the longitudinal. The lateral movement test was performed at different set velocities. It can be seen from the test results that the velocity error increased as the set Sensors 2020, 20, 75 21 of 37 velocity value increased, as shown in Table 7. The loss of velocity during the movement will have a great impact on the trajectory of the robot. Figure 15 shows the trajectory, displacement, and velocity curves of the robot in lateral translation motion test at a set velocity of 500 mm/s.  Table 7. When the moving velocity of the robot is set to 100 mm/s, the longitudinal velocity error is 2.18% and the lateral translation velocity error is 12.38%. It can be seen that the lateral movement velocity error is much larger than the longitudinal. The lateral movement test was performed at different set velocities. It can be seen from the test results that the velocity error increased as the set velocity value increased, as shown in Table 7. The loss of velocity during the movement will have a great impact on the trajectory of the robot. Figure 15 shows the trajectory, displacement, and velocity curves of the robot in lateral translation motion test at a set velocity of 500 mm/s.    Figure 15a shows the trajectory of the robot in the lateral translation test at a set velocity of 500 mm/s. The robot moves laterally 3800 mm in the x-axis direction and offset 100 mm in the y-axis direction. The first section of the x-direction velocity curve is an acceleration process; when the set value is reached, the velocity fluctuates at a certain value, as shown in Figure 15d. The lateral displacement curve in the x-direction is shown in Figure 15b. The test displacement and velocity curve of the robot in the y-direction in lateral translation motion is shown in Figure 15c,e, respectively. In the process of lateral translation, in the initial stage, the robot slides to one side, then changes direction, and the sliding displacement accumulates continuously. The velocity curve of the robot along the y-axis fluctuates around the zero line, resulting in a small fluctuation in the trajectory, shown in Figure 15a.
(2) Test of Oblique Motion at 45 • Figure 16a shows the test trajectory of the robot in oblique 45 • translation, and Figure 16b,c show the displacement and velocity curves of the robot in the x-and y-directions in oblique 45 • translation. In the first half of the 45 • oblique motion, the velocity in the y-axis direction is faster than that in the x-axis direction. In the second half of the motion, the velocity in the y-axis direction decreases gradually and is less than that in the x-axis direction. The trajectory of the robot deviates from the theoretical trajectory line in the positive direction of the y-axis, and the second half of the robot gradually approaches the theoretical trajectory line. The trajectory is not a straight line, and there is a certain radian. The maximum distance of the trajectory deviating from the theoretical straight line is 207.7 mm. In this test, the set velocity of the robot is 500 mm/s and the theoretical velocity of the x-and y-directions should be 353.6 mm/s. The average test velocities of the x-and y-axes in the steady state are 319.4 mm/s and 330.7 mm/s, respectively, and the relative errors are −9.67% and −6.47%, respectively.
Sensors 2020, 20, x FOR PEER REVIEW 22 of 37 direction. The first section of the x-direction velocity curve is an acceleration process; when the set value is reached, the velocity fluctuates at a certain value, as shown in Figure 15d. The lateral displacement curve in the x-direction is shown in Figure 15b. The test displacement and velocity curve of the robot in the y-direction in lateral translation motion is shown in Figures 15c,e, respectively. In the process of lateral translation, in the initial stage, the robot slides to one side, then changes direction, and the sliding displacement accumulates continuously. The velocity curve of the robot along the y-axis fluctuates around the zero line, resulting in a small fluctuation in the trajectory, shown in Figure 15a.
(2) Test of Oblique Motion at 45° Figure 16a shows the test trajectory of the robot in oblique 45° translation, and Figure 16b,c show the displacement and velocity curves of the robot in the x-and y-directions in oblique 45° translation. In the first half of the 45° oblique motion, the velocity in the y-axis direction is faster than that in the x-axis direction. In the second half of the motion, the velocity in the y-axis direction decreases gradually and is less than that in the x-axis direction. The trajectory of the robot deviates from the theoretical trajectory line in the positive direction of the y-axis, and the second half of the robot gradually approaches the theoretical trajectory line. The trajectory is not a straight line, and there is a certain radian. The maximum distance of the trajectory deviating from the theoretical straight line is 207.7 mm. In this test, the set velocity of the robot is 500 mm/s and the theoretical velocity of the xand y-directions should be 353.6 mm/s. The average test velocities of the x-and y-axes in the steady state are 319.4 mm/s and 330.7 mm/s, respectively, and the relative errors are -9.67% and -6.47%, respectively. (3) Test of Motion of Turning on the Spot Figure 17 shows the trajectories of points A and B on the robot in turning on the spot, where point A is very close to the center and point B is about 200 mm from the center. Figure 17a shows the trajectories of points A and B in the first cycle, and the trajectory of point B in five cycles. Figure 17b is an enlarged view of the trajectory of point A in the first cycle in Figure 17a. While the robot rotates around one circle, point B draws a perfect circular curve. There is local tortuosity in the circular curve, which indicates that the position of the geometric center of the robot slips during its rotation. The trajectories of point B of multiple rotation cycles do not coincide, and they constitute a trajectory band. By marking and drawing the trajectory of the center point of the robot, we can observe the change of the center point more directly in the course of rotation, but it is difficult to mark accurately. Marker point A can only be located at the geometric center of the robot as far as possible. Marker point B drifts in a small range, and the trajectory of point B in one cycle is shown in Figure 17b, which indicates the displacement variation of the center point of the robot in the process of rotation, and also reflects the trend of circular trajectory caused by the small deviation from the center point (about 3 mm). Table 8 shows the loss and relative errors of the angular velocity of the robot during turning on the spot. This angular velocity error cannot be ignored in precise motion control.
Sensors 2020, 20, x FOR PEER REVIEW 23 of 37 band. By marking and drawing the trajectory of the center point of the robot, we can observe the change of the center point more directly in the course of rotation, but it is difficult to mark accurately. Marker point A can only be located at the geometric center of the robot as far as possible. Marker point B drifts in a small range, and the trajectory of point B in one cycle is shown in Figure 17b, which indicates the displacement variation of the center point of the robot in the process of rotation, and also reflects the trend of circular trajectory caused by the small deviation from the center point (about 3 mm). Table 8 shows the loss and relative errors of the angular velocity of the robot during turning on the spot. This angular velocity error cannot be ignored in precise motion control.   According to Table 7, the compensation coefficients of longitudinal and lateral velocity are set to 1.04 and 1.145, respectively. Then, the compensation matrix of velocity is set to C3.   Figure 18 shows trajectories of translation motion along circle, 8-shaped, and simple harmonic curves. According to the curve equations in Table 2, in the tests, the parameters of the circle curve are r = 1000 mm, ω = π/10; the parameters of the 8-shaped curve are l = 1000 mm and ω = π/10; and the parameters of the simple harmonic curve are a = 300 mm, b = 1000 mm and ω = π/5. According to Table 7, the compensation coefficients of longitudinal and lateral velocity are set to 1.04 and 1.145, respectively. Then, the compensation matrix of velocity is set to C 3 .

Velocity Compensation in Motion Test of Complex Trajectories
Sensors 2020, 20, x FOR PEER REVIEW 24 of 37 Figure 18. Trajectories, displacements, and velocities, respectively, of the robot in translation motion along (a-c) a circle, (d-f) an 8-shaped curve, and (g-i) a simple harmonic curve. DX or DY, displacement in x-or y-direction; DXC or DYC, displacement in x-or ydirection with velocity compensation; VX or VY, velocity in x-or y-direction; VXC or VYC, velocity in x-or y-direction with velocity compensation.
The blue curve in Figure 18a is the test trajectory of the robot in translation motion along a circle, which shifts to the inside of the theoretical circle (dashed red line) and is a vertical approximate  Figure 18. Trajectories, displacements, and velocities, respectively, of the robot in translation motion along (a-c) a circle, (d-f) an 8-shaped curve, and (g-i) a simple harmonic curve. DX or DY, displacement in x-or y-direction; DXC or DYC, displacement in x-or y-direction with velocity compensation; VX or VY, velocity in x-or y-direction; VXC or VYC, velocity in x-or y-direction with velocity compensation.
The blue curve in Figure 18a is the test trajectory of the robot in translation motion along a circle, which shifts to the inside of the theoretical circle (dashed red line) and is a vertical approximate ellipse. During the movement, the robot slips or its heading angle changes, causing the trajectory to fail to coincide. After motion compensation, the trajectory of the robot (solid black line) is very consistent with the theoretical curve, as shown in Figure 18a. From the velocity comparison curves (Figure 18c) and displacement comparison curves (Figure 18b) in the x-and y-directions, we can see that velocity and displacement in the motion test with velocity compensation at the same time are greater than the uncompensated velocity and displacement. In the test of translation motion along an 8-shaped curve, the test motion trajectory also deviates greatly from the theoretical trajectory. In the lateral direction, the displacement of the robot is obviously less than the theoretical value, as shown in Figure 18d. In the process of motion, the robot slips and the heading angle changes. Therefore, the test trajectory of the robot deflects. In the first motion cycle, the initial point and termination point of the robot cannot coincide. The displacements and velocities in translation motion along an 8-shaped curve are shown in Figure 18d,f, respectively. In the test of translational motion of a simple harmonic curve, the error of longitudinal displacement (x-axis direction) is small, while that of transverse displacement (y-axis direction) is large. With the increased motion cycles of the robot, the deviation between transverse displacement and theoretical trajectory curve increases, as shown in Figure 18g. This is consistent with the trend of dynamic simulation. The robot maintains a uniform velocity in the x-direction, and the speed in the x-direction with compensation (solid blue line) is slightly larger, as shown in Figure 18i; the displacement curve in the x-direction with motion compensation is at the top, as shown in Figure 18h. In the y-direction, the velocity and displacement of the robot change periodically, and are compensated accordingly. Table 9 shows the relative errors of displacement in motion simulation along a circle, 8-shaped curve, and simple harmonic curve in Figure 18. From Table 9, improvement of the robot's displacement accuracy using the kinematic model with velocity compensation can be observed more clearly. From the results of the motion compensation test, we can get the following conclusions: (1) in the real environment, the wheel slippage of a Mecanum-wheeled robot has a great impact on the robot's motion accuracy, especially lateral motion; (2) for motion with small velocity change, the constant velocity compensation matrix can have a good compensation effect, and the actual trajectories of the robot are very consistent with the theoretical curves; (3) more accurate motion compensation requires the use of a variable velocity compensation matrix, which is dynamically set according to the robot's velocity; (4) in the process of motion, the robot will deviate from its heading, which will cause its trajectory to deviate; and (5) according to the robot's parameters and environment information detected by position, attitude, and speed sensors, the motion direction and velocity can be adjusted in real time by using this kinematic model with speed compensation, which is more conducive to the robot's motion accuracy.

Motion Simulation of Combined Configurations of Multiple Four-Mecanum-Wheeled Robots
In Section 3, the typical motions of the four-Mecanum-wheeled robot were simulated. In Section 4, motion tests were performed using the robot physical prototype. Although the geometry of the virtual prototype used for motion simulation and the physical prototype used for testing are different, the motion test results verify the feasibility and effectiveness of motion simulation of Mecanum-wheeled robots by RecurDyn software. Therefore, the motion model of combined configurations of multiple-Mecanum-wheeled robots can also be verified by this simulation method, which can also make up for the shortcomings that we only have one physical robot sample and cannot achieve a collaborative combined multiple-robot experiment. In this section, kinematic simulations of a variety of configurations are performed to verify the kinematic model of multiple robots. According to the motion simulations in Section 3, the relative error between the simulation result and the theoretical value is small. Even without motion compensation, the simulation trajectory of the robot is very close to the theoretical curve. Therefore, in this section, in order to carry out simulations efficiently, simulations without motion compensation are first carried out for various combination configurations, and then two robot combination configurations are selected for simulation with motion compensation.

Multiple Four-Mecanum-Wheeled Robot Combination Configurations for Simulation
In this section, seven multiple robot combination configurations are used for motion simulation. The configurations and their main dimensions are shown in Figure 19. Figure 19a,b show two four-Mecanum-wheeled robots arranged laterally. The center distances of the two robots are the same; the two robots in Figure 19a are connected in the middle, called side-by-side connected configuration of two robots (SCC-2), and the two in Figure 19b are separate, called side-by-side unconnected configuration of two robots (SUC-2). Figure 19c,d show longitudinally arranged dual four-Mecanum-wheeled robots with the same center-to-center distance. The two robots in Figure 19c are connected, called tandem connected configuration of two robots (TCC-2), and the two in Figure 19d are separate, called tandem unconnected configuration of two robots (TUC-2). Figure 19e shows a tandem connected configuration of three robots (TCC-3). Figure 19f shows an arbitrary unconnected configuration of two robots arranged at an angle of 30 • (AUC-2). If the two robots are connected, this is an arbitrary connected configuration (ACC-2). Mecanum-wheeled robots by RecurDyn software. Therefore, the motion model of combined configurations of multiple-Mecanum-wheeled robots can also be verified by this simulation method, which can also make up for the shortcomings that we only have one physical robot sample and cannot achieve a collaborative combined multiple-robot experiment. In this section, kinematic simulations of a variety of configurations are performed to verify the kinematic model of multiple robots. According to the motion simulations in Section 3, the relative error between the simulation result and the theoretical value is small. Even without motion compensation, the simulation trajectory of the robot is very close to the theoretical curve. Therefore, in this section, in order to carry out simulations efficiently, simulations without motion compensation are first carried out for various combination configurations, and then two robot combination configurations are selected for simulation with motion compensation.

Multiple Four-Mecanum-Wheeled Robot Combination Configurations for Simulation
In this section, seven multiple robot combination configurations are used for motion simulation. The configurations and their main dimensions are shown in Figure 19. Figure 19a,b show two four-Mecanum-wheeled robots arranged laterally. The center distances of the two robots are the same; the two robots in Figure 19a are connected in the middle, called side-by-side connected configuration of two robots (SCC-2), and the two in Figure 19b are separate, called side-by-side unconnected configuration of two robots (SUC-2). Figure 19c,d show longitudinally arranged dual four-Mecanumwheeled robots with the same center-to-center distance. The two robots in Figure 19c are connected, called tandem connected configuration of two robots (TCC-2), and the two in Figure 19d are separate, called tandem unconnected configuration of two robots (TUC-2). Figure 19e shows a tandem connected configuration of three robots (TCC-3). Figure 19f shows an arbitrary unconnected configuration of two robots arranged at an angle of 30° (AUC-2). If the two robots are connected, this is an arbitrary connected configuration (ACC-2).

Simulation of Turning on Spot of Multiple Four-Mecanum-Wheeled Robots
In situ rotation motion simulation of the five configurations of robots was carried out, as shown in Figure 20. Figure Figure 20c. We can see that after the robot rotates in place for many cycles, the trajectory of the marker point on the robot forms a circular trajectory range. In comparison, the trajectories of two unconnected robots (SUC-2) are relatively tortuous, as shown in Figure 21b. This is because the robots randomly move in the plane during in situ rotation. Figure 22 shows the trajectories of the center points of the combined robots and displacement curves in both x-and y-directions during the rotation. The center trajectory coverage of the unconnected robots is larger than that of the connected robots, as shown in Figure 22a,b, and their centers also have larger displacements in the x-or y-direction, as shown in Figure 22c

Simulation of Turning on Spot of Multiple Four-Mecanum-Wheeled Robots
In situ rotation motion simulation of the five configurations of robots was carried out, as shown in Figure 20. Figure Figure 20c. We can see that after the robot rotates in place for many cycles, the trajectory of the marker point on the robot forms a circular trajectory range. In comparison, the trajectories of two unconnected robots (SUC-2) are relatively tortuous, as shown in Figure 21b. This is because the robots randomly move in the plane during in situ rotation. Figure 22 shows the trajectories of the center points of the combined robots and displacement curves in both x-and y-directions during the rotation. The center trajectory coverage of the unconnected robots is larger than that of the connected robots, as shown in Figure  22a,b, and their centers also have larger displacements in the x-or y-direction, as shown in Figure  22c,d. With rotation, there are certain periodic variations in displacement, and displacement has an increasing tendency, as shown in Figure 22c,d.

Simulation of Turning on Spot of Multiple Four-Mecanum-Wheeled Robots
In situ rotation motion simulation of the five configurations of robots was carried out, as shown in Figure 20. Figure Figure 20c. We can see that after the robot rotates in place for many cycles, the trajectory of the marker point on the robot forms a circular trajectory range. In comparison, the trajectories of two unconnected robots (SUC-2) are relatively tortuous, as shown in Figure 21b. This is because the robots randomly move in the plane during in situ rotation. Figure 22 shows the trajectories of the center points of the combined robots and displacement curves in both x-and y-directions during the rotation. The center trajectory coverage of the unconnected robots is larger than that of the connected robots, as shown in Figure  22a,b, and their centers also have larger displacements in the x-or y-direction, as shown in Figure  22c,d. With rotation, there are certain periodic variations in displacement, and displacement has an increasing tendency, as shown in Figure 22c,d.

Translation Simulation along a Circle of Multiple Four-Mecanum-Wheeled Robots
Three robot configurations, SUC-2, TUC-2, and AUC-2, were simulated for translational motion around a circle with a radius of 2500 mm. The simulation process, shown in Figure 23, verifies the correctness of the motion model of multiple-robot combination configurations. Figure 24a shows the trajectories of the centers O2 of SUC-2 and O4 of TUC-2; these two trajectories are very close, but because the simulation does not perform velocity compensation, they are smaller than the theoretical circle in the x-direction. In the process of simulation, the distance between the two robots in two configurations, SUC-2 and TUC-2, changes. The change curves of the distance between the two robots in these two configurations in the simulation process are shown in Figure 24b, which shows that the relative positions of the two robots change. From Figure 24c, it can be seen that the trajectory of center O6 of configuration AUC-2 changes greatly during several motion cycles, and the change of distance between the two robots in this configuration is relatively large and has certain periodicity, as shown in Figure 24d.

Translation Simulation along a Circle of Multiple Four-Mecanum-Wheeled Robots
Three robot configurations, SUC-2, TUC-2, and AUC-2, were simulated for translational motion around a circle with a radius of 2500 mm. The simulation process, shown in Figure 23, verifies the correctness of the motion model of multiple-robot combination configurations. Figure 24a shows the trajectories of the centers O 2 of SUC-2 and O 4 of TUC-2; these two trajectories are very close, but because the simulation does not perform velocity compensation, they are smaller than the theoretical circle in the x-direction. In the process of simulation, the distance between the two robots in two configurations, SUC-2 and TUC-2, changes. The change curves of the distance between the two robots in these two configurations in the simulation process are shown in Figure 24b, which shows that the relative positions of the two robots change. From Figure 24c, it can be seen that the trajectory of center O 6 of configuration AUC-2 changes greatly during several motion cycles, and the change of distance between the two robots in this configuration is relatively large and has certain periodicity, as shown in Figure 24d.      Figure 25 shows the translation motion simulation process of the seven combination configurations in Figure 19 along an 8-shaped curve, which is defined in Table 2. Figure 26 shows the state after these robots have moved for several cycles. In Figure 27, the motion trajectories of these simulations are compared. According to Figure 27, since no velocity compensation is performed in the motion simulations, there is obvious displacement error in the x-direction for all simulation trajectories. For the combination of two robots arranged side-by-side or end-to-end, the error of their motion trajectories in the first cycle is similar whether they are connected or not, as shown in Figure 27a,b. For the combination of connected multiple robots, the error of the simulation track does not show an obvious difference, as shown in Figure 27c. After several cycles of motion, the trajectory of each moving robot will shift to a certain extent, and these trajectories form a trajectory band, as shown in Figure 26. Taking AUC-2 as an example, the trajectory band after five cycles of robot combination rotation is shown in Figure 27d.  Figure 25 shows the translation motion simulation process of the seven combination configurations in Figure 19 along an 8-shaped curve, which is defined in Table 2. Figure 26 shows the state after these robots have moved for several cycles. In Figure 27, the motion trajectories of these simulations are compared. According to Figure 27, since no velocity compensation is performed in the motion simulations, there is obvious displacement error in the x-direction for all simulation trajectories. For the combination of two robots arranged side-by-side or end-to-end, the error of their motion trajectories in the first cycle is similar whether they are connected or not, as shown in Figure  27a,b. For the combination of connected multiple robots, the error of the simulation track does not show an obvious difference, as shown in Figure 27c. After several cycles of motion, the trajectory of each moving robot will shift to a certain extent, and these trajectories form a trajectory band, as shown in Figure 26. Taking AUC-2 as an example, the trajectory band after five cycles of robot combination rotation is shown in Figure 27d.     Figure 25 shows the translation motion simulation process of the seven combination configurations in Figure 19 along an 8-shaped curve, which is defined in Table 2. Figure 26 shows the state after these robots have moved for several cycles. In Figure 27, the motion trajectories of these simulations are compared. According to Figure 27, since no velocity compensation is performed in the motion simulations, there is obvious displacement error in the x-direction for all simulation trajectories. For the combination of two robots arranged side-by-side or end-to-end, the error of their motion trajectories in the first cycle is similar whether they are connected or not, as shown in Figure  27a,b. For the combination of connected multiple robots, the error of the simulation track does not show an obvious difference, as shown in Figure 27c. After several cycles of motion, the trajectory of each moving robot will shift to a certain extent, and these trajectories form a trajectory band, as shown in Figure 26. Taking AUC-2 as an example, the trajectory band after five cycles of robot combination rotation is shown in Figure 27d.    Figure 28 shows the translation motion simulation process of the six combination configurations in Figure 19 along a simple harmonic motion curve, which is defined in Table 2. Figure 29 shows the motion simulation trajectories of the robots. It can be seen that the simulation trajectories are close to the theoretical trajectories in the first motion cycle, but as the motion cycles increase, the error of the robot in the x-direction accumulates, causing the simulation trajectories to lag behind the theoretical trajectories. In the y-direction, trajectory errors are not accumulated, and because the longitudinal velocity errors are small, the displacements of the robots in the y-direction can almost reach the theoretical value.

Translation Simulation along Simple Harmonic Motion Curve of the Configurations of Multiple
Four-Mecanum-Wheeled Robots Figure 28 shows the translation motion simulation process of the six combination configurations in Figure 19 along a simple harmonic motion curve, which is defined in Table 2. Figure 29 shows the motion simulation trajectories of the robots. It can be seen that the simulation trajectories are close to the theoretical trajectories in the first motion cycle, but as the motion cycles increase, the error of the robot in the x-direction accumulates, causing the simulation trajectories to lag behind the theoretical trajectories. In the y-direction, trajectory errors are not accumulated, and because the longitudinal velocity errors are small, the displacements of the robots in the y-direction can almost reach the theoretical value.

Motion Simulation after Velocity Compensation for TCC-3
In Section 5.2, the motion simulations of the combination configurations of robot were carried out, and the kinematic model of the multiple-robot combination was verified. However, in the simulation experiment, there was no velocity compensation, so there were obvious errors between the robot's trajectories and the theoretical curves, especially in the x-direction of the local coordinate system. In this section, the robot configuration of TCC-3 and ACC-2 are selected and motion compensation simulation is carried out to verify the applicability of the velocity compensation model of multiple robot combination configurations. The velocity compensation coefficients used in this simulation are also obtained through the longitudinal and lateral translation motion simulation. In Figure 30, compared with the simulation curves without velocity compensation drawn by solid blue lines, the trajectories of the robot after compensation, drawn by solid black lines, are significantly improved, and are very consistent with the theoretical curves drawn by the dashed red lines.
Sensors 2020, 20, x FOR PEER REVIEW 33 of 37 In Section 5.2, the motion simulations of the combination configurations of robot were carried out, and the kinematic model of the multiple-robot combination was verified. However, in the simulation experiment, there was no velocity compensation, so there were obvious errors between the robot's trajectories and the theoretical curves, especially in the x-direction of the local coordinate system. In this section, the robot configuration of TCC-3 and ACC-2 are selected and motion compensation simulation is carried out to verify the applicability of the velocity compensation model of multiple robot combination configurations. The velocity compensation coefficients used in this simulation are also obtained through the longitudinal and lateral translation motion simulation. In Figure 30, compared with the simulation curves without velocity compensation drawn by solid blue lines, the trajectories of the robot after compensation, drawn by solid black lines, are significantly improved, and are very consistent with the theoretical curves drawn by the dashed red lines.

Conclusions
In order to solve the problem of transporting large-scale goods or equipment in industry, a combination system of multiple Mecanum wheels or multiple-Mecanum-wheeled robots is usually adopted. The kinematics of this combined mobile system is the theoretical basis for the motion control, path planning, and navigation of the combined system. In a real environment, slippage of Mecanum wheels has a great impact on the robot's motion accuracy, especially lateral motion. In this work, based on studying the kinematic constraints of a single Mecanum wheel in a mobile system, a kinematic model of Mecanum-wheeled robots with multiple wheels is derived, then a kinematic model with velocity compensation of a combination mobile system composed of multiple Mecanum-wheeled robots is created. Taking four-Mecanum-wheeled robots as example, motion simulations of a virtual prototype on RecurDyn and motion tests of a physical prototype were carried out to verify the kinematic model of the robot and the motion compensation model. For motion with small velocity change, the constant velocity compensation matrix has a good compensation effect, and the actual trajectories of the robot are very consistent with the theoretical curves. Finally, the motions of a variety of combined mobile systems composed of multiple Mecanum-wheeled robots were simulated on RecurDyn. Motion simulations and tests prove that the kinematic model of single robot and multi-robot combined mobile systems is correct, and the inverse kinematic correction model with velocity compensation matrix is feasible. Through simulation or experiment, the velocity compensation coefficient of the robots can be obtained and the velocity compensation matrix can be created. This modified inverse kinematic model can effectively reduce the errors of motion caused by wheel slippage and improve the motion accuracy of the robot mobile system.
More accurate motion compensation requires the use of a variable velocity compensation matrix, which is dynamically set according to the robot's velocity. According to the robot's parameters and environment information detected by position, attitude, and speed sensors, the robot's motion direction and velocity can be adjusted in real time by using this kinematic model with speed compensation, which is more conducive to motion accuracy.