Analysis of Tread ICRs for Wheeled Skid-Steer Vehicles on Inclined Terrain

The instantaneous centers of rotation (ICRs) for the two treads of skid-steer vehicles moving with low inertia on hard horizontal terrain almost remain with constant local coordinates, which allows to establish an equivalence with differential-drive locomotion. However, this significant kinematic relationship has not been analyzed yet on sloped ground. One relevant difficulty of studying ICR behavior on inclined terrain, even on a flat surface, is the continuous variation of pitch and roll angles while turning. To overcome this problem, this paper analyzes a dynamic simulation of a skid-steer vehicle on horizontal ground where gravity is substituted by an equivalent external force in such a way that pitch and roll are kept constant. Relevant tread ICR variations on inclined ground have been deduced, which have a significant impact on skid-steer kinematics. These new findings have been corroborated experimentally with a four-wheeled mobile robot that turns on an inclined plane.


I. INTRODUCTION
Skid-steer locomotion is of common use for mobile manipulators [1] and mobile robots [2] because of its good maneuverability and mechanical simplicity. Consequently, it is employed by ground vehicles in several applications such as agriculture [3], mining [4], surveillance [5], manufacturing [6], planetary exploration [7] or search and rescue [8].
In skid-steer vehicles, wheels or tracks have a fixed direction, so tread slippage is required for turning, which implies an extra power consumption [9], specially for tracked vehicles due to their larger contact surface [10]. Thus, in contrast with differential-drive kinematics, the motion of the vehicle is heavily conditioned by dynamics-related aspects such as the position of the center of gravity, terrain characteristics, or inclination on the ground [11], [12].
Nonetheless, differential-drive kinematics can be still employed on skid-steered vehicles by introducing slippage ratios [13], [14] or by using the instantaneous centers of The associate editor coordinating the review of this manuscript and approving it for publication was Hassen Ouakad . rotation (ICRs) of the treads on the ground plane [15], [16]. Some works combine both approaches for high-speed motion [17], while others only employ the distance between tread ICRs as their unique kinematic parameter [18], [19].
The actual positions of tread ICRs vary continuously on account of vehicle-terrain dynamics, but contrary to the vehicle's ICR, they remain within bounded regions at both sides of the vehicle when moving on hard horizontal ground [15]. In fact, ICR variations have been attributed mostly to the inertial force [20]. Thus, constant approximations of tread ICR values for flat terrain under low inertia can be identified with offline experiments [16], [21] or online during navigation [22].
ICR kinematics has been employed as the basis for path tracking [23], for online power consumption estimation [24], for dynamic modeling [25], for online terrain classification [26] and, also, as constraints for simultaneous localization and mapping [27].
Nevertheless, the restriction of ICR-based kinematics to flat horizontal ground constitutes one major limitation for a wider application of this method. One critical difficulty to analyze the ICR behavior for a skid-steer vehicle on nonhorizontal surfaces is that different tread speeds provoke the simultaneous variation of pitch and roll angles.
To overcome the horizontality restriction, this paper explores the position of tread ICRs while turning on sloped terrain. With this purpose, a dynamic study is proposed where gravity is substituted by an equivalent external force so that pitch and roll keep constant. The impact on skid-steer kinematics is deduced by analyzing ICR variations, which have been confirmed experimentally with a four-wheeled skidsteer vehicle.
To the best of our knowledge, the behavior of tread ICRs on inclined hard ground has not been studied before. Related works include the estimation of the minimum turning radii on slopes with limited motor torque [28], slip compensation on loose soil inclines [29] and dynamics-based path tracking on tilted surfaces [30].
In this way, the main contributions of the paper are the following: • Novel conditions are applied to dynamic simulations to obtain ICR positions on sloped terrain.
• Relevant kinematics findings for skid-steer vehicles on inclined surfaces are deduced from ICR locations.
• These new considerations are corroborated experimentally with a real mobile robot on a tilted platform.
The rest of the paper is organized as follows. Next section obtains ICR positions on inclined terrain through special dynamic simulations. In Section III, the resulting ICRs are analyzed to deduce their effect on skid-steer kinematics. Then, these findings are checked experimentally in Section IV. Finally, the last section discusses the results obtained.

