In Situ Soil Property Estimation for Autonomous Earthmoving Using
  Physics-Infused Neural Networks
A novel, learning-based method for in situ estimation of soil properties
using a physics-infused neural network (PINN) is presented. The network is
trained to produce estimates of soil cohesion, angle of internal friction,
soil-tool friction, soil failure angle, and residual depth of cut which are
then passed through an earthmoving model based on the fundamental equation of
earthmoving (FEE) to produce an estimated force. The network ingests a short
history of kinematic observations along with past control commands and predicts
interaction forces accurately with average error of less than 2kN, 13% of the
measured force. To validate the approach, an earthmoving simulation of a bladed
vehicle is developed using Vortex Studio, enabling comparison of the estimated
parameters to pseudo-ground-truth values which is challenging in real-world
experiments. The proposed approach is shown to enable accurate estimation of
interaction forces and produces meaningful parameter estimates even when the
model and the environmental physics deviate substantially.
W. Jacob Wagner, Ahmet Soylemezoglu, Dustin Nottage, Katherine Driggs-Campbell
| In Situ Soil Property Estimation for Autonomous Earthmoving Using Physics-Infused Neural Networks

###### Abstract

A novel, learning-based method for in situ estimation of soil properties using a physics-infused neural network (PINN) is presented. The network is trained to produce estimates of soil cohesion, angle of internal friction, soil-tool friction, soil failure angle, and residual depth of cut which are then passed through an earthmoving model based on the fundamental equation of earthmoving (FEE) to produce an estimated force. The network ingests a short history of kinematic observations along with past control commands and predicts interaction forces accurately with average error of less than 2kN, 13% of the measured force. To validate the approach, an earthmoving simulation of a bladed vehicle is developed using Vortex Studio, enabling comparison of the estimated parameters to pseudo-ground-truth values which is challenging in real-world experiments. The proposed approach is shown to enable accurate estimation of interaction forces and produces meaningful parameter estimates even when the model and the environmental physics deviate substantially.

Soli Property Estimation Fundamental Equation of Earthmoving Physics-Infused Neural Networks Autonomous Earthmoving

## 1 Introduction

Autonomous earthmoving is important for enhancing the safety, productivity and efficiency of construction and mining operations, especially in hazardous and remote environments where human operators face significant risks. One big challenge in designing these systems is enabling the equipment to operate efficiently and robustly across a wide variety of soils. Soils are complex heterogeneous granular materials whose properties vary significantly, not just between disparate locations, but also locally. As soil strength properties vary, the force required to shear the soil changes as well, which has a large effect on the earthmoving process. For example, in low strength soils a bulldozer can fill its blade in a short distance by making a deep cut without the vehicle stalling or tracks slipping. In contrast, cut depth is limited for higher strength soils because force limits of the machine are exceeded by higher required shear forces that increase with depth.

#### 1.1. Autonomous Earthmoving

Early approaches to autonomous excavation are posed as trajectory control problems where a desired trajectory is generated, typically based on some heuristic such as ensuring the swept volume is equal to the capacity of the bucket, and tracked using a position controller (Bradley and Seward, 1998). However, these kinematic trajectory control approaches fail when the generated trajectory is infeasible due to force limits of the machine. The simplest approach to address this problem is to assume the soil is very high-strength and to conservatively select trajectories that ensure that the machine does not stall and tracks do not slip during execution of the motion. Unfortunately, this is very inefficient in lower-strength soils and limits the practicality of autonomous earthmoving.

To address these issues, alternative control strategies have been proposed that enable adaptation to changing soil parameters. In bulldozers, blade control technologies that adjust blade position to ensure the machine is exerting a force near the optimal load without producing track slip are commercially available (Hayashi et. al., 2013; Jackson, 2017). Others have proposed impedance control (Ha et. al., 2000) or power maximization (Sotiropoulos and Asada, 2019) as a way to balance the position control objective with force limitations. Instead of following a kinematic trajectory, a prototypical force-torque trajectory can be tracked which produces kinematic trajectories that vary depending on soil composition (Jud et. al., 2017). More recently, reinforcement learning methods have been used to train control policies that achieve greater fill-factors in excavation and scooping tasks (Azulay and Shapiro, 2021; Backman et. al., 2021; Egli et. al., 2022).

However, local adaptation to soil conditions does not ensure efficient completion of a global terrain shaping task, implying that adaptation should be considered at the planning level as well. For an excavation task, iterative learning control has been used to adjust the desired trajectory between dig cycles to generate more feasible paths that result in increased bucket fill (Meada and Rye, 2012). Another approach is to estimate soil properties in situ and use these estimated properties to inform planning and control (Singh, 1995a). At the planning level, these estimated properties can be used in combination with a model of earthmoving to ensure generated trajectories are feasible and efficient (Singh, 1995b). To improve tracking of the trajectory, the interaction
forces can be predicted and used with impedance control (Tan 2005).

#### 1.1.2 The Fundamental Equation of Earthmoving

The first step in estimating soil properties is to assume a model, as the properties are tautologically tied to a soil model. In most soil property identification for autonomous earthmoving work, the Mohr-Coulomb model of soil shear strength is assumed which relates the shear strength of the soil, \(\tau\), to the applied normal force, \(\sigma_{n}\)via

\[\tau=c+\sigma_{n}\tan(\phi) \tag{1}\]

where the parameters \(c\) and \(\phi\) are the soil cohesion and angle of internal friction respectively.

Autonomous earthmoving requires a model of how the soil properties affect this process. Using the method of trial wedges where soil is assumed to fail along a plane producing a wedge, (Mckyes 1989; Reece 1964), an equation can be derived that describes the applied force \(F\) required to shear soil with a flat blade or bucket moving horizontally through the soil, as depicted in Fig. 1. The remaining forces acting on this soil wedge include the force of the loose soil accumulated on the blade which is referred to as surcharge or \(Q\), the weight of the soil wedge \(W\), the frictional component of the soil shear force combined with the normal force \(R\), the soil cohesive force, and the soil-tool adhesive force. The forces are assumed to be in equilibrium, i.e. neither the soil nor the blade is accelerating, and are summed up along \(\bar{x}\) and \(\bar{x}\) directions to arrive at the fundamental equation of earthmoving (FEE):

