World Modeling for Autonomous Wheel Loaders

Abstract

This paper presents a method for learning world models for wheel loaders performing automatic loading actions on a pile of soil. Data-driven models were learned to output the resulting pile state, loaded mass, time, and work for a single loading cycle given inputs that include a heightmap of the initial pile shape and action parameters for an automatic bucket-filling controller. Long-horizon planning of sequential loading in a dynamically changing environment is thus enabled as repeated model inference. The models, consisting of deep neural networks, were trained on data from a 3D multibody dynamics simulation of over 10,000 random loading actions in gravel piles of different shapes. The accuracy and inference time for predicting the loading performance and the resulting pile state were, on average, 95% in 1.2 ms and 97% in 4.5 ms, respectively. Long-horizon predictions were found feasible over 40 sequential loading actions.

1. Introduction

The advances in artificial intelligence suggest that computerized control of construction and mining equipment has the potential to surpass that of human operators. Fully or semi-autonomous wheel loaders, not relying on experienced operators, can be an important solution to the increasing labor shortage. Wheel loaders typically operate on construction sites and quarries, repeatedly filling the bucket with soil and dumping it onto load receivers. They should move mass at a maximum rate with minimal operating cost without compromising safety. Recent research on autonomous loading control has focused on increasing performance and robustness using deep learning to adapt to soil properties [1,2,3,4,5,6,7]. However, these studies have been limited to a single bucket filling and have not considered the task of sequential loading from a pile. A challenge with sequential loading is that the pile state changes with every loading action [8,9,10]. The altered state affects the possible outcomes of the subsequent loading process and, ultimately, the total performance. A greedy strategy of always choosing the loading action that maximizes the performance for a single loading might be sub-optimal over a longer horizon. End-to-end optimization thus requires the ability to predict the cumulative effect of loading actions over a sequence of tasks. This involves having a model of the world and how it changes under the selected actions.

This paper introduces wheel loader world models for predicting the outcome of a loading action given the shape of the pile. The outcome includes the loaded mass, loading time, and work, as well as the resulting pile shape. The net outcome of a sequence of loading actions can then be predicted by repeated inferences of the model on the pile state, thus predicting its sequential evolution as well. The model aims at supporting optimal planning for autonomous wheel loaders, with the best sequence of loading actions computed from the initial pile surface only. We imagine that the plan would be updated with some regularity (e.g., daily, hourly, or after each loading cycle) from a new observation of the pile surface. Depending on the planning horizon and dimensionality of the action space, the optimizer requires numerous evaluations of the world model in a short time. Although a simulator based on multibody dynamics and the discrete element method (DEM) can predict the outcome of a loading [11], it is far too computationally intensive and slow for this optimization problem. Instead, we explore using a simulator to produce ground truth data and trained deep neural networks for “instantaneous” prediction of the outcome of particular loading actions on piles of different shapes. The learned model is informed specifically of the physics of the particular wheel loader and the type of soil used in the simulator.

The main contribution of this paper is a methodology for learning wheel loader world models and an investigation of how sensitive the model accuracy, inference speed, and required amount of training data are to the model complexity. The models predict the loading outcome, which includes the resulting pile state, loaded mass, time, and work, given the pile’s initial state and a selected loading action. The net performance of a sequence of loading actions can thus be predicted by repeated model inference. The world models are shown to be differentiable with respect to the action parameters. This enables the use of gradient-based optimization algorithms for the planning of action sequences. The paper ends with an analysis of the models’ computational properties, leaving the solution of the associated optimization problem to future work.

2. Related Work

Previous scientific work on data-driven models for predicting bucket-soil interactions have studied this from the perspective of automatic control [12,13,14,15] and bucket motion planning [16,17]. Model predictive control of earth-moving operations relies on a model to predict the dig force or soil displacement given a control signal acting over some time horizon and the current (and historic) system state [12,13]. A variational autoencoder (VAE) [15] and convolutional autoencoder [14] were used for learning a reduced representation of the terrain surface and combined with recurrent long short-term memory and mixture density neural networks to learn the time-evolution of the surface during an earth-moving task. Unfortunately, the errors grow during rollout, unless intermediate observations are added, making the resulting state unreliable and the method unstable for long horizons. Other researchers have developed models that predict the bucket fill factor and accumulated work from a planned bucket trajectory and initial soil pile shape [16,17]. This approach is infeasible in practical use-cases, where it is not possible to track a prescribed bucket trajectory [18]. In the present paper, we draw inspiration from the idea of using a VAE for representing the local soil surface but we focus on accurately predicting the end-state after breakout and on the accumulated mass, time, and work of the loading cycle. Instead of using a parameterized bucket trajectory as input, we use the control parameters of a force-based automatic bucket-filling controller.

3. Methodology

The method of developing a wheel loader world model includes creating a simulator [11,19] for the particular machine and its environment. The simulator is based on contacting 3D multibody dynamics with a real-time deformable terrain model that has been shown to produce digging forces and soil displacements with an accuracy close to that of resolved discrete elements and coupled multibody dynamics [20]. The simulated wheel loader is equipped with a type of admittance controller [21] for automatic bucket-filling, parameterized by four action parameters that determine how the boom and bucket actuation respond to the current dig force. An annotated dataset is created from simulated loading cycles carried out on soil piles of different shapes using different combinations of control parameters. Each simulation results in one data point, consisting of the pile’s local heightmap before and after loading, loaded mass, time and work, and the set control parameters that were used. Models are trained to predict how the loading outcome depends on the input values. The models are finally tested on validation data withheld during training. The validated models can then be used to predict the result of new loading operations or be embedded in optimization routines for finding sequences of loading actions that are optimal with respect to productivity or energy efficiency. The method is illustrated in Figure 1.