II. DYNAMIC SIMULATIONS
Let the origin of the local reference frame XYZ of a wheeled skid-steer vehicle be at the center of the polygon defined by the wheel contact points with the ground. The Y axis is aligned with the forward motion direction, whereas Z and X point upwards and to the right, respectively. Let roll α and pitch φ be the angles of X and Y with respect to a horizontal plane, respectively (see Fig. 1).
To overcome the difficulty caused by the continuous variation of roll and pitch angles when a vehicle turns on inclined ground, it is proposed to perform dynamics simulations under the following special conditions: • The vehicle with a weight W is placed stationary on a hard horizontal plane with its wheels always in contact with the ground.
• Constant left V l and right V r speeds are applied to the treads.
• The tyre model includes specific friction coefficients for lateral µ x and longitudinal µ y motions to take into account slippage [31].
• The effect of gravity is eliminated. Instead, an external force F, whose components in the local reference frame depend on the desired constant values for α and φ: is exerted to the center of gravity (COG) of the vehicle (see Fig. 1). These conditions allow to obtain the steady values for angular ω z , lateral v x and longitudinal v y speeds of the vehicle, if they exist. Then, the coordinates of ICRs for the left (x l ICR , y ICR ) and right (x r ICR , y ICR ) treads on the local XY plane can be computed for any combination of constant values of V l , V r , α and φ as: where both ICRs have the same Y coordinate in this asymmetric kinematic model [16].

A. CASE STUDY
The above presented conditions have been applied to perform dynamic simulations of the mobile robot Lazaro [32], which is a battery-powered vehicle with skid-steer locomotion (see Fig. 2). In this work, Lazaro is regarded as a pure skid-steer vehicle with four identical wheels, so the additional contact point provided by the caster wheel of the onboard leg is not employed. Lazaro weights W = −255.1 N. The separation between the rear and front wheels is 0.4 m, and the distance between the left and right wheels is 0.398 m. The local coordinates of the vehicle's COG when its leg is in the rest position are located at x 0 = −0.0047 m, y 0 = −0.0024 m, and z 0 = 0.167 m.
The left and right tread speeds are provided by two direct-current motors equipped with angular encoders. The top speed for V l and V r is ±0.28 m/s, which coincides with the maximum speed of the vehicle during straight line motion. The employed dynamic model is similar to the one for the tracked vehicle in [15]. Basically, it consists of force and moment equations that depend on soil mechanics and vehicle properties such as mass, inertia and geometry. A relevant parameter is the inertia with respect to Z , which is I z = 1.34 kg m 2 . A punctual pressure distribution under the four wheels and anisotropic friction have been considered to shape the interaction between the rubber wheels with hard smooth ground. Concretely, lateral and longitudinal friction coefficients are defined as µ x = 0.4 and µ y = 0.54, respectively.

B. TREAD ICR POSITIONS
All in all, 6561(9 4 ) simulations have been performed by using constant tread speeds V l and V r between ±0.28 m/s with steps of 0.07 m/s and constant angles α and φ between ±12 • with steps of 3 • . Simulations have been carried out in Matlab by integrating the dynamic equations to obtain the possible stationaries for v x , v y and ω z .
From the 81(9 2 ) simulations for a given combination of pitch and roll angles, 9 of them correspond to straight-line motion with v x = 0 m/s, ω z = 0 rad/s and v y = V l = V r . These results confirm that there is enough friction for every tested pair of pitch and roll angles, but the calculation of ICRs with (2) is prevented because the resulting divisions are indeterminate (0/0).
The remaining 72 simulations may or may not achieve steady values for v x , v y and ω z to which (2) can be applied. Fig. 3 shows the number of solutions for each α and φ combination with V l = V r , where the multiply and the asterisk signs indicate 72 and 36 outcomes, respectively. Red and blue asterisks correspond to results obtained only with negative and positive values of sign(V r − V l ), respectively. Fig. 4 shows the resultant distribution for the left and right ICRs on the local XY plane with 8784 dots obtained by applying (2) on the 4392 steady solutions. ICR positions cover a wide area, but they always lie further from the tread contact  lines with the ground [15] at x = ±0.199 m. In addition, it can be observed that y ICR coordinates are limited by the contact lines of the front and rear wheels with the surface at y = ±0.2 m.