\[F=f_{FEE}(\Theta)=\gamma d^{2}wN_{r}+cdwN_{c}+QN_{Q}+c_{a}wN_{a} \tag{2}\]

where the coefficients are given by

\[\begin{array}{l}N_{r}=\frac{[\cos(\phi)+\cos(\bar{\theta})]\sin(\alpha+\phi +\beta)}{2\sin(\eta)},\ \ N_{c}=\frac{\cos(\phi)}{\sin(\rho)\sin(\eta)}\\ N_{Q}=\frac{\sin(\alpha+\phi+\beta)}{\sin(\eta)},\ \ N_{a}=\frac{-\cos(\phi+\phi+\beta)}{ \sin(\rho)\sin(\eta)}\end{array}\] (3-6)

where \(\eta=\delta+\rho+\phi+\beta\) is defined for notational convenience. This particular formulation of the FEE was developed by Holz et. al. (2013) to enable consideration of inclined surfaces. Initially, the equation has three unknown variables \((F,R,\beta)\), but applying the equilibrium conditions reduces this down to the single unknown soil failure angle \(\beta\). Typically, it is assumed that the optimal soil failure angle should be chosen such that

\[\beta^{*}=\arg\min_{\beta}N_{r} \tag{7}\]

In some cases, it is possible to derive a closed form expression for \(\beta^{*}\), but often numerical methods are relied upon. The interaction force can be broken into its \(\bar{x}\bar{x}\) subcomponents by

\[\begin{array}{l}F=f_{FEE}(\Theta)=\\ \left[F\cos(90^{*}-\rho-\delta+\alpha),\ F\sin(90^{*}-\rho-\delta+\alpha)\right] \end{array} \tag{8}\]

During bulloozing, the blade may move vertically in addition to horizontally as the desired cut depth is adjusted. This causes a change in the relative movement between the soil wedge and the blade. In the case that the wedge is moving down with respect to the blade, this results in the soil-tool adhesive and frictional forces changing direction. In order to account for these effects, Holz et. al. (2013) make the following modification

\[\delta^{\prime}=\tanh(-C_{1}(\overline{\eta_{b}}\cdot\overline{\nu}))\,\delta,\ c_{a}^{\prime}=\tanh(-C_{1}(\overline{\eta_{b}}\cdot\overline{\nu}))\,c_{a} ^{\prime}\] (9-10)

where \(C_{1}\) is a positive constant, \(\overline{\eta_{b}}\) is the unit vector pointing upwards along the blade, and \(\overline{\nu}\) is the velocity vector of the blade represented in the \(\bar{x}\bar{x}\) coordinate system. The values \(\delta^{\prime}\) and \(c_{a}^{\prime}\) are used in place of \(\delta\) and \(c\) in equations 2-6.

### Soil Property Estimation

There has been considerable attention on developing methods for estimating and predicting vehicle traversability for the purposes of autonomous off-road navigation (Borges et. al. 2022). Of particular interest for this work is the characterization of soil properties from on-board vehicle sensing including vision and propriception, sometimes referred to as inverse terramechanics (IT) estimation. In-situ soil property estimation using these methods has been studied within the field of terramechanics and is often motivated by the need for autonomous navigation of planetary rovers (Lopez-Arreguin et. al. 2021). This problem is formulated as the estimation of parameters of a soil-wheel interaction model and its often solved using optimization-based approaches (Iagnemma et. al. 2004), although filtering methods have been developed (Dallas et. al. 2020) and learning based methods have been used (Lopez-Arreguin and Montenegro 2021).

For autonomous earthmoving tasks, a similar IT approach can be taken. However, as the vehicle is more forcefully interacting with the terrain, different models must be assumed. Typically, a subset of the model parameters is assumed to be known, and the remaining parameters are found by fitting the model to observed interaction forces (Singh 1995a, Luengo et. al. 1998, Tan et. al. 2005, Althoefer et. al. 2009). It is often difficult to determine the accuracy of the estimated parameters directly, but Tan et. al. (2005) report predictions of the soil failure force to within 10%-30%.

### Contributions

The estimation of soil properties is a key challenge for autonomous earthmoving operations, such as bulldozing, Existing IT methods have mainly focused on the soil-tire interaction for vehicle navigation applications and, while motivating, these methods do not directly apply to earthmoving tasks due to different underlying physics. The relatively limited prior work in IT methods for construction tasks leverage optimization approaches to find parameters of the earthmoving model given observations of the interaction force. In this work, we propose a novel IT method based on physics-infused neural networks

Figure 1: FEE soil wedge geometry and force diagram. This is a cross-sectional view depicting a blade, shown in gray, moving in the \(\bar{\nu}\) direction through the soil at a depth of cut \(d\) and failing the soil along the purple line at an angle \(\beta\) with respect to the surface. The soil surface, shown in green, is inclined at the angle \(a\) from the horizon and the flat blade is at an angle \(\rho\) with respect to the soil surface. The passive failure case is assumed where the soil wedge moves up with respect to the blade and forward with respect to the un-sheared soil. The cohesive forces resist this movement and are drawn accordingly. (Holz et. al. 2013)
(PINNs) that can estimate soil properties in situ from kinematic and control observations of a bulldozer blade. The main contributions of this work are:

* Development of a PINN-based in situ soil-property estimation method that incorporates physical laws and constraints into the learning process without requiring knowledge of true soil properties
* Evaluation of the accuracy and robustness of the proposed method using simulation data with ground-truth soil properties
* Introduction of a physical inconsistency loss function that enables enforcement of minimization or maximization of certain physical quantities
* Demonstration of a novel technique to estimate the blade-soil interaction forces directly from kinematic and control data