The learned models were investigated for accuracy, inference speed, and the required amount of training data. In general, model accuracy can be improved by increasing the input and output dimensions and the number of internal model parameters. However, larger models normally require more training data and may have lower inference speeds. To quantify the difference, we developed both a high-dimensional and a low-dimensional model using different representations of the pile state. The high-dimensional model takes a well-resolved heightmap of the pile state as input. This may capture many details of the pile surface at the cost of additional model parameters associated with convolutional layers. The low-dimensional model takes a heavily reduced representation of the pile surface, involving only four scalar parameters for its slope and curvature.

4. Wheel Loader World Models

4.1. Preface and Nomenclature

This paper defines the loading cycle as starting with a machine heading in direction $\mathit{t}$ to dig location $\mathit{x}$ in a pile, represented by the current pile state H. The cycle ends when the machine has returned to its initial location after digging in, filling the bucket, breaking out, and leaving the pile in a new state $H^{\prime }$. The wheel loader is equipped with an automatic bucket filling controller that is parametrized by some set of action parameters $\mathit{a}$. The controller is engaged when the machine reaches the dig location $\mathit{x}$. Each loading cycle can be assigned a performance $\mathcal{P}$, which in this study includes the loaded mass, loading time, and mechanical work. The performance is a consequence of the dynamics of the machine and the soil under the selected action $\mathit{a}$ and the initial state of the pile. The expected performance of a loading cycle, indexed by $n\in \mathbb{N}$, is therefore given by some unknown performance predictor function $\mathbf{\Psi }$. In other words,

$${\widehat{\mathcal{P}}}_n=\mathbf{\Psi }(H_n,{\mathit{x}}_n,{\mathit{t}}_n,{\mathit{a}}_n),$$

(1)

where we use the hat to distinguish predictions from actual values. The net effect of a sequence of N loading cycles is the accumulated outcome of a sequence of loading actions, ${\mathit{a}}_1,\dots ,{\mathit{a}}_N$, at pose $({\mathit{x}}_1,{\mathit{t}}_1),\dots ,({\mathit{x}}_N,{\mathit{t}}_N)$ applied on a sequence of pile states $H_1,\dots ,H_N$. Each loading cycle transforms the pile from its previous state to the next, $H_n\to H_n^{\prime }\equiv H_{n+1}$, according to some unknown pile state predictor function $\mathbf{\Phi }$:

$${\widehat{H}}_n^{\prime }=\mathbf{\Phi }(H_n,{\mathit{x}}_n,{\mathit{t}}_n,{\mathit{a}}_n).$$

(2)

The net outcome of N sequential loading actions is thus given by the sum

$${\widehat{\mathcal{P}}}_{1:N}\equiv \sum _{n=1}^N{\widehat{\mathcal{P}}}_n({\widehat{H}}_{n-1}^{\prime },{\mathit{x}}_n,{\mathit{t}}_n,{\mathit{a}}_n)$$

(3)

with initial pile state $H_0^{\prime }=H_1$. The predictions are associated with some error, that accumulates over repeated loadings from the evolving pile. The accumulated error in the pile state and loading performance over a horizon of N cycles is

$$\begin{array}{cc}\hfill {\mathcal{E}}_{1:N}^{{\widehat{H}}^{\prime }}\equiv & ∥{\widehat{H}}_N^{\prime }-H_N^{\prime }∥,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\mathcal{E}}_{1:N}^{\widehat{\mathcal{P}}}\equiv & \sum _{n=1}^N∥{\widehat{\mathcal{P}}}_n-{\mathcal{P}}_n∥.\hfill \end{array}$$

(5)

When solving the problem of finding the sequence of dig locations and loading actions that maximize the net performance ${\widehat{\mathcal{P}}}_{1:N}$, it is beneficial to use gradient-based optimization methods. These require the sought function $\mathbf{\Psi }$ to be differentiable with respect to $\mathit{a}$.

4.2. Global and Local Pile State

We represent the global pile state and the surrounding terrain by a height surface function $z=H(x,y)$ in Cartesian coordinates. When making sequential predictions, the aim is to accurately predict the global state of the pile, which may take any shape consistent with the physics of the soil. To simplify matters, the predictor models take only the local pile state as input, assuming that the outcome of a loading depends only on the local state and not on the global shape of the pile. A dig location, $\mathit{x}=[x,y]$ and heading $\mathit{t}$ defines a local frame with basis vectors $\{{\mathit{e}}_x,{\mathit{e}}_y,{\mathit{e}}_z\}$, where ${\mathit{e}}_x$ is aligned with the dig direction, and ${\mathit{e}}_z$ is the vertical direction aligned with the gravitational field. We represent the local pile state with a discrete heightmap of $I\times J$ grid cells with the side length $\Delta l$. The heightmap, in the local frame, is

$$h_{ij}=H(\mathit{x}-\delta {\mathit{e}}_x+i\Delta l{\mathit{e}}_x+j\Delta l{\mathit{e}}_y)-H(\mathit{x}-\delta {\mathit{e}}_x)$$

(6)

with integers $i\in [0,I]$ and $j\in [-J/2,J/2]$ and constant displacement $\delta$ that ensures a certain fraction of the ground in front of the pile is present in the heightmap. The length and width of the heighmap, $I\Delta l$ and $J\Delta l$, must be large enough to cover the interaction region; see Section 6.3.

The simplified problem of predicting the loading outcome, Equations (1)–(3), using the local heightmap is

