This paper presents a generalized, multibody dynamics model for a tracked vehicle equipped with a towing winch and control strategies that enhance vehicle mobility by regulating track slip based on real-time terrain characterization and automating winch use. The vehicle model is validated under conditions where no action is taken by the winch. Thereafter, two mobility enhancing control strategies are outlined. The first strategy regulates track slip to a real-time estimated value that generates maximum net traction. This is done by computing state-force estimates from a Kalman filter that are compared to terrain traction models using a Bayesian hypothesis selection approach. If the vehicle is traction limited and the first strategy fails, a second strategy that automates winch use is activated. Simulation results are shown for both scenarios.

## Introduction

Tracked vehicles have good performance in off-road locomotion over soft terrain when compared to wheeled vehicles. This is because the vehicle weight is distributed over a large area, which reduces sinkage and increases traction [1,2]. The demand for reliable autonomy for heavy class tracked vehicles is growing for operations involving agriculture, construction, and logistics. These tasks require towing heavy loads where immobilization can lead to long down times.

To maintain mobility autonomously, two approaches are proposed in this paper. The first is a traction control architecture outlined in Fig. 1. Here, a discrete time Kalman filter (DTKF) computes state-force estimates and compares these to terrain traction models stored in memory prior to operation. The mostly likely models are selected in a Bayesian hypothesis selection method. Based on the models selected, a “peak slip” point is approximated where the vehicle can operate with a track slip ratio that generates maximum net traction. This is the reference input to a traction control system that uses a proportional-integral-derivative (PID) controller to regulate track slip. Individual pieces of such an architecture have been proposed and implemented via simulation, but none implement traction control using a validated model that incorporates vehicle-terrain dynamics and limitations of the vehicle powertrain as presented in this paper. Related work in terrain characterization has focused on lighter weight, wheeled platforms [37] where the underlying terrain models assumed do not hold well compared with heavy-duty vehicles considered here. The approach of computing state-force estimates and comparing them to traction models stored offline to characterize terrain is similar to Ref. [6]. However, the formulation here uses a DTKF instead of an extended Kalman filter (EKF) so that the filter can run with fixed gain to reduce the real-time computation cost. Furthermore, we evaluate 320 terrain traction models in the Bayesian hypothesis selection where only 21 were evaluated in Ref. [6]. Previous work in traction controllers for tracked vehicles is limited. Fan derives a sliding mode controller and proves convergence via Lyaponuv stability theory [8,9]. However, it is implied that the track slip ratio can be measured perfectly for state feedback and does not take into account realistic limitations of the vehicle powertrain as an actuator.

The second approach employed to maintain mobility autonomously uses a winch control architecture. A picture of a CASE Quadtrac equipped with a winch is shown in Fig. 2. Hydraulic winches allow heavy-duty vehicles to reduce loading at the tractor's hitch or drawbar temporarily in soft terrain by releasing cable. Once firmer terrain is reached, the winch reels cable in to recover the payload. The winch control approach here uses a low bandwidth logic controller based on slip estimates from the DTKF. This gives the traction control mode time to respond when the two control modes are active simultaneously so that it can respond to varying drawbar loading. Previous work in cable towing systems has been investigated for marine vehicles where one vehicle tows another [1013]. In these studies, the path of the towed vehicle is the control variable of interest. Therefore, knowing the cable position and trajectory is essential for the desired control performance. Modeling of the cable includes finite-segment, rigid link lumped parameter elements. The control objective of interest motivating the modeling presented in this paper is mobility of the towing vehicle. Therefore, the simplified model of the cable here treats it as a rigid body that can only support tensile loads.

The work is motivated by the National Science Foundation South Pole Overland Traverse (SPoT) and Greenland Inland Traverse (GrIT) which use AGCO MT865 and CASE Quadtrac tractors, respectively. The traverses resupply research stations in Antarctica and Greenland, respectively, and optimizing the mobility of the tractor fleets has the potential to increase the efficiency of operations. The model and control development is for an MT865 tractor but the same approaches can be used for other vehicles of the same class.

The outline of this paper is as follows: Section 2 summarizes the vehicle-terrain model and validation of the model using experimental data from an AGCO MT865 tractor in Antarctica. Section 3 describes the tractor-winch-sled model, which incorporates hydraulic dynamics and hybrid dynamics of the winch system along with rigid-body sled dynamics. Section 4 describes an algorithm for real-time terrain characterization. Section 5 explains the design of the feedback traction controller. Section 6 describes criteria for entering and exiting the traction control modes based on the operating terrain parameter space. Section 7 shows simulation results of the traction control mode for five different terrain types to show the robustness of the control approach. Section 8 describes the design and justification for the rule set used in the logic-based winch controller. Sections 9 shows simulation results using both the traction and winch control modes, where the use of modes provides a robust mobility solution for heavy-tracked vehicles. Section 10 concludes the paper.

## Vehicle-Terrain Model Validation

Experimental data for an AGCO MT865 tractor were collected in Antarctica using the vehicle's CAN network, GPS, and a load pin to measure the drawbar force. Data were acquired while the vehicle towed a sled containing 80,000 kg of fuel bladders. The measured values include engine rpm ΩE, normalized driver input throttle ΠI, the transmission gear number or ratio gGR, the driver speeds $φ˙$, tractor ground speed vT, and drawbar load at the hitch DB. These values allow longitudinal motion validation of the model presented in Ref. [1]. The equations of motion reduce to
$v˙T=F−R−DBmT$
(1)

$φ¨=τ−Fr−ζφ˙JS$
(2)
where F is the gross traction generated by both tracks, R is the resistance from the terrain imparted on the track, DB is the drawbar load, mT = 25,000 kg is the tractor's nominal mass, τ is the torque of the driver sprocket in the track, ζ = 10 (kg/m2)/s a damping coefficient, and JS = 300 kg/m2 is the track's moment of inertia. The equations for gross traction and terrain resistance come from Wong–Bekker terramechanics theory [1,14] and are given by
$F=(Ac+W tan Φ)(1−Kil(1−e−ilK))$
(3)