## 2 Methods

The creation of a simulation to study the in-situ soil-property estimation problem is first discussed, and the PINN approach to soil property estimation is then developed.

### Simulation

In this work, the proposed soil property estimation method is evaluated to determine how capable the system is at extracting characteristics of the soil. Real world evaluation a challenging task, as evaluation of the full system would require significant instrumentation of a bulldozer, time consuming real-world experimentation across various soil types and levels of compaction, and careful setup (e.g., compaction, moisture content) and characterization of these soils to establish a soil property ground-truth to enable validation of the approach. Instead, using a simulation as an analog for experiments on a real bulldozer is proposed. To perform these experiments, a simulator that exhibits similar behavior to a real system is required, but faithful representation of the dynamics of a specific piece of earthmoving equipment for a specific soil is not strictly necessary.

#### 2.1.1 Vortex Studio

Vortex Studio1 is a simulation software which enables modelling of multi-body dynamics and supports accurate real-time soil-tool interaction physics. In Vortex, the soil is modeled using a hybrid heightfield-particle representation. An FEE based model is used to determine the reaction force of the soil prior to shearing. As the relative density \(I_{\alpha}\) of the soil increases, \(\phi\) and \(c\) are varied to increase the resistance of the soil to shearing due to stronger particle interlocking, and \(\gamma\) is varied to adjust the weight of the soil wedge and surcharge (Holz, 2009).