$$\begin{array}{cc}\hfill {\widehat{\mathcal{P}}}_n& =\psi ({\mathit{h}}_n,{\mathit{a}}_n),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\widehat{\mathit{h}}}_n^{\prime }& =\varphi ({\mathit{h}}_n,{\mathit{a}}_n),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\widehat{\mathcal{P}}}_{1:N}& \equiv \sum _{n=1}^N{\widehat{\mathcal{P}}}_n({\widehat{\mathit{h}}}_{n-1}^{\prime },{\mathit{a}}_n).\hfill \end{array}$$

(9)

with local predictor functions $\psi$ and $\varphi$ for the performance and pile state. The computational process is described in Algorithm 1 and illustrated in Figure 1.

Algorithm 1: Long-horizon prediction using world models

4.3. Local Pile Characteristics

As an alternative to representing the local pile state in terms of a heightmap, we introduce a low-dimensional characterization, $\tilde{\mathit{h}}\equiv [\theta ,\alpha ,{\kappa }_x,{\kappa }_y]$, in terms of four scalar quantities: slope angle $\theta$, incidence angle $\alpha$, longitudinal curvature ${\kappa }_x$, and lateral curvature ${\kappa }_y$. These aim to capture the essence of the local pile shape from the perspective of bucket-filling [22,23]. First, the mean unit normal $\mathit{n}$ is computed from the local heightmap. The slope angle relative to the horizontal plane is then computed $\theta \equiv {cos}^{-1}\left[\mathit{n}·{\mathit{e}}_z\right]$. The incidence angle is the angle between the attack direction ${\mathit{e}}_x=\mathit{t}/∥\mathit{t}∥$ and the pile normal projected on the horizontal plane ${\mathit{n}}_{\perp }=\mathit{n}-(\mathit{n}·{\mathit{e}}_z){\mathit{e}}_z$, that is, $\alpha \equiv {cos}^{-1}(\left[{\mathit{n}}_{\perp }·{\mathit{e}}_x\right]/∥{\mathit{n}}_{\perp }∥)$. See Figure 2 for an illustration.

Taking inspiration from [24], we calculate the local mean curvature in the ${\mathit{e}}_x$ and ${\mathit{e}}_y$ directions by fitting a quadratic surface $h\left(\mathit{x}\right)\approx b+{\mathit{c}}^{\mathrm{T}}\mathit{x}+\frac{1}{2}{\mathit{x}}^{\mathrm{T}}\mathit{Q}\mathit{x}$ with surface parameters $\mathit{Q}=\mathrm{diag}[-{\kappa }_x,-{\kappa }_y]$, $\mathit{c}\in {\mathbb{R}}^2$, and $b\in \mathbb{R}$. This sign convention makes the curvatures positive for convex pile shape, which is the recommended shape for high-performance loading [23].

4.4. Performance Predictor Model

We developed two different models to predict the loading performance, referred to as the high-dimensional model and the low-dimensional model, respectively. The high-dimensional model, ${\psi }^{\mathrm{high}}(\mathit{h},\mathit{a})$, takes the local heightmap $\mathit{h}$ as input, while the low-dimensional model, ${\psi }^{\mathrm{low}}(\tilde{\mathit{h}},\mathit{a})$, uses the local pile characteristics $\tilde{\mathit{h}}$. Both models take the same action parameters $\mathit{a}$ as input, and they both output a performance vector $\widehat{\mathcal{P}}$. The name distinction comes from $\dim \left(\mathit{h}\right)\gg \dim \left(\tilde{\mathit{h}}\right)$. The difference is reflected in the different neural network architectures of the models (Figure 3). The high-dimensional model uses three convolutional layers and a fully connected linear layer to encode $\mathit{h}$ into a latent vector of length 32 before it is concatenated with $\mathit{a}$. $\mathit{h}$ is interpolated to a size of $32\times 32$ before encoding. The convolutions use ten filters of $3\times 3$ kernels and unit stride with zero-padding. They are followed by batch normalization and max pooling with window size $2\times 2$. The activation function is subject to hyperparameter tuning. The concatenation steps in both ${\psi }^{\mathrm{high}}$ and ${\psi }^{\mathrm{low}}$ are followed by multilayer perceptrons (MLPs) of identical architectures that are also subject to hyperparameter tuning.

4.5. Pile State Predictor Model

The pile state predictor model combines a VAE architecture with an MLP to make predictions via three steps: First, an encoder network compresses the initial pile state $\mathit{h}\in {\mathbb{R}}^{52\times 52}$ into a lower-dimensional, regularized latent representation $\mathit{z}\in {\mathbb{R}}^{64}$. Next, given $\mathit{z}$ along with the loading action $\mathit{a}$ and a scale factor $\Delta h\equiv max\left(\mathit{h}\right)-min\left(\mathit{h}\right)$, an MLP predicts the regularized latent representation of the new pile state ${\mathit{z}}^{\prime }$. Last, a decoder network constructs the resulting pile state ${\widehat{\mathit{h}}}^{\prime }\in {\mathbb{R}}^{52\times 52}$ from ${\mathit{z}}^{\prime }$. Figure 4 illustrates the inference process. For the VAE encoder/decoder blocks, we use the same CNN architecture as in [25].

The MLP block consists of two hidden layers with 1024 nodes and uses Leaky ReLU activation. Note that $\mathit{h}$ is interpolated to size $64\times 64$ before encoding and back to $52\times 52$ after decoding to fit the off-the-shelf architecture.

4.5.1. Post-Processing

