# Machine Learning Based Path Planning for Improved Rover Navigation (Pre-
Print Version)

Neil Abcouwer  Shreyansh Daftry  Tyler del Sesto
Olivier Toupet  Masahiro Ono
Jet Propulsion Laboratory
California Institute of Technology
4800 Oak Grove Drive
Pasadena  CA 91109  USA
<EMAIL_ADDRESS>Siddarth Venkatraman
Manipal Institute of Technology
Udupi - Karkala Rd  Eshwar Nagar  Manipal
Karnataka 576104  India
<EMAIL_ADDRESS>Ravi Lanka
Rakuten Inc.
19 Primrose Road  Bengaluru
Karnataka 560025  India
<EMAIL_ADDRESS>Jialin Song  Yisong Yue
Computing and Mathematical Sciences
California Institute of Technology
1200 E California Blvd
Pasadena  CA 91125  USA
{jssong<EMAIL_ADDRESS>©2020. All rights reserved. This work has been
submitted to the IEEE for possible publication. Copyright may be transferred
without notice, after which this version may no longer be accessible.

###### Abstract

Enhanced AutoNav (ENav), the baseline surface navigation software for NASA’s
Perseverance rover, sorts a list of candidate paths for the rover to traverse,
then uses the Approximate Clearance Evaluation (ACE) algorithm to evaluate
whether the most highly ranked paths are safe. ACE is crucial for maintaining
the safety of the rover, but is computationally expensive. If the most
promising candidates in the list of paths are all found to be infeasible, ENav
must continue to search the list and run time-consuming ACE evaluations until
a feasible path is found. In this paper, we present two heuristics that, given
a terrain heightmap around the rover, produce cost estimates that more
effectively rank the candidate paths before ACE evaluation. The first
heuristic uses Sobel operators and convolution to incorporate the cost of
traversing high-gradient terrain. The second heuristic uses a machine learning
(ML) model to predict areas that will be deemed untraversable by ACE. We used
physics simulations to collect training data for the ML model and to run Monte
Carlo trials to quantify navigation performance across a variety of terrains
with various slopes and rock distributions. Compared to ENav’s baseline
performance, integrating the heuristics can lead to a significant reduction in
ACE evaluations and average computation time per planning cycle, increase path
efficiency, and maintain or improve the rate of successful traverses. This
strategy of targeting specific bottlenecks with ML while maintaining the
original ACE safety checks provides an example of how ML can be infused into
planetary science missions and other safety-critical software.

Figure 1:  A view of the ENav simulation environment. The green-yellow terrain
shows how the Gradient Convolution heuristic, developed in this work, has
assessed the cost of traversing the terrain (yellow regions are higher cost)
and steers the rover toward safer regions.

###### Contents

  1. 1 Introduction
  2. 2 Baseline Navigation
    1. 2.1 Process Disparity
    2. 2.2 Analyze Terrain
    3. 2.3 Select Path
      1. 2.3.1 Candidate Paths
      2. 2.3.2 Initial Path Ranking
      3. 2.3.3 ACE on Paths
    4. 2.4 Actuation
    5. 2.5 Performance
  3. 3 Gradient Convolution Heuristic
    1. 3.1 Integration
  4. 4 Learned Heuristic
    1. 4.1 Model Architecture and Training
    2. 4.2 Integration
  5. 5 Experiments
    1. 5.1 Experiment 1: The effects of heuristics
    2. 5.2 Experiment 2: No minimum evaluations if solution found
    3. 5.3 Experiment 3: Using a broader tree of paths
    4. 5.4 Experiment 4: Using a deeper tree of paths
  6. 6 Conclusions

##  1 Introduction

The communications latency between Mars and Earth has driven Mars rover
missions to develop significant surface autonomy. On-board planning [1]
generates schedules and models resource usage over time. Autonomous target
selection [2] is used to opportunistically gather science data on desirable
specimens. Autonomous rover navigation [3] allows rovers to traverse unknown
terrain of increasing complexity and at greater rates without human drivers in
the loop.

Future missions will demand even more capable autonomous mobility. The
proposed Mars Sample Return mission’s fetch rover will be more focused on fast
and reliable mobility than gathering its own science data [4]. The stated
goals of the Mars Exploration Program Analysis Group (MEPAG) call for a
greater spatial coverage of the Martian surface [5]; this would necessitate a
greater number of rovers and landers, which could not be operated without
either a vastly increased amount of ground operators or more significant on-
board autonomy. Missions to ocean worlds in the outer solar system are of
great interest, but the increased communication delays may require rovers to
operate autonomously for several earth days without human input[6].

The Mars 2020 mission [7] and its Perseverance Rover will use the Enhanced
Navigation (ENav) library [8] to plan paths on the Martian surface. ENav takes
as input stereo imagery, maintains a 2.5D heightmap describing the terrain,
and chooses the best maneuver to safely move the rover toward the global goal.
ENav uses the Approximate Clearance Evaluation (ACE) algorithm [9] to evaluate
a sorted list of paths for safe traversal. Running the ACE algorithm on dozens
of rover poses along hundreds of candidate rover paths represents a
significant computational burden, especially if the list is sorted poorly and
many paths fail the ACE check, which is more likely in complex and challenging
terrain.

In this paper, we describe ENav and its baseline performance (Section 2). We
introduce two heuristics, one hand-designed and one machine-learned, that can
be used to more effectively sort the list of candidate paths. Because the
heightmap is represented as a two-dimensional array, a multitude of efficient
techniques from computer vision can be used for analysis. We designed the
Gradient Convolution heuristic to use convolution and Sobel operators [10] to
estimate the cost of terrain traversal (Section 3). We also used a data set of
heightmaps and corresponding ACE evaluations to train a machine learning (ML)
classifier to infer ACE values (Section 4). These two heuristics were
integrated into ENav (Figure 1) and used to more effectively sort the rover
paths before the ACE evaluation step. By incorporating the ML classifier into
the ranking process, but still checking paths for safety with ACE, we show how
ML can be integrated into ENav without sacrificing safety requirements. We
present our results for various experiments and describe how each heuristic
affected the performance of ENav in Monte Carlo simulations across multiple
terrains (Section 5). We show that integrating the heuristics improved path
efficiency, greatly reduced ACE evaluations, computation time, and the
likelihood of ”overthinking” each planning cycle, and maintained or improved
success rates compared to the baseline performance. Finally we discuss the
implications and limitations of our results (Section 6).

##  2 Baseline Navigation

The ENav library is the core of the Mars 2020 autonomous rover navigation and
constitutes the starting point and baseline performance level for this work.

At a high level, the ENav planning cycle involves following steps:

###  2.1 Process Disparity

As the rover drives, stereo cameras acquire images. Stereo correlation is used
to generate a disparity image, which is then converted into a 3D point cloud.
The point cloud is passed to ENav, which updates an internal 2.5D heightmap
centered on the rover’s position. Each point is added to the appropriate X-Y
cell, and each cell with new points is assigned a height equal to the average
height of the enclosed 3D points.

###  2.2 Analyze Terrain

Once the heightmap is updated, the terrain is analyzed and the results are
saved in a costmap. Compared to the heightmap, the costmap is lower
resolution, but covers a larger area. If a cell in the costmap corresponds to
an area that was just imaged, the tilt angle and roughness (the average
squared distance between a fitted plane and the points in the plane) are
updated. The cost of the cell may be a weighted sum of the tilt, the
roughness, and a minimum time needed to traverse a cell. The cell will be
given an infinite cost if the tilt is extreme, plane-fitting fails, or if
designated keep-in or keep-out zones (an area set by rover operators) are
violated. Costmap cells representing areas that have not yet been imaged are
given a finite cost value that is greater than that of known flat terrain, but
more than that of known obstacles.

###  2.3 Select Path

With an updated costmap, ENav can select the next path for the rover to take
to make progress towards the goal:

####  2.3.1 Candidate Paths

ENav considers a parameterized tree of candidate paths the rover can drive.
The hardware of the Perseverance rover and its progenitors does not allow
steering while driving. Thus the rover can move in fixed-curvature arcs (i.e.
constant turning radius). Special cases of arcs include driving straight (zero
curvature) and turning in place (infinite curvature). The first branch in the
tree of paths generally consists of various turns (or no turning at all) and
subsequent branches consist of fixed-length arcs of various curvatures. The
default tree of paths for our experiments is composed of 14 candidate turns-
in-place (including the option of not turning), followed by 11 3-meter arcs of
various curvatures, followed by another set of 11 3-meter arcs. This means the
rover is planning about 6 meters ahead of its current position, and
considering 1694 potential paths.

####  2.3.2 Initial Path Ranking

Each path in the tree is analyzed and assigned a total cost, which is a
weighted sum of various cost factors. Firstly, each maneuver takes a certain
amount of time to actuate; paths with less steering changes or no initial turn
have a smaller cost. Secondly, the costmap (which was updated during the
Analyze Terrain step) is sampled at regular intervals along the maneuvers in
each path. Finally, Dijkstra’s algorithm [11] is used to estimate the minimum
cost to travel from the endpoint of the considered path to the goal. After all
paths have been assessed, the list of paths is sorted by cost.

####  2.3.3 ACE on Paths

Starting with the lowest-cost path, the ACE algorithm is applied at regular
intervals along the path. Given the rover pose and heightmap, ACE finds the
min and max heights within the potential footprint of each wheel, computes
bounds on the rover attitude, suspension angles, and clearance between the
rover belly pan and the terrain. If the bounds violate parameterized limits,
ACE can return an infinite cost, rendering that path infeasible. Otherwise, a
finite cost is returned based on proximity to the limits, which accumulates
over the length of the path. Once a path with finite cost is found, ENav
continues to analyze other candidate paths until reaching a threshold, and
returns the lowest-cost ACE-validated path found.

###  2.4 Actuation

Given the lowest-cost path, only the first maneuver (either the first meter of
an arc or the first 30 degrees of a turn in place) of that path is actuated,
stereo imagery is acquired and integrated, and the planning cycle repeats.

###  2.5 Performance

This planning cycle starts when a disparity image is available for ENav, and
ends when the solution is found. If the new maneuver is found before the rover
finishes driving its current maneuver, the rover can immediately transition to
the new maneuver. Otherwise, the rover must stop driving until a solution is
found. Failing to find a solution before the rover stops is called
”overthinking”. Overthinking is undesirable, as it reduces the average
traverse rate of the rover, increases wear and tear on the wheels and brakes,
and increases mission risk overall.

Overthinking is most likely to occur due to a poor ranking of candidate paths
and excessive ACE evaluations taking significant computation time. In our
experiments, we track ENav’s overthink rate, defined as the percentage of
planning cycles where the number of ACE evaluations exceeds 275. Each instance
of ACE is estimated to take between 10 and 20 ms on the RAD750 computer used
in the Perseverance rover, so this represents around 3 to 4 seconds of
computation. In the theoretical scenario where the entire baseline tree of
1694 candidate paths were to be evaluated by ACE at 25 cm intervals, ACE would
be called over 22 thousand times, taking over 3 minutes of computation time.

In benign terrain (terrain with little slope or few rocks, fully defined in
section 5), the computational budget is rarely exceeded, but in complex
terrain where unsafe obstacles are likely to occur between the rover and the
goal, simulations predict an overthink rate of 20%. The goal of this work is
to reduce the planning computation time and increase the success rate and path
efficiency in complex terrain. To achieve this, we added heuristics rank the
potential rover paths more effectively with regards to the likelihood of
passing ACE, reducing computation time and the overthink rate.

##  3 Gradient Convolution Heuristic

(a) Height Map

(b) Gradient Map

(c) Rover Kernel

(d) Cost Map

Figure 2: Given a height map (a), Sobel operators can be used to create a
gradient or gradient squared map (b). The gradient map can be convolved with a
kernel representing the orientation-agnostic footprint of the rover (c) to
form the gradient cost map (d), which estimates the cost to traverse a
location within the map.

In order to more effectively sort ENav’s list of candidate rover paths such
that the most highly ranked paths are more likely to pass ACE, we designed the
Gradient Convolution heuristic to estimate costs that would be correlated with
ACE evaluations, but with a lower computation time.

The ACE algorithm considers how uneven the terrain is under the rover wheels.
The Gradient Convolution heuristic was designed to assess terrain roughness at
the points where the rover wheels might be. Given a heightmap, the heuristic
is computed as follows:

  1. 1.

Convolve the heightmap (Figure 2(a)) with normalized 3x3 Sobel operators to
find the local x and y gradient:

${\bf G_{x}}=\frac{1}{2r}\frac{1}{4}\begin{bmatrix}+1&0&-1\\\ +2&0&-2\\\
+1&0&-1\end{bmatrix}*{\bf A}$ (1)

${\bf G_{y}}=\frac{1}{2r}\frac{1}{4}\begin{bmatrix}+1&+2&+1\\\ 0&0&0\\\
-1&-2&-1\end{bmatrix}*{\bf A}$ (2)

where ${\bf A}$ is the heightmap, $r$ is the width of one square cell in the
heightmap, and $*$ is the convolution operator.

  2. 2.

Find the squared gradient magnitude map (Figure 2(b)) as:

${\bf G_{sq}}={\bf G_{x}}\circ{\bf G_{x}}+{\bf G_{y}}\circ{\bf G_{y}}$ (3)

where $\circ$ is element-wise multiplication. Using the squared gradient
magnitude is less computationally expensive (as no square root is required)
and creates a quadratic cost, penalizing extreme values.

  3. 3.

Convolve the squared gradient map with a kernel representing the heading-
agnostic footprint of the rover (Figure 2(c)). The footprint is an annulus,
with an outer radius corresponding to the furthest extent of any wheel and an
inner radius corresponding the closest any wheel gets to the pivot point. If
elements of the kernel are within the annulus, they have a value of 1,
otherwise 0. Normalize by dividing by the sum of the non-zero elements in the
footprint kernel, and multiply by a parameterized cost factor.

${\bf G_{c}}=\frac{k}{\sum{\bf R}_{i}}{\bf R}*{\bf G_{sq}}$ (4)

where ${\bf R}$ is the rover kernel, and $k$ is the cost factor. This provides
the gradient convolution cost map (Figure 2(d)).

Convolution and the Sobel operators are very commonly used in computer vision
contexts. As computer vision becomes increasingly important for rover
navigation, and dedicated vision-processing computers are added to Mars
missions [3], the likelihood of these methods being computed quickly and
efficiently is high.

###  3.1 Integration

During the Analyze Terrain step of the ENav planning cycle (discussed in
Section 2.2), the gradient convolution cost is calculated for each cell in the
costmap and included in the overall cost of the cell. This causes the
heuristic to contribute not only to the cost of terrain traversal along each
path, but also to the Dijkstra-based estimate of cost from the end of each
path to the overall goal.

##  4 Learned Heuristic

Figure 3: An example of a Learned Heuristic. Sets of terrain heightmaps (left)
and maps generated by the ACE algorithm (center) were used to train a neural
network to generate an inferred ACE probability map (right).

In order to find a more accurate heuristic for ACE cost, we trained a model to
predict ACE values based on heightmap data (Figure 3). In contrast to the
Gradient Convolution heuristic, which was hand-coded by domain experts, the
learned heuristic is automatically encoded using a data-driven framework. More
specifically, we developed a deep convolutional neural network (DCNN) based
model that can directly predict the outcome of the ACE algorithm for a given
terrain heightmap. Using this prediction, ENav can more optimally sort its
initial list of potential paths and hence reduce the average number of ACE
evaluations required until finding a safe path.

###  4.1 Model Architecture and Training

We formulate this problem as a supervised-learning based classification. Our
DCNN model is based on a modified encoder-decoder style U-Net architecture
[12], and implemented using the Tensorflow framework [13]. The encoder
consists of a series of convolutional layers that down-samples the input to a
low-dimensional feature map, and a decoder that consists of up-sampling layers
with convolutions that then take this feature map and increase their
resolution to that of the original input. U-Net also has a series of residual
connections from the encoder to the decoder feature maps that helps restore
the high-resolution details lost during down-sampling and also prevents
vanishing gradients during training.

The input to our model is a heightmap and the output is an ACE map, such that
the value for each pixel in the ACE map corresponds to the expected ACE cost
for the corresponding terrain parameters. However, ACE cost depends not only
on the terrain but also on the rover heading. We encode the rover heading as
part of the learning problem itself by extending the output to have a multi-
channel representation such that each channel represents a cardinal heading
angle for the rover. In our experiments, we have found a discretization of 8
heading angles (at 45 degree intervals) to be sufficient. Sigmoid activation
is applied to each channel to give a value in the range [0, 1] corresponding
to the probability of a cell being infinite ACE cost or not.

Training data was gathered by running a Monte Carlo simulation of the baseline
ENav algorithm on 1500 terrains, randomly sampling 8 heightmaps from each
trial. For each cell in each sampled heightmap, the ACE algorithm was run with
the eight fixed rover heading values, resulting in an “ACEmap” where each cell
has eight heading-specific values. Of the 12000 total heightmap, ACE map
pairs, 9500 were used as a training set, and 2500 were used as the validation
set. The learned heuristic model achieved 97.8% training accuracy and 95.3%
validation accuracy.

Note the result of this prediction is a probability of ACE returning a safety
violation, which is different from the output of the ACE algorithm itself. ACE
can return finite or infinite costs, where finite costs represent how close
the rover is to safety violations, and infinite costs represent safety
violations. Predicting the actual ACE costs is a dual problem of
classification and regression, and is more difficult that segmentation alone.
Efforts to predict these values have not yet yielded results.

###  4.2 Integration

To integrate the model inference into the ENav algorithm, the heightmaps are
passed to a TensorFlow process, which returns an ACE prediction map to ENav.
When ENav does the Initial Path Ranking step, the ACE map is sampled at
regular intervals along the paths, proportionately weighting the predictions
for the two fixed headings closest to the current heading. The ACE predction
is multiplied by a cost factor, and adds to the accumulated cost of each path.
The average ACE map value also contributes to the costmap and is used for the
Dijkstra-based estimate of cost from the end of each oath to the overall goal.

Importantly, the sorted list of paths is still evaluated by the “true” ACE
algorithm. This approach allows us to leverage the benefits of ML while
maintaining the same safety guarantees.

##  5 Experiments

Described at length in [8], a Monte Carlo simulation environment was built for
testing ENav and Mars 2020 navigation (Figure 4). The Robotics Operating
System (ROS) [14] was used for inter-process communication. One ROS node wraps
ENav, and another wraps the HyperDrive Simulator (HDSim), which simulates
rover motion, terrain settling and slipping, and disparity images for JPL
rover missions. Simulated terrains, representing various slopes and rock
densities, are loaded into HDSim. Rock densities are classified by the
cumulative fractional area covered by rocks (CFA) [15].

Figure 4: ROS software architecture for testing ENav against the HDSim physics
simulator

The rover starts at one side of the map and ENav is given a global goal 80m
away at the other end of the map. The trial is run until the rover either
successfully reaches the goal, or until a failure condition is found, such as
when no feasible path can be found, when the safety limits of the rover are
violated, or when the duration of the trial exceeds a time limit.

We leverage this existing simulation setup to run our experiments. In the case
of simulations employing the machine learning model, we use an additional ROS
node running TensorFlow, which receives heightmaps and publishes ACE
estimates. Because HDSim currently only runs on 32-bit systems, and TensorFlow
only runs on 64-bit systems, we use a ROS multi-master system to communicate
between two computers.

The discussion of the experiments notes potential gains in computation time.
These gains are theoretical and predicated on the notion that each call of the
ACE algorithm takes 10 to 20 ms on the RAD750 flight processor, these calls
make up a very large portion of ENav’s computation time, and reduction of
these calls translates to a gain in computation time. The use of non-flight-
like computers, the need for multi-master ROS, and lack of optimization in the
added algorithms for this research prevent meaningful comparisons of wall
clock time.

Figure 5: A summary of key rover path planning performance metrics across
various experiments. Incorporating heuristics produced more efficient paths,
reduced the number of costly ACE evaluations, and maintained or slightly
increased the rate of successfully reaching the goal.

The results of the following experiments are summarized in Figure 5. Each
simulation tracks results on two subsets of terrains: Benign terrains have a
CFA value of 7% or less, and a slope of 15∘ or less. Complex terrains have
greater slope or CFA. Each performance metric is calculated as a weighted
average across each subset of terrains (Benign or Complex), with terrains of
greater complexity being less likely to occur on Mars, and accordingly given
less weight. We are more concerned with tracking performance on complex
terrain. Within each Monte Carlo simulation, 780 trials are run on complex
terrain, so $n=780$ will be used when calculating 95% confidence margins of
error (MOE).

The following performance metrics are tracked:

  * •

Success Rate is the percentage of trials for each terrain that result in the
rover reaching the global goal without timing out, reaching a point with no
feasible paths, or violating safety constraints. Higher values are better.

  * •

Average Path Inefficiency is defined as the average length of the path taken
by the rover divided by the Euclidean distance from the start to the goal,
minus 1, expressed as a percentage. For example, if the rover goal was 100 m
away, and the rover needed to travel a circuitous route with a length of 125 m
to avoid obstacles and reach the goal, the path inefficiency was 25%. Lower
values are better.

  * •

Average ACE Evaluations is the average number of ACE evaluations conducted per
planning cycle. Lower values are better. Each evaluation takes an estimated 10
to 20 ms on a RAD750 processor like that on Perseverance, so a conversion to
average cycle time can easily be made.

  * •

Overthink Rate is the average percent of ENav planning cycles that required
more ACE evaluations than a threshold value (275 by default), which indicates
that the highest-ranked candidate paths were all deemed unsafe by ACE, and
therefore the initial ranking of paths was unsuitable. When the number of ACE
evaluations exceeds the threshold, it indicates that ENav is ”overthinking”
and the rover may need to stop driving until a solution is found. Lower values
are better.

###  5.1 Experiment 1: The effects of heuristics

The first experiment tested whether adding the designed and learned heuristics
to the baseline ENav software yields gains in the tracked performance metrics.
Monte Carlo simulations were run with the same terrains and parameters as the
baseline simulation. Note that the terrains used for these experiments are a
separate set from the set of terrains used to train the learned model, to
prevent unrepresentative performance due to over-fitting.

The simulation with the gradient convolution heuristic showed mixed results
(Figure 5). Benign terrain performance was fairly similar to the baseline.
Success rate dropped slightly, but the overthink rate decreased. Complex
performance showed more obvious changes. Success rate slightly worsened, from
69.9% to 67.1%, but other metrics were promising. There were significant
decreases in path inefficiency, from 25.4% to 19.9% (MOE: 2.8%). In other
words, given a goal 100 m away, the new method would be expected to drive a
120 m path, compared to 125 m for the baseline, saving time and energy. The
ACE metrics also improved, with the overthink rate reducing significantly,
from 20.0% to 14.2% (MOE: 2.5%). In summary, the designed heuristic did show
promise for improving ENav performance in complex terrains.

More impressively, the simulation with the learned heuristic showed
improvements across almost every metric. Complex success rate improved
slightly, from 69.9% to 72.5% (within the MOE). Path inefficiency
significantly improved, especially for complex terrains, going from 25.4% to
20.4%. The number of ACE evaluations also reduced, especially in terms of the
overthink rate, which plummeted to 7.1% compared to 20.0% for the baseline.
This would result in far fewer cases of the rover needing to stop driving
before the next path can be found.

We asserted that ML could be added without sacrificing the safety guarantee of
the ACE algorithm, and this assertion holds. No trial failures were caused by
violated safety constraints; all failures are due to either timeouts or
failures to find paths to the goal. This was true for all experiments.

This experiment answers the most fundamental question of this work: can
heuristics, designed or learned, improve the performance of the baseline rover
navigation algorithm? The answer is yes. Heuristics can more effectively rank
the set of candidate paths, reduce the average computation time needed to find
a safe path, choose maneuvers that increase path efficiency, and increase the
likelihood of successfully reaching the goal, while maintaining rover safety.

###  5.2 Experiment 2: No minimum evaluations if solution found

As discussed in Section 2, by default when a path that passes ACE is found,
ENav will continue to evaluate paths up to a parameterized limit. This helps
if the first safe path has a poor (but finite) ACE cost; a path with a lower
cost relative to ACE evaluations can be found.

But if these heuristics rank candidate paths more effectively, then this
continued evaluation may no longer be necessary and computational resources
might be better used elsewhere. This second experiment quantifies how much
computation time can be saved, and a what cost to other metrics. We repeated
the simulation with each heuristic, this time setting up ENav parameters such
that the first ACE-feasible path from the sorted list is chosen for the next
maneuver.

The simulation with the designed heuristic showed similar results to
Experiment 1, but with a vast reduction in the number of ACE evaluations. The
average number of ACE evaluations in complex terrain was only 156, compared to
377 for the baseline. On the Perseverance RAD750 processor, this would
translate to an average ACE computation time of around 1.9 seconds, rather
than 4.5 seconds for the baseline. However the complex success rate reduced to
63.1%.

The simulation with the ML heuristic showed promising results. Success rates
remained similar to the baseline, but the complex success rate of 69.3% does
represent a small decrease from learned heuristic success rate (72.5%) in
Experiment 1. Path efficiencies improved even more than Experiment 1, but the
biggest improvement was in the average number of ACE evaluations. In benign
terrain, the average number of ACE evaluations went from 275 to 39. In Complex
terrain, ACE Evaluations went from 377 to 90. For a RAD750 processor, this
would represent an improvement from about 4.5 seconds to 1.1 seconds, saving
3.4 seconds of computation time per planning cycle.

This experiment showed that the inclusion of learned heuristics needed less
than one fourth the number of ACE evaluations, and therefore significantly
reduced time, to achieve similar success rates as the baseline software.

###  5.3 Experiment 3: Using a broader tree of paths

Experiment 2 showed that using the ML heuristic could significantly reduce
computation time in the ENav planning cycle. Reduction of computation time is
only useful if that freed computational budget is used for methods that might
increase the success rate. One potential approach is to increase the number of
candidate paths, which can increase the probability of finding a safe path, or
provide a wider range from which to choose a more efficient path through the
terrain.

As discussed in Section 2, the baseline tree of paths has 14 potential tuns-
in-place (including no turn) at depth 1, followed by two levels of 11
potential arcs, for a total of 1694 paths. For Experiment 3, the possible
choices at each level of the tree was increased by 4, giving 4050 possible
paths, more densely covering the terrain. This means that analyzing all paths
could take around 2-3 times as much time, if not effectively sorted by the
heuristics.

Experiment 3 was conducted with the machine learning model in the loop, no
minimum number of ACE evaluations, and the broader tree. The experiment showed
the best success rates and path efficiencies so far, improving complex success
rate from 69.9% to 73.9% (MOE: 3.1%) and inefficiency from 25.4% to 17.6%.
While timing results showed more average ACE evaluations and a higher chance
of overthinking than Experiment 2, these metrics were still significantly
superior to the baseline. However, there would be an increased amount of time
required to perform the initial ranking of the larger set of paths, which is
not captured by our current experimental setup.

Because prior experiments showed that the learned heuristic was superior to
the designed heuristic, further experiments with the designed heuristic were
not conducted.

###  5.4 Experiment 4: Using a deeper tree of paths

Experiment 3 showed how assessment of a broader tree of paths becomes viable
and shows high performance with the addition of the learned heuristic. The
next experiment tests whether a deeper tree can also yield promising results.
We add an additional set of 11 arcs to the previous leaf nodes in the tree,
extending the planning horizon of the rover to nine meters ahead.

Merely adding another depth of 11 arcs to the baseline tree ($14\textrm{
turns}\times 11\textrm{ arcs}\times 11\textrm{ arcs}=1694\textrm{ paths}$),
would cause the time required to analyze the tree to multiply by 11. Instead,
we evaluate the tree up to the depth of two arcs, prune away all but the
lowest cost paths within each set of 11 depth-3 arcs, making a tree of
($14\textrm{ turns}\times 11\textrm{ arcs}\times 1\textrm{ arc}=$) 154 paths,
then extending and evaluating 11 arcs from the remaining paths. This combats
the branching factor, keeping the number of possible paths the same and
roughly doubling the amount of computation required to analyze the paths, in
line with the twice-as-large tree from Experiment 3.

The results from this approach were poor. The success rate in complex terrain
tumbled to 59.9%, and no other success metrics improved relative to
Experiments 2 or 3. Additional complexity may have caused trial timeouts to be
more likely. These losses were not offset by gains in planning; stereo
estimation of terrain nine meters away may be too inaccurate or extending the
paths to nine meters may not be useful when slip is accumulated over that
distance. Better pruning strategies may yield better results in future
experiments.

##  6 Conclusions

Our experiments show that heuristics for terrain safety can improve the
performance of the baseline ENav software.

All the heuristic experiments showed significantly increased path efficiency
in complex terrain. This may be because the ACE algorithm is likely to find
narrow safe paths that disappear once slip occurs, causing backtracking or
tight turns. In contrast, the heuristics form a smoother estimate and steer
the rover away from such problematic areas.

The experiments with no minimum number of evaluations per planning cycle
showed that the computation time associated with ACE evaluations can be
greatly reduced with only small decreases in success rate. The success rate
can be increased by using the freed computation time to evaluate a broader set
of paths, yielding results superior to the baseline in every performance
metric. Other methods might leverage the computation for even better results.

We found that the ML model outperformed the hand-designed heuristic in most of
the performance metrics in most of the experiments, but the designed heuristic
showed potential improvements compared to the baseline ENav performance. The
designed heuristic could be considered for further testing in representative
hardware and for architectures where ML inference can not be deployed.

An important caveat to our claims of improved performance is that gains in
computation time are measured purely in terms of reduced ACE evaluations. If
the time savings are completely consumed by increased time used to calculate
the heuristics, these claims are no longer valid. We hope for opportunities to
more accurately test the performance in a flight-like computing environment.

More broadly, these results offer an example of how ML can be used to improve
the performance of a safety-critical flight system without sacrificing safety
guarantees. While the ACE-estimating model is a black box, the result is used
only for ranking, not for safety-critical decisions. Any path returned from
ENav has passed a safety evaluation by the true ACE algorithm. This method
could be applied to various critical robotic or spaceflight tasks, such as
scheduling or kinematics, where the range of possible solutions is vast and
calls for an inference-based approach, but the result can be evaluated by a
known traditional algorithm that can be verified and validated.

###### Acknowledgements.

The research described in this paper was performed at the Jet Propulsion
Laboratory, California Institute of Technology, under a contract with the
National Aeronautics and Space Administration (80NM0018D0004). The authors
would like to thank the JPL Research and Technology Development (R&TD) program
for supporting this research.

## References

  * [1] W. Chi, S. A. Chien, J. Agrawal, G. Rabideau, E. Benowitz, D. M. Gaines, E. Fosse, S. Kuhn, and J. Biehl, “Embedding a scheduler in execution for a planetary rover.” in _ICAPS_ , 2018, pp. 312–320.
  * [2] R. Francis, T. Estlin, G. Doran, S. Johnstone, D. Gaines, V. Verma, M. Burl, J. Frydenvang, S. Montaño, R. C. Wiens, S. Schaffer, O. Gasnault, L. DeFlores, D. Blaney, and B. Bornstein, “Aegis autonomous targeting for chemcam on mars science laboratory: Deployment and results of initial science team use,” _Science Robotics_ , vol. 2, no. 7, 2017. [Online]. Available: https://robotics.sciencemag.org/content/2/7/eaan4582
  * [3] M. McHenry, N. Abcouwer, J. Biesiadecki, J. Chang, T. D. Sesto, A. Johnson, T. Litwin, M. Maimone, J. Morrison, R. Rieber, O. Toupet, and P. Twu, “Mars 2020 autonomous rover navigation,” in _AAS_ , 2020.
  * [4] B. K. Muirhead, A. K. Nicholas, J. Umland, O. Sutherland, and S. Vijendran, “Mars sample return campaign concept status,” _Acta Astronautica_ , vol. 176, pp. 131 – 138, 2020. [Online]. Available: http://www.sciencedirect.com/science/article/pii/ S0094576520303933
  * [5] “Mars science goals, objectives, investigations, and priorities: 2020 version,” Mars Exploration Program Analysis Group (MEPAG), Tech. Rep., March 2020\. [Online]. Available: https://mepag.jpl.nasa.gov/reports.cfm
  * [6] B. Sherwood, J. Lunine, C. Sotin, T. Cwik, and F. Naderi, “Program options to explore ocean worlds,” _Acta Astronautica_ , vol. 143, pp. 285–296, 2018\.
  * [7] K. H. Williford, K. A. Farley, K. M. Stack, A. C. Allwood, D. Beaty, L. W. Beegle, R. Bhartia, A. J. Brown, M. de la Torre Juarez, S.-E. Hamran _et al._ , “The nasa mars 2020 rover mission and the search for extraterrestrial life,” in _From Habitability to Life on Mars_. Elsevier, 2018, pp. 275–308.
  * [8] O. Toupet, T. D. Sesto, M. Ono, S. Myint, J. V. Hook, and M. McHenry, “A ros-based simulator for testing the enhanced autonomous navigation of the mars 2020 rover,” in _IEEE Aerospace_ , 2020.
  * [9] K. Otsu, G. Matheron, S. Ghosh, O. Toupet, and M. Ono, “Fast approximate clearance evaluation for rovers with articulated suspension systems,” _Journal of Field Robotics_ , vol. 37, no. 5, pp. 768–785, 2020. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/ rob.21892
  * [10] I. Sobel, “An isotropic 3x3 image gradient operator,” _Presentation at Stanford A.I. Project 1968_ , 02 2014.
  * [11] E. W. Dijkstra _et al._ , “A note on two problems in connexion with graphs,” _Numerische mathematik_ , vol. 1, no. 1, pp. 269–271, 1959.
  * [12] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [13] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard _et al._ , “Tensorflow: A system for large-scale machine learning,” in _12th $\\{$USENIX$\\}$ symposium on operating systems design and implementation ($\\{$OSDI$\\}$ 16)_, 2016, pp. 265–283.
  * [14] M. Quigley, K. Conley, B. Gerkey, J. Faust, T. Foote, J. Leibs, R. Wheeler, and A. Y. Ng, “Ros: an open-source robot operating system,” in _ICRA workshop on open source software_ , vol. 3, no. 3.2. Kobe, Japan, 2009, p. 5.
  * [15] M. Golombek and D. Rapp, “Size-frequency distributions of rocks on mars and earth analog sites: Implications for future landed missions,” _Journal of Geophysical Research: Planets_ , vol. 102, no. E2, pp. 4117–4129, 1997.

Neil Abcouwerabcouwer.jpg received his M.S. in Robotics from Carnegie Mellon
University. He is a Robotics Electrical Engineer at the Jet Propulsion
Laboratory. He has worked on various robotic software projects, including
mobility, pointing, pose estimation, and compression modules for the
Perseverance Mars rover.

Shreyansh Daftrydaftry.jpg received his M.S. in Robotics from Carnegie Mellon
University. He is a Robotics Technologist at Jet Propulsion Laboratory. At
JPL, he has worked on mission formulation for Mars Sample Return, and
technology development for autonomous navigation of ground, airborne and
subterranean robots.

Siddarth Venkatramanvenkatraman_resized.jpg is working toward his B.E. in
Computer Science and Computational Mathematics at Manipal Institute of
Technology. He was an Intern in the Perception Systems group at the Jet
Propulsion Laboratory. He has worked on perception and control implementations
with Machine Learning for various robotics projects, including a semi-
autonomous car.

Tyler del SestoTylerDelSestoBioPicture.jpg received his M.S. in Mechanical
Engineering from Carnegie Mellon University. He is a Robotic Systems Engineer
at the Jet Propulsion Laboratory. He has worked robotic operations for the
Curiosity rover as well as software and systems level testing of mobility and
autonav for the Perseverance rover.

Olivier Toupetotoupet received his M.S. in Aeronautics and Astronautics from
the Massachusetts Institute of Technology. He is currently the supervisor of
the Robotic Aerial Mobility Group at the Jet Propulsion Laboratory. His
research also includes the development of novel robotic technologies for Mars
rovers, such as the path planner for the Perseverance rover mentioned in this
paper, and the traction control algorithm for the Curiosity rover.

Ravi Lankaravi.png received his M.S. in Computer Science from New York
University. He is currently a Research Scientist at Rakuten, India. He was
previously a Research Technologist in the Machine Learning and Instrument
Autonomy group at the Jet Propulsion Laboratory. His research interests lie in
developing scalable approaches to combinatorial optimization problems, by
leveraging the interplay between machine learning and optimization.

Jialin Songsong_resized.jpg received his B.S. in Computer Science and
Mathematics from the University of Toronto. He is currently a Ph.D. candidate
at California Institute of Technology. He has worked on various projects on
using machine learning to speed up solvers for optimization problems such as
planning and structural optimization.

Yisong Yueyue_headshot.jpg received his Ph.D. in Computer Science from Cornell
University. He is a Professor of Computing and Mathematical Sciences at
Caltech. His research interests span both the theory and application of
machine learning. His research has been applied to recommender systems,
robotics, planning, automation, and experiment design in science.

Masahiro Onoono.png received his Ph.D. in Aeronautics and Astronautics from
the Massachusetts Institute of Technology. He is a Research Technologist at
the Jet Propulsion Laboratory. His broad interest is centered around the
application of robotic autonomy to space exploration, with an emphasis on
machine learning applications to perception, data interpretation, and decision
making.


# Discovery of sparse hysteresis models for piezoelectric materials

Abhishek Chandra  Authors to whom correspondence should be addressed: Abhishek
Chandra<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Bram Daniels
Electromechanics and Power Electronics Group, Department of Electrical
Engineering, Eindhoven University of Technology, 5600 MB Eindhoven, The
Netherlands  Mitrofan Curti  Electromechanics and Power Electronics Group,
Department of Electrical Engineering, Eindhoven University of Technology, 5600
MB Eindhoven, The Netherlands  Koen Tiels  Control Systems Technology Group,
Department of Mechanical Engineering, Eindhoven University of Technology, 5600
MB Eindhoven, The Netherlands  Elena A. Lomonova  Electromechanics and Power
Electronics Group, Department of Electrical Engineering, Eindhoven University
of Technology, 5600 MB Eindhoven, The Netherlands  Daniel M. Tartakovsky
Department of Energy Science and Engineering, Stanford University, 367 Panama
Street, Stanford, California 94305, United States

###### Abstract

This article presents an approach for modelling hysteresis in piezoelectric
materials, that leverages recent advancements in machine learning,
particularly in sparse-regression techniques. While sparse regression has
previously been used to model various scientific and engineering phenomena,
its application to nonlinear hysteresis modelling in piezoelectric materials
has yet to be explored. The study employs the least-squares algorithm with a
sequential threshold to model the dynamic system responsible for hysteresis,
resulting in a concise model that accurately predicts hysteresis for both
simulated and experimental piezoelectric material data. Several numerical
experiments are performed, including learning butterfly-shaped hysteresis and
modelling real-world hysteresis data for a piezoelectric actuator. The
presented approach is compared to traditional regression-based and neural
network methods, demonstrating its efficiency and robustness. Source code is
available at https://github.com/chandratue/SmartHysteresis.

††preprint: AIP/123-QED

Constitutive relationships are unique to specific materials and cannot be
universally applied across different materials. Ottosen and Ristinmaa (2005);
Bertotti (1998) Despite this, these materials’ properties and characteristics,
such as hysteresis, can be similar. For example, hysteresis in piezoelectric
materials (piezoelectric actuators) used in micro and nano-positioning devices
can reduce their positioning precision. Lin and Lin (2012); Gan and Zhang
(2019) Accurate calculations of the electric field-strain relationship between
voltage and displacement require a constitutive relationship specific to the
material.

In practice, constitutive relationships exhibit nonlinear hysteresis
behaviour, and accurate hysteresis models are crucial for tasks such as high-
resolution positioning. Conventional methods for modelling hysteresis in
piezoelectric materials include experimental or phenomenological models based
on established theories. Examples of phenomenological models include the Duhem
model, Bouc-Wen model, Jiles-Atherton model and Preisach model. Bertotti and
Mayergoyz (2005) While these models have proven effective, they also have
limitations. For example, the Preisach model requires a weight function often
estimated through empirical methods or polynomial and spline fitting, adding
complexity to the modelling process. Zeinali, Krop, and Lomonova (2019);
Daniels  _et al._ (2023) Additionally, the Preisach model does not provide
physical insights into the underlying mechanisms of hysteresis, making it
difficult to interpret the results and make predictions for new
situations.Bertotti (1998)

With the availability of increased measurement data, several data-driven
approaches for modelling and predicting hysteresis have emerged. Deep Neural
Networks (DNNs), known for their exceptional performance in data-science
applications, have been utilized to model hysteresis. Kim, Kwon, and Song
(2019); Wei and Sun (2000) Although DNNs excel at prediction, they lack
interpretability due to the absence of explicit governing equations.Montáns
_et al._ (2019) This drawback can be overcome by evolutionary algorithms,
which may uncover the structure of nonlinear dynamic systems from data and are
thus a viable alternative for modelling hysteresis.Kyprianou, Worden, and
Panet (2001) Despite this, their implementation is computationally intensive
and difficult to scale to larger situations. To address this, machine learning
algorithms that promote sparsity offer a computationally efficient option,
blending the capabilities of symbolic regression with computational
feasibility. Consequently, such techniques may result in parsimonious and
interpretable models. Montáns  _et al._ (2019)

We build on the work by Chandra et al. Chandra  _et al._ (2022) and present a
white-box technique to model hysteresis in piezoelectric materials using the
Sparse Identification of Nonlinear Dynamics (SINDy) framework. This framework,
introduced in Brunton, Proctor, and Kutz (2016) has demonstrated its potential
in learning nonlinear dynamical systems. Lai and Nagarajaiah (2019) The
proposed technique does not require assumptions such as rate dependency,
monotonicity, and congruency during the modelling phase. However, the
resulting model exhibits these properties for a piezoelectric material due to
the accurate input-output mapping. As a result, the white-box model makes
efficient system analysis possible for piezoelectric materials.

The major contributions of the manuscript are as follows,

  * •

A sparse regression-based methodology for modelling hysteresis in
piezoelectric materials is introduced.

  * •

It is shown that the presented methodology could be used to model constitutive
hysteresis relationships in a white-box dynamical form.

  * •

Several numerical experiments are presented for piezoelectric datasets and
compared with traditional regression and neural network algorithms.

The structure of the letter is as follows: First, we present our methodology
for modelling hysteresis as a sparse dynamical equation, leveraging the
inherent sparsity and dynamics of hysteresis. To validate the effectiveness of
our methodology, we conduct two numerical experiments employing the Duhem and
Bouc-Wen models. These experiments serve as a proof of concept, demonstrating
the capability of our approach to capture and replicate hysteresis behaviour
accurately. Moreover, we conduct a comprehensive comparative analysis,
contrasting our methodology with conventional regression-based and deep neural
network approaches. The evaluation considers training data size, noisy data,
and library selection.

Additionally, we apply our proposed methodology to model a real-world
piezoelectric dataset, assessing the resulting model’s accuracy and training
time compared to previously published models. Lastly, we address the modelling
of butterfly-shaped hysteresis data and present a straightforward extension of
the Duhem model to encompass butterfly-type hysteresis. Through these
investigations, we emphasize our approach’s superior performance and
versatility in accurately modelling hysteresis phenomena for piezoelectric
materials. The letter concludes with a summary of findings and implications.

In the following, we detail the methodology employed to model sparse
hysteresis relationships. To obtain a simple hysteresis model, the approach
involves obtaining time-series data of input and output ($u\in\mathbb{R}^{n}$
and $w\in\mathbb{R}^{n}$), with time derivatives either measured directly or
computed numerically. The data are then split into training and testing sets
and organized into matrices $X$ and $\dot{X}$ of size $\mathbb{R}^{n\times
2}$, with $1\leq i\leq n$.

$X=\begin{pmatrix}\vline&\vline\\\ w(t_{i})&u(t_{i})\\\
\vline&\vline&\end{pmatrix}$, $\dot{X}=\begin{pmatrix}\vline&\vline\\\
\dot{w}(t_{i})&\dot{u}(t_{i})\\\ \vline&\vline&\end{pmatrix}$,
where (.) refers to the time derivative of the quantity. In order to obtain a
simple model, a set of candidate functions ($\Theta\in\mathbb{R}^{n\times m}$)
is selected, with the assumption that a sparse combination of these functions
dictates the hysteresis dynamics. The columns of $\Theta$ represent the
potential space of basis functions, which include nonlinear functions.

$\Theta(X,\dot{X})=\begin{pmatrix}\vline&\vline&\vline&\vline&\vline\\\
\theta_{1}&\vline&\theta_{j}&\vline&\theta_{m}\\\
\vline&\vline&\vline&\vline&\vline\end{pmatrix}$,
$\Xi=\begin{pmatrix}\vline&\vline\\\ \xi_{1_{m}}&\xi_{2_{m}}\\\
\vline&\vline\end{pmatrix}$.

The objective is to obtain a set of sparse coefficient vectors
$\Xi\in\mathbb{R}^{m\times 2}$ that represent the coefficients for the linear
combination of basis functions in the library $\Theta$. The problem of
approximating hysteresis dynamics can be expressed succinctly as

$\dot{X}=\Theta(X,\dot{X})\Xi.$ (1)

The study solves the hysteresis modelling problem using the Sequential
Threshold Least-Squares (STLSQ) algorithm. Brunton and Kutz (2019) The STLSQ
algorithm is applied to the approximation problem, which is defined using the
time-series input and output data and a candidate library of basis functions.
The algorithm returns a sparse coefficient vector $\Xi$. Most of the
coefficients become zero, and the basis functions corresponding to the non-
zero elements of $\Xi$ represent the sparse dynamic relationship that dictates
the hysteresis behaviour. The accuracy of the obtained models is evaluated
using the metric relative percent error. The library is prepared using terms
from the Bouc-Wen hysteresis model, Duhem hysteresis model, and higher-order
polynomials and their combinations. The threshold ($\lambda$) for the STLSQ
algorithm is chosen based on the dataset to ensure that the obtained model is
accurate and sparse.

(a)  (b)

Figure 1: Learning hysteresis model from the data generated by simulating
Duhem model (top:) Input voltage and output displacement varying with respect
to time. (bottom:) Learned voltage vs displacement relationship compared
against the simulated data acting as ground truth.

(a)  (b)

Figure 2: Learning hysteresis model from the noisy data generated by
corrupting the Bouc-Wen output with $5\%$ Gaussian noise (top:) Input voltage
and noisy output displacement varying with respect to time. (bottom:) Learned
voltage vs displacement relationship compared against the noisy simulated data
acting as ground truth.

Next, we present numerical experiments to validate the methodology. Two
numerical experiments have been carried out on computer-simulated data to
validate the method with threshold $\lambda=0.1$. Source code is available at
https://github.com/chandratue/SmartHysteresis. Comprehensive instructions for
accessing and executing the code are provided therein. For the first
experiment, the data are produced utilizing the Duhem hysteresis model, which
is appropriate for explaining hysteresis in piezoelectric materials. Lin and
Lin (2012) PySINDy,de Silva  _et al._ (2020) the SINDy library implemented in
Python was used to perform time differentiation of input-output data. The
Duhem model characterizes the relationship between the input and output as
$\dot{w}(t)=\alpha|\dot{u}(t)|u(t)-\beta|\dot{u}(t)|w(t)+\gamma\dot{u}(t)$,
where $\alpha$, $\beta$, and $\gamma$ are parameters controlling the
hysteresis loop’s shape. Taking $\alpha=0.4$, $\beta=0,5$, and $\gamma=0.25$
with input as presented in Fig. 1. Once the STLSQ algorithm is applied, the
resulting hysteresis model is expressed as

$\dot{w}(t)=0.4|\dot{u}(t)|u(t)-0.5|\dot{u}(t)|w(t)+0.251\dot{u}(t)$ (2)

accounting for $5.8\times 10^{-6}\%$ error. The relative percentage error
$\mathcal{R}$ is chosen for the error metric.

$\displaystyle\mathcal{R}=\frac{||w^{*}-w||_{2}}{||w||_{2}}\times 100$ (3)

where $w^{*}$ is the simulated solution of the approximated model.

For the second experiment, the data are generated from the Bouc-Wen model
Ottosen and Ristinmaa (2005) given by,
${\dot{w}(t)=\alpha\dot{u}(t)-\beta|\dot{u}(t)||w(t)|^{n-1}w(t)-\gamma\dot{u}(t)|w(t)|^{n}}$,
with $\alpha=5$, $\beta=0.25$, $\gamma=0.5$, and $n=1$. The learned model has
$\mathcal{R}=4.6\times 10^{-6}\%$ error and the learned model is

$\dot{w}(t)=4.99\dot{u}(t)-0.25|\dot{u}(t)||w(t)|^{n-1}w(t)-0.49\dot{u}(t)|w(t)|^{n}$
(4)

The presented experiments on the simulated dataset for the Bouc-Wen and Duhem
models show that the method is applicable to a hysteresis dataset. However,
given the input-output dataset, numerous traditional regression-based and
neural network-based methods could be used to establish the underlying
relationship. We compare our method with three traditional methods: linear
regression, ridge regression and feedforward deep neural network. However, the
methods mentioned above need a relatively smaller feature library. For the
next two experiments, we simulate the Bouc-Wen model and use its features for
these methods for a fair comparison.

First, we test the performance of sparse discovery when the number of data
points is limited and generated from the Bouc-Wen model with parameters
$\alpha=1$, $\beta=0.5$, $\gamma=2$, and $n=1$. We perform experiments with
$1000$, $5000$, and $10000$ training points and observe that as the number of
training points increases, the method performs better and even achieves a
near-perfect model for the simulated smooth dataset, which other methods could
not achieve, as presented in Table 1. Furthermore, for as less as $1000$
training points, SINDy outperforms the traditional methods and learns the
relationship with the least error.

Next, we test our method for a noisy dataset. The simulated dataset generated
from the Bouc-Wen model with parameters $\alpha=0.4$, $\beta=0.5$,
$\gamma=0.25$, is corrupted with $1$ percent and $5$ percent Gaussian noise,
respectively. The results for SINDy and traditional methods are presented in
Table 2 and Fig. 2. $1$ percent Gaussian noise is common for raw experimental
data under controlled conditions. Our method learns the model with
$\mathcal{R}=0.01$ percent accuracy by learning the exact form of the
dynamical system with errors in the coefficients. Additionally, in the case of
$5$ percent Gaussian noise, which could happen under less controlled
experiments,Lai and Nagarajaiah (2019) the method still achieves less than
$0.3$ percent accuracy and outperforms the traditional methods.

(a)  (b)  (c)

Figure 3: Learning hysteresis model for real-world piezoelectric actuator
Ayala, Rakotondrabe, and dos Santos Coelho (2020) (top:) Input voltage and
output displacement varying with respect to time. (middle:) Learned voltage vs
displacement relationship compared against the experimental data acting as
ground truth. (bottom:) Absolute error in the learned hysteresis model varying
with respect to time.

However, measurement data may be corrupted in ways other than Gaussian noise.
An alternative form of data corruption is the presence of outliers within the
dataset, which are individual data points that exhibit substantial deviations
from the standard pattern of the data. Outliers can emerge due to various
factors, such as measurement errors, data corruption, or other forms of noise.
When outliers are present in the dataset, they can distort the identification
of the governing equations of the system, leading to inaccurate or unreliable
results. SINDy works by minimizing a loss function that penalizes the
complexity of the identified model, encouraging the selection of only the most
important terms that capture the underlying dynamics. However, outliers can
significantly increase the loss function and bias the selection of terms
towards capturing the outliers, leading to overfitting and poor
generalization. Therefore, it is important to preprocess the data and remove
outliers before applying SINDy. Several outlier detection techniques can be
used to identify and remove outliers from the data. Additionally, robust
variants of SINDy Tran and Ward (2017); Schaeffer and McCalla (2017) that are
designed to handle outliers can be used to improve the accuracy and
reliability of the identified models.

Table 1: Examining the performance of the proposed methodology to model hysteresis given a number of training points and comparing it with traditional regression methods and deep neural network using relative percent error as a metric. Method $\backslash$ Training size | 1000 | 5000 | 10000
---|---|---|---
SINDy |  $2.8\times 10^{-4}$% |  $1\times 10^{-5}$% |  $\sim$0%
Linear Regression |  $3.3\times 10^{-4}$% |  $2.3\times 10^{-5}$% |  $1\times 10^{-3}$%
Ridge Regression |  $9\times 10^{-2}$% |  $4.1\times 10^{-3}$% |  $9\times 10^{-3}$%
Neural Network |  $1.8\times 10^{-3}$ % |  $4.4\times 10^{-4}$% |  $3.9\times 10^{-4}$%
Table 2: Examining the performance of the proposed methodology to model hysteresis for noisy training data in the form of Gaussian noise and comparing it with traditional regression methods and deep neural network using relative percent error as a metric. Method | Error for 1% noise | Error for 5% noise
---|---|---
SINDy | 0.0101% | 0.25%
Linear Regression | 0.0104% | 0.26%
Ridge Regression | 1.46% | 1.24%
Neural Network | 1.54% | 33.28%

The presented two comparisons on limited and noisy datasets assumed we know
the feature library for linear regression, ridge regression and feedforward
neural network in advance. However, in practice, this is not feasible. One
could only assume the features to employ these methods in a real-world
situation. For the next comparison, we take the output generated from the
Duhem model but the features from the Bouc-Wen model, as the Bouc-Wen model is
used more frequently to model hysteresis. Ottosen and Ristinmaa (2005)
However, the library for SINDy remains the same. In such a scenario, SINDy
outperforms the traditional methods exceptionally by achieving similar
accuracy as the last case; however, the traditional methods approximate the
dataset with more than a five percent error. Hence, SINDy proves suitable when
no prior information about the model is known.

Next, we model a real-world piezoelectric actuator hysteresis dataset. An
open-source voltage-displacement dataset Ayala, Rakotondrabe, and dos Santos
Coelho (2020) was used to model a real-world dataset. The dataset is generated
for a cantilever beam-shaped piezoelectric actuator used for positioning. The
actuator was subjected to harmonic excitations, and $45000$ data points were
available. The current study used only $15000$ data points as in Ayala,
Rakotondrabe, and dos Santos Coelho (2020) without preprocessing with
threshold $\lambda=0.01$. However, several black-box treatments Ayala  _et
al._ (2015, 2020) have been carried out for the dataset. This work presents
the white-box dynamical system for the piezoelectric actuator given by,

${\dot{w}(t)=-0.17-2.38u+0.58\dot{u}+0.12|\dot{u}|-0.07w|w|}$ (5)

and the learned loop is presented in Fig. 3 along with the experimental data.
Fig. 3 also presents the voltage-displacement data varying with time used for
modelling. Also, the absolute error obtained in the learned model $|w-w^{*}|$
is presented against time in Fig. 3. The modelled equation represents the
relationship between the displacement $w$ of a physical system and the voltage
$u$ applied to it. Here, $\dot{w}$ represents the derivative of displacement
with respect to time $t$, and $\dot{u}$ represents the derivative of voltage
with respect to time. The physical interpretation of this equation can be
broken down as follows:

  * •

The first term, -$0.17$, represents a constant bias or offset that affects the
rate of change of displacement, regardless of the voltage applied.

  * •

The second term, -$2.38u$, represents the effect of the voltage u on the rate
of change of displacement. This term suggests that increasing the voltage will
decrease the rate of change of displacement while decreasing the voltage will
increase the rate of change of displacement.

  * •

The third term, $0.58\dot{u}$, represents the effect of the rate of change of
voltage on the rate of change of displacement. This term suggests that
increasing the rate of change of voltage will increase the rate of change of
displacement.

  * •

The fourth term, $0.12|\dot{u}|$, represents the effect of the absolute value
of the rate of change of voltage on the rate of change of displacement. This
term suggests that an increase in the absolute value of the rate of change of
voltage will increase the rate of change of displacement.

  * •

The fifth term, -$0.07w|w|$, represents the effect of the displacement on the
rate of change of displacement. This term suggests that an increase in the
absolute value of displacement will decrease the rate of change of
displacement, while a decrease in the absolute value of displacement will
increase the rate of change of displacement.

Overall, this equation describes how the displacement of the piezoelectric
actuator responds to changes in the voltage applied to it and the voltage
change rate. The specific coefficients in the equation determine the strength
and direction of these relationships.

In addition, we compare our approximation results with the methods that have
modelled this dataset in literature. An evolutionary algorithm Ayala  _et al._
(2015) and a deep neural network-based method Ayala  _et al._ (2020) are used
to model the dataset with an $R2$ score of $0.99$. Our obtained $R2$ score is
$0.99$ as well. However, it is well known that an evolutionary algorithm and
DNN-based methods could take minutes and even hours to train, whereas our
method takes seconds for training and testing. In another work, Nguyen and
Chauhan (2021) the authors modelled the dataset using four deep learning-based
methods, namely, Long Short-Term Memory (LSTM), Gated Recurrent Unit (GRU),
Temporal Convolutional Network (TCN) and Fully Convolutional Neural Network
(FCN). Table 3 presents the accuracy and computational time results.
Normalized root mean square error (NRMSE) was chosen as the error metric as it
was used in the previous study Nguyen and Chauhan (2021) defined as

$\text{NRMSE}=\frac{\sqrt{\frac{1}{n}\sum_{i=1}^{n}(w_{i}-w^{*}_{i})^{2}}}{w_{\max}-w_{\min}}$

where $w^{*}$ is the simulated solution of the approximated model. $w_{\max}$
and $w_{\min}$ are the maximum and minimum values for the experimental data
$w$. SINDy attained comparable accuracy to deep learning but could not beat
TCN. This behavior is because of TCN’s superior performance in modelling and
predicting time series data that are memory dependent. The hysteresis data are
itself memory dependent and hence is a suitable dataset for TCN. However,
SINDy learned the model four orders of magnitude faster in terms of
computational time. Computational time in modelling is important when the
interest is in modelling and simulation simultaneously, where the overall
computation time could be decreased using SINDy.

Table 3: Performance of the proposed methodology to model hysteresis for the dataset of piezoelectric actuator presented by Ayala et al. Ayala, Rakotondrabe, and dos Santos Coelho (2020) and comparison with the results obtained by deep learning-based methods presented by Nguyen et al. Nguyen and Chauhan (2021) using normalized root mean square error (NRMSE) as the metric. Method | NRMSE | Training Time | Testing Time
---|---|---|---
SINDy | 0.23 |  $1.8\times 10^{-2}$ s |  $4.5\times 10^{-4}$ s
LSTM Nguyen and Chauhan (2021) | 0.25 |  $104.83$ s |  $7.1\times 10^{-3}$ s
GRU Nguyen and Chauhan (2021) | 0.3 |  $101.31$ s |  $5.4\times 10^{-3}$ s
TCN Nguyen and Chauhan (2021) | 0.19 |  $61.57$ s |  $3.0\times 10^{-3}$ s
FCN Nguyen and Chauhan (2021) | 0.6 |  $74.92$ s |  $3.4\times 10^{-3}$ s

The presented experiments so far dealt with single-loop hysteresis. However,
hysteresis loops fall into two categories: single and multiple. Single-loop
hysteresis is common in literature for piezoelectric materials, while
multiple-loop hysteresis exhibits complex and nonlinear behaviour, with two
branches resembling butterfly wings. The butterfly hysteresis loop comprises
two single hysteresis loops with both clockwise and counter-clockwise
directions. However, constructing a suitable model to explain the butterfly
hysteresis effect is challenging due to the lack of available hysteresis
models for this type of behaviour. Unlike single-loop hysteresis, which only
exhibits one direction, the butterfly hysteresis effect has both directions,
making it challenging to create mathematical models that can accurately
describe it. Following, we aim to model butterfly-shaped hysteresis through
the proposed methodology.

The dataset exhibits a butterfly-shaped pattern, learned using a threshold of
$\lambda=0.1$. Fig. 4 displays the input voltage and output displacement
variation over time and compares the learned voltage-displacement relationship
and the data employed for modelling. The learned butterfly hysteresis model is
given by

$\begin{split}\dot{w}(t)=6.4|\dot{u}(t)|u(t)y(t)-3.4|\dot{u}(t)|w(t)y(t)\\\
+0.8\dot{u}(t)y(t)\end{split}$ (6)

where $w(t)$ is the output displacement and $u(t)$ is the input voltage.
$y(t)$ is the output obtained from the Bouc-Wen model

$\dot{y}(t)=3.2|\dot{u}(t)|u(t)-1.7|\dot{u}(t)|y(t)+0.4\dot{u}(t)$ (7)

The hysteresis is well-learned with the proposed methodology, as indicated by
$\mathcal{R}=3.1\times 10^{-4}$ for the butterfly-shaped hysteresis. One could
infer that $\dot{w}=2y\dot{y}$ by examining the model obtained for the
butterfly-shaped loop. Hence, a butterfly-shaped loop is a generalization of
the Duhem model itself and could be found through an intermediate step of the
first modelling single hysteresis loop. Also, this suggests that a general
butterfly-shaped hysteresis loop could be expressed as $w=y^{2}+c$ where $c$
is a constant. In literature, complicated analytical relationships have been
derived to model this relationship.Li  _et al._ (2021) However, this physical
insight could potentially help tackle more complicated butterfly responses in
the future. Also, as the considered Duhem model is known to model symmetrical
loops, similar generalizations of asymmetric Duhem type models could help
model asymmetric butterfly-shaped hysteresis loops.

(a)  (b)

Figure 4: Learning butterfly-shaped hysteresis model from the simulated data
(top:) Input voltage and output displacement varying with respect to time.
(bottom:) Learned voltage vs displacement relationship compared against the
simulated data acting as ground truth.

Finally, we present the conclusions and insights obtained from the study. The
following conclusions could be drawn from this study. This paper introduced an
approach for modelling hysteresis in piezoelectric materials using a sparsity-
promoting machine learning technique. The STLSQ algorithm was applied to
assist in sparsity during the learning process. The algorithm only requires
input-output data and accurately discovers the white-box hysteresis model with
minimal error. The relationship was expressed as a dynamic system, and the
resulting ordinary differential equation (ODE) could be simulated using any
standard ODE solver to predict hysteresis. Additionally, the method did not
require prior assumptions or knowledge about the underlying model governing
hysteresis during the modelling phase.

The presented algorithm was tested on a small training set with just $1000$
data points with an error of less than $1\times 10^{-3}$ percent. Furthermore,
the method was validated on raw experimental and noisy simulated data with up
to $5$ percent Gaussian noise. The algorithm proved accurate for an open-
source piezoelectric actuator dataset. It achieved the fastest computational
time compared to state-of-the-art methods for the same dataset in the
literature. The presented methodology also modelled the butterfly-shaped
hysteresis efficiently. Overall, the presented numerical experiments indicate
the sparse discovery of hysteresis models for piezoelectric materials to be
accurate, robust and computationally inexpensive.

## Author Declarations

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Abhishek Chandra: Conceptualization; Data curation; Formal analysis;
Investigation; Validation; Software; Writing – original draft; Writing –
review $\&$ editing. Bram Daniels: Writing – review $\&$ editing. Mitrofan
Curti: Supervision; Funding acquisition; Writing – review $\&$ editing. Koen
Tiels: Supervision; Funding acquisition; Writing – review $\&$ editing. Elena
A. Lomonova: Resources; Supervision; Funding acquisition; Writing – review
$\&$ editing. Daniel M. Tartakovsky: Supervision; Writing – review $\&$
editing.

### Data Availability Statement

The data that support the findings of this study are openly available at
https://github.com/chandratue/SmartHysteresis.

## References

  * Ottosen and Ristinmaa (2005) N. S. Ottosen and M. Ristinmaa,  _The mechanics of constitutive modeling_ (Elsevier, 2005).
  * Bertotti (1998) G. Bertotti,  _Hysteresis in magnetism: for physicists, materials scientists, and engineers_ (Gulf Professional Publishing, 1998).
  * Lin and Lin (2012) C.-J. Lin and P.-T. Lin, “Tracking control of a biaxial piezo-actuated positioning stage using generalized duhem model,” Computers & Mathematics with Applications 64, 766–787 (2012).
  * Gan and Zhang (2019) J. Gan and X. Zhang, “A review of nonlinear hysteresis modeling and control of piezoelectric actuators,” AIP Advances 9, 040702 (2019).
  * Bertotti and Mayergoyz (2005) G. Bertotti and I. D. Mayergoyz,  _The Science of Hysteresis: 3-volume set_ (Elsevier, 2005).
  * Zeinali, Krop, and Lomonova (2019) R. Zeinali, D. C. Krop, and E. A. Lomonova, “Comparison of preisach and congruency-based static hysteresis models applied to non-oriented steels,” IEEE Transactions on Magnetics 56, 1–4 (2019).
  * Daniels  _et al._ (2023) B. Daniels, T. Overboom, M. Curti, and E. Lomonova, “Everett map construction for modeling static hysteresis: Delaunay based interpolant versus b-spline surface,” IEEE Transactions on Magnetics , 1–1 (2023).
  * Kim, Kwon, and Song (2019) T. Kim, O.-S. Kwon, and J. Song, “Response prediction of nonlinear hysteretic systems by deep neural networks,” Neural Networks 111, 1–10 (2019).
  * Wei and Sun (2000) J.-D. Wei and C.-T. Sun, “Constructing hysteretic memory in neural networks,” IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 30, 601–609 (2000).
  * Montáns  _et al._ (2019) F. J. Montáns, F. Chinesta, R. Gómez-Bombarelli, and J. N. Kutz, “Data-driven modeling and learning in science and engineering,” Comptes Rendus Mécanique 347, 845–855 (2019).
  * Kyprianou, Worden, and Panet (2001) A. Kyprianou, K. Worden, and M. Panet, “Identification of hysteretic systems using the differential evolution algorithm,” Journal of Sound and vibration 248, 289–314 (2001).
  * Chandra  _et al._ (2022) A. Chandra, M. Curti, K. Tiels, E. A. Lomonova, and D. M. Tartakovsky, “Data-driven sparse discovery of hysteresis models for piezoelectric actuators,” in  _2022 IEEE 20th Biennial Conference on Electromagnetic Field Computation (CEFC)_ (IEEE, 2022) pp. 1–2.
  * Brunton, Proctor, and Kutz (2016) S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proceedings of the national academy of sciences 113, 3932–3937 (2016).
  * Lai and Nagarajaiah (2019) Z. Lai and S. Nagarajaiah, “Sparse structural system identification method for nonlinear dynamic systems with hysteresis/inelastic behavior,” Mechanical Systems and Signal Processing 117, 813–842 (2019).
  * Brunton and Kutz (2019) S. L. Brunton and J. N. Kutz,  _Data-driven science and engineering: Machine learning, dynamical systems, and control_ (Cambridge University Press, 2019).
  * de Silva  _et al._ (2020) B. M. de Silva, K. Champion, M. Quade, J.-C. Loiseau, J. N. Kutz, and S. L. Brunton, “PySINDy: a Python package for the sparse identification of nonlinear dynamics from data,” arXiv preprint arXiv:2004.08424 (2020).
  * Ayala, Rakotondrabe, and dos Santos Coelho (2020) H. V. H. Ayala, M. Rakotondrabe, and L. dos Santos Coelho, “Piezoelectric micromanipulator dataset for hysteresis identification,” Data in brief 29, 105175 (2020).
  * Tran and Ward (2017) G. Tran and R. Ward, “Exact recovery of chaotic systems from highly corrupted data,” Multiscale Modeling & Simulation 15, 1108–1129 (2017).
  * Schaeffer and McCalla (2017) H. Schaeffer and S. G. McCalla, “Sparse model selection via integral terms,” Physical Review E 96, 023302 (2017).
  * Ayala  _et al._ (2015) H. V. H. Ayala, D. Habineza, M. Rakotondrabe, C. E. Klein, and L. S. Coelho, “Nonlinear black-box system identification through neural networks of a hysteretic piezoelectric robotic micromanipulator,” IFAC-PapersOnLine 48, 409–414 (2015).
  * Ayala  _et al._ (2020) H. V. H. Ayala, D. Habineza, M. Rakotondrabe, and L. dos Santos Coelho, “Nonlinear black-box system identification through coevolutionary algorithms and radial basis function artificial neural networks,” Applied Soft Computing 87, 105990 (2020).
  * Nguyen and Chauhan (2021) X. A. Nguyen and S. Chauhan, “Characterization of flexible and stretchable sensors using neural networks,” Measurement Science and Technology 32, 075004 (2021).
  * Li  _et al._ (2021) Z. Li, H. Xu, X. Zhang, F. Shu, and Z. Li, “Development of a modified Bouc-Wen model for butterfly hysteresis behaviors,” in  _2021 IEEE 7th International Conference on Control Science and Systems Engineering (ICCSSE)_ (IEEE, 2021) pp. 159–163.