Footnote 1: [https://www.cm-labs.com/en/vortex-studio/](https://www.cm-labs.com/en/vortex-studio/)

When the force exerted by the tool on the soil reaches this threshold, the portion of the heightfield in contact with the blade is converted into particles whose behavior is governed by a discrete element model (DEM) simulation which more accurately characterizes the behavior of the disturbed soil. As particles come to rest, they are then reintegrated back into the heightfield representation. The initial relative density is defined as a fixed value but may vary throughout the simulation as the soil is compacted by the tool or as sheared particles are merged back into the heightfield (Holz, 2009; Holz et. al., 2013). The remainder of the soil properties remain constant for a given soil type and do not vary within a simulation episode.

Haeri et al. (2020) showed that Vortex Studio is capable of predicting the soil-tool interaction forces for a lunar simulant material within 20%-30% of the measured force obtained from real-world experimentation. They introduce some modifications to the original Vortex model, including a term that reduces the contribution of the surcharge by a factor of 10 in equation 2. Additionally, they find that selecting \(\beta\) such that the top surface of the soil wedge matches extent of the surcharge pile, better matches the experimental data as compared to minimization of \(N_{\gamma}\) as in Equation 7.

One benefit of using a simulation to study soil-property estimation is that it can provide ground-truth knowledge of the soil parameters for validation of a given approach. However, in Vortex, this is not straightforward as the simulation is not a purely physics-based model, but rather a combination of physics and heuristics that have been added to improve the realism, stability, and performance of the system. In addition to the modifications discussed by Haeri et. al. (2020), these heuristics include various adjustments to the FEE force calculation, such as limiting the force based on the submerged surface area, the addition of penetration forces, scaling of the FEE force, inclusion of additional frictional, elastic, and damping contact forces between the heightmap and tool, etc. These heuristics are not directly derived from the soil properties, but rather must be tuned to match experimental data or may be modified to incorporate operator feedback, improving the realism of the simulation.

Vortex provides four default soil configurations (Clay, Loam, Sand, and Gravel), but further tuning is recommended for specific applications. However, it is not obvious how to configure these parameters to achieve realistic performance for different soil types and blade geometries. Despite these limitations, Vortex meets the needs for a system that exhibits realistic earthmoving phenomenon and is sufficiently complex to act as an analog for experiments on a real machine.

#### 2.1.2 Simple Bladed Vehicle Model

A simplified simulation of a bladed earthmoving vehicle is developed, see Figure 2, and consists of two dynamic

Figure 2: Rendering of simple bladed vehicle simulation. The chassis is shown in green and the flat blade collision geometry is shown in blue. The soil particles generated due to the blade failing the soil constitute the surcharge.
components: a chassis and a blade. The chassis is modeled as a rectangular prism (3m, 2m, 2m) with a mass of 5,000kg. The blade is flat 3.164m wide, 0.660m tall, 0.001m thick2 with a mass of 400kg.

Footnote 2: Originally, the blade was made to have a more reasonable thickness, but initial experimentation revealed problems with the vehicle being able to penetrate the soil surface even with large forces for denser soil configurations. The authors believe that the is an artifact of the way contact detection is used to compute the soil blade angle \(\rho\) internal to Vortex. The hypothesis is that the system computes the angle using the bottom surface of the blade instead of the leading edge causing the simulator to produce unrealistically large reaction forces. Reducing the blade width alleviated most of these problems.

Vortex constraints, which enable restricting the motion between two or more parts, are used to control the motion of the chassis and the blade. Motorized constraints enable control of a part's velocity and exert a force proportional to the velocity error. The inverse of the loss parameter specifies this constant. Locked constraints enable control of a part's position and exert a viscoelastic force on the part equivalent to a spring-damper system (CM Labs Simulations, 2016). A motorized constraint is used to control the forward velocity of the chassis and a locked constraint is used control its vertical position. The lateral position and orientation of the chassis are locked, limiting the motion of the vehicle to the \(\bar{xz}\) plane. The blade is attached to the chassis body using a locked constraint enabling control of the vertical position relative to the chassis.

To emulate limits on a real machine, the relative position of the blade with respect to the chassis is limited to enable the bottom edge of the blade to reach 0.3m. below and 1.0m above the surface when the chassis is on level ground. The soil height is measured at four locations underneath the chassis and the average of these height is used to command the chassis vertical position to emulate a simple vehicle suspension. Additionally, the vertical and horizontal constraint forces are limited to 30kN and 20kN respectively to emulate the tractive and penetrative force limits of a bladed earthmoving machine.

On a bulldozer, hydraulic pressures in the vehicle hydrostatic drive circuit and in the cylinders controlling the blade position can be measured to derive estimates of the tractive force given knowledge of the machine's dynamics (Yamamoto et. al. 1997). The simulation developed in this effort has a simpler drive system, therefore, the simulator chassis-blade constraint forces are used to measure the interaction force instead. Significant effort was given to tune the loss, stiffness, and damping coefficients of the constraints to yield a stable simulation with minimal oscillations in the constraint forces and part positions. However, oscillations on the order of 5kN and 10kN in the x and z forces respectively are commonly observed in the experiments, see Fig. 3. These oscillations arise from the interaction of the constraints used to control the vehicle and the constraints imposed by the soil-tool interaction model. In particular, when a portion of the heightfield is converted into particles, the blade momentarily loses contact with the heightfield resulting in a significant disturbance in the force applied to the blade. This amounts to a step input to the control constraints and can lead to transient violations of the force limits, as Vortex constraints are essentially dynamical systems.

During simulation, at a rate of 60Hz, observations are collected of the blade \(x\) and \(z\) positions and chassis \(z\) positions, \(\mathbf{p}=[p_{x}^{b},p_{y}^{b},p_{z}^{c}]\), the blade \(z\) velocity and chassis \(x\) and \(z\) velocities, \(\mathbf{v}=[v_{x}^{b},v_{x}^{c},v_{z}^{c}]\), and the \(x\) and \(z\) cutting forces, \(\mathbf{F}=[f_{x}^{b},f_{y}^{b}]\). Additionally, the commanded blade relative and absolute \(z\) position, and the chassis \(x\) velocity \(\mathbf{u}=[u_{x}^{b},u_{x}^{b},u_{x}^{c}]\) are collected. At each timestep \(t\), the observations and actions are concatenated into the observation vector \(\mathbf{o}_{t}=[\mathbf{p}_{x},\mathbf{v}_{x},\mathbf{u}_{t}]\) which is consumed by the PINN. Pseudo-ground truth soil-propense are collected from the simulation as \(\mathbf{\theta}=[\phi,c,\delta,\alpha_{x},\gamma,\rho,a,w,d,d,\partial,\vec{v}]\), where \(Q\) is obtained using the Vortex soil mass sensor plugin. The simulation does not provide access to the soil failure angle \(\beta\), but for notational simplicity it will be included in \(\mathbf{\theta}\) for the remainder of the text.

#### 2.1.3 Data Collection Controller

To ensure clean data is collected to train an effective model, a simple data collection controller was designed. In the case that the force required to shear the soil is larger than the machine limits, the blade becomes stalled, and the soil is not sheared. Therefore, an anti-stall blade controller is developed which enables automatic data collection.

For a given desired forward velocity \(\vec{v}_{x}\) and target absolute cut depth \(\bar{d}_{x}\), the controller adjusts the blade relative position command \(u_{x_{r}}\) to ensure that the velocity tracking error \(e_{v_{x}}\) does not drop below a target velocity dependent threshold \(\vec{e}_{v_{\text{min}}}=\vec{v}_{x}e_{v_{x}}\). This is accomplished using a gain-scheduled integral controller which reduces the desired depth of cut by a depth offset \(\Delta\bar{d}_{x}^{k}\) at each timestep \(k\); this quantity is limited to enable raising of the blade only a small amount above the surface. The controller is defined as follows

\[\Delta\bar{d}_{x}^{k}=\text{K}_{\text{App}}\left(\text{e}_{v_{x}}^{k}-\text{e }_{v_{\text{min}}}\right)+\Delta\bar{d}_{x}^{k-1}\quad\in\left[0,\bar{d}_{x}+ \Delta\bar{d}_{x_{\text{max}}}\right] \tag{11}\] \[\text{u}_{x_{r}}=-(\bar{d}_{x}-\Delta\bar{d}_{x}^{k}) \tag{12}\]

The scheduled depth offset gain term, \(K_{\text{App}}=K_{v}/\bar{v}_{x}\), modulates the depth offset integral gain \(K_{v}\) to enable improved forward velocity tracking for low desired velocities \(\vec{v}_{x}\). Overall, the effect of this controller is that vehicle maintains forward motion with some low-frequency oscillation in the depth of cut and forward velocity around the targets. As the soil becomes

Figure 3: Measured forces and cut depth for a typical episode within the \(D_{\text{default}}\) dataset. This episode was recorded for a loan soil with \(\bar{d}_{x}=60\%\) and illustrates the large amplitude noise in force measurements observed particularly for denser soils. Also, note how the surcharge grows with time indicating the accumulation of soil on the blade from the cutting operation.
harder to shear and vehicle force limits are reached, the desired depth of cut becomes harder to track, but stalling is prevented ensuring the soil failure condition assumed by the FEE is met.

### Physics-Infused Neural Network (PINN)

Broadly, the approach taken in this effort is to train a neural network to produce estimates of the soil parameters which are then fed through the FEE model to produce an estimated force. The network ingests a \(T\) step history of observations and control commands \(\mathbf{\theta}_{[tx+T]}\) along with known FEE parameters \(\mathbf{\theta}_{\mathbf{k}}=[c_{a},\mathbf{r},\rho,a,w,d,Q,\overline{\mathbf{v}}]\), and outputs an estimate of the unknown parameters \(\widehat{\mathbf{\theta}}_{\mathbf{u}}=[\widehat{\phi},\hat{c},\delta,\widehat{\theta },\underline{\mathbf{\theta}},\underline{\mathbf{\theta}},\underline{\mathbf{\theta}}]\) and the estimated interaction force \(\widehat{\mathbf{F}}\). It is important to include control actions in addition to kinematic information, as tracking error can be thought of as a measure of the interaction force. For example, Bradley and Seward (1998) use the position tracking error to achieve "software force feedback" that can be used to determine when the bucket comes into contact with the ground.

By properly constraining the network, the hypothesis is that the model will converge to estimating reasonable soil parameters because this representation is the most parsimonious representation of the data. In other words, the model should learn to estimate meaningful soil properties as they explain the observations in the simplest form.

#### 2.2.1. Network Architecture

While it may be possible to estimate FEE parameters from observations and commands at a single time-step, this is difficult given the noise present. Therefore, observations and actions across \(T=60\) timesteps (1 sec.) are incorporated, and a temporal convolution is leveraged to compress the \(T\times 9\) dimensional observation \(\mathbf{\theta}_{[tx+T]}\) into a \(l\times N\) vector, where \(l\) and \(N\) are model hyperparameters. This vector is passed through a stacked transformer style encoder network with multi-headed attention, modeled after the network presented by Zerveas et. al. (2021), producing another \(l\times N\) vector. This vector is then combined with the known FEE parameters \(\mathbf{\theta}_{\mathbf{k}}\) and passed through a dense network with 2 hidden layers of 20 neurons each that is referred to as the integration network. The role of this network is to incorporate the knowledge of the existing soil parameters with the latent state produced by the transformer from the observation history. This enables estimation of the unknown soil parameters. Without this network, the system cannot learn to account for variation in \(\mathbf{\theta}_{\mathbf{k}}\) by changing \(\widehat{\mathbf{\theta}}_{\mathbf{u}}\). A linear output layer produces the estimated unknown FEE parameters \(\widehat{\mathbf{\theta}}_{\mathbf{u}}\) and the residual end to end force \(\mathbf{F}_{r}\). The FEE parameters are combined \(\widehat{\mathbf{\theta}}=[\mathbf{\theta}_{\mathbf{k}},\widehat{\mathbf{\theta}}_{\mathbf{u}}]\), clipped to remain within the limits defined in Table 1, and passed to the FEE network to produce the FEE reaction force \(\widehat{\mathbf{F}}^{\text{FEE}}=\mathbf{f}_{\text{FEE}}(\widehat{\mathbf{\theta}})\). The residual and FEE reaction forces are then summed to produce the estimated interaction force \(\widehat{\mathbf{F}}=\mathbf{F}^{\text{FEE}}+\mathbf{F}^{\prime}\). See Fig. 4 for a visual depiction of this architecture.

The residual depth of cut which augments the depth of cut by \(d^{\prime}=d+\Delta d\), was added to the model after initial experimentation revealed significant sensitivity of the FEE model to the depth of cut; major improvements in reconstructing the observed force entirely from the ground-truth parameters was shown to be possible with hand tuning of \(\Delta d\). Similarly, the residual interaction force \(\mathbf{F}^{\prime}\) was added to the network to avoid non-FEE components of the simulation (e.g. noise in the measured force and non-zero measured force when \(d<0\)) corrupting the parameter estimates.

#### 2.2.2. Network Training

To allow for stable training, observations \(\mathbf{o}\) are normalized to have unit mean and variance where the normalizing mean and variances are computed from the training set observations. Values for the residual force are normalized \(\mathbf{F}^{\text{F}}\) to 10% of the value of \(\mathbf{F}\). Both known and unknown FEE parameters are scaled using min-max normalization with ranges provided in Table 1. The ranges for \(\phi,c,\delta\), and \(c_{a}\) are obtained from multiple sources listing geotechnical parameters (Geotechnata.info 2023, MATLAB 2023, and Fine Software 2023) while the remaining parameters are scaled using knowledge of machine limits and geometry of the FEE. During tuning of the model, some liberty was taken to adjust parameter ranges for \(\beta\) and \(c\) from the initially assumed values in order to avoid the model converging to a poor local optimum.

The loss function for the network is comprised five types of losses: a mean average error force prediction loss (\(\mathcal{L}_{MAE}^{F}\)), a MAE regularization of the residual force (\(\mathcal{L}_{MAE}^{F}\)), ReLU-based regularization of FEE parameters (\(\mathcal{L}_{\text{ReLU}}^{T}\)), a mean square error (MSE) based regularization of the residual depth of cut (\(\mathcal{L}_{MSE}^{\Delta d}\)), and a MAE regularization of the gradient \(\partial N_{\gamma}/\partial\beta\) (\(\mathcal{L}_{MAE}^{\Delta N_{\gamma}/\partial\beta}\)). It is defined as follows

(13) \[\mathcal{L}_{MAE}^{F}(\mathbf{\bar{\theta}},\widehat{\mathbf{F}},\mathbf{F}^{ \prime})=\lambda_{F}\mathcal{L}_{MAE}^{F}(\mathbf{F},\widehat{\mathbf{F}},\mathbf{w}_{xx}) +\lambda_{F}\mathcal{L}_{MAE}^{F}(\mathbf{F}^{\prime},\mathbf{0},\mathbf{w}_{xx})+\] \[\sum_{x\in\mathbb{Z}}[\lambda_{x}\mathcal{L}_{\text{ReLU}}^{F} (\lambda,\mathcal{L}_{\text{ReLU}}^{F}(\lambda,\mathcal{L})]+\lambda_{AdMSE} \mathcal{L}_{MSE}^{Ad}(\Delta d,\underline{\mathbf{\theta}})+\] \[\lambda_{0N_{\gamma}/\partial\beta}\mathcal{L}_{MAE}^{\Delta N_{ \gamma}/\partial\beta}\left(\frac{\Delta N_{\gamma}}{\partial\beta},\mathbf{0}, \mathbf{w}_{\partial N_{\gamma}/\partial\beta}\right)\] \[\mathcal{L}_{MAE}^{F}(\mathbf{\bar{\beta}},\mathbf{y},\mathbf{w})=\sum_{i=0}^{ \left|\mathbf{y}\right|}\mathbf{w}_{i}|\overline{\mathbf{y}}_{i}-\mathbf{y}_{i}|\] (14) \[\mathcal{L}_{MSE}^{MSE}(\Delta d,\underline{\mathbf{\theta}})=\left( \underline{\mathbf{\theta}}\ -\ \underline{\mathbf{\theta}}\ \right)^{2}\] (15) \[\mathcal{L}_{\text{ReLU}}^{T}(\lambda,\mathcal{L}^{T})=max(0, \mathcal{L}^{T}_{i}-\chi)\ +\ max(0,\chi-l_{u}^{T})\] (16)

The FEE parameters regularized to remain within the limits listed in Table 1 are \(\mathbf{X}=[\phi,c,\delta,\beta,\eta,\zeta,\underline{\mathbf{\zeta}},\underline{ \mathbf{\zeta}},\underline{\mathbf{\zeta}},d,d^{\prime}]\) and the weights \(\lambda_{\chi}=1\) for all of these parameters except for \(\lambda_{u},\lambda_{\beta}=10\). These values are chosen to be higher because when violated, they lead to singularities in Eq. 2 causing \(\widehat{\mathbf{\Gamma}}\) to become unstable, derailing training of the model. The other loss weights are tuned to achieve low force estimation error and low residual contributions. Their values are: \(\lambda_{F}=0.5\), \(\lambda_{F_{r}}=3e-2\), \(\lambda_{AdMSE}=5e-3\), \(\lambda_{0N_{\gamma}/\partial\beta}=0.2\).

\begin{table}
\begin{tabular}{l c c c} Parameter & Norm. Range & Limits \([l_{L},l_{u}]\) & Units \\ \hline \(\phi\) & \([17,45]\) & \([0,90]\) & - \\ \(c\) & \([0,10e3]\) & \([0,\cdot]\) & Pa \\ \(\delta\) & \([11,35]\) & \([0,90]\) & deg. \\ \(c_{a}\) & \([0,10e3]\) & \([0,\cdot]\) & Pa \\ \(\gamma\) & \([14e3,22e3]\) & \([0,\cdot]\) & N/m3  \\ \(\rho\) & \([2,179]\) & \([2,178]\) & - \\ \(a\) & \([10,10]\) & \([-30,30]\) & - \\ \(w\) & \([0,3.164]\) & \([0,3.164]\) & m \\ \(d^{\prime}\) & \([0,0.3]\) & \([0,0.660]\) & m \\ \(\Delta d\) & \([-5e-2,5e-2]\) & \([-5e-2,5e-2]\) & m \\ \(Q\) & \([0,10e3]\) & \([0,\cdot]\) & N \\ \(v_{x}\) & \([0,1]\) & \([-2,2]\) & m/s \\ \(v_{x}\) & \([-1,1]\) & \([-2,2]\) & m/s \\ \(g\) & \([11.5,34.5]\) & \([11.5,34.5]\) & - \\ \(\beta N_{\gamma}/\partial\beta\) & \([-10,10]\) & \([-,\cdot]\) & unitless \\ \(\eta\) & \([2,178]\) & \([2,178]\) & - \\ \(\zeta=\phi-\delta\) & \([11,35]\) & \([0,\cdot]\) & - \\ \hline \end{tabular}
\end{table}
Table 1FEE parameter normalization ranges and limits.
Although not explicitly expressed in Eq. 13, averaging is performed over all samples in the batch for \(L^{X}_{\text{ReLU}}\), \(L^{\partial N_{Y}/\partial B}_{MAE}\), and \(L^{\Delta g}_{MSE}\). The \(L^{\nabla}_{MAE}\) loss is only averaged over samples for which the limit constraints for \(\eta\) and \(\beta\) are not violated to avoid the gradient step causing large changes to the network parameters when near a singularity of the FEE model. Similarly, \(L^{\nabla}_{MAE}\) is only averaged over samples for which the residual augmented depth of cut \(d^{\prime}>0\) to avoid penalizing the residual force from compensating for when the FEE model produces a zero-output force. All of the losses are computed using the normalized values, which helps to expedite hyperparameter tuning.

To encourage the model to estimate \(\beta\) so that it is consistent with the assumptions of the FEE model, it is limited to lie within an interval for which the minimum is expected to exist using ReLU-based physical inconsistency regularization, Equation 16. Additionally, \(\partial N_{Y}/\partial B\) is regularized to be small. It is important to note that only gradients of this loss \(L^{\partial N_{Y}/\partial B}_{MAE}\) with respect to \(\beta\) are allowed to flow, i.e. gradients with respect to the other model parameters of \(\partial N_{Y}/\partial B\) are zero. Combining these losses encourages the network to converge to the critical point \(\beta^{*}\) that is the global optimum of Eq. 7 in a _soft_ fashion, allowing for deviations from this optimum when significant force prediction accuracy improvements can be achieved. This approach is novel and expands the type of physical inconsistency loss functions that can be accounted for when developing a PINN beyond equality and inequality constraints, as outlined by Karpatne et. al. (2022).

## 3 Results & Discussion

To evaluate the approach proposed in this effort, the PINN derived parameters are compared with the ground truth parameters used in the simulation.

### Experimental Setup

Two sets of experiments using two different datasets collected from different configurations of the Vortex soil model are performed. The first dataset, \(D_{\text{FEE}}\), is collected by configuring the soil model where the heuristics discussed in Section 2.1.1 are disabled such that only the reaction forces generated by the FEE model are exerted on the blade. Additionally, particle generation is disabled, i.e. \(Q=0\), to ensure that the forces generated by particles do not confound the parameter estimation. The goal of making these changes is to collect a dataset for which accurate ground-truth knowledge of the soil parameters is available. This enables evaluation of the accuracy of the parameter estimation approach in an ideal scenario, where the physics that governs the system matches the model that is assumed.

The second dataset, \(D_{\text{default}}\), uses the default configuration for the soils, leaving the generation of particles and all of the heuristics enabled. In this case, there exists a significant mismatch between the physics governing the system and the assumed model. As discussed previously, while the simulation does not necessarily reflect the real-world physics to some quantifiable accuracy, the system exhibits similar phenomena as observed in real-world earthmoving operations. The \(D_{\text{default}}\) dataset can therefore be treated as an analog to experimentation on a real-world system.

For all experiments, the ground is assumed to be flat and level, meaning that \(\alpha=0\) and \(d=p_{\theta}^{b}\). The mass of the particles in front of the blade is obtained using a soil mass sensor yielding the surcharge force \(Q\). A simple windowed average is used to filter the velocity measurement \(\vec{v}_{t}=\text{avg}(\left[v_{t}^{x},v_{t}^{b}\right]_{t=10:t})\) to reduce oscillations in the velocity-based scaling of the soil-tool parameters. The blade angle is fixed for all experiments at \(\rho=80^{\circ}\) as is the tool width at \(w=3.164\)m. The initial values of \(\phi,c\), and \(\gamma\) are obtained at the beginning of an episode based on the initial relative density \(I_{d}\), as it is not supported to obtain the changing parameter values while running the simulation. This is not problematic in practice as the system only moves forward and does not interact with any previously disturbed terrain. The adhesion is fixed for all soil types, \(c_{\alpha}\)=200 Pa, as is the soil-tool friction angle, \(\delta=10^{\circ}\). Since the blade angle is fixed at \(\rho=80^{\circ}\), the FEE predicts zero force in the \(z\) direction for larger forward velocities. To enable observation of some non-zero forces in the \(z\) direction for \(D_{\text{FEE}}\), the default soil-tool angle is overridden to \(\delta=15^{\circ}\). All of these values are collected to form the pseudo-ground-truth FEE parameters \(\mathbf{\theta}\).

The \(D_{\text{FEE}}\) dataset consists of 4,320 sequences collected from 432 separate 10-second-long episodes using the default simulation update frequency of 60 Hz. Data is collected across all 4 default soil types and the relative density of the soil is varied across episodes, \(I_{d}\in[0,\!100]\), as are target velocity and cut depth commands, \(\hat{v}_{x}\in[0.3,\!1.0]\)m/s, \(\hat{d}_{x}\in[0.05,\!0.3]\)m. The \(D_{\text{default}}\) dataset contains 5,760 sequences from 576 episodes and is collected similarly. The PINN is implemented in PyTorch, training performed on a laptop equipped with a GPU, and the network loss in Eq. 13 is minimized using the Adam optimizer.

### Results

The model does well predicting the interaction force, achieving an average magnitude error \(|\mathbf{F}-\mathbf{\bar{F}}|\) of approximately 730N, or 11% of the measured force, on \(D_{\text{FEE}}\) and 1940N, 13% of measured, on \(D_{\text{default}}\). To put this in perspective, normalizing the error by machine limits reveals 2% error on \(D_{\text{FEE}}\) and 5% error on \(D_{\text{default}}\). This is compelling and shows that the network is doing a good job of learning the dynamics of both the simplified and complex system. To evaluate the ability of the model to identify soil parameters, parameter estimates \(\mathbf{\bar{\theta}}\) and estimated forces \(\mathbf{\bar{F}}^{\text{FEE}}\) are compared to their pseudo-ground-truth counterparts as a function of relative density. Parameter estimates from each episode within the dataset where relative density \(I_{d}\) and soil type are the same are aggregated.

\[\mathbf{\bar{\theta}}_{i}^{j}=\left[\mathbf{\bar{\theta}}^{j}\in\mathbf{D}\text{ s.t. }I_{d}=i\right]\text{for}\,j\in\mathbf{\bar{\theta}} \tag{17}\]

The mean and variance of the aggregated parameters are then computed at each \(I_{d}\) for each soil type and visually depicted in the upper half of Figs. 5 and 6. This enables direct comparison between individual parameter estimates and their corresponding pseudo-ground-truth.

Figure 4: FEE parameter estimation PINN diagram.
[https://doi.org/10.56884/ABCD1234](https://doi.org/10.56884/ABCD1234) (_ed. note: ISTVS will update your paper's DOI_) ISBN 978-1-942112:55-6 0 2023 International Society for Terrain-Vehicle Systems. All rights reserved.

Fig. 5: \(D_{\text{free}}\) dataset (top) FEE parameters estimated by the model (blue) and the pseudo-ground-truth values obtained from Vortex. (green), (bottom) Soil failure force predicted by the FEE component of the model. Units for all angles are in radians,
Figure 6: \(D_{\rm default}\) dataset (top) FEE parameters estimated by the model (blue) and the pseudo-ground-truth values obtained from Vortex. (green), (bottom) Soil failure force predicted by the FEE component of the model. Units for all angles are in radians,
For the \(D_{\text{FEE}}\), we can see that the network latches on to these true parameter values. It struggles a bit in correctly estimating \(\phi\), but the error is relatively small for smaller \(I_{d}\). However, the magnitude of the residual force \(\mathbf{F}^{r}\) also increases correspondingly, hinting that this magnitude may be a predictor of degrading parameter estimation performance. This is intuitive because as the larger \(\mathbf{F}^{r}\) becomes, the less of the force is accounted for by the FEE model. Additionally, as \(I_{d}\) increases the average depth of cut \(d\) across the dataset also decreases. This is because the vehicle is experiencing higher resistance to movement and increased stalling, which forces the blade controller to raise the blade to account for poor tracking of the forward velocity command. This indicates that the assumption of soil failure made by the FEE is increasingly being violated as the soil becomes more compact. When this assumption does not hold, then the model becomes invalid, and it becomes impossible to correctly estimate parameters from the observations alone.

This analysis is useful but does not inform the effect errors in the parameters have on the estimated failure force. Small parameter error is a sufficient but not necessary condition to achieve low error of the failure force. This is because the FEE model, Eqs. 2-10, is under constrained, meaning that there are a multitude of parameter combinations that will produce the same failure force. This presents a significant challenge, particularly as the number of unknown model parameters increases. In large part, this is why the loss function, Eqs. 13-16, contains so many regularization components. These components all constrain the estimation to encourage the network to learn the true underlying parameter values.

To complement the direct parameter error analysis and provide insight regarding the effect these parameters have on the force, the mean of the parameters at each relative density are used to produce predictions of a hypothetical soil failure force for which the depth of cut is fixed, and the surcharge is zero. This means of visualizing the results provides an alternative for evaluating the quality of the parameter estimation and illustrates the effect of the ambiguity present in the FEE model (see Figs. 5 and 6).

In general, this approach enables reasoning about the effect the estimated parameters have on the expected soil failure force, which is the relevant information autonomous soil-property aware earthmoving planners require. For the \(D_{\text{FEE}}\) dataset, the predicted force obtained from the averaged estimated parameters \(\mathbf{\hat{P}^{EE}}=f_{FEE}(\mathbf{\bar{\theta}}|\mathbf{Q}=\mathbf{0},\mathbf{d}=\mathbf{0}.\mathbf{1})\) agrees well with the predicted force obtained from the ground-truth parameters \(\mathbf{F^{FE}}=f_{FEE}(\mathbf{\theta}|\mathbf{Q}=\mathbf{0},\mathbf{d}=\mathbf{0}.\mathbf{1})\). This indicates that while the parameter estimates do not match exactly the ground-truth, the resulting predicted force is still accurate. This is important from an autonomy perspective, because if the parameter estimates are consistently biased, meaning that the system estimates the same parameter values for the same soil conditions, and the system is able to accurately predict the cutting force, then the model can be used to map soil conditions during earthmoving and utilize this map to plan efficiently.

Moving to the real-world analog dataset \(D_{\text{default}}\), the analysis becomes a bit more challenging. The parameter estimates are noisier than in the simplified dataset, and largely only \(\$\) matches with the pseudo-ground-truth values. This is expected, due to the increased complexity of the full Vortex model that results in interaction forces only partially derived from the FEE and the pseudo-ground-truth parameters. Broadly, the forces observed in this dataset are larger than those in \(D_{\text{FEE}}\), in particular \(\mathbf{F_{x}}\) is much larger.

That being said, for the non-cohesive soils, \(\phi\) is estimated well for lower compaction. For the cohesive soils, the general trend of both \(\mathbf{\hat{\phi}}\) and \(\ell\) to increase with relative density is observed. This indicates that the network has uncovered a relationship between the observational history \(\mathbf{o}_{[t:t+\tau]}\) and these properties that can be explained by \(I_{d}\). Essentially, this means that the model has implicitly learned to estimate \(I_{d}\). The predicted force plots reveal a more compelling result. As the compaction increases, a gradual rise in the predicted force is observed, mirroring the smaller changes of the pseudo-ground-truth derived forces, until \(I_{d}>60\), where the force rises rapidly and exceed the machine tractive force limits. Note that the mean depth of cut \(d\) also drops to nearly zero at \(I_{d}=80\), particularly for the non-cohesive soil, again indicating that soil failure assumption is not being met.

## 4 Nomenclature

\begin{tabular}{l l l} c & Soil cohesion & [Pa] \\ c\({}_{\text{a}}\) & Soil-tool adhesion & [Pa] \\ \(d\) & Depth of cut & [m] \\ \(F\) & Cutting force & [N] \\ N\({}_{\text{a}}\) & FEE Soil-tool adhesion coeff. & [unitless] \\ N\({}_{\text{c}}\) & FEE Soil cohesion coeff. & [unitless] \\ N\({}_{\text{r}}\) & FEE Soil wedge weight coeff. & [unitless] \\ N\({}_{\text{q}}\) & FEE Surcharge coeff. & [unitless] \\ \(Q\) & Surcharge force & [N] \\ \(R\) & Soil friction and normal force & [N] \\ \(\ddot{x}\) & Longitudinal direction & [unit vector] \\ \(W\) & Soil wedge weight & [N] \\ \(w\) & Tool width & [m] \\ \(\dot{Z}\) & Vertical direction & [unit vector] \\ \(\alpha\) & Soil surface inclination & [rad] \\ \(\beta\) & Soil failure angle & [rad] \\ \(\gamma\) & Soil moist unit weight & [N/m\({}^{2}\)] \\ \(\delta\) & Soil-tool friction angle & [rad] \\ \(\rho\) & Soil-tool angle & [rad] \\ \(\alpha_{\text{n}}\) & Normal force & [N] \\ \(\tau\) & Shear force & [N] \\ \(\phi\) & Soil internal friction angle & [rad] \\ \end{tabular}

## 5 Conclusions

In this work, a novel physics-infused neural network approach to soil property estimation approach is introduced, a simulation is developed to enable attaining ground truth values of soil parameters, and the system is evaluated across two different datasets representing two environments with significantly different physics. The approach is shown to work well when the assumed model aligns with the environment's physics and is shown to produce meaningful parameter estimates when the model and the environmental physics deviate substantially. For both datasets, the parameter estimates enable prediction of interaction forces that are informative for the planning of autonomous earthmoving operations.

While this work provides a novel framework for soil property estimation, a number of future improvements are being considered including: estimation of parameter uncertainty, development of a soil-property mapping system, extending the model to account for changes in the surface profile by varying \(\alpha\), estimation of additional parameters of the FEE model (\(\gamma,Q\), and \(\alpha_{\text{a}}\)), and modifications to the network to enable supervision of the model with only knowledge of kinematic history and without
interaction force measurements to enable more rapid integration on existing equipment.

## 6 Acknowledgements

The authors are grateful to Xavier Trudeau-Morin, Marek Teichmann, and Laszlo Kovacs from CM Labs for their support on the Vortex soil simulation and for Xavier's assistance in obtaining force measurements in the simulator.