VAEs are known to produce somewhat blurry images [26]. In our case, this led to a loss of detail in the decoded pile states. In our preliminary tests, we found that we could improve the prediction accuracy on average by reusing some information from the initial pile state. By interpolating between $\mathit{h}$ and ${\widehat{\mathit{h}}}^{\prime }$ along the edges, the details of the outer area are retained. This makes the predicted local pile states blend seamlessly into the global pile state after the substitution. Figure 5 illustrates the post-processing.

4.6. Low-Dimensional Pile State Prediction Using Cellular Automata

The low-dimensional predictor model did not have direct access to the local heightmap. Instead, we constructed a pile state predictor method that acts directly on the global pile state. It works in two stages. First, the predicted load mass (part of the performance vector ${\widehat{\mathcal{P}}}_n$) is removed from the global pile state at the dig location. This is performed by computing the corresponding volume and subtracting this from the global pile along a strip as wide as the bucket, starting at ${\mathit{x}}_n$ and stretching along ${\mathit{t}}_n$. In the second stage, soil mass is redistributed to eliminate any slopes steeper than the set angle of repose. The mass-preserving cellular automata algorithm in [27] was employed, with the “velocity of flowing matter” set to $z_{+}=0.2\Delta l$.

4.7. Delimitations

The predictor models developed in this paper are limited to a single type of soil, gravel, which was assumed to be homogeneous. However, the process of learning world models is applicable to any other soil type supported by the simulator. This paper does not consider soil spillage on the ground after breaking out with an overfilled bucket. Spillage would affect subsequent performance, either by causing a loss in control precision or traction or by requiring the ground to be cleared. In this paper, “loading” refers to the bucket-filling phase, including approaching the pile straight ahead and reversing in the opposite direction after the breakout, following the definition in [16]. To account for the full loading performance, one should incorporate the V-cycle maneuver, emptying the bucket into the load receiver and clearing spillage from the ground.

5. Simulator and Loading Controller

5.1. Simulator

We collected data using a simulator developed [11] and validated [19] in a previous work. Images and videos from the simulator are found in Figure 1, Figure 6 and Supplementary Video S1. The simulator combines a deformable terrain model with a wheel loader model and runs approximately in real time. The terrain model, introduced in [20], combines the representation of soil as a continuous solid, distinct particles, and rigid multibodies. When a digging tool comes in contact with the terrain, the active zone of moving soil is predicted and resolved in terms of particles using DEM. When particles come to rest on the terrain surface, they merge back into the solid soil. The computational processes preserve the total mass. The wheel loader is modeled as a rigid multibody system, matching roughly the dimensions and physical properties of a Komatsu WA320-7. The model has five actuated joints, with hinge motors powering the driveline and steering and linear motors for the boom and bucket cylinders. The motors are controlled by specifying a momentaneous joint target speed and force limits. The actuators will run at the set target speed only if the system dynamics and required constraint force do not exceed their respective force limits. The limits are set to match the strength of the real machine [28]. The transmission driveline model includes the main, front, and rear differential couplings to the wheels. The force interaction between the vehicle and the terrain occurs through the tire-terrain and bucket-terrain contacts. The resulting forces on the bucket depend on the shape and amount of active soil, its dynamic state, and mechanical properties, including bulk density, internal friction angle, cohesion, and dilatancy angle. As in [19], these are set to the values 1727 kg/m³, ${32}^{\circ }$, 0 kPa, and $8^{\circ }$, which are intended to reflect the properties of dry gravel. The simulations were run with the physics engine AGX Dynamics [29] using a $0.01$ s time-step and terrain grid cell-size of 0.1 m × 0.1 m.

5.2. Loading Controller

The loading scenario starts with the machine driving at 8 km/h towards a pile 5 m from the target dig location. The bucket is lowered and held level to the ground during the approach. The target drive velocity is kept constant during the bucket-filling phase, although the actual velocity will vary due to digging resistance. Throughout the cycle, the machine maintains the same heading. When the bucket reaches the set dig location, the automatic bucket-filling controller is engaged. This controller was inspired by the admittance controller from [21], which regulates the bucket actuator velocity using the measured boom cylinder force. Our admittance controller uses the same method but applies it to both the boom and bucket actuators. The controller determines the target velocities of the boom and bucket cylinders, $v_{\mathrm{bm}}$ and $v_{\mathrm{bk}}$, as follows:

$$\begin{array}{cc}\hfill v_{\mathrm{bm}}& =\mathrm{clip}\left(k_{\mathrm{bm}}\left[f_{\mathrm{bm}}/f_{\mathrm{b}0}-{\delta }_{\mathrm{bm}}\right],0,1\right)v_{\mathrm{bm}}^{\max },\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill v_{\mathrm{bk}}& =\mathrm{clip}\left(k_{\mathrm{bk}}\left[f_{\mathrm{bm}}/f_{\mathrm{b}0}-{\delta }_{\mathrm{bk}}\right],0,1\right)v_{\mathrm{bk}}^{\max },\hfill \end{array}$$

(11)

where $f_{\mathrm{bm}}$ is the momentaneous boom cylinder force and $\mathrm{clip}(value,min,max)$ limits $value$ to the maximum and minimum values. The ramp function has four parameters (velocity gains $k_{\mathrm{bm}}$ and $k_{\mathrm{bk}}$ and force thresholds ${\delta }_{\mathrm{bm}}$ and ${\delta }_{\mathrm{bk}}$) that parameterize the behavior of the controller. The actuator max speeds, $v_{\mathrm{bm}}^{\max }$ and $v_{\mathrm{bk}}^{\max }$, and the normalizing boom force $f_{\mathrm{b}0}$ are set using specifications from the manufacturer. The control parameters are collected in the action vector

