World Modeling for Autonomous Wheel Loaders

: 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.2ms and 97% in 4.5ms, respectively. Long-horizon predictions were found feasible over 40 sequential loading actions.


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 1 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 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.

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 earthmoving 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 [13,12].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 earthmoving 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 [17,16].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.

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 Fig. 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

Preface and nomenclature
This paper defines the loading cycle as starting with a machine heading in direction t to dig location 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 ′ .The wheel loader is equipped with an automatic bucket filling controller that is parametrized by some set of action parameters a.The controller is engaged when the machine reaches the dig location x.Each loading cycle can be assigned a performance 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 a and the initial state of the pile.The expected performance of a loading cycle, indexed by n ∈ N, is therefore given by some unknown performance predictor function Ψ.In other words, 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, a 1 , . . ., a N , at pose (x 1 , t 1 ), . . ., (x N , t N ) applied on a sequence of pile states H 1 , . . ., H N .Each loading cycle transforms the pile from its previous state to the next, H n → H ′ n ≡ H n+1 , according to some unknown pile state predictor function Φ: The net outcome of N sequential loading actions is thus given by the sum with initial pile state H ′ 0 = 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 When solving the problem of finding the sequence of dig locations and loading actions that maximize the net performance P 1:N , it is beneficial to use gradient-based optimization methods.These require the sought function Ψ to be differentiable with respect to a.

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, x = [x, y] and heading t defines a local frame with basis vectors {e x , e y , e z }, where e x is aligned with the dig direction, and e z is the vertical direction aligned with the gravitational field.We represent the local pile state with a discrete heightmap of I × J grid cells with the side length ∆l.The heightmap, in the local frame, is with integers i ∈ [0, I] and j ∈ [−J/2, J/2] and constant displacement δ 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∆l and J∆l, must be large enough to cover the interaction region, see Sec. 6.3.
The simplified problem of predicting the loading outcome, Eqs.(1)-(3), using the local heightmap is with local predictor functions ψ and ϕ for the performance and pile state.The computational process is described in Algorithm 1 and illustrated in Fig. 1.

Local pile characteristics
As an alternative to representing the local pile state in terms of a heightmap, we introduce a lowdimensional characterization, h ≡ [θ, α, κ x , κ y ] , in terms of four scalar quantities: slope angle θ, incidence angle α, longitudinal curvature κ x , and lateral curvature κ 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 n is computed from the local heightmap.The slope angle relative to the horizontal plane is then computed The incidence angle is the angle between the attack direction e x = t/ ∥t∥ and the pile normal projected on the horizontal plane Taking inspiration from [24], we calculate the local mean curvature and b ∈ R.This sign convention makes the curvatures positive for convex pile shape, which is the recommended shape for high-performance loading [23].
Algorithm 1: Long-horizon prediction using world models Initialization:

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

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 h ∈ R 52×52 into a lowerdimensional, regularized latent representation z ∈ R 64 .Next, given z along with the loading action a and a scale factor ∆h ≡ max(h) − min(h), an MLP predicts the regularized latent representation of the new pile state z ′ .Last, a decoder network constructs the resulting pile state ĥ′ ∈ R 52×52 from z ′ .Fig. 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 1,024 nodes and uses Leaky ReLU activation.Note that h is interpolated to size 64 × 64 before encoding and back to 52 × 52 after decoding to fit the off-the-shelf architecture.

MLP Encoder Decoder
Figure 4: The model architecture of the pile state predictor model.A VAE is paired with an MLP that learns state transition in the latent space.

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 h and ĥ′ 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.Fig. 5 illustrates the post-processing.

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 P n ) is removed from the global pile state at the dig location.This is done by computing the corresponding volume and subtracting this from the global pile along a strip as wide as the bucket, starting at x n and stretching along 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∆l.

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 manoeuver, emptying the bucket into the load receiver, and clearing spillage from the ground.

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 Fig. 1, Fig. 6 and Supplementary Video 1.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 bucketterrain 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 3 , 32 • , 0 kPa, and 8 • , which are intended to reflect the properties of dry gravel.The simulations were run with the physics engine AGX Dynamics [29]

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 bm and v bk , as follows where f bm is the momentaneous boom cylinder force and clip(value, min, max) limits value to the maximum and minimum values.The ramp function has four parameters (velocity gains k bm and k bk , and force thresholds δ bm and δ bk ) that parameterize the behavior of the controller.The actuator max speeds, v max bm and v max bk , and the normalizing boom force f b0 are set using specifications from the manufacturer.The control parameters are collected in the action vector Parameters δ bm and δ bk regulate what force magnitude is required to trigger the lift and tilt reactions while parameters k bm and k bk control how rapid the reaction is.Different values of the control parameter thus render different scooping motions, as illustrated by the examples in Fig. 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 bm = 0.6v max bm and tilt v bk = 0.6v max bk to reach a boom angle of −20 • 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.

Data collection and model training
This section describes how data was collected from the simulator and used for learning and evaluating the predictor models.

Loading outcome
The models were developed to predict the outcome of a loading cycle in terms of the performance P and the resulting local pile state h ′ given the initial pile h and action a.The performance was measured using the three essential scalar metrics, namely load mass (tonne), loading time (s), and work (kJ).Hence, P ∈ 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.