III. ICR ANALYSIS ON INCLINED TERRAIN
The ICRs for the 72 speed combinations of the α = φ = 0 • case, that are shown with black dots in Fig. 5, concentrate at the center of the ICR shapes of Fig. 4. Their mean values arex l ICR = −0.4151 m,x r ICR = 0.4157 m andȳ ICR = −0.0023 m, which can be used to approximate kinematics on horizontal ground [16].
Another combinations of pitch and roll angles present 36 stationary solutions separated from the other 36 (if they exist), according to sign(V r − V l ) and maintain small ICR variations regarding |V r − V l |. This has been illustrated in Fig. 5 with red dots for the (α = −3 • ∧ φ = −9 • ) case, where lower and higher y ICR coordinates correspond to V r > V l and V l > V r , respectively.
Having made evident the importance of sign(V r − V l ), the ICR dots have been represented separately by this factor in Fig. 6. The resulting diamond shapes with 2196 dots each are rotated to the left and to the right for V r > V l and for V l > V r , respectively.
Besides, Fig. 7 shows the contour lines of |φ − α| and |φ + α| on the ICR diamond shapes for V r > V l and V l > V r , respectively, where larger norms involve ICRs farther away from the shortest diagonal.
In fact, the diamond shapes of Fig. 6 can be divided into two parts, as shown in Fig. 8. The separation plotted with black dots corresponds to φ = α and φ = −α for V r > V l and for V l > V r , respectively. The yellow dots of the upper parts occur with φ < α and φ < −α for positive and negative signs of V r − V l , respectively. Conversely, the green dots of the bottom parts happen with φ > sign(V r − V l ) α. Furthermore, by plotting with black dots both the φ = α and the φ = −α cases, four separated ICR zones with 432 dots each can be distinguished in Fig. 9 for every diamond shape: I: III: (φ < α ∧ φ < −α) ∨ −φ > |α|, in blue color. IV: (φ < α ∧ φ > −α) ∨ α > |φ|, in yellow color.

A. KINEMATIC EFFECTS
The following results can be deduced for the center of each ICR zone of Fig. 9 with respect to the horizontal surface case given by (x l ICR ,ȳ ICR ) and (x r ICR ,ȳ ICR ). For V r > V l , it can be stated that: I: In the red case the ICRs shift to the right: the left ICR (x l ICR ) approaches its tread contact line, whereas the right ICR (x r ICR ) moves away from its tread contact line. This means that the left tread prevails over the right tread for traction. II: In the green case the ICRs move backward chiefly (y ICR decreases) and the rear wheels become predominant for traction. III: In the blue case occurs the opposite of case I: the ICRs shift to the left and the traction of the right tread prevails over the left tread traction. IV: In the yellow case occurs the opposite of case II: the ICRs move forward mainly (y ICR increases), which means that the front wheels become predominant for traction. 550 VOLUME 11, 2023 For V l > V r , it can be similarly stated that: I: The ICRs shift to the left and the right tread prevails over the left tread for traction. II: The ICRs move forward and the front wheels become predominant for traction. III: The ICRs shift to the right and the left tread prevails over the right tread for traction. IV: The ICRs move backward and the rear wheels become predominant for traction. This roughly happens for the center of each zone, but mixed behaviors can be observed in Fig. 9 along the transition between adjacent zones. Moreover, these effects grow with higher |φ| and |α|, that cause greater values for |φ − α| or for |φ + α| (see Fig. 7).
Cases II and IV do not affect much steering on slopes as the skid-steer vehicle maintains traction wheels at both sides with front or rear traction. However, cases I and III clearly imply a loss of vehicle steerability as one of the treads predominates over the other.

B. INCLINED PLANE
To illustrate the kinematic effects, the case of turning on spot on an inclined plane will be studied in detail. Let γ be the angle of the Y axis of the vehicle with respect to the direction Y g of the maximum slope β > 0 in the plane. Thus, the roll and pitch angles of the vehicle can be calculated as: Then, the following equivalence between the γ quadrants (see Fig. 10) and the four ICR zones of Fig. 9    When turning counterclockwise with V r = −V l > 0, the following will happen at the center of the quadrants: I: The traction of the right tread that points upwards is worse than the left tread that points downwards. II: The rear traction becomes predominant over the front traction. III: The traction of the left tread that points upwards is worse than the right tread that points downwards. IV: The front traction becomes predominant over the rear traction. Consequently, the vehicle will slide down while it tries to turn on spot. In addition, this effect will become more evident with greater values of β.
Similarly, when turning clockwise with V l = −V r > 0, it will happen: I: The traction of the left tread that points upwards is worse than the right tread that points downwards. II: The front traction becomes predominant over the rear traction. III: The traction of the right tread that points upwards is worse than the left tread that points downwards. IV: The rear traction becomes predominant over the front traction. This also implies that the vehicle will reduce its angular speed and it will descend the slope, mainly along quadrants I and III. The transition between adjacent quadrants will also entail a blend of their traction profiles.
All these results are counterintuitive with respect to the projection of the vehicle's COG moving with low inertia on an inclined plane, where it is expected that, independently of the sign(V r − V l ), the rear, left, forward and right tractions predominate in the I, II, III and IV quadrants, respectively. Contrary to the analysis performed, this approach does not predict any sliding down at all.