$$\mathit{a}=\left[{\delta }_{\mathrm{bm}},k_{\mathrm{bm}},{\delta }_{\mathrm{bk}},k_{\mathrm{bk}}\right].$$

(12)

Parameters ${\delta }_{\mathrm{bm}}$ and ${\delta }_{\mathrm{bk}}$ regulate what force magnitude is required to trigger the lift and tilt reactions, while parameters $k_{\mathrm{bm}}$ and $k_{\mathrm{bk}}$ control how rapid the reaction is. Different values of the control parameter thus render different scooping motions, as illustrated by the examples in Figure 7. The challenge is to select the parameters most appropriate for the local pile state. Note that the controller operates only on the lift and tilt actuators. The vehicle keeps thrusting forward to reach the set target drive velocity. If there is insufficient traction, the wheels may slip.

The bucket-filling controller is stopped if the bucket achieves the tilt end position, reaches a penetration depth of $3.2$ m, breaks out of the initial surface, or if the loading duration exceeds 15 s. The bucket is then tilted with maximum speed to its end position (if not there already). The brake is applied for at least $0.5$ s to let the agitated soil settle. After that, the vehicle is driven in reverse with a target speed of 8 km/h, lift $v_{\mathrm{bm}}=0.6v_{\mathrm{bm}}^{\max }$ and tilt $v_{\mathrm{bk}}=0.6v_{\mathrm{bk}}^{\max }$ to reach a boom angle of $-{20}^{\circ }$ relative to the horizontal axis, and the bucket tilt ends. The scenario ends when the vehicle has reversed 5 m from the dig location. To simplify time and energy comparison between simulations, they all use the same start and end distance from the dig location and all end with identical boom and bucket angles. The target speeds and force measurements are smoothed using a $0.1$ s moving average to avoid a jerky motion.

6. Data Collection and Model Training

This section describes how data were collected from the simulator and used for learning and evaluating the predictor models.

6.1. Loading Outcome

The models were developed to predict the outcome of a loading cycle in terms of the performance $\mathcal{P}$ and the resulting local pile state ${\mathit{h}}^{\prime }$ given the initial pile $\mathit{h}$ and action $\mathit{a}$. The performance was measured using the three essential scalar metrics, namely load mass (ton), loading time (s), and work (kJ). Hence, $\mathcal{P}\in {\mathbb{R}}^3$. The load mass is measured as the amount of soil the bucket carries at the end of each loading cycle. The loading time is the time elapsed between each loading cycle’s start and end. The work is the energy consumed by the boom and bucket actuators and the forward drive. It includes the energy required to fill the bucket, break out, raise the bucket, and accelerate the vehicle and soil. Much of the work is lost to frictional dissipation internally in the soil and at the bucket–soil interface. The physics-based simulator accounts for this dissipation. Energy dissipation in the vehicle engine and hydraulics is not included and should be added if needed.

6.2. Data Collection

We collected a total of $10,718$ samples, ${\{{\mathit{h}}_n,{\mathit{a}}_n,{\mathit{h}}_n^{\prime },{\mathcal{P}}_n\}}_{n=1}^{\mathrm{10,718}}$, by repeating the loading scenario (Section 5.2) in the simulation. First, six different initial piles were prepared: two triangular, two conical, and two wedged with respective slope angles of ${20}^{\circ }$ and ${30}^{\circ }$. These are illustrated in Figure 8. Perlin noise [30] was applied to the surfaces to increase the variability in pile shape. For each of the six initial piles, 30 consecutive loading cycles were simulated. Each loading cycle produced one data point $({\mathit{h}}_n,{\mathit{a}}_n,{\mathit{h}}_n^{\prime },{\mathcal{P}}_n)$. Since the resulting pile state was used as the initial state in the next loading cycle, a variety of piles of different shapes was achieved. This process was repeated 60 times for each of the six seed piles.

In each simulation, a random dig location was selected, and the wheel loader was positioned at ${\mathit{x}}_n$, a distance of 5 m from the pile, and the heading was given a random disturbance of $\pm {10}^{\circ }$; see Figure 8. Each loading used a set of action parameters, ${\mathit{a}}_n$, randomly sampled by Latin Hypercube Sampling (LHS) with the ranges $0.0\le {\delta }_{\mathrm{bm}},{\delta }_{\mathrm{bk}}\le 0.7$ and $0.0\le k_{\mathrm{bm}},k_{\mathrm{bk}}\le 5.0$. Note that the actual slope angle $\theta$ and incidence angle $\alpha$ varied according to the dig location because of the Perlin noise.

6.3. Local Heightmap Settings

We found that the performance and pile state predictor models benefited from using different sizes of the local heightmap and displacement $\delta$ in the dig direction. We used a $3.6$ m sided quadratic heightmap discretized by $36\times 36$ grid cells for the performance predictor. The size leaves a margin of $0.5$ m on the sides of the bucket. Due to avalanching, the pile state predictor needed a larger heightmap to capture the state that the pile finally settles into after the breakout. For this, we used a heightmap with 5.2 m side lengths, discretized in $52\times 52$ cells, and displaced the heightmap 1.0 m towards the approach direction.

6.4. Model Training

This section describes the training processes of the performance predictor (Section 6.4.1) and the pile state predictor model (Section 6.4.2). The models were implemented using PyTorch.

6.4.1. Performance Predictor Model