Data collection
We collected a total of 10, 718 samples, {h n , a n , h ′ n , P n } 10718 n=1 , by repeating the loading scenario (Sec.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 • and 30 • .These are illustrated in Fig. 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 (h n , a n , h ′ n , 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 Figure 8: Six initial seed piles are used, triangular, conical, and wedged, with slope angle θ.The dig location x and heading are randomized at each loading.positioned at x n , a distance of 5 m from the pile, and the heading was given a random disturbance of ±10 • , see Fig. 8.Each loading used a set of action parameters, a n , randomly sampled by Latin Hypercube Sampling (LHS) with the ranges 0.0 ≤ δ bm , δ bk ≤ 0.7 and 0.0 ≤ k bm , k bk ≤ 5.0.Note that the actual slope angle θ and incidence angle α varied according to the dig location because of the Perlin noise.

Local heightmap settings
We found that the performance and pile state predictor models benefited from using different sizes of the local heightmap and displacement δ in the dig direction.We used a 3.6 m sided quadratic heightmap discretized by 36 × 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 × 52 cells, and displaced the heightmap 1.0 m towards the approach direction.

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

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 done via grid search over the number of hidden layers (1, 2, or 3), the number of units (2 3 , 2 4 , . . ., 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 9,646 samples, while the test set size was fixed at 1,072 samples.

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, The standardization in Eq. ( 13) centers each height distribution around 0 across the entire data set before applying the min-max scaling in Eq. ( 14).We found this helps to preserve the pile-ground boundary in our reconstructions.The normalization in Eq. ( 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), L VAE = L MSE + 0.1L 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 (h, h ′ ) into latent pairs (z, z ′ ).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.

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.

Performance predictor model
In total, we developed 480 model instances during the hyperparameter tuning for both the high and low-dimensional performance predictors, ψ high and ψ 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 2, 000), calculating the median over the entire training.The inference speed is measured as the model execution time, taking the average of 1,000 model executions.
We selected the MLPs with two hidden layers with 2, 048 units as models of special interest.These are denoted ψ high ⋄ and ψ 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, ranges between 3.5 and 7 % with insignificant dependence on the activation function.For the low-dimensional model, the selected models of special interest, ψ low ⋄ and ψ low ⋆ , are smaller with two hidden layers and 512 units per layer.This amounts to roughly 5 × 10 5 parameters (no convolutional network is involved).It was trained on a smaller dataset with only 3,000 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 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 (Fig. 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 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 Fig. 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.Figure 10: The ten worst load mass predictions of ψ low ⋄ .The listed values are the errors for the low and high-dimensional models and the ground truth load mass.The heightmap of the initial pile state is shown with the grid cells color-coded by the decrease (gold) or increase (cyan) after loading.Some of the heightfields have been mirrored for illustrative purposes.

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 ĥ′ compared to the simulated ground truth heightmap h ′ .MAE and MRE are calculated in terms of the volume difference of the enclosing surfaces.In detail, the MRE is computed as 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.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 (Sec.4.5.1)reduces the MAE from 0.84 to 0.75 m 3 , with insignificant overhead in inference time.The improvement by the post-processing can also be seen from the distribution of the errors in Fig. 11.Fig. 12 visualizes five selected test predictions, marked by diamonds in Fig. 11(b).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 Fig. 12. Post-processing is most beneficial in cases where the initial state is irregular and when there is little change at the boundary (Fig. ??).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 (Fig. ??).

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, {x n , t n , a n } 40 n=1 .The initial pile is the one shown in Fig. 1.The selected headings and dig locations are shown in Fig. 13.The actions were picked randomly from the training data set.The simulated sequential loading can be seen in Supplementary Video 2. The performance predictor ψ high ⋄ is used together with the pile state predictor ϕ for high-dimensional model predictions.The low-dimensional performance predictor ψ low ⋄ is combined with the cellular automata model to predict the next pile state.
The evolution of the pile shape is shown in Fig. 13 and in Supplementary Video 2. The accumulated load mass, time, work, and residual pile volume, evolving with the number of loadings n, are summarized in Table 3.The error accumulation over time is shown in Fig. 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 Φ for the local pile state prediction, making no change outside the local heightmap as discussed in Sec.4.5.The accumulated error eventually leads to a predicted pile state with a slope steeper than the angle of repose (Fig. 13).
The high-and low-dimensional models' accumulated inference time was also measured over the 40 loading cycles.The high-dimensional model took 0.  0.01 s.When the cellular automata is 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 lowdimensional 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.

Diggability map
The predictor model can be used to search for the most favourable 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.Fig. 15 shows a map created by using Ψ high ⋄ to predict the performance using a fixed action 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.

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 (Fig. 14) might then become an obstacle.If that is the case, the pile state predictor can be replaced with mass removal and cellular automata, as described in Sec.4.6, making the cellular automata the computational bottleneck.In our implementation, the cellular automata was 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.

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 movement of both wheel loaders and haul trucks in a way that is optimal for coordinating multiple loading and hauling vehicles with the multi-objective goal of executing individual tasks efficiently while not obstructing the work of the other machines.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.

Conclusion
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.The parameters are changed by the number of units/layer and the number of hidden layers (H).The activation function used is either Leaky ReLU (L) or Swish (S).Note that the number of units is fixed at 512 for Ψ low and 2048 for Ψ high in the figures of #training data vs.MRE and time/epoch.The number of training/validation data is fixed at 3,000 for Ψ low and 9,646 for Ψ high in the other figures.The filled symbols ⋄ and ⋆ identify the models of Swish and Leaky ReLU as the extreme examples of higher inference speed Ψ low , and higher accuracy model Ψ high .Table 1 shows these specific results.

Figure 1 :
Figure 1: Overview of the wheel loader world model, its development process, and intended usage.First, a dataset of the outcome of parametrized wheel loading actions a on a pile with local shape h is collected using a simulator.Two models are trained.One model predicts the expected loading performance P, in terms of loaded mass, time, and work, given a and h.Another model predicts the new shape h ′ that the pile transitions into after the completed loading.The outcome of a sequence of loading actions on a pile with global state H can be predicted by repeated model inference, given the dig location x and heading t for each loading.

Figure 2 :
Figure 2: The mean normal n of the local heightmap h defines a slope angle θ relative to the horizontal plane.The attack angle α is the angle between the dig direction t and the normal projected onto the horizontal plane.in the e x and e y directions by fitting a quadratic surface h(x) ≈ b + c T x + 1 2 x T Qx with surface parameters Q = diag[−κ x , −κ y ], c ∈ R 2 , and b ∈ R.This sign convention makes the curvatures positive for convex pile shape, which is the recommended shape for high-performance loading[23].

Figure 3 :
Figure 3: Illustration of the model architectures.In the high-dimensional model (a), a convolutional neural network encodes the heightmap before being fed with the action vector to an MLP.In the low-dimensional model (b), the pile characteristics and action vector are inputs directly to the MLP.

Figure 5 :
Figure 5: The post-processing interpolates between the interior of the predicted ĥ′ and the edges of the initial heightmap h to preserve the details along the boundary.The pile in the heightmaps, color-coded by height, was dug from the left.
using a 0.01 s time-step and terrain grid cell-size of 0.1 m × 0.1 m.

Figure 6 :
Figure 6: Image sequence of a simulated loading cycle.A video version is found in Supplementary Video 1.

Figure 9 :
Figure 9: The distribution of relative errors for the performance predictor visualized by scatter plots and moving averages (curves).The gray histograms show the data distribution.

3 ]Figure 11 :
Figure 11: The distribution of the prediction test error (a) and AE correlation plot with and without post-processing (p-p) (b).The five points of interest, marked by diamonds (⋄) in (b), are visualized in Fig. 12.

Figure 12 :
Figure12: Five selected samples of initial pile state h and the corresponding resulting ground truth pile state h ′ and model predictions ĥ′ , with and without post-processing.The VAE reconstruction capability is included in the second column.Samples (a) have small errors in the test data set, (b) have large errors, and (c) are close to the median.In samples (d) and (e), the post-processing increases and decreases the errors, respectively.The surface under the pile shows the absolute difference compared to the ground truth final pile state, except that only in the second column is it compared to the initial pile state.

Figure 13 : 20 Figure 14 :
Figure 13: Pile state evolution after every five sequential loadings for (a) ground truth simulation, (b) high-dimensional model, and (c) low-dimensional model combined with cellular automata.Identical dig poses (indicated by arrows) and action parameters are used in the three cases.20 (a) Load mass [tonne].(b) Time [s]. (c) Work [kJ].

Figure 15 :
Figure 15: Diggability maps that show the predicted loaded mass, time, and work at 150 locations with normal headings, indicated by the lines around the edge.The color coding is explained by the histograms for each of the predicted quantities.

Figure 16 :
Figure 16: All bucket-tip trajectories in the collected dataset, randomly colored.

Figure 17 :
Figure 17: The distribution of the collected performance measurements and how they are correlated with the characteristic pile slope θ, incidence angle α, and curvature κ x and κ z .The line is the performance moving average of each characteristics variable.

Figure 18 :
Figure 18: The trend of generalization error (MRE), inference time, and training time (time/epoch) with the number of training/validation data (#training data) and parameters (#params).The parameters are changed by the number of units/layer and the number of hidden layers (H).The activation function used is either Leaky ReLU (L) or Swish (S).Note that the number of units is fixed at 512 for Ψ low and 2048 for Ψ high in the figures of #training data vs.MRE and time/epoch.The number of training/validation data is fixed at 3,000 for Ψ low and 9,646 for Ψ high in the other figures.The filled symbols ⋄ and ⋆ identify the models of Swish and Leaky ReLU as the extreme examples of higher inference speed Ψ low , and higher accuracy model Ψ high .Table1shows these specific results.

Figure 19 :
Figure 19: The function values (left) and the Autograd gradients (right) of Leaky ReLU and Swish models with respect to the action a along a direction in the action space.The prediction values are normalized.

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

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

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.
3 s in total.The low-dimensional model took