$zS=(W/A(kc/b)+kΦ)1/n$
(4)

$z=zS+iSzS$
(5)

$R=b(kcb+kΦ)zn+1n+1$
(6)

where A =2.97 m2 is the nominal contact area of one track, c is the terrain cohesion in kPa, W is the vehicle nominal weight in kN, Φ is the terrain friction angle, l =3.25 m is the nominal contact length a track, i is the track slip ratio, K is the terrain shear deformation modulus in cm, kc is the cohesion modulus in kN/mn+1, kΦ is the friction angle modulus in kN/mn+2, n is the unitless exponent modulus, z is the total track sinkage in m, S is the unitless slip-sinkage coefficient, and zS is the static track sinkage in m. Equation (5) extends the original Wong–Bekker theory to account for slip-sinkage effects. zS is the static sinkage and the iSzS term predicts additional sinkage due to track slippage and associated terrain excavation. S values used in this study are based on experimental data in Ref. [15].

Engine speed estimates based on model predictions are derived from track speeds are given by
$ωE=gGRgFDφ˙$
(7)
where gGR is the gear number, gFD is the final drive gear ratio and ωE is the engine speed in rad/s [16]. The torque delivered to the tracks is a function of the engine speed and is calculated by
$τ=gGRgFD(ΠτE,max−τload)$
(8)
where Π is the normalized throttle input, τE,max is the maximum available torque at the engine's current operating speed that is calculated by interpolation of tabulated torque-speed engine data provided in Ref. [16], and τload is the load for the winch hydraulics. An overriding engine controller, ΠC is modeled to prevent stalling. Π and ΠC are given by MathWorks [17] and
$Π=max(ΠI,ΠC)$
(9)

$Π˙C=1δE((0.5−0.5tanhωE−ωIωT)−ΠC)$
(10)

For a more detailed discussion of the coupling of the powertrain and vehicle terrain dynamics, refer to Refs. [18] and [19].

The model is validated by simulation where a tractor starts from rest or zero initial conditions at t =0 using the inputs of the driver during data collection of ΠI and gGR and measured DB. Other inputs to the model for simulation are the estimated terrain parameters for compacted firn on which experiments were conducted: K =2 cm, c =10 kPa, and Φ = 19 deg and the total resistance Rtotal = R + DB. In Fig. 3, model predictions or simulation results are plotted against experimental data as blue and black lines, respectively, for 4 min of data during which the tractors starts from rest and ramps up to a steady velocity while shifting through gears. Measured and simulated responses for vT, $φ˙$, i, ΩE, and engine power PE are provided in Fig. 3. The blue line on the Π plot is the output of the overriding engine controller, ΠC.

The overriding engine controller overrides the driver input for the first 25–30 s in third gear and at several other points during gear shifts. This highlights the model's simplifying assumptions as presented, which assume constant track resistance forces. It is plausible that different terrain was encountered during the beginning of data collection such that this assumption does not hold for all time. In addition, it is postulated that the discrepancies during gear shifts can be accounted for by the use of oversimplified dynamics that do not account for control mechanisms. However, predictions at all other points prove to be accurate and show that the model sufficiently couples the powertrain with the vehicle-terrain dynamics.

## Tractor–Winch–Sled Model

Two important model elements are presented in this section to simulate the tractor–winch–sled system: the multibody vehicle dynamics that incorporate the winch and the hydraulics components that power it.

### Rigid Body Dynamics.

The rigid-body dynamics are derived using Newton–Euler methods. Figure 4 shows the free body diagram of the tractor–winch–sled system. A one-dimensional model is used since the dynamics of a heavy, tracked vehicle do not involve aggressive maneuvers and the height of the drawbar hitch is below the vehicle's center of gravity so as to prevent pitching. Because the cable is incapable of sustaining a compression load and maximum and minimum cable length limitations exist, this model has different sets of differential equations (phases) that govern its dynamics depending on its operating condition. A system of this type is called a hybrid system and is said to have hybrid dynamics. The two phases considered in this work capture the multibody dynamics under the following conditions: (1) Cable in tension, winch is unlocked and (2) Cable in tension, winch is locked. The first phase governs dynamics when tension in the cable, the load is being let out or reeled back in, and the minimum or maximum cable length limits have not been exceeded. The state vector and control input for this phase are $x1=[vT,ψ˙,vSD,DB,φ˙]T$ and u = [τ, τW]T, where ψ is the winch angular position, $ψ˙$ is the winch angular velocity, vSD is sled's velocity, τW is the winch torque, and the remaining variables have been defined previously. The governing differential equations and constraint are given by
$x˙1=M1−1Γ1(x1,u) s.t. DB≥0,ψ≥0,ψrW≤LC$
(11)
with mass matrix M1
$M1=[JS00000mT00100JW0−rW000mSD−10−1rW10]$
(12)
and force vector $Γ1$
$Γ1=[τ−Fr−ζφ˙F−R−τW−sign(ψ˙)ζWψ˙−RSD0]$
(13)

where JW = 0.82 kg2 is the winch rotational inertia, rW = 0.19 m is the winch radius, ζW = 3 (kg/m2)/s is the winch damping coefficient, LC is the total cable available, and ψrW indicates the amount of cable let out by the winch.