The models were trained until convergence with a mean squared error (MSE) loss using the Adam optimizer using a learning rate of ${10}^{-5}$. Hyperparameter tuning was performed via grid search over the number of hidden layers (1, 2, or 3), the number of units ($2^3,2^4,\dots ,2^{12}$) in each hidden layer with 0.1 dropout rate, as well as the activation function in fully connected and convolutional layers (Leaky ReLU or Swish [31]). The dataset was first min-max normalized. The training and validation set (split 90/10) was sequentially increased in size from 100 to 9646 samples, while the test set size was fixed at 1072 samples.

6.4.2. Pile State Predictor Models

The pile state predictor model was trained in a two-stage process. First, we trained the VAE to perform heightmap reconstruction. Before training, the heightmaps were re-scaled by subtracting the average slope and applying min-max scaling to them individually:

$$\begin{array}{cc}\hfill {\mathit{h}}_n& :={\mathit{h}}_n-\frac{1}{N}\sum _{m=1}^N{\mathit{h}}_m\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\mathit{h}}_n& :=\frac{{\mathit{h}}_n-min\left({\mathit{h}}_n\right)}{max\left({\mathit{h}}_n\right)-min{\mathit{h}}_n}.\hfill \end{array}$$

(14)

The standardization in Equation (13) centers each height distribution around 0 across the entire data set before applying the min-max scaling in Equation (14). We found this helps to preserve the pile–ground boundary in our reconstructions. The normalization in Equation (14) makes the VAE agnostic to scale, which we found increases its overall performance. The VAE was trained using the Adam optimizer with a learning rate of ${10}^{-3}$. We used a weighted sum of the element-wise MSE and the Kullback–Leibler Divergence (KLD), ${\mathcal{L}}_{\mathrm{VAE}}={\mathcal{L}}_{\mathrm{MSE}}+0.1{\mathcal{L}}_{\mathrm{KLD}}$ for the loss function.

In the second training phase, we used the encoder trained in the first phase to encode sampled pile state pairs $(\mathit{h},{\mathit{h}}^{\prime })$ into latent pairs $(\mathit{z},{\mathit{z}}^{\prime })$. We then trained an MLP on the latent pairs with a 0.1 dropout rate for the hidden layers and used the Adam optimizer with a learning rate ${10}^{-5}$.

The full dataset was used with a split ratio of 80/10/10 in training, validation, and test data for the pile state predictor. Since the wheel loader geometry and actions are left-right symmetric, we applied random reflections to the heightmaps to augment the dataset.

7. Results

The models for predicting the loading performance and resulting pile state were evaluated separately on the hold-out test data. The best models were selected and used in combination to test their ability to predict the outcome of sequential loadings. An overview of the full dataset is given in Appendix A.

7.1. Performance Predictor Model

In total, we developed 480 model instances during the hyperparameter tuning for both the high- and low-dimensional performance predictors, ${\psi }^{\mathrm{high}}$ and ${\psi }^{\mathrm{low}}$. The models were evaluated using three metrics: mean relative error (MRE), training time, and inference speed. The MRE is relative to the simulation ground truth and is calculated from the average of five distinct training results. The training time was counted as the time per epoch (number of epochs ranging up to about 2000), calculating the median over the entire training. The inference speed is measured as the model execution time, taking the average of 1000 model executions. We selected the MLPs with two hidden layers with 2048 units as models of special interest. These are denoted ${\psi }_{\diamond }^{\mathrm{high}}$ and ${\psi }_{★}^{\mathrm{high}}$, where the difference is the use of Swish and Leaky ReLU for activation functions. These models have roughly ${10}^7$ model parameters, including the convolutional layers. The loading performance MRE, listed in

Table 1. Key properties of the selected performance predictor models.

	${\mathit{\psi }}^{\mathbf{high}}$		${\mathit{\psi }}^{\mathbf{low}}$
	L ReLU	Swish	L ReLU	Swish
mass MRE [%]	$4.47$	$5.70$	7.77	7.58
time MRE [%]	$3.61$	$4.43$	5.66	5.41
work MRE [%]	$5.90$	$6.98$	8.51	8.47
training time [s]	1.9	1.92	$0.27$	$0.27$
inference time [ms]	1.13	1.27	$0.31$	$0.36$

, ranges between 3.5 and 7% with insignificant dependence on the activation function. For the low-dimensional model, the selected models of special interest, ${\psi }_{\diamond }^{\mathrm{low}}$ and ${\psi }_{★}^{\mathrm{low}}$, are smaller with two hidden layers and 512 units per layer. This amounts to roughly $5\times {10}^5$ parameters (no convolutional network is involved). It was trained on a smaller dataset with only 3000 samples, as adding additional samples did not significantly improve the performance. The loading performance MRE for the low-dimensional model saturates in the range of about 5.5–8.5%, in other words, with 20–70% larger errors than for the high-dimensional model. The training time and inference time are found in Table 1. The low-dimensional model trains about seven times faster, and inference runs three times faster than the selected high-dimensional model. The computational overhead of using Swish over Leaky ReLU is marginal. The high- and low-dimensional models have a memory footprint of 16.7 Mb and 0.027 MB, respectively. The effect of hyperparameters on the model performance is found in Appendix B.

It is interesting to see how the errors are distributed in the space of the predictions (Figure 9). We observe that the relative error (RE) remains comparatively small for high load masses, while it increases for the smallest load masses. This suggests that the model is more reliable for high-performing loading actions, which is the general goal, than for low-performing ones.

It is also interesting to understand why the model fails occasionally. Therefore, we identified the test samples with the ten worst load mass predictions. These are displayed in Figure 10 with the local heightmaps. The common factors are (a) that the loaded mass is small in the ground truth samples while overestimated by the model and (b) that the heightmap is skewed, with most mass distributed on either the left or right side. It is understandable that the low-dimensional model, with only four parameters characterizing the pile state, has more difficulty with these piles as it cannot distinguish between uniform and irregular pile surfaces. That makes accurate load mass prediction difficult for complex pile surfaces.