IV. EXPERIMENTAL RESULTS
To confirm the predicted behavior on an inclined plane, the mobile robot Lazaro has been placed on a rectangular melamine board 2.4 m length and 1.8 m width. Different VOLUME 11, 2023  angles β have been tested by tilting the platform along its longitudinal direction (see Fig. 11).
Let the origin of the global frame on the plane be at the center of the board with its Y g axis aligned with the length of the rectangle and pointing upwards. The X g axis, in turn, coincides with the width direction and points to the right.
The global pose of the vehicle on the platform, i.e., its Cartesian coordinates x g , y g and heading γ , has been calculated by matching the laser scans obtained by the onboard laser scanner every T = 0.1 s with a previous twodimensional (2D) map of the environment. For this purpose, the 2D Iterative Closest Point (ICP) method is employed [33] and the board is surrounded with wooden planks with distinctive corners to avoid left-right and up-down symmetries. In addition, v x , v y and ω z have been estimated by differentiating local variations in x, y and γ along time: The employed 2D range-finder is a Hokuyo URG-04LX-UG01, which has a maximum range of 4 m, an accuracy of ±3% of the measurement, a 240 • field of view and an angular resolution of 0.352 • . The sensor has been placed 0.095 m behind the origin of the local frame of the vehicle looking back (see Fig. 11).
In all the experiments, Lazaro has been commanded with V l = −V r = 0.21 m/s to turn clockwise. On horizontal plane (i.e., β = 0 • ), the 2D laser scanner traces a circular path when the vehicle approximately turns on spot (see Fig. 12) at a constant angular speed (see Fig. 13).
Even with a small inclination of β = 2.5 • , it is clearly visible in Fig. 14 that the mobile robot slides down while it rotates with almost constant ω z (see Fig. 15). The peaks to the 552 VOLUME 11, 2023  left with x g negative of Fig. 14 correspond to γ ≈ −π/2 rad in the IV quadrant. Conversely, the peaks to the right with x g positive occur near γ ≈ π/2 rad in the II quadrant.
Sliding is even more evident in Fig. 16 with a greater slope of β = 4.4 • , but small variations in angular speed, apart from measurement noise, are perceptible now (see Fig. 17).
The experiment with β = 7.6 • lasts much less time than the others to avoid the collision of the vehicle with the bounds at the bottom of the platform while it slides down (see Fig. 18). The maximum of |ω z | (see Fig. 19) occurs at γ ≈ ±π/2 rad, i.e., the centers of the II and IV quadrants. On the contrary, the minimum of |ω z | happens with γ = ±π rad, i.e., the center of the III quadrant. Table 1 contains the mean values of angularω z and slidinḡ v s speeds during the experiments. It can be observed that when the β angle increases the rotation speed decreases and the sliding speed increases, as predicted. Fig. 20 shows the ICRs calculated using (2) with the estimated values of v x , v y and ω z obtained with (4) for all the   experiments. The ICRs have been colored according to the quadrant where they belong. It can be observed that most of the ICRs for β = 0 • have a positive y ICR component, meaning that actual COG must also have a positive y 0 local coordinate, unlike it was assumed in the simulations. It is also visible in Fig. 20 that ICRs move away from the horizontal case as β increases. Nevertheless, ICRs always lie away from the tread contact lines and inside the front and rear wheel contact VOLUME 11, 2023  lines. For β > 0 • , the distribution of ICRs zones corresponds with the one obtained in the simulations (see Fig. 9-bottom).

V. DISCUSSION
Many interesting results have come from the qualitative analysis of tread ICRs of wheeled skid-steer vehicles on inclines performed in this paper. The most relevant ones are the following: • Significant changes in the local positions of tread ICRs on sloped terrain with respect to horizontal surfaces have been observed.
• The location of ICRs heavily depends on pitch and roll angles, as well as the sign of the difference between tread speeds. However, |V r −V l | has not a significant influence under low inertia.
• Four ICR zones, that depends on φ, α and sign(V r − V l ), have been identified with different traction behaviors. In two of these zones, the front or rear traction prevails.
In the other two zones one tread predominates over the other, which make steering difficult.
• Precisely, the traction of the tread with higher absolute speed tends to worsen or to improve when the vehicle is heading upwards and downwards, respectively, which provokes that the vehicle slides down while turning.
• This effect, which has also been observed experimentally in [30], is noticiable with small inclines and become more relevant with greater values of |φ| and |α|, i.e., greater β angles on an inclined plane.
• Sliding down while turning on moderated slopes is not caused by the lack of friction [34], but by the slippage inherent to skid-steering.
As a conclusion, we hope that all these findings can serve to formulate a parametric model for ICRs on inclined terrains capable of extending horizontal kinematics.