The second operating phase governs the dynamics of the tractor sled vehicle when there is tension in the cable and the winch is locked. This condition is met when the winch has let out the payload to its maximum length, reeled in the towed load, or has braked to stop the release of cable length. The state vector and control input in this phase are $x2=[vT,vSD,DB,φ˙]T$ and u = τ. The governing differential equations are
$x˙2=M2−1Γ2(x2,u) s.t. DB≥0,ψ≥0,ψrW≤LC,ψ˙=0$
(14)
with mass matrix M2
$M2=[JS0000mT0100mSD−10−110]$
(15)
and force vector $Γ2$
$Γ2=[τ−Fr−ζφ˙F−R−RSD0]$
(16)

### Winch Hydraulics.

Modeling of a hydraulic system allows limitations of the winch as an actuator and the load it puts on the engine to be realized. The system components contain a pressure-compensated, variable displacement piston pump connected to the engine through gear ratio gP = 1, a fixed displacement piston motor connected to the winch through the gear ratio gM = 18, a four-way two-position valve in the center, a pressure relief valve B1, and a braking valve B2 [2022]. The schematic for the proposed system is shown in Fig. 5 in two modes of operation: reeling in the winch and braking. In both modes, the kinematic constraints for the pump and motor speeds are
$ωP=ωEgP$
(17)

$ωM=ψ˙gM$
(18)
In Fig. 5(b), the position selected is for braking (II) and connects all four ways of the four-way two-position valve. This eliminates any pressure differential between PH and PO and activates the braking valve denoted as B2. The equations describing the pressure-rise rate within the section of hose PO in this mode of operation are given by [23]
$PH=0$
(19)

$P˙O=βVOQO$
(20)

$QO=DMωM−qB2$
(21)