7.2. Pile State Predictor Model

The developed pile state predictor model is evaluated by the mean absolute error (MAE) and mean relative error (MRE) of a prediction ${\widehat{\mathit{h}}}^{\prime }$ compared to the simulated ground truth heightmap ${\mathit{h}}^{\prime }$. MAE and MRE are calculated in terms of the volume difference of the enclosing surfaces. In detail, the MRE is computed as

$$\mathrm{MRE}=\frac{1}{N}\sum _n^N{\epsilon }_n^{\mathrm{rel}}=\frac{1}{N}\sum _n^N\frac{{\sum }_{i,j}|{\widehat{h}}_{ij}^{\prime n}-h_{ij}^{\prime n}|\Delta l^2}{{\sum }_{i,j}\left|h_{ij}^{\prime n}\right|\Delta l^2},$$

(15)

where the volume of each cell, indexed by $ij$, is calculated between the surface and the extended ground plane in the grid. The MAE is calculated the same way, excluding the volume normalization. The results are summarized in

Table 2. Test results of the pile state predictor model, with and without the post-processing (p-p).

p-p	MAE [m3]	MRE [%]	Inference Time [ms]
on	0.75	3.04	4.51
off	0.84	3.41	4.48

. The MAE of misplaced volume constitutes roughly 25% of the bucket capacity but only 3% of the volume under the local heightmap. The post-processing step (Section 4.5.1) reduces the MAE from 0.84 to 0.75 m³, with insignificant overhead in inference time. The memory footprint of the pile state predictor model (Section 4.5) was 10.7 MB, with 6.1 MB and 4.6 MB for the VAE and MLP parts, respectively.

The improvement by the post-processing can also be seen from the distribution of the errors in Figure 11.

Figure 12 visualizes five selected test predictions, marked by diamonds in Figure 11b. Overall, the predicted pile surfaces look reasonable compared to the ground truth. We find that the model performance tends to be better when the initial pile has a smooth, regular shape, which is to be expected as the VAE decoder tends to produce smooth outputs. Examples of the smoothness issue can be seen in the reconstruction column in Figure 12.

Post-processing is most beneficial in cases where the initial state is irregular and when there is little change at the boundary (Figure 12e). When the initial pile is smooth with a uniform slope close to the angle of repose, loading will induce avalanches that affect the pile state at the rear boundary. This explains why the AE is sometimes increased by post-processing (Figure 12d).

7.3. Sequential Loading Predictions

This section compares the evolution of 40 sequential loading cycles using the high and low-dimensional world models to the simulation ground truth. The simulation and the predictor models start with the same initial pile state, $H_1$, and use identical sequences of the dig location, heading, and action parameters, ${\{{\mathit{x}}_n,{\mathit{t}}_n,{\mathit{a}}_n\}}_{n=1}^{40}$. The initial pile is the one shown in Figure 1. The selected headings and dig locations are shown in Figure 13. The actions were picked randomly from the training data set. The simulated sequential loading can be seen in Supplementary Video S2. The performance predictor ${\psi }_{\diamond }^{\mathrm{high}}$ is used together with the pile state predictor $\varphi$ for high-dimensional model predictions. The low-dimensional performance predictor ${\psi }_{\diamond }^{\mathrm{low}}$ is combined with the cellular automata model to predict the next pile state.

The evolution of the pile shape is shown in Figure 13 and in Supplementary Video S2. The accumulated load mass, time, work, and residual pile volume, evolving with the number of loadings n, are summarized in

Table 3. The accumulated load mass, time, work, and residual pile volume during sequential loading. The number of loading cycles is denoted by n, and GT stands for ground truth.

	${\widehat{\mathcal{P}}}_{1:\mathit{N}}$									${\widehat{\mathit{H}}}_{\mathit{N}}^{\prime }$
	Load Mass [tonne]			Loading Time [s]			Work [MJ]			Pile Volume [m3]
N	GT	Ψhigh	Ψlow	GT	Ψhigh	Ψlow	GT	Ψhigh	Ψlow	GT	Φ	cell.aut.
5	14.7	15.7	16.1	55.8	53.2	54.5	2.0	2.0	2.0	768.6	765.9	767.7
10	28.2	31.5	31.3	106.1	101.0	101.9	4.1	3.9	3.8	760.7	755.9	758.8
15	40.8	45.3	47.0	161.1	155.4	154.5	6.5	6.1	6.0	753.3	745.6	749.7
20	56.8	62.0	64.4	218.0	210.8	208.1	8.9	8.5	8.3	744.0	728.7	739.5
25	71.0	75.4	80.6	273.7	265.7	260.2	11.1	10.6	10.4	735.7	714.7	730.1
30	86.1	90.5	96.6	327.3	318.7	308.9	13.2	12.8	12.4	726.9	700.1	720.7
35	101.4	106.3	114.3	389.4	385.1	369.9	15.7	15.3	14.6	718.0	679.1	710.4
40	116.0	118.6	130.9	443.0	441.9	418.8	17.9	17.4	16.6	709.5	664.2	700.6

. The error accumulation over time is shown in Figure 14. While the pile state prediction error for the low-dimensional model grows at a nearly constant rate, the growth rate of the high-dimensional model is more irregular, starting at a lower rate but growing more rapidly after 15 sequential loadings. The error in loaded mass grows nearly linearly and roughly at the same rate for both models. Also, the errors in loading time and work grow mainly linearly. Both models follow the same trend.