$qB2={qleak PO≤Psetqleak+(qB2,Max−qleak)PMax(PO−Pset) Pset
(22)

$DP=0$
(23)

$τP=0$
(24)

$τW=DMPOgM$
(25)

The model used for the braking valve B2 is the same as for a pressure relief valve B1 to regulate hydraulic pressure. The valve set relief pressure, Pset is an input that modulates how much braking torque to use for the hydraulic motor in the winch.

The other position of the four-way two-position valve (I) isolates the sections of hose, PH and PO to create a pressure differential driven by the piston pump connected to the engine. The equations in this mode of operation are given by
$P˙H=βVHQH$
(26)

$PO=0$
(27)

$QH=DPωP+DMωM−qB1$
(28)

$qB1={qleak PH≤PRqleak+(qB1,Max−qleak)PMax(PH−PR) PR
(29)

$DP={DP,Max PH≤PsetDP,Max−DP,Max(PH−Pset)(PMax−Pset) Pset
(30)

$τP=DPPHgP$
(31)

$τW=DMPHgM$
(32)

This mode of operation for the hydraulics allows the winch to pull in the towed load. The value Pset in Eq. (30) is an input to the variable displacement, pressure-compensated pump to maintain a pressure set point for PH. This modulates the amount of torque the motor uses to pull in the towed load.

## Real-Time Terrain Characterization

Characterizing terrain mobility conditions in real-time requires knowledge of track slip and terrain-vehicle forces. However, neither of these can be directly measured, and therefore, force and slip are estimated from predicted and measured trajectory motion using a Kalman filter. These estimates are then compared with terrain slip–force curve models stored prior to operation where the most likely terrain and its associated peak slip are identified recursively using Bayes rule.

The Kalman filter design used here requires algebraic manipulations of Eqs. (1) and (2) since the only assumed measurements are $v˙T$, vT, $φ˙$, and DB. By defining the net traction as Fnet = FR and the resistance torque at the driver sprocket as $τres=F+ζφ˙$, Eqs. (1) and (2) can be rewritten as
$v˙T=Fnet−DBmT$
(33)

$φ¨=τ−τresJS$
(34)
where Fnet and τres and DB are considered part of the filter's state vector. In order to estimate these forces and torque, a second-order random walk model is used to augment the state equations
$[F̂˙F̂¨]=[0100][F̂F̂˙]$
(35)
$F̂$ is the estimated force or torque, and $F̂˙$ and $F̂¨$ are its first and second time derivatives. The state vector of the DTKF can now be written as
$x=[vTφ˙FnetF˙netτresτ˙resDBDB˙]T$
Since the process model is linear in the state vector x, the model is written in a continuous time state space form
$x˙=Ax+Bu$
(36)
where
$A=[001mT000−1mT00000−1JS000000100000000000000000100000000000000000100000000]$
(37)

$B=[0000000001JS000000]T$
(38)

$u=[0τ]T$
(39)
This set of linear, time-invariant, first-order differential equations is discretized using a zero-order hold filter for the system input u [24,25]
$xk=Fxk−1+Guk−1+wk−1$
(40)
Equation (40) yields a discrete time prediction model at an integration time-step of 0.05 s using additive Gaussian process noise $w∼N(08×1,R8×8)$ to model system disturbances with the covariance matrix tuned to
$R=diag([1×10−3m/s, 1×10−5rad/s, 2×106N, 2×105N/s, 2×107N⋅m, 2×108N⋅m/s, 1×106N, 2×105N])$
Trajectory motion is measured using the linear model
$yk=Hxk+vk$
(41)
where the measurement vector is $y=[v˙TvTφ˙DB]T$, and the measurement matrix is
$H=[001mT000−1mT0100000000100000000000010]$
and vk is additive white Gaussian noise defined by $vk∼N(04×1,Q4×4)$ where
$Q=diag([1×10−1m/s2, 3×10−3m/s, 2.8521×10−4rad/s,2.6795×107N])$
is estimated from the experimental data.
Since the derived Kalman filter uses LTI prediction and measurement models, the standard filter equations [25,26] converge to a steady-state filter gain Kss. Therefore, the filter reduces to
$x̂k−=Fx̂k−1++Gu$
(42)

$x̂k+=x̂k−+Kss(y−Hx̂k−)$
(43)
for online implementation where $xk−$ and $xk+$ are the prior and posterior state estimates. The posterior state estimate is used to approximate terrain conditions. Specifically, a smoothed slip estimate $ĩ$ and force vector $F̂≡[F̂net,τ̂res]T$ compute likelihoods of a fixed set of 320 terrain models parameterized by the vector $p=[cΦnkeqKS]T$. At each time-step $ĩ$ is given by
$ĩ=1−ṽT/(rφ˙̂)$
(44)
where $ṽT$ is a smoothed vehicle speed estimate from a first-order IIR filter running at 20 Hz
$ṽT[k]=0.8187ṽT[k−1]+0.1813v̂T[k]$
(45)
Terrain models are stored prior to operation as sets of [i, Fnet] and [i, τnet] curves computed from Eqs. (3) to (6) on a mesh of slip values i = [0.01% 0.11%…100%]T. To match a mesh value online, $ĩ$ is updated by
$ĩ=argimin(i−ĩ)2$
(46)
Using the updated smoothed slip estimate, the likelihood Gaussian probability density function (PDF) of $F̂k$ for the jth terrain parameter vector is
$pr(F̂k|pj)=1(2π)m|Σ|1/2erT(pj,ĩk)Σ−1r(pj,ĩk)$
(47)
where the covariance matrix is tuned to $Σ=diag([5×108N2, 5×108N2m2)$ and $r(pj,ĩk)$ represents the residual vector defined as
$r(pj,ĩk)=Fk(pj,ĩk)−F̂k$
(48)
The probability of each terrain model can then be computed recursively using Bayes rule by
$Pr(pj)k=Pr(pj|F̂k)=pr(F̂k|pj)Pr(pj|F̂k−1)∑i=1Npr(F̂k|pi)Pr(pi|F̂k−1)$
(49)
It is assumed that the actual operating terrain will not perfectly match a hypothesis. Therefore, the operating terrain is estimated by a weighted sum of the hypotheses vectors p1 to pN by
$p̂k=[p1p2…pN][Pr(p1)kPr(p2)k⋮Pr(pN)k]$
(50)

The parameter estimate $p̂k$ along with the slip mesh previously defined are used to recompute Eqs. (3)(6) where the slip value that corresponds to the maximum Fnet value is referred to as the estimated peak slip $îpk$ and is the reference input for the traction controller outlined in Sec. 5.

This approach requires a priori knowledge about the range of terrain parameters expected since only a finite number of models can be considered. The a priori knowledge assumed here are bounds on the six-dimensional terrain parameter space defined by the terrain parameter vector $p=[cΦnkeqKS]T$. Discretization of the terrain parameter space is as follows: $c∈{1,2,…,8}, Φ≡20 deg, n≡1, (kc/b)+kΦ≡keq∈{100,200,…,500}, K∈{0.5,1.5,…,7.5}, S≡(60%/33%)z$.

The justification for the values chosen here comes from Refs. [1], [2], [15], and [27]. Estimated values for Antarctic snow terrain cohesion are given in Ref. [27]. Data from Wong [1,2] show values slightly above and below 20 deg for Φ, and Lever estimates the range to be 15 deg ≤ Φ ≤ 25 deg [27]. Furthermore, c and Φ cannot be uniquely estimated or defined to produce the same traction forces due to the parameterization in Eq. (3), so there is no added benefit to discretizing both parameters and determining the likelihood of different combinations. The Bekker sinkage exponent n is defined to be constant as [27] estimates n =1 or close to this value for all Antarctic snow. keq is discretized to match the variety of different sinkage conditions experienced by SPoT defined by the Bekker relationship in Eq. (4). The discretization of K is larger due to its effect on the slip point at which the maximum net traction occurs referred to as the peak slip point, ipk. In order to allow for a range of ipk, K between 0.5 cm and 7.5 cm is considered. The last parameter S is defined as a constant due to the large amount of experimental data in Ref. [15] that shows consistent dynamic slip-sinkage behavior.

## Gain-Scheduled Traction Control

A block diagram of the closed-loop longitudinal tractor dynamics under traction control is shown in Fig. 6. Here, the throttle Π in the input to the engine and gear ratio selection gGR is the input to the transmission. The tractor speed vT and driver speed $φ˙$ are the states of the rigid body dynamics. The engine is considered as a time-varying gain depending on the engine operating speed, and the transmission is a piecewise constant gain value that is the product of gGR and gFD. The approximate continuous range for the engine gain is 1750 ≤ KE ≤ 2400 and the product of gGRgFD takes on discrete values dependent on the gear number (1–16).

The traction controller takes the form of a PIDF controller. The limited range of continuous engine gains can be handled by a single PIDF controller with fixed gains. However, the range of the 16 different transmission gains gGRgFD is too large for a fixed-gain controller, and the controller is a gain scheduled based on the gear ratio gGR.

The block diagram of the closed-loop system for regulating slip is shown in Fig. 6. The slip reference $îpk$ is the output from the Bayes' estimator in Sec. 4. In this control loop, slip is indirectly controlled by regulating the driver speed at the track. To calculate driver speed reference $φ˙ref$ the smoothed velocity estimate $ṽT$ is used
$φ˙ref=ṽTr(1−îpk)$
(51)
The lower limit at each time-step is $φ˙refgGRgFD≥136 rad/s(∼1300 rpm)$ so that the traction controller cannot stall the engine. The driver speed error $φ˙e$ is the input to the PIDF controller and is given by
$φ˙e=φ˙ref−φ˙̂$
(52)
The output of the PIDF controller ΠPID is combined with a feed-forward throttle term ΠFF and is the input to the engine
$Π=ΠPID+ΠFF$
(53)

This feed-forward term is the throttle input that is in effect just prior to entering the traction control mode and provides a stable transition between the leader–follower and traction control modes.

The structure of the PIDF controller is given by
$ΠPID(s)φ˙e(s)=Kp,t+Ki,ts+Kd,tstfs+1$
(54)
where the derivative term is low pass filtered to limit its bandwidth. The gains selected for the controller are Kp,t = 0.28, Ki,t = 0.17, Kd,t = 0.03, and the filter time constant is tf = 0.3 s. These gains are divided by the gear ratio gGR in each gear number to have a PIDF controller gain schedule for each gear. This way, the controller uses smaller gain values in gear numbers that have higher gear ratios and larger gain values in higher gear numbers with lower gear ratios. The aim here is that controller performance should be similar across different gear selections. This controller is discretized to operate the traction control loop in Fig. 6 at 20 Hz using the difference equation given by
$ΠPID[k]=−a1a0ΠPID[k−1]−a2a0ΠPID[k−2]+b0a0φ˙e[k]+b1a0φ˙e[k−1]+b2a0φ˙e[k−2]$
(55)

The final result is 16 sets of controller coefficients a0, a1, and a2 and b0, b1, and b2 for all 16 gear numbers, which are provided in Ref. [28].

Gear selections in the traction control mode force the engine to operate within a bounded engine speed range of ΩE ∈ [1350, 1950] so that the engine cannot stall or operate against the upper rpm limit (2100 rpm in gears 1–15 and 2300 rpm in gear 16). If the engine operating speed drifts below the lower bound of 1350 rpm, the transmission shifts down in gear number and if the engine speed goes above the upper bound of 1950 rpm, it shifts up in gear number. Gear shifts are limited to one change in gear number every 2 s, while detection of needed gear shifts occurs every 0.05 s so as to avoid limiting gear changes to every 2 s. For example, if a gear shift occurs at a time of t =12 s in the simulation, and the engine operates outside desired bounds at t =14.5 s but not t =14 s, a gear shift can occur before t =16 s. These changes in gear number are required for operation in the traction control mode, if there are changes in the tractor's ground speed so that the controller can track updated $φ˙ref$ values.

## Traction Control Mode Activation

The activation and deactivation of the traction control mode requires criteria for detecting soft terrain. There are a wide variety of possible approaches to this problem. Here, slip thresholds or limits are used to enter and exit the control mode. These thresholds are based on the bounds and discretization of the terrain parameter space in Sec. 4 and simulation experiments. Figure 7 shows a histogram of the peak slip ipk for the maximum values of Fnet for each model or hypothesis considered for estimation. The histogram shows that the upper and lower bounds for peak slip values are approximately ipk ∈ [5,40]. Since most of the terrain models have ipk ≤ 25% traction control is initiated once $î>25%$. Higher and lower threshold values were evaluated. Higher values led to increased failure rates of the control mode since the controller turns on too late to prevent immobilization. Lower values risk the traction controller turning on when there is sufficient traction to pull the payload, and the controller may not be able to converge to $îpk$. The slip threshold used for exiting the traction control mode is $î<5%$. This is chosen since data collected in Antarctica in Fig. 3 show slip values below this threshold on firm terrain, and it is below almost all ipk values in Fig. 7.

## Traction Control Mode Simulation

The results in this section simulate five different tractors colored magenta, cyan, green, blue, and red. Each one tows an 80,000 kg payload and transitions into the traction control mode at $t=16.25 s, t=18.05 s, t=18.1 s, t=18.5 s$, and t =19.45 s, respectively, from open-loop throttle and gear commands. These are Π = [0.24, 0.26, 0.3, 0.4, 0.5, 0.6]T and gGR = [2, 3, 4, 5, 6, 7]T for the time vector t = [0, 4, 8, 12, 16, 18]T. All tractors start on a firm terrain and travel into lower traction terrains stretching from 30 m to 130 m whose terrain parameters are in Table 1. Terrain parameters do not match any of the hypotheses of the 320 models from Sec. 4. Plots of slip i versus Fnet,DB for these lower traction terrains are provided in Fig. 8 where Fnet,DB is defined as
$Fnet,DB=Fnet−DB$
and DB is calculated across a mesh of slip values i by
$DB=(FnetmT+RSDmSD)(mTmSDmT+mSD)$
(56)

This is done by rearranging the equations from Sec. 3 for the dynamics in phase 2. The zoomed in plot in Fig. 8(b) highlights the narrow domain of slip values where tractors must operate with positive Fnet,DB values in order to maintain mobility.

A two-dimensional plot of lateral and longitudinal tractor positions taken at snapshots in time is shown in Fig. 9 for the five tractors colored magenta, cyan, green, blue, and red. All tractors are able to maintain tractor mobility using the traction control mode architecture except the red tractor, which operates on a terrain with the lowest cohesion c and shear deformation modulus K and becomes immobilized ∼53 m into the 100 m soft terrain stretch. Figure 10 shows the trajectories of all five tractors for the east or lateral position X, tractor speed vT, engine speed in rpm ΩE, selected gear ratio gGR, slip ratio i, and throttle input Π. The cyan and green tractors maintain mobility across the soft terrain stretch and continuously increase their ground speeds and shift up in gear. The magenta and blue tractors also maintain tractor mobility but are only able to maintain and not increase their ground speeds. The red tractor, however, continuously shifts down in its gear selection to try and regulate to its $îpk$ slip ratio.

The performance of the traction control mode for each tractor depends on the accuracy of the peak slip point estimate $îpk$ from the Bayesian hypothesis selection. Figure 11 plots the estimated peak slip point $îpk$ and maximum net traction estimate $F̂net,max$ for each tractor. Solid colored lines correspond to true values and dots to estimated values. Since none of the terrain parameter vectors given in Table 1 match any of the hypotheses, the peak slip point ipk cannot be perfectly estimated. In fact, each of the five terrain conditions differs in all of its terrain parameters from the 320 hypotheses except in the parameter n. However, the cyan and green tractors $îpk$ estimates converge close to their true ipk values of 24.5% and 21.5% so the traction control mode is able to maintain mobility and increase tractor ground speeds across the soft terrain patches. The $îpk$ estimates for the magenta and blue tractors converge to ∼24% and ∼22.5% when the true values for their terrains are 29% and 17%. Even though both of these estimates are off by 5%, the tractors to maintain mobility. The red tractor, however, converges to $îpk$ value of ∼21% when the true value is 8.5%. This error is much larger and leads to tractor immobilization.

The value of ipk for the red tractor's soft terrain is significantly different than the other soft terrains due to the very small shear deformation modulus K =0.7 cm. This small K value makes it difficult to accurately estimate the correct ipk value for two reasons. First, once the tractor operates under closed-loop traction control at a higher slip, determining whether or not the underlying terrain has a high or low K value is difficult since this parameter implies a steep slope of force-slip curves at lower slip values.

## Winch Control Mode

Section 7 presented the architecture of the traction control mode and demonstrated its effectiveness in simulation results. However, if there is no domain of slip ratios that provide positive Fnet,DB values, this mode will not prevent the tractor from being immobilized and the drawbar load on the vehicle must be reduced, which requires a towing winch to release the payload. An approach to automate winch use is proposed in this section and is referred to as the winch control mode. The hybrid rigid body dynamics and hydraulic actuator modeling for simulating this mode was discussed in Sec. 3.

The winch control mode uses a low bandwidth 4 Hz logic-based controller. This allows winch torques to be varied iteratively based on the vehicle's response to drawbar loading conditions. The low bandwidth of the controller provides the traction control mode time to respond to new loading conditions when the two are active simultaneously.

There are two stages during winch use. The first stage starts by activating the winch control mode. At this point, the tractor detects that the traction control mode alone cannot maintain vehicle mobility when $î≥45%$ and lets out the payload. Subsequently, the hydraulic winch slowly increases its hydraulic brake torque so as to bring the sled back into motion and reduce the rate at which the cable is being let out until it reaches zero. This is done by slowly increasing the hydraulic pressure of the brake valve. However, if the brake torque is increased so as to overload the drawbar, the hydraulic pressure is reduced, and then, the pressure is slowly increased again at a reduced rate. If the tractor reaches firm enough terrain, this will eventually result in the sled payload being brought back into motion and the winch will stop once the tractor and sleds speed match. The second stage of winch use pulls the payload back in. If stage 1 reaches completion, the terrain has sufficient shear strength to fully recover the payload and stage 2 will recover the payload more aggressively. The specifics of each stage are discussed in this section.

### Stage 1: Winch Activation and Initializing Payload Recovery.

Since recovering lost cable length can be difficult and time consuming to recover, there should be a strong indication that the traction control architecture alone cannot maintain tractor mobility before using the winch. At the same time, waiting too long to activate the winch could result in an immobilized tractor that has excavated itself into the terrain. To strike a balance between both concerns, the bounds on the terrain parameter space and their resulting slip–force curves are used as a guide. Figure 7 shows a histogram of the peak slip points ipk for all terrain hypothesis models. The largest value of ipk in Fig. 7 is ipk ≃ 41%. Therefore, once the estimates of $î≥45%$, this indicates that traction control alone cannot maintain tractor mobility, and the winch is activated to let out payload and reduce drawbar load.

The starting position of the four-way two-position valve in the winch hydraulics is shown in Fig. 5(b) before winch activation. This position is maintained once the control mode has been activated and the Pset value transitions from Pset = 2700 psi to Pset = 0 psi in Eq. (22) for the braking valve B2. This allows more hydraulic fluid to flow out of the PO node to reduce fluid pressure and subsequently the braking torque τW.

Once the winch has been activated, attempts to recover the payload begin immediately. This is due to the fact that soft terrain patch lengths are unknown to the tractor. Furthermore, recovering all the cable let out once the winch has been activated can be a time-consuming process.

Algorithm 1 runs at 4 Hz and summarizes the logic used to initialize payload recovery once the winch control mode has been activated. Line 4 gathers the latest original, unsmoothed slip estimate $î$ from the DTKF. Once this has been obtained, the hydraulic braking valve regulates the pressure set point Pset value based on current and past slip values. When the tractor first enters the control mode, Pset = 0 psi and increments at levels of incPSI = 15 psi in line 6 until the track slip reaches 20%. This threshold is chosen so that if the current drawbar load is at its limit and slip increases, the traction control mode has a chance to maintain mobility without operating under continuous load increases. The limiting factor of performance here is that tractors must initialize payload recovery on terrains having positive Fnet with i 20% as is the case with all curves in Fig. 8 except the magenta curve. Otherwise the algorithm cannot recover the payload. It should be noted, however, that the magenta curve has a K on the edge of the bounded terrain parameter space. If the pressure reference Pset that has been maintained in line 8 causes the tractor to detect $î>35%$, the Pset value is significantly reduced until the drawbar load is reduced to a point where $î<35%$ in line 10. Each time Pset is significantly reduced, incPSI is reduced at a rate determined by aCoeff where incPSI ≥ 5 psi. Each time this occurs, the tractor is allowed to make another attempt to pull in the load without having to reset Pset = 0 psi and can increase Pset to reload the drawbar at a slower rate so that the tractor operates at the appropriate drawbar load for the terrain. While Algorithm 1 runs at 4 Hz, $î$ is calculated at 20 Hz to make sure the tractor is not operating above 45% track slip. If this occurs, the braking pressure is reset to Pset = 0

Stage 1 is considered complete once the winch has stopped rotary motion and is no longer releasing cable length. At this point, the tractor is ready to enter stage 2 of the winch control mode to make a full recovery of the payload.

### Stage 2: Full Payload Recovery.

Stage 1 uses the braking valve B2 to reduce the rate of cable loss $ψ˙rW$ until $ψ˙rW=0$. At this point, stage 2 is initiated to make a full recovery of the payload and cable. When stage 2 is initiated, the four-way two-position valve is put into position I. Figure 5(a) illustrates the valve position change from II to I. This allows the variable displacement, pressure compensated pump to push hydraulic fluid in the H node, increase the pressure, and apply torque to the fixed displacement hydraulic motor to pull in the payload. Once this occurs, the Pset value in Eq. 30 is set to Pset = 15 psi [20].

This small initial Pset value is used since PH = 0 psi at the instant the valve position change occurs and the variable displacement, pressure compensated pump model puts the pump to its maximum displacement value DP,Max if PHPset. This minimizes the amount of time the pump fully loads the engine to provide a stable transition and to avoid stalling out the tractor's engine. After the initial transition, the winch follows a similar procedure to Algorithm 1. This brief rule-based procedure is summarized in Algorithm 2 and is identical to the procedure used in stage 1 but omits the cautionary step of greatly reducing drawbar pull if the vehicle detects $î≥35%$. There are several reasons for this. First, if the tractor has progressed from stages 1 to 2, the immediate terrain should provide enough traction to reel in the payload. Second, if the tractor is not able to recover the entire length of the cable due to another immediate soft terrain encounter, it is likely that the soft terrain patches that occur in the operating environment are too frequent and/or too long. Third, the drawbar requirements in the second stage are slightly different than in the first stage. In the first stage, the winch is attempting to bring the speed of the winch $ψ˙rW$ to zero. To make progress, the second derivative $ψ¨rW$ must always be negative. This acceleration requires a drawbar load constantly in excess of what is required to pull the payload in order to complete the first stage. However, in stage 2, a constant winch speed results in progress toward full recovery of the payload. Therefore, additional drawbar load in excess of what is required to tow the payload or a nonzero value of $ψ¨rW$ is only needed briefly instead of constantly as in stage 1.

## Winch Control Mode Simulation

In Sec. 7, results showed that all five tractors except the red one were able to maintain mobility across 100 m soft terrain stretch. However, the red tractor did traverse ∼53 of 100 m. Now, the red tractor will use its winch to traverse the remaining ∼47 m and beyond. Furthermore, it should be noted that the combined capability of the traction control mode and the winch control mode allows for terrain traversal without having to forfeit cable across the entire 100 m to minimize the amount of cable length required.

Figures 1214 show trajectory data of the red tractor using both the traction control mode and winch control mode. Figure 12 shows a two-dimensional plot of the tractor position at snapshots in time. The terrain operating conditions for this red tractor are identical to those in Sec. 7 and can be referenced in Table 1 and Fig. 8. Here, the tractor operates as before until ∼60 s into the simulation when the winch control mode becomes activated in stage 1 as triggered by a slip estimate $î≥45%$. This is evident in Figs. 13 and 14 where the slip i, pressure in the hydraulics PO, and drawbar load DB all drop to values close to zero to let the payload out at ∼60  s. Subsequently, Pset is slowly incremented as outlined in Algorithm 1. This can be seen in the slow increase in PO (Fig. 13) and in the drawbar load DB (Fig. 14). At ∼83 s, the tractor has detected that the drawbar is overloaded and decreases Pset and incPSI in lines 10 and 11 of Algorithm 1. This rapidly decreases the load on the tractor and increases the drawbar load more slowly for another attempt to initialize the payload recovery. At ∼90 s, the tractor exits the soft terrain stretch and is able to bring the winch to a locked state where no more cable is being let out. This can be seen in Fig. 14 where ψrW is the amount of cable length let out and $ψ˙rW$ is the speed at which cable is being let out. This process takes ∼15 s and stage 2 is entered ∼105 s into the simulation. The initial kickback from switching the discrete four-way two-position valve to position I creates a spike in the hydraulic pressure PH which highlights the reason for initializing the Pset value to a small value of 15 psi upon transition into stage 2. The Pset value is then slowly increased as outlined in Algorithm 2 to slowly increase the load on the engine τload and the speed at which the payload is pulled in $ψ˙rW$.

## Conclusion

A comprehensive tracked vehicle dynamics model equipped with a towing winch has been developed to include the effects of the power-train, terrain–vehicle interaction, and payload on deformable terrain. This allows for an iterative, model-based design approach for a traction and winch control modes. The traction control mode implements a complete control architecture of terrain identification and slip setpoint control based on realistic sensor data, actuator limitations, and a validated model. Simulation results validate the controller in different terrain conditions. Since the traction control mode cannot perfectly estimate terrain conditions due to noisy data and the nonlinear nature of slip traction curves, a logic-based feedback controller is proposed for the winch so that vehicles can extend the domain under which mobility is retained. Simulation results showed the effectiveness of the controller in allowing a tractor to maintain mobility by reducing drawbar load, where it previously could not use only the traction control mode.

Future work will investigate more sophisticated approaches for detecting transitions between the activation and deactivation of control modes. In this work, threshold values were used based on the terrain parameter space. This would allow for a more general applicability of the control architectures presented here by eliminating the need for tuned, heuristic thresholds.

## Acknowledgment

The authors thank Jeffery Zimmerman, Darin Motz, Ravi Godbole, and others from AGCO corporation who provided intellectual support in accessing tractor MT865 CAN bus data in Antarctica. Furthermore, the collection of these data was also made possible through coordination efforts by South Pole Traverse manager Kory Hotel and tractor operators Brad Johnson and Robert Shaw who drove tractors during experiments.

## References

References
1.
Wong
,
J.
,
2008
,
Theory of Ground Vehicles
,
1st and 4th ed.
,
Wiley
, Hoboken, NJ.
2.
Wong
,
J. Y.
,
2009
,
Terramechanics and Off-Road Vehicle Engineering: Terrain Behaviour, Off-Road Vehicle Performance and Design
,
Butterworth-Heinemann
, Oxford, UK.
3.
Iagnemma
,
K.
,
Kang
,
S.
,
Shibly
,
H.
, and
Dubowsky
,
S.
,
2004
, “
Online Terrain Parameter Estimation for Wheeled Mobile Robots With Application to Planetary Rovers
,”
IEEE Trans. Rob.
,
20
(
5
), pp.
921
927
.
4.
Hutangkabodee
,
S.
,
Zweiri
,
Y. H.
,
Seneviratne
,
L. D.
, and
Althoefer
,
K.
,
2006
, “
Soil Parameter Identification for Wheel-Terrain Interaction Dynamics and Traversability Prediction
,”
Int. J. Autom. Comput.
,
3
(
3
), pp.
244
251
.
5.
Ojeda
,
L.
,
Borenstein
,
J.
,
Witus
,
G.
, and
Karlsen
,
R.
,
2006
, “
Terrain Characterization and Classification With a Mobile Robot
,”
J. Field Rob.
,
23
(
2
), pp.
103
122
.
6.
Ray
,
L. E.
,
2009
, “
Estimation of Terrain Forces and Parameters for Rigid-Wheeled Vehicles
,”
IEEE Trans. Rob.
,
25
(
3
), pp.
717
726
.
7.
Ray
,
L. R.
,
Brande
,
D. C.
, and
Lever
,
J. H.
,
2009
, “
Estimation of Net Traction for Differential-Steered Wheeled Robots
,”
J. Terramech.
,
46
(
3
), pp.
75
87
.
8.
Fan
,
Z.
,
Koren
,
Y.
, and
Wehe
,
D.
,
1995
, “
A Simple Traction Control for Tracked Vehicles
,” American Control Conference (
ACC'95
), Seattle, WA, June 21–23, pp. 1176–1177.
9.
Slotine
,
J.-J. E.
, and
Li
,
W.
,
1991
,
Applied Nonlinear Control
, Vol.
199
,
Prentice-Hall
,
Englewood Cliffs, NJ
.
10.
Sanders
,
J. V.
,
1982
, “
A Three-Dimensional Dynamic Analysis of a Towed System
,”
Ocean Eng.
,
9
(
5
), pp.
483
499
.
11.
Kamman
,
J. W.
, and
Huston
,
R. L.
,
2001
, “
Multibody Dynamics Modeling of Variable Length Cable Systems
,”
Multibody Syst. Dyn.
,
5
(
3
), pp.
211
221
.
12.
Kamman
,
J. W.
, and
Huston
,
R. L.
,
1999
, “
Modeling of Variable Length Towed and Tethered Cable Systems
,”
J. Guid., Control, Dyn.
,
22
(
4
), pp.
602
608
.
13.
Milinazzo
,
F.
,
Wilkie
,
M.
, and
Latchman
,
S.
,
1987
, “
An Efficient Algorithm for Simulating the Dynamics of Towed Cable Systems
,”
Ocean Eng.
,
14
(
6
), pp.
513
526
.
14.
Bekker
,
M.
,
1956
,
Theory of Land Locomotion
,
University of Michigan Press
,
Ann Arbor, MI
.
15.
,
M.
,
2010
, “
Slip Sinkage Effect in Soil–Vehicle Mechanics
,”
J. Terramech.
,
47
(
1
), pp.
21
31
.
16.
Caterpillar
,
2002
, “
Challenger MT800 Series
,” Caterpillar, Peoria, IL.
17.
MathWorks
,
2015
, “
Generic Engine
,” MathWorks, Natick, MA.
18.
Cook
,
J. T.
,
Ray
,
L. E.
, and
Lever
,
J. H.
,
2016
, “
Multi-Body Dynamics Model of a Tracked Vehicle Using a Towing Winch for Optimal Mobility Control and Terrain Identification
,”
ASME
Paper No. DSCC2016-9626.
19.
Cook
,
J. T.
,
Ray
,
L. E.
, and
Lever
,
J. H.
,
2017
, “
Dynamics Modeling and Robotic-Assist, Leader-Follower Control of Tractor Convoys
,”
J. Terramech.
,
75
, pp.
57
72
.https://www.sciencedirect.com/science/article/pii/S0022489816301136
20.
MathWorks
,
2015
, “
Variable-Displacement Pressure-Compensated Pump
,” MathWorks, Natick, MA.
21.
MathWorks
,
2015
, “
Hydraulic Motor
,” MathWorks, Natick, MA.
22.
MathWorks
,
2015
, “
Pressure Relief Valve
,” MathWorks, Natick, MA.
23.
Manring
,
N.
, and
Luecke
,
G. R.
,
1998
, “
Modeling and Designing a Hydrostatic Transmission With a Fixed-Displacement Motor
,”
ASME J. Dyn. Syst., Meas., Control
,
120
(
1
), pp.
45
49
.
24.
Juang
,
J.-N.
, and
Phan
,
M. Q.
,
2001
,
Identification and Control of Mechanical Systems
,
Cambridge University Press
, Cambridge, UK.
25.
Stengel
,
R. F.
,
2012
,
Optimal Control and Estimation
,
Courier Corporation
, North Chelmsford, MA.
26.
Simon
,
D.
,
2006
,
Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches
,
Wiley
, Hoboken, NJ.
27.
Lever
,
J.
,
Ray
,
L.
,
Streeter
,
A.
, and
Price
,
A.
,
2006
, “
Solar Power for an Antarctic Rover
,”
Hydrol. Process.
,
20
(
4
), pp.
629
644
.
28.
Cook
,
J. T.
,
2017
,
Dynamics Modeling and Control of Autonomous Tractors for Polar Traverses
,
Dartmouth College
, Hanover, NM.