The high-dimensional model uses $\mathbf{\Phi }$ for the local pile state prediction, making no change outside the local heightmap, as discussed in Section 4.5. The accumulated error eventually leads to a predicted pile state with a slope steeper than the angle of repose (Figure 13).

The high- and low-dimensional models’ accumulated inference time was also measured over the 40 loading cycles. The high-dimensional model took 0.3 s in total. The low-dimensional model took 0.01 s. When the cellular automata are included, this amounts to 44.3 s, but it should be noted that the implementation of the cellular automata was not optimized or adapted for GPU computing.

To summarize, the high-dimensional model was better for shorter time horizons, but the low-dimensional model was more stable over longer horizons. The stability might be the same if the high-dimensional model used the cellular automata in place of the pile state predictor, but this would come with an additional computational cost.

7.4. Diggability Map

The predictor model can be used to search for the most favorable dig location around the pile. We demonstrate this by creating diggability maps, where the model has been inferred around the entire edge of a pile. Figure 15 shows a map created by using ${\mathbf{\Psi }}_{\diamond }^{\mathrm{high}}$ to predict the performance using a fixed action $\mathit{a}=[0.68,4.51,0.17,4.46]$ at 150 dig locations along the local normal direction. The map highlights the regions with higher pile angles where the loaded volume would be relatively high and would be dug with lower energy cost. Low-performance regions are visible where the pile bulges. The diggability maps can be further improved by searching also for the locally optimal heading and loading action.

8. Discussion

8.1. Feasibility

The relatively small errors of single loading prediction, 3–8.5%, are encouraging. There are, however, several practical limitations with the the method in its current form. The model is particular to the simulator’s specific wheel loader and soil, which was homogeneous gravel. Consequently, a new model needs to be trained for each combination of wheel loader model and soil type unless this is made part of the model input. There is a greater need for automation solutions for more complex soils, such as blasted rock or highly cohesive media, than for homogeneous gravel. Although these materials can be numerically simulated, the problem of learning predictor models may be harder and require different data, as discussed for the loading process in [7]. For blasted rock, the local rock fragmentation probably needs to be included in the model input. On the other hand, the level of accuracy needed ultimately depends on the specific application and the prediction horizon required.

It is an attractive option to train on actual field data (rather than simulated data) to eliminate any simulation bias. That would entail collecting data from several thousand loadings with variations in pile shape and action parameters. This would be demanding but not unrealistic given the number of vehicles (of the same model) that are in operation worldwide. If all were equipped with the same bucket-filling controller, it would ultimately be a question of instrumenting these machines and sites to scan the pile shape before and after loading and tracking the vehicle’s motion and force measurements.

Long-horizon planning of sequential loading requires long rollouts of repeated model inference. The error accumulation in the pile state predictor model (Figure 14) might then become an obstacle. Although the prediction accuracy at each loading cycle is improved by the post-processing (Figure 11), it does not involve avalanching and the pile may therefore evolve steep slopes during long rollouts (Figure 12). If that is the case, the pile state predictor can be replaced with mass removal and cellular automata, as described in Section 4.6, making the cellular automata the computational bottleneck. In our implementation, the cellular automata were about 200 times slower than model inference, about 1 s versus 5 ms per loading. An optimized implementation can be an order of magnitude faster at least. An alternative is to use the cellular automata as a post-processor, forcing the predicted pile state to be consistent with the soil’s angle of repose. If the post-processing is only occasionally needed, instead of after each pile state prediction, the computational overhead may be marginal.

8.2. Applications

We envision that the prediction models can be used in several different ways. They can be used to select the next loading action of an autonomous wheel loader or to plan the coordinated movement of both wheel loaders and haul trucks in a way that is optimal for multiple loading and hauling vehicles with the multi-objective goal of executing individual tasks efficiently while not obstructing the work of the other machines. This assumes adding a model for the V-cycle maneuver and emptying of the bucket into the load receiver. If the problem of optimal sequential planning is computationally intractable, the prediction models can be useful in developing good planning policies, for instance, using model-based reinforcement learning. Optimal dig location for repeated loading considering the characteristics of the current pile shape has been explored previously [32,33]. In contrast, our model supports optimization with respect also to future pile states. That means that it can be used to realize a specific pile shape that, for some reason, is beneficial for work site productivity and safety. If it is desirable to embed the model on a vehicle, it is worth noting that no extraordinary hardware is needed. The memory requirement, for instance, is less than 100 MB.

8.3. Implications and Future Research

This study demonstrates the feasibility of learning models for predicting the outcome of loading actions, including the future states of the pile. When these models are sufficiently precise and robust, it becomes possible to optimize wheel loading operations with respect to time and fuel consumption over long time horizons. Such a system can potentially achieve an efficiency that exceeds that of skilled human operators. This would unlock new opportunities for reducing fuel and emissions in mining and construction. This implies several research questions that first need to be addressed. Firstly, the stability of the model during long rollouts needs to be further improved. Secondly, models for more complex materials must be trained, e.g., including cohesive and fragmented rock. Finally, methods for solving the long-horizon optimization problem efficiently need to be researched.

9. Conclusions

This paper shows the feasibility of learning wheel loader world models that predict the outcome of single loading cycles given the local shape of the pile and the choice of control parameters for automatic bucket-filling. The proposed models can be used for automatic planning and control to maximize the net performance of a sequence of loading cycles predicted through repeated model inference. Topics left to explore in future work include handling more complex materials and how the proposed methods can be used for optimal planning. The latter would also provide better insight into model accuracy and inference speed requirements.