11institutetext: Jose Guadalupe Hernandez 22institutetext: BEACON Center for
the Study of Evolution in Action and Department of Computer Science and
Ecology, Evolutionary Biology, and Behavior Program, Michigan State
University, East Lansing, MI, USA 22email<EMAIL_ADDRESS>33institutetext:
Alexander Lalejini 44institutetext: Department of Ecology and Evolutionary
Biology, University of Michigan, Ann Arbor, MI, USA 44email:
<EMAIL_ADDRESS>55institutetext: Emily Dolson 66institutetext: BEACON
Center for the Study of Evolution in Action and Department of Computer Science
and Ecology, Evolutionary Biology, and Behavior Program, Michigan State
University, East Lansing, MI, USA 66email<EMAIL_ADDRESS>
# What can phylogenetic metrics tell us about useful diversity in evolutionary
algorithms?

Jose Guadalupe Hernandez  Alexander Lalejini  and Emily Dolson

###### Abstract

It is generally accepted that “diversity” is associated with success in
evolutionary algorithms. However, diversity is a broad concept that can be
measured and defined in a multitude of ways. To date, most evolutionary
computation research has measured diversity using the richness and/or evenness
of a particular genotypic or phenotypic property. While these metrics are
informative, we hypothesize that other diversity metrics are more strongly
predictive of success. Phylogenetic diversity metrics are a class of metrics
popularly used in biology, which take into account the evolutionary history of
a population. Here, we investigate the extent to which 1) these metrics
provide different information than those traditionally used in evolutionary
computation, and 2) these metrics better predict the long-term success of a
run of evolutionary computation. We find that, in most cases, phylogenetic
metrics behave meaningfully differently from other diversity metrics.
Moreover, our results suggest that phylogenetic diversity is indeed a better
predictor of success.

##  1 Introduction

Maintaining a sufficiently diverse population to successfully solve
challenging problems is a central challenge in all branches of evolutionary
computation. If the population’s diversity collapses, an evolutionary
algorithm can prematurely converge on a sub-optimal solution from which it is
unable to escape (Goldberg et al.,, 1987). While many diversity maintenance
techniques have been designed to combat this challenge, we currently lack a
clear understanding of what factors contribute to their success or failure in
any given situation. Broadly speaking, diversity maintenance techniques can
fail in two ways: 1) failure to maintain a diverse population at all, and 2)
failure to maintain diversity that is actually helpful to solving the problem.
Because more effort has historically been paid to the former category, here we
will focus on the latter.

Figure 1: Example populations with different kinds of diversity. This figure
shows two different phylogenies. Arrows show parent-child relationships. Each
node is a taxonomically unique phenotype (i.e., a phenotype with a unique
evolutionary origin). For simplicity, leaf nodes in these diagrams are assumed
to be the current set of taxa in the population; in reality, there could be
non-leaf nodes corresponding to extant taxa. A) A population with high
phenotypic diversity (phenotypic richness = 5) and low phylogenetic diversity
(mean pairwise distance = 2). B) A population with low phenotypic diversity
(phenotypic richness = 2) and high phylogenetic diversity (mean pairwise
distance = 6).

In evolutionary computation, diversity is usually evaluated by counting the
number of unique “types” in the population. These “types” may be genotypes,
phenotypes, ecotypes, species, or other descriptors of a group of solutions.
The most commonly used types are phenotypes, which are often described using
error vectors, behaviors, or output values. Sometimes a population’s diversity
is measured as the raw count of unique types (in biology, such metrics are
called “richness”). Alternatively, sometimes different metrics are used that
take into account how evenly the population is distributed across the types
(as in Shannon diversity/entropy). Occasionally more nuanced metrics are used
that consider the level of similarity of the types in the population (e.g.,
calculating Euclidean distance between error vectors, or including some sort
of clustering step). However, the metrics commonly used in evolutionary
computation are a small subset of the broader range of ways one could measure
diversity. Prior work suggests that some diversity metrics are more predictive
of success than others (Mouret and Doncieux,, 2009; Jackson,, 2010), so
investigating the implications of a wider variety of diversity metrics is
worthwhile.

A variety of diversity metrics from ecology and evolutionary biology have yet
to be examined in the context of evolutionary computation. One particularly
interesting class of diversity metrics, called phylodiversity metrics (Tucker
et al.,, 2017), takes into account the evolutionary history of a population to
calculate its diversity. For intuition behind how these metrics are different
from other diversity metrics, see Figure 1. They do so by measuring the
topology of the phylogenetic tree for the population (i.e., the ancestry
tree). Intuitively, types that are more evolutionarily distant from each other
(i.e., share a more distant common ancestor) are likely farther apart from
each other in the fitness landscape. As such, the ability to maintain
evolutionarily distinct types may be important for effective problem-solving
in evolutionary computing, due to the likely difficulty of re-evolving such
distinct types from each other. Notably, biological simulations suggest that
phylogenetic diversity provides different information about a population than
a specific type of phenotypic richness (functional diversity) does; moreover,
the extent of this difference varies across different scenarios (Tucker et
al.,, 2018).

Preliminary data supports the hypotheses that 1) phylogenetic diversity
captures information that other diversity metrics do not (Dolson et al.,
2018b, ; Tucker et al.,, 2018), and 2) phylogenetic diversity may be a better
predictor of success in genetic programming than other diversity metrics
(Dolson et al., 2018b, ). Here, we investigate these hypotheses in the context
of a range of problems and diversity maintenance techniques. Our findings
provide further support for these two hypotheses. In particular, using a
technique called causality analysis, we show that phylogenetic diversity is a
stronger predictor of future fitness than phenotypic diversity is.

##  2 Methods

To identify phylogenetic metrics that correlate with success in evolutionary
computation, we ran an evolutionary algorithm on a variety of problems using a
range of parent selection methods in a full-factorial design. Specifically, we
applied each of tournament, random, fitness sharing, lexicase, and Eco-EA
selection to the following problems: an exploration diagnostic, Count Odds,
Logic-9, sorting networks, and NK Landscapes.

###  2.1 Selection methods

#### Tournament

In tournament selection, $T$ individuals are randomly selected from the
population to form a “tournament”. The individual in the tournament that has
the highest fitness is selected to reproduce. Tournament selection is included
here as a control, as it is simple and introduces no pressure for diversity.

In the exploration diagnostic experiments, we used a tournament size of $T=8$
for consistency with (Hernandez et al.,, 2021). In the other experiments, we
used a tournament size of $T=2$.

#### Random

Random selection is also included in this paper as a control. In this
selection method, individuals are selected at random to reproduce.

#### Fitness sharing

Fitness sharing is one of the earliest diversity maintenance techniques
(Goldberg et al.,, 1987). It is very similar to tournament selection, except
that before the tournament occurs, the fitness of each individual is
discounted in proportion to the number of similar individuals and the degree
of their similarity. As a result, fitness sharing creates pressure for
increased population diversity.

Fitness sharing requires two parameters: $\alpha$ (which controls the shape of
the sharing function) and $\sigma_{s}h$ (which controls the niche width).
Here, we used $\alpha=1$ and $\sigma_{s}h=2$.

#### Lexicase selection

Lexicase selection is designed to work in situations where fitness can be
broken down in to multiple constituent parts (Spector,, 2012). Conventionally,
these constituent parts, or fitness criteria, are individual test cases in a
genetic programming problem, but lexicase selection has been shown to work
well in other scenarios too (Dolson et al., 2018a, ; Metevier et al.,, 2019).
To select a parent with lexicase selection, the fitness criteria are randomly
ordered. Then the criteria are stepped through in order, with all but the best
performers on each criterion being eliminated from consideration for this
selection event. This process continues either until only one individual
remains to be selected or until there are no criteria left, in which case an
individual is randomly selected from those remaining. Lexicase selection is
known to generally maintain high levels of diversity (Helmuth et al.,, 2015,
2016), and preliminary evidence suggests that lexicase selection is also
particularly effective at maintaining phylogenetic diversity (Dolson et al.,
2018b, ).

#### Eco-EA

Like lexicase selection, Eco-EA is designed to work on problems that have
various sub-components (Goings and Ofria,, 2009; Goings et al.,, 2012). The
original inspiration for Eco-EA came out of problems in which solutions to
multiple simpler tasks served as helpful building blocks for solving a more
complex overall task. These building block tasks are analogous to fitness
criteria in lexicase selection (Dolson et al., 2018a, ; Dolson et al., 2018b,
). Eco-EA associates each sub-task with a limited resource. Individuals that
perform a task gain fitness while resources associated with that task are
available; however, this decreases that task’s associated resources. When the
resource is depleted, there is no fitness benefit to performing the task. New
resources flow in and out at a constant rate over the course of each update.
As such, there is pressure for diversity.

The specific fitness gain associated with performing a task is described by
the following function:

$\text{total\\_fitness}=\text{base\\_fitness}*2^{\text{min(score}^{\alpha}\text{
}*\text{ }C_{f}\text{ }*\text{ }\text{resource, max\\_bonus})}$ (1)

where base fitness is the individual’s fitness according to the global fitness
function, score is the individual’s score on the current sub-task, $\alpha$ is
a variable that tunes the shape of the relationship between sub-task
performance and fitness gain, $C_{f}$ is the fraction of available resources
that can be consumed, resource is the amount of the resource available, and
max bonus is the maximum allowed bonus. To strengthen the negative frequency
dependence, a cost can be applied for attempting to use more resource than is
available.

For the exploration diagnostic, we used $C_{f}=1$, $\alpha=.25$, cost $=,1$,
and max bonus $=5$. Resources flow in at a rate of 250 units per generation.
Every generation, .01% of the total quantity of each resource flows out. For
the other experiments, we used $C_{f}=.01$, $\alpha=2$, cost $=3$, and max
bonus $=5$. Resources flow in at a rate of 50 units per generation and 5% of
each resource flows out each generation.

###  2.2 Problems

#### Exploration diagnostic

Recently, Hernandez et al. proposed a suite of “diagnostic” fitness landscapes
designed to test evolutionary algorithms in a controlled environment
(Hernandez et al.,, 2021). Among these is an exploration diagnostic, which
tests an algorithm’s ability to simultaneously explore many different paths
through a fitness landscape. In this diagnostic fitness landscape, genomes are
vectors of floating point numbers. An individual’s fitness is determined by
finding the highest value in its genome, which represents the start of the
values that contribute to fitness. This value and all subsequent monotonically
decreasing values are summed together to produce the total fitness. For
lexicase selection and Eco-EA, the fitness contribution of each individual
site was used as the selection criterion or sub-task respectively.

Our hypotheses are predicated on the idea that phylogenetic diversity is a
better measure of fitness landscape exploration than more conventionally used
diversity metrics. As such, the exploration diagnostic represents an ideal
context in which to test our hypotheses. Moreover, this fitness landscape
provides a smooth and controlled environment with more intuitive diversity
dynamics. In contrast, more complex fitness landscapes can yield idiosyncratic
crashes in diversity that vary based on the path(s) the population takes.
Thus, we will use the exploration diagnostic for our in-depth analysis, and
provide results in the more complex landscapes as a follow-up investigation of
the robustness of our findings.

#### Count Odds

Count Odds is part of a benchmark suite of program-synthesis problems derived
from a set of real programming problems given to introductory computer science
students (Helmuth and Spector,, 2015). These problems are specified using sets
of testcases where each testcase contains inputs and expected outputs.
Solutions to these problems consist of evolved code that produces the correct
output for a given set of inputs. In the Count Odds problem, programs are
given a list of numbers as input and must output an integer indicating how
many of those numbers are odd.

For lexicase selection and Eco-EA, each test case is used as one selection
criterion or sub-task, respectively. For fitness sharing, similarity is
calculated based on the euclidean distance of the vector of performances on
all test cases.

#### Logic-9

In Logic-9, a solution comprises code that, when executed, performs the nine
one- and two-input bitwise Boolean logic tasks (not, and, or, nand, orn, andn,
nor, xor, and equ) (Ofria and Wilke,, 2004; Dolson et al.,, 2019). To evaluate
a program, we provide two numeric inputs that the program may use during its
execution. Over the course of evaluation, a program may output numbers. We
check whether each output is the result of performing one of the nine logic
operations. If so, the program is counted as having completed that operation.
Every logic operation completed increases the program’s score by one. In
addition to the Boolean logic tasks, programs get credit for solving the echo
task, in which the program must output one of the numbers it received as
input. This task is included because it is known to be helpful in scaffolding
the evolution of the more complex logic tasks.

For lexicase selection and Eco-EA, each logic task is used as one selection
criterion or sub-task, respectively. For fitness sharing, similarity is
calculated based on the euclidean distance of the vector of which tasks were
completed.

#### NK Landscapes

In NK Landscapes, individuals are bitstrings of length $N$ (Kauffman and
Levin,, 1987). Each position in the bitstring has an associated table used to
look up the fitness contribution of that position, which is based on the value
(0 or 1) at that position and the values of $K$ neighboring positions. Thus,
each lookup table comprises $2^{K+1}$ entries, one for each possible sequence
of relevant bits. Total solution quality for a given bitstring individual is
calculated by summing up the fitness contributions from each site.

An NK Landscape is generated by randomly selecting a value between 0 and 1 to
be the fitness contribution for each entry in each site’s fitness contribution
lookup table. $N$ determines the size of the fitness landscape, and $K$
determines the landscape’s level of epistasis or ruggedness (i.e., how
interdependent the fitness contribution of each site is on other sites). For
example, if $K$ is 0, each site has two possible fitness contributions, one
for if that position is set to 1 and one for if it is set to 0. If $K$ is set
to 1, each site has four possible fitness contributions based on the value of
that site and the value of its neighbor (one each for 00, 01, 10, and 11).
Every increase in $K$ doubles the possible different fitness contributions.

For lexicase selection and Eco-EA, the fitness contribution from each site is
used as one selection criterion or sub-task, respectively. For fitness
sharing, similarity is based on the Hamming distance between the solutions.
For the experiments presented here, we used $K=3$ and $N=20$.

#### Sorting Networks

Sorting networks are computational units designed to sort fixed-length
sequences of numbers via a set of comparisons between pre-specified positions
in the sequence (Sekanina and Bidlo,, 2005). When the network compares two
positions, the numbers in them are swapped if they are out of order. In this
problem, individuals are represented as sequences of comparators between
positions. As with the Count Odds problem, fitness is assessed via a sequence
of test cases. Once all test cases are solved correctly, individuals can
receive additional fitness bonuses for having as few comparators as possible.
For lexicase selection and Eco-EA, each test case is used as a single
selection criterion or sub-task respectively.

Here, we evolved sorting networks to sort 30 values and test them on 100 test
cases. The maximum allowed number of comparators per network was 128.

###  2.3 Computational Substrates

For the genetic programming problems (Count Odds and Logic-9), we evolved
linear genetic programs where each genome is a sequence of simple
computational instructions. Most notably, the instruction set is designed to
support the evolution of modularity by supporting the encapsulation of
subroutines into “scopes”. Programs have a set of read-only memory spaces used
to provide input and a set of write-only memory spaces to be used as output.
Each instruction in the genome is executed in sequence. If execution reaches
the end of the genome before the program runs out of evaluation time,
execution will loop back around to the front of the genome. We propagated
programs asexually and applied the following mutations to offspring, each with
a probability of 0.005: 1) instruction substitutions, 2) point insertions and
deletions, and 3) substituting the argument being supplied to an instruction.
This linear genetic programming representation is described in more detail in
(Dolson et al.,, 2019).

For NK Landscapes, we evolved bitstrings. We propagated individuals asexually,
and we applied bit flip mutations at a per-bit rate of 0.01.

In the exploration diagnostic, we evolved vectors of 50 floating point values
where each value ranged between 0.0 and 25.0. We reduced the size of these
ranges from those used by Hernandez et al., (2021) to ensure that Eco-EA had
sufficient time to solve the problem within our allocated computational
budget. We used asexual reproduction, and each position had a 0.007
probability of mutating. Mutations modify a value by a number drawn from a
normal distribution with mean 0 and standard deviation 1.

Sorting networks are represented as sequences of pairwise comparators. They
can mutate via insertions of new comparators (0.0005 probability),
duplications of existing comparators (0.0005 probability), deletions of
existing comparators (0.001 probability), swapping pairs of existing
comparators (0.001 probability), and substituting different indices in
existing comparators (0.001 probability).

###  2.4 Other Parameters

For the exploration diagnostic, we used a population size of 500 and allowed
runs to evolve for 500,000 generations. This length was selected to ensure
that all selection schemes (most notably Eco-EA) had adequate time to find a
good solution. For the other fitness landscapes, we used a population size of
1000 and allowed runs to evolve to 1,000 generations.

###  2.5 Phylogenetic Diversity Metrics

A wide variety of phylogenetic diversity metrics have been developed, and the
extent to which they capture different information from each other is an area
of active research (Tucker et al.,, 2017). All of them require that you keep
track of the full phylogeny (ancestry tree) of a population. For further
discussion of building phylogenies in the context of evolutionary computation,
see (Dolson et al.,, 2020). Here, we focus on two classes of metrics: pairwise
distance metrics and evolutionary distinctiveness metrics.

Pairwise distance metrics calculate the number of edges111Weighted edges can
also be used, in which case the weights along the path should be summed. Here,
we use unweighted edges. in the shortest path between each pair of nodes
associated with extant taxa (i.e., taxonomic units corresponding to
individuals in the current population) (Webb et al.,, 2002). The resulting set
of distances can then be summarized by calculating its minimum, maximum, mean,
or variance. Each of these statistics produces a different phylogenetic
diversity metric with different properties.

Evolutionary distinctiveness metrics assign an evolutionary distinctiveness
score to each extant taxon (Isaac et al.,, 2007). This score takes into
account each branch’s age, and represents how evolutionarily distant each
taxon is from all other extant taxa. To calculate evolutionary
distinctiveness, the age of each branch is calculated and divided by the
number of extant taxa the branch ultimately leads to. A taxon’s evolutionary
distinctiveness is the sum of the values calculated for all branches between
that taxon and the tree’s root. As with pairwise distances, the set of
evolutionary distinctiveness scores can be summarized by taking its minimum,
maximum, mean, or variance.

###  2.6 Analysis techniques

#### Statistics

Correlations among variables at a fixed time point were measured using
Spearman correlations. We used Spearman correlations rather than Pearson
correlations due to the fact that many of the relationships being measured
were non-linear (but still monotonically increasing). To compare different
conditions, we used Kruskal-Wallis tests with subsequent pairwise Wilcoxon
rank-sum tests and a Bonferonni correction for multiple comparisons. To ensure
statistical rigor, we decided what set of statistical comparisons to perform
based on analysis of an initial pilot data-set. We then re-ran all experiments
with different random seeds and performed the pre-determined analysis on this
new data to generate the results presented here.

#### Transfer entropy

As alluded to previously, there is a positive feedback loop between diversity
and evolutionary success during the initial phase of evolution. To further
complicate matters, fully solving a problem can lead diversity to crash as
that solution sweeps through the population. For these reasons, looking at
correlations between fitness and diversity at any given time point gives us
only limited information about their relationship. One approach to getting
around this problem is to take all measurements at the time step when a
perfect solution is first discovered (Dolson et al.,, 2017). While this
approach helps, it does not address the question of what is driving the
feedback loop. Moreover, even before a perfect solution evolves, the evolution
of partial solutions may initiate partial selective sweeps.

In order to understand the specific role that different types of diversity
play in driving the feedback loop between diversity and success, we turn to an
analytical approach called causality analysis. As the nature of causality can
quickly drift into murky philosophical territory, here we specifically use a
notion of causality called Granger causality (Granger,, 1969; Bressler and
Seth,, 2011). We say that $X$ Granger-causes $Y$ if past values of $X$ contain
information about the current value of $Y$ above and beyond the information
that past values of $Y$ contain about the current value of $Y$. This
definition comes from the insights that 1) the past causes the future, not the
other way around, and 2) if $X$ and $Y$ are jointly caused by an external
process, that process will be captured by the information that past values of
$Y$ have about the current value of $Y$.

Granger causality is normally measured in the context of vector auto-
regressive of models. However, our data do not match the assumptions of such
models (particularly stationarity). Thus, we measure Granger-causality with an
information theoretic metric called Transfer Entropy (Schreiber,, 2000; Yao
and Li,, 2020). In information theoretic terms, transfer entropy is
$I(Y_{t};X_{t-k}|Y_{t-k})$, the conditional mutual information between $Y_{t}$
and $X_{t-k}|y_{t-k}$. Here, $Y$ is the variable being predicted, $t$ is the
time point it is being predicted at, $X$ is the variable we are using the
predict $Y$, and $k$ is the “lag”. The lag specifies the time scale on which
we are interested in measuring Granger causality. In other words, it indicates
which past value of $X$ we are attempting to use to predict the current value
of $Y$.

Often, the goal of these measurements is simply to establish the direction of
Granger-causality. In those cases, very short lags are often used. In this
case, however, it would also be valuable to know whether predictive capability
is maintained over large lags. If it were, there would be a variety of useful
practical implications. For example, we could potentially use the diversity at
a relatively early time point to predict whether a given run of evolutionary
computation will be successful. For this reason, we measure Transfer Entropy
using lags ranging from 10 generations to 100,000 generations.

Here, we measure both the Transfer Entropy between fitness and phylogenetic
diversity and the Transfer Entropy between fitness and phenotypic diversity.
Note that, because Transfer Entropy is a value calculated based on two time
series, we cannot meaningfully lump multiple replicates into the same
calculation. Thus, for each condition we will end up with a distribution of
Transfer Entropy values.

###  2.7 Code Availability

All code used in this paper is open source and freely available in the
supplemental material (Dolson,, 2021). Research code was written in C++ using
the Empirical library (Ofria et al.,, 2018). Data analysis was performed using
the R statistical computing language, version 4.0.4 (R Core Team,, 2021), the
ggplot2 (Wickham,, 2016), ggpubr (Kassambara,, 2020), and infotheo (Meyer,,
2014) libraries.

A C++ implementation of phylogeny tracking and all phylodiversity metrics used
here is available in the Empirical library (Ofria et al.,, 2018). This
implementation is designed to plug into any computational evolution code
written in C++.

##  3 Results and Discussion

###  3.1 Do phylogenetic metrics provide novel information?

If phylogenetic diversity is to tell us anything useful, a necessary first
step is that it provide information that more commonly-used metrics do not. It
is particularly important to establish this distinction in light of the fact
that phylogenetic metrics are typically more computationally expensive to
calculate.

Figure 2: Relationship between phenotypic diversity and phylogenetic diversity
at a single time point for the exploration diagnostic. Values were measured at
the final time point (generation 500,000). Spearman correlation coefficients
are shown for each selection scheme, along with a 95% confidence interval
around regression lines. Note that axis scales vary between panels.

First, we investigated the relationships between different metrics of
phenotypic and phylogenetic diversity. The two measurements of phenotypic
diversity that we analyzed were phenotypic richness and phenotypic Shannon
diversity, the two most commonly used diversity metrics (see supplemental
material (Dolson,, 2021)). Unsurprisingly, phenotypic richness and phenotypic
Shannon diversity are closely correlated across all conditions. We performed
the rest of the analyses in this paper using both richness and Shannon
diversity, but observed qualitatively the same results for both. For
simplicity, we only present the results using richness (for Shannon diversity
results, see supplemental material (Dolson,, 2021)).

There are many ways of measuring phylogenetic diversity (Tucker et al.,,
2017). Here, we focus on the pairwise distance and evolutionary
distinctiveness metrics (see 2.5 for more information). Minimum pairwise
distance was not informative, as there are nearly always at least two taxa a
distance of 1 away from each other. In the conditions observed here, the other
pairwise distance metrics tended to correlate fairly closely with each other.
Correlations among evolutionary distinctiveness metrics were weaker and less
consistent, but still present in most cases. In contrast, we did not observe
consistent relationships between pairwise distance metrics and evolutionary
distinctiveness metrics. While there were sometimes strong correlations within
a condition, the direction of these correlations varied. For this reason, we
performed subsequent analyses using both mean pairwise distance and mean
evolutionary distinctiveness. As we observed qualitatively the same results
for both, here we present only the results using mean pairwise distance (for
evolutionary distinctiveness results, see supplemental material (Dolson,,
2021)).

Next, we compared the phenotypic metrics to the phylodiversity metrics.
Intuitively, we might expect them to be highly correlated. In practice,
however, we see that the correlation between phenotypic and phylogenetic
diversity is not consistently222The correlation for tournament selection in
the exploration diagnostic is incredibly high, however 1) the observed range
of mean pairwise distance is so low that the correlation is almost certainly
an artifact, and 2) this correlation is not observed for other fitness
landscapes. significantly different from 0 (see Figure 2). In some cases, the
correlation is even negative.

(a) Phenotypic richness over time for each selection scheme on the exploration
diagnostic. Shaded areas represent 95% confidence interval around the mean for
all 50 replicates for each selection scheme. Note that the x-axis is on a log
scale.

(b) Phylogenetic diversity (mean pairwise distance) over time for each
selection scheme on the exploration diagnostic. Shaded areas represent 95%
confidence interval around the mean for all 50 replicates for each selection
scheme. Note that the x-axis is on a log scale.

As previously noted, diversity metrics can be sensitive to the exact time at
which they are measured. To confirm that this lack of instantaneous
correlation is indicative of a consistent lack of relationship, we plotted
phenotypic and phylogenetic diversity over time (see Figures 3(a) and 3(b)).
Indeed, phenotypic and phylogenetic diversity behave differently over long
temporal scales as well.

Figure 4: Transfer entropy between phenotypic and phylogenetic diversity on
the exploration diagnostic. Each boxplot shows the distribution of observed
transfer entropies for each direction of transfer, lag, and selection scheme.
Results shown here are for 50 replicate runs of each selection scheme on the
exploration diagnostic. Note that the x-axis is on a log scale.

Having established that phenotypic and phylogenetic diversity are not reliably
correlated, we next asked whether either diversity metric can predict the
other. Based on our hypothesis that phylogenetic diversity is useful because
it more directly indicates a population’s spread across the fitness landscape,
we might expect current phylogenetic diversity to predict future phenotypic
diversity. We tested this hypothesis by measuring the transfer entropy from
phylogenetic diversity to phenotypic diversity. Consistent with our
hypothesis, transfer entropy from phylogenetic diversity to phenotypic
diversity was generally higher than transfer entropy from phenotypic diversity
to phylogenetic diversity (see Figure 4). Thus, in the Granger sense of
causality, we can say that phylogenetic diversity produces phenotypic
diversity to a greater extent than phenotypic diversity produces phylogenetic
diversity.

In the exploration diagnostic fitness landscape, the only exceptions to this
observation are our two control selection schemes: tournament selection and
random selection. In tournament selection, neither form of diversity is
particularly predictive of the other. This behavior is unsurprising, as
tournament selection generally maintains minimal levels of both forms of
diversity. In random selection, both forms of diversity are somewhat
predictive of each other. In the other fitness landscapes we observe more
variation in transfer entropy (see supplemental material (Dolson,, 2021)).

Taken together, these results provide strong evidence that, in an evolutionary
computation context, phylogenetic diversity metrics provide information that
phenotypic diversity metrics do not. We base this conclusion on the lack of
consistent correlation between these metrics at a fixed point in time, the
differences in their long-term behavior, and that fact that the transfer
entropy from phylogenetic diversity to phenotypic diversity is higher than the
other way around (implying that phylogenetic diversity contains information
about future phenotypic diversity that current phenotypic diversity does not
contain about future phylogenetic diversity). We may now proceed to ask
whether that information is actually useful for understanding problem-solving
success in evolutionary computation.

###  3.2 Do phylogenetic metrics predict problem-solving success?

Having established that phylogenetic metrics provide different information
than phenotypic diversity metrics, the next question to ask is what that
information can tell us. In the exploration diagnostic landscape, we can see
some intuitive connections between both types of diversity and fitness (see
Figure 5). Excluding random selection, the final performance of a selection
scheme appears to be correlated with the final level of phylogenetic diversity
maintained by that selection scheme (but not with the level of phenotypic
diversity). In the other fitness landscapes, however, the connection between a
selection scheme’s ability to maintain diversity of either type and its
ability to succeed is less obvious (see supplemental material (Dolson,,
2021)). This increased complexity is unsurprising, as success in the
exploration diagnostic landscape is based primarily on the ability to explore;
succeeding on the other fitness landscapes is more complicated.

Figure 5: Fitness over time for each selection scheme on the exploration
diagnostic. Shaded areas represent 95% confidence interval around the mean for
all 50 replicates for each selection scheme. Note that the x-axis is on a log
scale. The y-axis shows the proportion of the maximum possible fitness
achieved.

It should also be noted that, while there is a correlation between
phylogenetic diversity and fitness in the aggregate on the exploration
diagnostic landscape, there is not a consistent correlation among individual
runs within a selection scheme (see supplemental material (Dolson,, 2021)).
Eco-EA, for instance, appears to lose phylogenetic diversity as it approaches
a good solution, leading to a negative correlation between fitness and
phylogenetic diversity. For most of the other selection schemes, there is no
significant correlation 333In the pilot data set, we observed a strong
positive correlation between phylogenetic diversity and fitness for lexicase
selection. However, this correlation disappeared when we re-ran the
experiments to generate the final data set..

To understand the precise dynamics driving the relationship between diversity
and success, we measured the transfer entropy of phylogenetic diversity to
fitness and the transfer entropy of phenotypic diversity to fitness. In the
exploration diagnostic landscape, we see that, for the three non-control
selection schemes, phylogenetic diversity is substantially more predictive of
future fitness than phenotypic diversity (see Figure 6). The predictive power
of both types of diversity weakens substantially in the other fitness
landscapes (see Figure 7). However, when there is a discernible difference,
phylogenetic diversity has higher transfer entropy than phenotypic diversity.

Figure 6: Transfer entropy from diversity to fitness for the exploration
diagnostic. Each boxplot shows the distribution of observed transfer entropies
for each type of diversity, lag, and selection scheme. Results shown here are
for 50 replicate runs of each selection scheme on the exploration diagnostic.
Note that the x-axis is on a log scale. Figure 7: Transfer entropy from
diversity to fitness for complex fitness landscapes. Each boxplot shows the
distribution of observed transfer entropies for each type of diversity, lag,
and selection scheme. Results shown here are for 50 replicate runs of each
selection scheme on each landscape. Note that the x-axis is on a log scale.

From these results, we can conclude that, for selection schemes that maintain
diversity, phylogenetic diversity “Granger-causes” success to a greater extent
than phenotypic diversity does. We interpret this as strong evidence that
phylogenetic diversity is, in general, more predictive of success than
phenotypic diversity (as measured by phenotypic richness and phenotypic
Shannon diversity).

##  4 Conclusion

We have demonstrated that, in the context of evolutionary computation,
phylogenetic diversity metrics capture information information that is
substantially different from the information captured by conventionally used
phenotypic diversity metrics (phenotypic richness and phenotypic Shannon
diversity). The extent of this difference appears to vary by problem and by
selection scheme, but it is evident in 1) the lack of consistent strong
correlation between phenotypic and phylogenetic diversity at a fixed point in
time, 2) the lack of similarity in long-term trends in phenotypic and
phylogenetic diversity, and 3) the fact that phylogenetic diversity is better
able to predict future phenotypic diversity than the other way around.

Moreover, our results also suggest that phylogenetic diversity is generally a
stronger driver of success than phenotypic diversity (in an evolutionary
computing context). This finding also varies by selection scheme and by
problem, with selection schemes that maintain diversity showing a stronger
relative effect of phylogenetic diversity than tournament selection does.
Impressively, this is true even at very long time lags; phylogenetic diversity
provides predictive information about fitness tens of thousands of generations
in the future.

Taken together, these results suggest that it may be worthwhile for
researchers studying diversity in evolutionary computation to measure
phylogenetic diversity in addition to or instead of phenotypic diversity.
Doing so will take us a step closer to identifying diversity that is helpful
to solving a given problem. Additionally, these results hint at the
possibility of using phylogenetic diversity early in a run of evolutionary
computation as a predictor of which runs will go on to be most successful.
Anecdotally, we have found phylogenetic diversity in just the first few
generations to be a useful indicator of whether we have correctly selected
parameters for fitness sharing and Eco-EA.

In this paper, we have barely scratched the surface of phylogenetic metrics
that may be relevant to evolutionary computation. Given our findings thus far,
a more thorough investigation of other phylogeny metrics in the context of
evolutionary computation is warranted. Additionally, future work should
evaluate the relationships between other conventionally used approaches to
measuring diversity (e.g., genetic diversity). As we improve our understanding
of these relationships, we may even be able to use them to make inferences
about fitness landscapes (Dolson et al.,, 2020). We hope that the results
presented here will inspire others to incorporate a phylogenetic perspective
into their evolutionary computation research.

##  5 Author contributions

ED conceptualized the questions and experiments in this chapter, wrote the
code for the non-exploration-diagnostic experiments, ran all experiments,
analyzed the data, and wrote the first draft of this chapter. JGH and AL wrote
the code for the exploration diagnostic experiments and assisted with the data
analysis and writing.

##  6 Acknowledgements

We thank members of the MSU ECODE lab, the MSU Digital Evolution lab, and the
Cleveland Clinic Theory Division for the conversations that inspired this
work. This research was supported by the National Science Foundation (NSF)
through the BEACON Center (Cooperative Agreement DBI-0939454). Michigan State
University provided computational resources through the Institute for Cyber-
Enabled Research. Any opinions, findings, and conclusions or recommendations
expressed in this material are those of the author(s) and do not necessarily
reflect the views of the NSF, UM, or MSU.

## References

  * Bressler and Seth, (2011) Bressler, S. L. and Seth, A. K. (2011).  Wiener–Granger Causality: A well established methodology.  NeuroImage, 58(2):323–329.
  * Dolson, (2021) Dolson, E. (2021).  Supplemental material for ”What can phylogenetic metrics tell us about useful diversity in evolutionary algorithms?” at GPTP 2021.  DOI: 10.5281/zenodo.4733407. URL: https://github.com/emilydolson/ecology_in_EC_parameter_sweep.
  * (3) Dolson, E., Banzhaf, W., and Ofria, C. (2018a).  Applying Ecological Principles to Genetic Programming.  In Banzhaf, W., Olson, R. S., Tozier, W., and Riolo, R., editors, Genetic Programming Theory and Practice XV, pages 73–88. Springer International Publishing, Cham.
  * Dolson et al., (2020) Dolson, E., Lalejini, A., Jorgensen, S., and Ofria, C. (2020).  Interpreting the Tape of Life: Ancestry-based Analyses Provide Insights and Intuition about Evolutionary Dynamics.  Artificial Life, 26(1):1–22.
  * Dolson et al., (2019) Dolson, E., Lalejini, A., and Ofria, C. (2019).  Exploring Genetic Programming Systems with MAP-Elites.  In Banzhaf, W., Spector, L., and Sheneman, L., editors, Genetic Programming Theory and Practice XVI, pages 1–16. Springer International Publishing, Cham.
  * Dolson et al., (2017) Dolson, E., Perez, S., Olson, R., and Ofria, C. (2017).  Spatial resource heterogeneity increases diversity and evolutionary potential.  bioRxiv.
  * (7) Dolson, E. L., Banzhaf, W., and Ofria, C. (2018b).  Ecological theory provides insights about evolutionary computation.  PeerJ Preprints, 6:e27315v1.
  * Goings et al., (2012) Goings, S., Goldsby, H. J., Cheng, B. H., and Ofria, C. (2012).  An ecology-based evolutionary algorithm to evolve solutions to complex problems.  Artificial Life, 13:171–177.
  * Goings and Ofria, (2009) Goings, S. and Ofria, C. (2009).  Ecological approaches to diversity maintenance in evolutionary algorithms.  In IEEE Symposium on Artificial Life, 2009. ALife ’09, pages 124–130.
  * Goldberg et al., (1987) Goldberg, D. E., Richardson, J., and Grefenstette, J. J. (1987).  Genetic algorithms with sharing for multimodal function optimization.  In Genetic algorithms and their applications: Proceedings of the Second International Conference on Genetic Algorithms, pages 41–49, Hillsdale, NJ. Lawrence Erlbaum.
  * Granger, (1969) Granger, C. W. J. (1969).  Investigating Causal Relations by Econometric Models and Cross-spectral Methods.  Econometrica, 37(3):424–438.  Publisher: [Wiley, Econometric Society].
  * Helmuth et al., (2016) Helmuth, T., McPhee, N. F., and Spector, L. (2016).  Lexicase Selection for Program Synthesis: A Diversity Analysis.  In Riolo, R., Worzel, W. P., Kotanchek, M., and Kordon, A., editors, Genetic Programming Theory and Practice XIII, Genetic and Evolutionary Computation, pages 151–167. Springer International Publishing.
  * Helmuth and Spector, (2015) Helmuth, T. and Spector, L. (2015).  General Program Synthesis Benchmark Suite.  In Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation, GECCO ’15, pages 1039–1046, New York, NY, USA. ACM.
  * Helmuth et al., (2015) Helmuth, T., Spector, L., and Matheson, J. (2015).  Solving Uncompromising Problems With Lexicase Selection.  IEEE Transactions on Evolutionary Computation, 19(5):630–643.
  * Hernandez et al., (2021) Hernandez, J. G., Lalejini, A., and Ofria, C. (2021).  An Exploration of Exploration: Measuring the ability of lexicase selection to find obscure pathways to optimality.  arXiv:2107.09760 [cs].  arXiv: 2107.09760.
  * Isaac et al., (2007) Isaac, N. J. B., Turvey, S. T., Collen, B., Waterman, C., and Baillie, J. E. M. (2007).  Mammals on the EDGE: Conservation Priorities Based on Threat and Phylogeny.  PLOS ONE, 2(3):e296.
  * Jackson, (2010) Jackson, D. (2010).  Promoting Phenotypic Diversity in Genetic Programming.  In Schaefer, R., Cotta, C., Kołodziej, J., and Rudolph, G., editors, Parallel Problem Solving from Nature, PPSN XI, Lecture Notes in Computer Science, pages 472–481, Berlin, Heidelberg. Springer.
  * Kassambara, (2020) Kassambara, A. (2020).  ggpubr: ’ggplot2’ Based Publication Ready Plots.  R package version 0.4.0.
  * Kauffman and Levin, (1987) Kauffman, S. and Levin, S. (1987).  Towards a general theory of adaptive walks on rugged landscapes.  Journal of Theoretical Biology, 128(1):11–45.
  * Metevier et al., (2019) Metevier, B., Saini, A. K., and Spector, L. (2019).  Lexicase Selection Beyond Genetic Programming.  In Banzhaf, W., Spector, L., and Sheneman, L., editors, Genetic Programming Theory and Practice XVI, Genetic and Evolutionary Computation, pages 123–136. Springer International Publishing, Cham.
  * Meyer, (2014) Meyer, P. E. (2014).  infotheo: Information-Theoretic Measures.  R package version 1.2.0.
  * Mouret and Doncieux, (2009) Mouret, J. and Doncieux, S. (2009).  Overcoming the bootstrap problem in evolutionary robotics using behavioral diversity.  In Evolutionary Computation, 2009. CEC’09. IEEE Congress on, pages 1161–1168. IEEE.
  * Ofria et al., (2018) Ofria, C., Dolson, E., Lalejini, A., Fenton, J., Jorgensen, S., Miller, R., Moreno, M. A., Stredwick, J., Zaman, L., Schossau, J., Gillespie, L., G, N. C., and Vostinar, A. (2018).  Empirical.  DOI: 10.5281/zenodo.1439475.
  * Ofria and Wilke, (2004) Ofria, C. and Wilke, C. O. (2004).  Avida: A software platform for research in computational evolutionary biology.  Artificial Life, 10(2):191–229.
  * R Core Team, (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Schreiber, (2000) Schreiber, T. (2000).  Measuring Information Transfer.  Physical Review Letters, 85(2):461–464.  Publisher: American Physical Society.
  * Sekanina and Bidlo, (2005) Sekanina, L. and Bidlo, M. (2005).  Evolutionary Design of Arbitrarily Large Sorting Networks Using Development.  Genetic Programming and Evolvable Machines, 6(3):319–347.
  * Spector, (2012) Spector, L. (2012).  Assessment of problem modality by differential performance of lexicase selection in genetic programming: a preliminary report.  In Proceedings of the 14th annual conference companion on Genetic and evolutionary computation, pages 401–408. ACM.
  * Tucker et al., (2017) Tucker, C. M., Cadotte, M. W., Carvalho, S. B., Davies, T. J., Ferrier, S., Fritz, S. A., Grenyer, R., Helmus, M. R., Jin, L. S., Mooers, A. O., Pavoine, S., Purschke, O., Redding, D. W., Rosauer, D. F., Winter, M., and Mazel, F. (2017).  A guide to phylogenetic metrics for conservation, community ecology and macroecology.  Biological Reviews, 92(2):698–715.
  * Tucker et al., (2018) Tucker, C. M., Davies, T. J., Cadotte, M. W., and Pearse, W. D. (2018).  On the relationship between phylogenetic diversity and trait diversity.  Ecology, 99(6):1473–1479.
  * Webb et al., (2002) Webb, C. O., Ackerly, D. D., McPeek, M. A., and Donoghue, M. J. (2002).  Phylogenies and Community Ecology.  Annual Review of Ecology and Systematics, 33(1):475–505.
  * Wickham, (2016) Wickham, H. (2016).  ggplot2: Elegant graphics for data analysis.  Springer-Verlag New York.
  * Yao and Li, (2020) Yao, C.-Z. and Li, H.-Y. (2020).  Effective Transfer Entropy Approach to Information Flow Among EPU, Investor Sentiment and Stock Market.  Frontiers in Physics, 8:206.


# Energy Optimal Point-to-Point Motion Profile Optimization

Nick Van Oosterwyck<EMAIL_ADDRESS>Foeke Vanbecelaere
111Equal contribution<EMAIL_ADDRESS>Ferre Knaepkens
<EMAIL_ADDRESS>Michael Monte<EMAIL_ADDRESS>Kurt
Stockman<EMAIL_ADDRESS>Annie Cuyt<EMAIL_ADDRESS>Stijn
Derammelaere<EMAIL_ADDRESS>Department of Electromechanics,
CoSys-Lab, University of Antwerp, Groenenborgerlaan 171, 2020 Antwerp, Belgium
Department of Electrical Energy, Metals, Mechanical Constructions and Systems,
Ghent University campus Kortrijk, 8500 Kortrijk, Belgium  Department of
Mathematics and Computer Science, University of Antwerp, Middelheimlaan 1,
2020 Antwerpen, Belgium  AnSyMo/CoSys, Flanders Make, the strategic research
center for the manufacturing industry  FlandersMake@UGent - Core Lab EEDT-MP
College of Mathematics and Statistics, Shenzhen University Shenzhen, Guangdong
518060, China

###### Abstract

Position-controlled systems driving repetitive tasks are of significant
importance in industrial machinery. The electric actuators used in these
systems are responsible for a large part of the global energy consumption,
indicating that major savings can be made in this field. In this context,
motion profile optimization is a very cost-effective solution as it allows for
more energy-efficient machines without additional hardware investments or
adaptions. In particular, mono-actuated mechanisms with position-dependent
system properties have received considerable attention in literature. However,
the current state-of-the-art methods often use unbounded design parameters to
describe the motion profile. This both increases the computational complexity
and hampers the search for a global optimum. In this paper, Chebyshev
polynomials are used to describe the motion profile. Moreover, the exact
bounds on the Chebyshev design parameters are derived. This both seriously
reduces the computational complexity and limits the design space, allowing the
application of a global optimizer such as the genetic algorithm. Experiments
validate the added value of the chosen approach. In this study, it is found
that the energy consumption can be reduced by 62.9% compared to a standard
reference motion profile.

###### keywords:

Motion profile optimization , Point-to-point motion , Energy efficiency ,
Validation , CAD model

###### MSC:

[2020] 49

##  1 Introduction

In the last decades, economic considerations and stricter government
regulations have driven engineers to come up with new techniques to reduce the
energy consumption of industrial machinery. Statistics indicate that electric
motors are generally responsible for about 2/3 of the industrial electricity
consumption, which indicates that major savings are to be made in this field
[1].

In this context, several technologies and methods have been developed to
reduce the electrical energy consumption of mechatronic systems. So-called,
motion profile optimization starts from the idea that in many industrial
applications, only part of the motion is constrained by the process
requirements. Hence, an optimization potential rises in the non-constrained
part of the position function $\theta(t)$, in between the start and endpoint
of the point-to-point (PTP) motions (Fig 1).

Figure 1: Motion profile of a PTP movement with constraints $\theta_{A}$,
$\theta_{B}$, $t_{A}$ and $t_{B}$.

By modifying the position function $\theta(t)$ between a fixed start-
$\theta(t_{A})=\theta_{A}$ and endpoint $\theta(t_{B})=\theta_{B}$, it is
possible to minimize one or several design objectives such as energy [2],
motion time [3, 4], machine cost [5] and vibrations [6]. Moreover, since these
industrial applications often involve repetitive movements, the motion profile
optimization effect will be perceptible every machine cycle, thus, making it
an indispensable step in modern energy-efficient machine design.

###  1.1 Related work

Dedicated machines with single-axis servo mechanisms are omnipresent in
industrial production processes, especially given the tendency to evolve
towards dedicated actuators for each machine movement [7]. In this context,
several approaches have been presented in the past literature to optimize the
motion profiles of these one DOF systems.

In [3, 6, 8, 9, 10], optimal motion profiles are obtained for mechanical
systems with constant load parameters such as the inertia $J$. Nevertheless,
as indicated in literature [11, 7, 12], it is essential to consider varying
loads to cover the majority of machine applications.

However, in those cases obtaining an analytic description of the position-
dependent system properties is a challenging task for industrial machine
designers. For instance, [12, 13, 14, 15] use Hamilton’s principle and
Lagrange multipliers to obtain the differential-algebraic equations of the
system. In [16], a method of virtual work is described to obtain the system
matrix while [17] determines the inertia profile using the method of kinetic
energy. Such approaches are cumbersome, complex and error-prone, and are not
easily applicable. Especially given the trend indicated in [18] that there is
a demand for methods that take into account the ease of implementation.

Fortunately, machine builders often already design their machines in 3D CAD
multibody software, which can be used to extract crucial information. Hence,
in [19] and [7], the authors of this paper describe a technique to derive the
position dependency of critical parameters inertia $J(\theta)$ and load torque
$\tau_{l}(\theta)$, based on only three CAD motion simulations. If the load
model is known, optimization algorithms can iterate on it to minimize energy
usage.

For what concerns these optimization algorithms, several approaches are
possible. On the one hand, [8] and [20] use an indirect approach such as
Pontryagin’s Maximum Principle to obtain the best possible control. However,
this method tends to be abandoned recently due to the small convergence area
and difficulties incorporating constraints [21].

On the other hand, direct approaches recast the optimization into a nonlinear
optimization problem, which can be solved with different numerical methods. In
particular, [2] and [11] use gradient-based methods such as Sequential
Quadratic Programming (SQP) or quasi-Newton methods that are known to have
very low solve times and good scalability. However, these algorithms can only
deliver local optimal solutions and are not suited for problems with multiple
minima. As indicated in [22], the optimum obtained with gradient-based methods
is greatly influenced by the selected starting points, which are to be chosen
arbitrarily.

To avoid this problem, heuristic optimization algorithms [19, 14, 12, 7] such
as generalized pattern search (GPS) or genetic algorithms (GA) are of
interest. In contrast to gradient-based algorithms that do not search the
entire design space, derivative-free algorithms like GA often sample a wide
part of the design space in order to be successful [23]. Nevertheless, because
these heuristic solvers do not exploit gradient information, they are not
computationally competitive with gradient-based methods [24].

Regarding the motion profile function, several papers rely on piecewise
position functions ([11, 7]), where either cubic, quintic, or trigonometric
splines are used. However, the objective functions in these works, are
characterized by many local minima, causing the risk of getting stuck in a
suboptimal solution. For instance, in [25], the usage of cubic splines
resulted in a savings difference of $18\%$ between the global and local
optimum.

On the other hand, continuous motion profile functions such as classic
polynomials ([19, 14, 12, 9, 6]) are also popular because they do not
introduce high jerk peaks into the system, which increases the wear of the
components. However, the resulting optimization problem is known to be badly
conditioned. For example in [19], the coefficients reached values up to
$1.8\,10^{20}$. Therefore, the authors of this paper proposed Chebyshev
polynomials and rescalings in [2] to obtain a more numerically stable problem
formulation and to increase the robustness against getting stuck in local
minima. Moreover, in contrast to the classical polynomial descriptions, the
design variables of Chebyshev polynomials can be bounded as is shown in this
paper.

###  1.2 Method

In light of the considerations mentioned above, a CAD-based method for
computing energy-optimal PTP (Point-to-Point) motion profiles of single DOF
mechanisms using Chebyshev polynomials has been previously presented by the
authors in [2]. However, due to the numerous symbolic calculations involved in
constructing the objective function, solve times of almost 2 hours were
reported. In addition, the solutions in [2] were obtained using gradient-based
solvers which have a high risk of getting stuck in local minima. Finally, only
theoretical reductions were reported, thus, leaving the feasibility of the
proposed motion profiles undetermined. This paper builds upon these previous
results by providing five critical improvements:

  * 1.

In order to reduce the computational burden, a discrete approach is presented
which eliminates the use of symbolic operations. To do so, the discrete system
property data which originates from the CAD motion simulations have to be
properly rescaled and interpolated.

  * 2.

As an accurate model of the system dynamics is crucial for a correct
optimization, the dynamics of the mechanism are extended by including damping
and friction into the optimization routine. In addition, a new identification
procedure is described which is validated on an industrial case.

  * 3.

A derivation for exact bounds of the Chebyshev polynomial coefficients is
introduced. This allows limiting the feasible design space. The latter is
essential for heuristic optimizers to reduce their computation time and the
chances that the global optimum remains unidentified.

  * 4.

To check the robustness of the proposed method against getting stuck in local
optima, the resulting optimization problem is solved with both a fast
gradient-based and a global heuristic solver (i.e. GA).

  * 5.

Experimental tests have been carried out on an industrial pick-and-place unit
to quantify the actual measured energy savings.

##  2 System Modeling

The complete mechatronic system can be divided into two subsystems (Fig. 2).
On the one hand, there is the mechanical subsystem which describes the
dynamics of a generic single-axis system. For high dynamical applications,
these systems usually consist of slider-crank mechanisms and four-bar linkages
[7]. Nevertheless, the approach is valid to any position-controlled system
where the mechanism is driven by a single actuator.

On the other hand, there is the actuation subsystem which converts the
electrical energy into mechanical energy and drives the mechanism. For the
envisaged position-controlled systems, PMSMs are becoming the industry
standard for rotary applications, whereas linear motors are used for fast and
precise linear movements [26]. In Fig. 2, the PMSM actuator is represented by
an equivalent DC model.

Figure 2: Schematic of the q-axis of a single axis mechanism.

###  2.1 Mechanical subsystem

The dynamics of a single axis DOF mechanism can be described by means of the
torque equation [11, 27] :

$\tau_{m}(t)=\tau_{l}(\theta)+\underbrace{J(\theta)\ddot{\theta}}_{\tau_{a}}+\underbrace{\frac{1}{2}\frac{\mathrm{d}J(\theta)}{\mathrm{d}\theta}(\dot{\theta})^{2}}_{\tau_{v}}+\tau_{f}(\dot{\theta})\,.$
(1)

With reference to Fig. 2 and equation (1), let us define $\theta=\theta(t)$ as
the Lagrangian coordinate which describes the angular position of the main
driving axis as a function of time $t$. The motor torque $\tau_{m}(\theta)$ is
defined as the driving torque generated by the motor. The load torque
$\tau_{l}(\theta)$ contains both gravitational forces as well as external
process powers that act on the mechanism.

Furthermore, all inertias of the mechanism’s components are related to the
main driving axis resorting to the concept of reduced moment of inertia.
Therefore, the reduced inertia of the complete system $J(\theta)$ is defined
as a combination of the reduced load inertia $J_{l}(\theta)$ and inertia of
the motor shaft itself $J_{m}$. Note that the position-dependent inertia of
the system $J(\theta)$, results in two torque components when it is reduced to
the motor side. The acceleration torque $\tau_{a}$ represents the part of the
motor torque responsible for the motor acceleration forces that arise during
the movement, while the variation torque $\tau_{v}$ compensates for the
variation of inertia in the system.

Finally, the frictional torque $\tau_{f}(\dot{\theta})$ is defined as the
result of frictional forces such as, for instance, viscous brush friction or
dry bearing friction in the motor bearings and mechanical system. A commonly
used model of friction shows three components of force: Coulomb (sliding)
friction, viscous damping, and static friction [28]. Regarding the PMSM, as
indicated in [29, p. 175], the only appreciable friction effect in operation
is viscous friction. Thus, coulomb and static frictions can be neglected in
the PMSM model. For what concerns the mechanical model, only the viscous
damping is modeled since the other friction components are constant and will
not have an effect on the optimal motion profile:

$\tau_{f}(\dot{\theta})=\mu_{v}\dot{\theta}\,,$ (2)

with $\mu_{v}$ the equivalent viscous friction coefficient.

The key benefit of the formulation in (1) is that it permits to model every
possible mechanism with a known geometry and allows to define a generic
optimization approach.

###  2.2 Actuation subsystem

Concerning the dynamics of the PMSM as depicted in Fig. 2 (represented by an
equivalent DC model), the electromechanical behavior can be easily described
by the following basic laws [30, p. 843]:

$\tau_{m}=k_{t}i\,,$ (3)

$u=Ri+L\frac{\mathrm{d}i}{\mathrm{d}t}+e=Ri+L\frac{\mathrm{d}i}{\mathrm{d}t}+pk_{v}\dot{\theta}\,,$
(4)

with electric back emf $e$, resistance $R$, back emf constant $k_{v}$, motor
torque constant $k_{t}$, and number of pole pairs $p$, which can be found in
the motor data sheet.

In equation (4), the voltage drop $L\frac{\mathrm{d}i}{\mathrm{d}t}$ due to
the armature inductance is omitted as the mean value of its reactive power
will be zero and therefore does not contribute to the system’s energy need
[11].

Depending on whether the electric power flows from the drive unit to the
PMSM’s or vice versa, the PMSM operates in respectively motor or generator
mode. In this latter condition, depending on the capabilities of the drive
unit, the generated electric power can be either stored in a capacitor,
dissipated as heat on a braking resistance, or transferred back to the energy
source. Recent commercial PMSM drives are sized so that no electric power is
actually dissipated during normal functioning so that the braking resistance
is actually activated only under emergency conditions [7]. Therefore, in what
follows, it is assumed that all the generated energy is returned to the grid
and no losses occur in the process.

For a correct model of the actuation subsystem and prediction of the energy
usage, it is important to model other losses such as cooling fans and drive
circuitry as well. Nevertheless, the power consumption of these devices is
generally considered constant and is therefore not affected by the motion
profile [31].

In order to minimize the total energy need $E$ of the application, it is
crucial to quantify the input energy of the complete system. Therefore,
similar to [7], a formulation of the input electrical energy $E$ is derived
and a torque-based design objective is obtained which allows minimizing the
energy solely based on the mechanical parameters.

Starting from equations (3) and (4), the instantaneous power $P_{e}$ is
defined as

$P_{e}=u\,i=\frac{R}{k_{t}^{2}}\tau_{m}^{2}+\frac{pk_{v}}{k_{t}}\tau_{m}\dot{\theta}\,.$
(5)

The motion profile is defined on the time interval $t\in[t_{A},t_{B}]$ and
must have zero initial and final speed and acceleration, i.e.
$\dot{\theta}(t_{A})=\dot{\theta}(t_{B})=\ddot{\theta}(t_{A})=\ddot{\theta}(t_{B})=0$.
The total energy can be expressed as

$\displaystyle E=\int_{t_{A}}^{t_{B}}P_{e}\,\mathrm{d}t$
$\displaystyle=\int_{t_{A}}^{t_{B}}\left[\frac{R}{k_{t}^{2}}\tau_{m}^{2}+\frac{pk_{v}}{k_{t}}\tau_{m}\dot{\theta}\right]\,\mathrm{d}t\,.$
(6)

Then, by incorporating the torque equation from (1), the total energy of the
motion is given by:

$\displaystyle
E=\frac{pk_{v}}{k_{t}}\underbrace{\int_{t_{A}}^{t_{B}}(\tau_{a}+\tau_{v})\dot{\theta}\,\mathrm{d}t}_{E_{k}}+\frac{pk_{v}}{k_{t}}\underbrace{\int_{t_{A}}^{t_{B}}\tau_{l}\,\dot{\theta}\,\mathrm{d}t}_{E_{p}}$
(7)
$\displaystyle+\underbrace{\int_{t_{A}}^{t_{B}}\left[\frac{R}{k_{t}^{2}}\tau_{m}^{2}+\frac{pk_{v}}{k_{t}}\tau_{f}\dot{\theta}\right]\mathrm{d}t}_{E_{l}}\,.$

Here, the first term $E_{k}$ represents the kinetic energy of the moving
masses in the system. Due to the rest-to-rest motion of the envisaged
applications, this term reduces to zero. Further, the term $E_{p}$ represents
the potential energy stored in the system. As this term $E_{p}$ only depends
on the fixed start $\theta_{A}$ and end position $\theta_{B}$, it is
disregarded in the optimization routine [7]. The final term $E_{l}$ represents
the energy that is lost due to the coil resistance and frictional forces and
is the only term that is affected by optimizing the motion profile
$\theta(t)$. Nevertheless, in many industrial applications, the frictional
forces $\tau_{f}$ are negligible [8], especially if the inertial loads are
predominant. Thus, the energy losses $E_{l}$ can be expressed as:

$E_{l}=\int_{t_{A}}^{t_{B}}\frac{R}{k_{t}^{2}}\tau_{m}^{2}\mathrm{d}t=\frac{R\Delta
t}{k_{t}^{2}}\tau_{rms}^{2}\,,$ (8)

where $\tau_{rms}$ is the RMS value of the motor Torque $\tau_{m}$. This
proves that the RMS torque $\tau_{rms}$ can be effectively used as an
optimization objective to minimize the total energy usage of the system. This
is very useful in situations where the motor coil properties are unknown or
where parameters are missing [7].

##  3 Identification

###  3.1 Inertia and Load Torque

Identification of all the position varying parameters in the highly nonlinear
differential torque equation (1) is not straightforward. Fortunately, machine
builders design their machines in 3D CAD multibody software. For this reason,
[7] and [19] describe a technique to derive the position dependency of
critical parameters inertia $J(\theta)$ and load torque $\tau_{l}(\theta)$,
based on three CAD motion simulations. (Fig. 3).

Figure 3: Schematic overview of the procedure for extracting position-
dependent properties $J(\theta)$ and $\tau_{l}(\theta)$ based on three
different CAD motion simulations [19].

In this paper, the identification routine is illustrated by applying it to an
industrial pick-and-place unit (Fig. 4) that performs repetitive movements
between start point $A$ with angular position $\theta_{A}=0$ and endpoint $B$
with angular position $\theta_{B}=173.6^{\circ}$. The resulting inertia
$J(\theta)$ and load torque $\tau_{l}(\theta)$ profiles are presented in Fig.
5. Because of the machine position limits $\theta_{A}$ and $\theta_{B}$, only
the green shaded part of the system properties is relevant during operation.

Figure 4: Experimental set-up (left) and schematic overview (right) of the
pick- and place unit. Figure 5: Values of system properties inertia
$J(\theta)$ and load torque $\tau_{l}(\theta)$.

###  3.2 Viscous Friction Coefficient

Once the system properties $J(\theta)$ and $\tau_{l}(\theta)$ are determined,
the only indefinite term in Eq. (1) is the friction torque $\tau_{f}$, and
more specifically $\mu_{k}$. In the previous description of the energy flows,
the friction torque $\tau_{f}$ was neglected, leading to a simple objective
(i.e. $\tau_{rms}$) to quantify the energy consumption. However, it is
important to verify this statement for the intended setup. Therefore this
section describes a method to quantify the frictional forces $\tau_{f}$.

Since the viscous friction coefficient $\mu_{k}$ parameter is highly dependent
on the practical setup, it is often only possible to determine this parameter
experimentally. Therefore, a first measurement is carried out by using an
arbitrary motion profile $\theta^{*}(t)$ as a set point and recording the
resulting actual motor torque $\tau_{m}^{e}(t)$ and position $\theta^{e}(t)$.
The arbitrary motion profile $\theta^{*}(t)$ can be determined by using a
default motion law such as a trapezoidal or s-curve profile.

After this measurement, a least squares fit can be used to determine the
experimental value of $\mu_{k}$, by fitting the torque model
$\tau_{m}(\theta^{e}(t),\mu_{k})$ with the measured torque $\tau_{m}^{e}$.

However, using the measured position $\theta^{e}$ and its time derivatives
$\dot{\theta}^{e}$, $\ddot{\theta}^{e}$ in the torque equation (1) leads to
unfeasible results since the derivatives amplify any noise that is present in
the measurement. Therefore, the measured position $\theta^{e}(t)$ is fitted
with an $n$-th degree polynomial $\theta^{p}(t)=\sum_{i=1}^{n}a_{i}t^{i}$ and
differentiated symbolically to smooth out any noise.

The friction parameter $\mu_{v}$ is thus determined by comparing the measured
torque $\tau_{m}^{e}(t)$ with virtual model and fitted motion profile
$\tau_{m}(\theta^{p}(t))$:

$\displaystyle\underset{\mu_{v}\,\in\,\mathbb{R}}{\text{minimize}}$
$\displaystyle||\tau_{m}^{e}(t)-\tau_{m}(\theta^{p}(t),\mu_{v})||_{2}\,.$ (9)

For the pick-and-place unit, a viscous damping coefficient of $0.0157Nms/rad$
was found. In Fig. 6, a comparison of the measured $\tau_{m}^{e}(t)$ and
virtual $\tau_{m}(\theta^{p}(t))$ torque is presented. The difference between
the virtual torques with and without friction is minimal, which indicates that
the friction can be neglected for the present case. The graph also shows a
close correlation between the virtual and measured torque, which indicates
that the virtual model can be effectively used to minimize the RMS torque
$\tau_{rms}$ and, by extension, the energy consumption $E$.

Figure 6: Comparison of the virtual $\tau_{m}(\theta^{p}(t))$ and measured
torque $\tau_{m}^{e}(t)$ (with and without friction).

##  4 Optimization Approach

###  4.1 Motion Profile Definition & Rescaling

In this paper, a Chebyshev polynomial $\sum_{i=0}^{n}p_{i}T_{i}(x)$ is used to
define the position profile $\theta(t)$, where $t\in[t_{A},t_{B}]$, in between
the start- ($\theta(t_{A})=\theta_{A}$) and endpoint
($\theta(t_{B})=\theta_{B}$) of the motion task. The sequence of orthogonal
Chebyshev polynomials $T_{k}(x)=T_{k}(\cos(\vartheta))$, defined on the
interval $x\in[-1,1]$, is obtained from the recurrence relation:

$\displaystyle T_{0}(x)$ $\displaystyle=1,\quad T_{1}(x)=x,$ (10)
$\displaystyle T_{k+1}(x)$ $\displaystyle=2xT_{k}(x)-T_{k-1}(x),$

Alternatively, the polynomials can be derived from the trigonometric
definition, which gives exactly the same results:

$T_{k}(x)=T_{k}(\cos(\vartheta))=\cos(k\vartheta).$ (11)

To use $T_{n}(x)$ as a representation for the position profile, a linear
transformation from $t$ into the range $[-1,1]$ of $x$ is required [32]:

$\displaystyle t$
$\displaystyle=\frac{1}{2}(t_{B}-t_{A})x+\frac{1}{2}(t_{B}+t_{A})=ax+b,$ (12)

where scale factors $a$ and $b$ are defined for the purpose of the following
paragraphs. In addition, the position $\theta\in[\theta_{A},\theta_{B}]$ is
also rescaled to the interval $\phi\in[-1,1]$, which makes it possible to
obtain strict bounds on the design space in (32). Thus, the rescaled motion
profile description $\phi(x)$ of degree $n$ with optimizable coefficients
$\mathbf{p}=[p_{0},p_{1},\ldots,p_{n}]^{T}$ is obtained.

$\phi(x)=\sum_{i=0}^{n}p_{i}T_{i}(x),\quad x\in[-1,1].$ (13)

The output of the motion simulations in the previous section deliver $n_{s}$
samples of inertia $\mathbf{J}=[J_{1},\ldots,J_{n_{s}}]^{T}$, load torque
$\boldsymbol{\uptau}_{l}=[\tau_{l,1},\ldots,\tau_{l,n_{s}}]^{T}$ and
corresponding angle query points
$\boldsymbol{\uptheta}=[\theta_{1},\ldots,\theta_{n_{s}}]^{T}$. Due to the
position rescaling of the motion profile $\phi(x)$, the angle query points
$\boldsymbol{\uptheta}$ have to be rescaled accordingly:

$\displaystyle\boldsymbol{\upphi}$
$\displaystyle=\frac{2}{(\theta_{B}-\theta_{A})}\,\boldsymbol{\uptheta}-\frac{(\theta_{B}+\theta_{A})}{(\theta_{B}-\theta_{A})}=c\,\boldsymbol{\uptheta}+d.$
(14)

Moreover, as the property description is now defined on the rescaled interval
$\phi\in[-1,1]$, the following relationship holds with regard to the
derivative properties such inertia variation
$\frac{\mathrm{d}J(\phi)}{\mathrm{d}\phi}$:

$\frac{\mathrm{d}J(\phi)}{\mathrm{d}\phi}=\frac{1}{2}(\theta_{B}-\theta_{A})\frac{\mathrm{d}J(\theta)}{\mathrm{d}\theta}=e\,\frac{\mathrm{d}J(\theta)}{\mathrm{d}\theta}.$
(15)

When using the rescaled position profile $\phi(x)$, it is important to rescale
the torque equation (1) as well. Otherwise, the resulting values of the torque
profile $\tau(x)$ are distorted which results in different objective values
(i.e. $\tau_{rms}$) and solutions. To preserve the motor torque’s absolute
values, the following rescaled torque equation is introduced:

$\tau_{m}(x)=\tau_{l}(\phi)+\frac{1}{2}\frac{\mathrm{d}J(\phi)}{\mathrm{d}\phi}\frac{1}{e}\left(\frac{\dot{\phi}}{a.c}\right)^{2}+J(\theta)\frac{\ddot{\phi}}{a^{2}.c}+\mu_{k}\frac{\dot{\phi}}{a.c}.$
(16)

An overview of the position and torque rescalings is presented in Fig. 7. The
new system equation (16) ensures the system dynamics are equally scaled and
the minima are not altered.

Figure 7: Original $\theta(t)$ and rescaled position profiles $\theta(x)$,
$\phi(x)$ with their corresponding torque equations.

For what concerns the constraints, the rest-to-rest motion requires zero speed
$\dot{\phi}$ and acceleration $\ddot{\phi}$ in the start and endpoint:

$\begin{array}[]{ccccc}\phi(-1)=-1&,&\dot{\phi}(-1)=0&,&\ddot{\phi}(-1)=0,\\\
\phi(1)=1&,&\dot{\phi}(1)=0&,&\ddot{\phi}(1)=0.\end{array}{}$ (17)

Referring to (13), and by incorporating the motion profile constraints (17),
the lower degree coefficients $[p_{0},...,p_{5}]^{T}$ can be written as a
function of the remaining coefficients $[p_{6},...,p_{n}]^{T}$, such that
$n-5$ degrees of freedom (DOF) are kept available for the optimization
algorithm [12]. Thus, the energy optimal motion profile problem is formulated
as the following minimization problem with design variable vector
$\mathbf{o}=[p_{6},...,p_{n}]^{T}$:

$\displaystyle\underset{\mathbf{o}\,\in\,\mathbb{R}^{n-5}}{\text{minimize}}$
$\displaystyle\tau_{rms}=\sqrt{\frac{1}{2}\int_{-1}^{1}{\tau_{m}(\phi(x,\mathbf{o}))}^{2}\,\mathrm{d}x}.$
(18)

In some applications, an additional constraint of zero jerk in the begin and
endpoint can be imposed to limit the vibrations:

$\dddot{\phi}(-1)=0\quad;\quad\dddot{\phi}(1)=0.$ (19)

Because of these two extra equations, the DOF is reduced to $n-7$ and the
design variable vector can be expressed as $\mathbf{o}=[p_{8},...,p_{n}]^{T}$.

###  4.2 Initialization & Design Space

In this paper, the resulting optimization problem is solved with both a fast
gradient-based solver, the BFGS (Broyden–Fletcher–Goldfarb–Shanno) quasi-
Newton method [33], and a global heuristic solver, the genetic algorithm [34].

For gradient-based optimization, a starting point needs to be defined. The use
of the Chebyshev basis $T_{i}(x)$ in representation (13) allows initializing
the optimization parameter vector at zero since the coefficients in a
convergent Chebyshev series development of the motion profile function
$\phi(x)$ would converge to zero [35]. Here, we can safely assume some similar
behavior for the coefficients $p_{i}$ in (13).

For what concerns the genetic algorithm, a similar approach is used for the
initialization of the population. However, because a GA often samples a wide
part of the design space [23], it is beneficial to determine the exact bounds
on the design vector $\mathbf{o}$. By doing so, the solver can cover a large
part of the design space and reveal the global optimal solution. In the
following paragraphs, thanks to the rescaled Chebyshev motion profile
$\phi(x)$, strict bounds on the design vector $\mathbf{o}$ can be derived.

To define these bounds, we take a look at the projection of the position
profile $\phi(x)$ onto the orthogonal Chebyshev polynomial basis $T_{l}(x)$.
Given that $x=\cos(\theta)$, we introduce the inner product $F$:

$\displaystyle F=\langle\phi(x),T_{l}(x)\rangle$
$\displaystyle=\int\limits_{-1}^{1}\frac{\phi(x)T_{l}(x)}{\sqrt{1-x^{2}}}\,\mathrm{d}x$
(20)
$\displaystyle=\int\limits_{0}^{2\pi}\phi(\cos\theta)T_{l}(\cos\theta)\,\mathrm{d}\theta.$

Then, by taking into account the position function definition (13), we find
the following result:

$\displaystyle F$
$\displaystyle=\int\limits_{0}^{2\pi}\left(\sum\limits_{k=0}^{n}p_{k}T_{k}(\cos\theta)\right)T_{l}(\cos\theta)\,\mathrm{d}\theta$
(21)
$\displaystyle=\sum\limits_{k=0}^{n}p_{k}\int\limits_{0}^{2\pi}T_{k}(\cos\theta)T_{l}(\cos\theta)\,\mathrm{d}\theta.$

Here, the integral
$I=\int\limits_{0}^{2\pi}T_{k}(\cos\theta)T_{l}(\cos\theta)\,\mathrm{d}\theta$
can be further simplified by using the Chebyshev polynomial orthogonality
properties, which are rederived here for the sake of readability. Because of
Eq. (11) and by using the inverse Simpson rule of trigonometry, the integral
$I$ can be written as:

$\displaystyle I$
$\displaystyle=\int\limits_{0}^{2\pi}\cos(k\theta)\cos(\ell\theta)\,\mathrm{d}\theta$
(22)
$\displaystyle=\frac{1}{2}\int\limits_{0}^{2\pi}\cos\big{(}(k+\ell)\theta\big{)}\,\mathrm{d}\theta+\frac{1}{2}\int\limits_{0}^{2\pi}\cos\big{(}(k-\ell)\theta\big{)}\,\mathrm{d}\theta.$

This integral can be split into three cases:

  1. 1.

$k=\ell=0$

$I=2\pi,$ (23)

  2. 2.

$k=\ell\neq 0$

$\displaystyle I=\pi,$ (24)

  3. 3.

$k\neq\ell$

$\displaystyle I=0.$ (25)

Thus, by taking into account (25), only the term for which $k=l$ remains in
the summation $F$:

$\displaystyle
F=p_{\ell}\int\limits_{0}^{2\pi}\cos^{2}(\ell\theta)\,\mathrm{d}\theta.$ (26)

This can be split into two cases. For $\ell=0$ and by making use of (23) and
(20) we find:

$p_{0}=\dfrac{1}{2\pi}\int\limits_{0}^{2\pi}\phi(\cos\theta)\,\mathrm{d}\theta,$
(27)

and for $\ell>0$, by making use of (24) and (20):

$p_{\ell}=\dfrac{1}{\pi}\int\limits_{0}^{2\pi}\phi(\cos\theta)\cos(\ell\theta)\,\mathrm{d}\theta.$
(28)

For $\theta\in[0,2\pi]$, $\cos\theta$ lies in interval $[-1,1]$. Because of
the position rescalings of the motion profile $\phi(x)$, the image
$\phi(\cos\theta)$ also lies in the interval $[-1,1]$. Thus, we find:

$|p_{0}|\leq\dfrac{1}{2\pi}\int\limits_{0}^{2\pi}|\phi(\cos\theta)|\,\mathrm{d}\theta\leq\dfrac{1}{2\pi}\int\limits_{0}^{2\pi}\,\mathrm{d}\theta=1.$
(29)

and

$|p_{\ell}|\leq\dfrac{1}{\pi}\int\limits_{0}^{2\pi}|\phi(\cos\theta)||\cos(\ell\theta)|\,\mathrm{d}\theta\leq\dfrac{1}{\pi}\int\limits_{0}^{2\pi}|\cos(\ell\theta)|\,\mathrm{d}\theta.$
(30)

To calculate this last integral, we use the periodicity of the function
$\cos(\ell\theta)$. This function has a period of $2\pi/\ell$, so goes $\ell$
times up and down on the interval $[0,2\pi]$. So, after taking the absolute
value of this function, we find $2\ell$ times the integral over the positive
part of a period, for example, the interval $[-\pi/2\ell,\pi/2\ell]$:

$\displaystyle\dfrac{1}{\pi}\int\limits_{0}^{2\pi}|\cos(\ell\theta)|\,\mathrm{d}\theta=\dfrac{2\ell}{\pi}\int\limits_{-\pi/2\ell}^{\pi/2\ell}\cos(\ell\theta)\,\mathrm{d}\theta=\dfrac{4}{\pi}.$
(31)

Thus, the following bounds for the coefficients $p_{i}$ are obtained:

$|p_{0}|\leq 1\hskip 14.22636pt\text{ and }\hskip
14.22636pt|p_{\ell}|\leq\dfrac{4}{\pi},\hskip 8.5359pt\ell=1,\dots,n.$ (32)

These constraints on the design space simplify the subsequent optimization.

##  5 Results

###  5.1 Motion Profile Optimization

In order to assess the performance of the proposed method, a set of
optimizations has been performed on the industrial pick-and-place unit
depicted in Fig. 4. The mechanism is required to move between its start
position $\theta_{A}$ of $0^{\circ}$ and end position $\theta_{B}$ of
$173.6^{\circ}$ and has a motion time $\Delta t$ of $73.5ms$. As for the
constraint, two different cases are considered, namely

  * 1.

Jerk Free (JF): Only the boundary constraints of (17) are taken into account.
The corresponding rescaled Chebyshev position profile $\phi(x)$ of degree $n$
is hereafter referred to as cheb”n”. A 5th-degree polynomial, hereafter
indicated as poly5, is taken as the reference motion profile for comparison
purposes. This is the smallest degree polynomial that satisfies the
constraints.

  * 2.

Jerk Zero (J0): In addition to the constraint of a jerk-free optimization, a
zero-jerk constraint is added in the start, and endpoint (19) is added to the
motion profile definition. The resulting $n$-th degree position profile
$\phi(x)$ is referred to as cheb”n”J0. The reference motion profile is in this
case a 7th-degree polynomial, hereafter referred to as poly7J0.

For every case, the resulting optimization problem is solved in a MATLAB
environment for degrees $n=7,9,11,$ and $13$. The results are presented in
Fig. 8 and Tables 1 & 2 where for every motion profile, the corresponding RMS
torque $\tau_{rms}$ and solve time $t_{sol}$ are displayed. Savings up to
$54.4\%$ are obtained in under $0.77$ s. The results clearly converge towards
a minimal value for increasing degree $n$. In general, the motion profiles
which include the jerk constraint (19), have slightly bigger $\tau_{rms}$
values, which is to be expected due to the fact that this extra constraint
limits the acceleration near the endpoints while it is desirable to have high
accelerations here since the inertia is low.

In Table 1, the $\tau_{rms}$ values of a conventional trapezoidal 1/3 motion
profile are presented as well, which accelerates during 1/3rd of the time,
moves at a constant speed during 1/3rd, and decelerates at the last 1/3rd [8].
What is interesting in this table is that the torque demand can already be
significantly reduced by selecting an adequate default motion law.
Notwithstanding that the greatest savings are realized after optimization.

It is worth noting that for the jerk-free motion profiles, the same solution
was found for both the genetic algorithm and gradient-based solver. However,
the calculation times with GA are considerably higher. When including the jerk
constraint, the GA comes close but does not completely reveal the full
optimization potential. Therefore, for what concerns the present study,
gradient-based optimizations algorithms are preferable. Since the GA did not
obtain a better solution for any motion profile in the bounded search space,
we can expect that the results obtained with the gradient-based method are
global optimal solutions.

Although only the forward motion is considered here, similar results can be
obtained for the return motion by simply changing the position constraints.

Table 1: Results of the motion profile optimization (Jerk Free). | Gradient-Based | Genetic Algorithm
---|---|---
JF | $\tau_{rms}\,[Nm]$ | $t_{sol}\,[s]$ | $\tau_{rms}\,[Nm]$ | $t_{sol}\,[s]$
poly5 (ref.) |  | 22.48
---
- |  | 22.48
---
-
trap |  | 17.16
---
-23.7%
- |  | 17.16
---
-23.7%
-
cheb7 |  | 13.78
---
-38.7%
0.21 |  | 13.78
---
-38.7%
3.28
cheb9 |  | 12.47
---
-44.5%
0.32 |  | 12.47
---
-44.5%
40.33
cheb11 |  | 12.33
---
-45.2%
0.51 |  | 12.33
---
-45.2%
67.05
cheb13 |  | 12.29
---
-45.4%
1.06 |  | 12.29
---
-45.4%
142.34
Table 2: Results of the motion profile optimization (Jerk 0). | Gradient-Based | Genetic Algorithm
---|---|---
J0 | $\tau_{rms}\,[Nm]$ | $t_{sol}\,[s]$ | $\tau_{rms}\,[Nm]$ | $t_{sol}\,[s]$
poly7J0 (ref.) |  | 28.44
---
- |  | 28.44
---
-
cheb9J0 |  | 16.12
---
-43.3%
0.27 |  | 16.12
---
-43.3%
6.15
cheb11J0 |  | 13.61
---
-52.2%
0.38 |  | 14.11
---
-50.4%
175.23
cheb13J0 |  | 12.98
---
-54.4%
0.77 |  | 13.15
---
-53.8%
195.02
Figure 8: Results of the motion profile optimization for different degrees
$n$.

###  5.2 Measurements

The theoretical results are validated against experimental measurements on the
pick-and-place unit (Fig. 4). The setup comprises a Beckhoff CX5140 PLC, a
Beckhoff AX5901 motor drive, and a Beckhoff AM3064 PMSM, which is connected to
the shaft of the mechanism. In order to measure the input electrical energy, a
Tektron PA4000 power analyzer is used to analyze the power supply (Fig. 9).

Figure 9: Schematic overview of the experimental setup.

The theoretical savings potential of the motion profile optimization is only
fulfilled when the motor is capable of following the optimized position
setpoint. Therefore, a performant motion controller needs to be designed in
order to keep the tracking error as low as possible. Here, similar to [19], a
cascade controller with torque and speed feedforward is employed as it has
proven to be successful for high dynamic systems. The look-up table for the
feedforward torque is determined using the torque equation (16).

Figure 10: Schematic overview of the cascade motion controller with
feedforward [19].

In Tables 3 and 4, the results of both the measured RMS torque $\tau_{rms}$
and measured input electrical energy $E$ for different motion profiles are
presented. As expected from the simulations, the lowest absolute energy
consumption is obtained when using jerk-free motion profiles. When the jerk
constraint is active, a decrease of 62.9% in energy consumption can be
achieved by optimizing the motion profile, while a relative saving of 52.5% is
possible if no extra constraint on the jerk is imposed.

The measured $\tau_{rms;meas}$ and calculated RMS motor torque $\tau_{rms}$
show a very high similarity, which confirms that the present system model is
valid.

Table 3: Experimental results with energy measurement (Jerk Free). JF | $\tau_{rms}\,[Nm]$ | $\tau_{rms;meas}\,[Nm]$ | $E_{meas}\,[Wh]$
---|---|---|---
poly5 | 22.48 | 19.59 | 312.2
trap |  | 17.16
---
-23.7%
| 15.88
---
-18.98%
| 215.1
---
-31.1%
cheb7 |  | 13.78
---
-38.7%
| 13.40
---
-31.6%
| 181.7
---
-41.8%
cheb9 |  | 12.47
---
-44.5%
| 12.07
---
-38.4%
| 152.3
---
-51.2%
cheb11 |  | 12.33
---
-45.2%
| 11.93
---
-39.1%
| 150.1
---
-51.9%
cheb13 |  | 12.29
---
-45.4%
| 11.83
---
-39.6%
| 148.2
---
-52.5%
Table 4: Experimental results with energy measurement (Jerk Zero). J0 | $\tau_{rms}\,[Nm]$ | $\tau_{rms;meas}\,[Nm]$ | $E_{meas}\,[Wh]$
---|---|---|---
poly7J0 | 28.44 | 25.30 | 458.5
cheb9J0 |  | 16.12
---
-43.3%
| 15.81
---
-37.5%
| 222.9
---
-51.4%
cheb11J0 |  | 13.61
---
-52.2%
| 13.08
---
-48.3%
| 170.3
---
-62.9%
cheb13J0 |  | 12.98
---
-54.4%
| 12.72
---
-49.7%
| 170.8
---
-62.7%

##  6 Conclusion

This study proposes a novel approach for motion profile optimization of PTP
motions with Chebyshev polynomials. At first, system properties have been
extracted from both CAD motion simulations and measurements to obtain an
accurate virtual twin of the system. A Chebyshev motion profile with scaling
laws is presented. Especially novel in this paper is the derivation of the
boundary conditions of this profile which enables to define bounds for the
design variables. The latter allows to use an optimizer that is designed to
obtain globally optimal solutions, i.e. Genetic Algorithm. In addition, the
solutions are validated with fast gradient-based optimization algorithms.
Finally, experimental optimization results have been considered to verify the
feasibility of the proposed solutions.

The numerical results, achieved on an exemplary model, clearly show that large
$\tau_{rms}$ savings of up to 53.8% can be achieved. In addition, it is shown
that by employing Chebyshev polynomials for the motion profile, a fast
gradient-based optimization can be effectively employed with solve times under
$0.8$s. At last, the validation measurements show that similar savings are
obtained on the real machine with a maximum energy reduction of $62.9\%$.

Due to the straightforward implementation of both the optimization itself and
integration of the resulting motion profiles in the motor drive, the proposed
method can be easily adopted in any existing configuration where the CAD is
data available. Therefore, the proposed method is expected to have a
beneficial impact on the energy usage of the envisaged PTP applications.

## Acknowledgements

Research funded by a PhD grant of the Research Foundation Flanders (FWO)
[1S88120N].

## References

  * [1] H. Bo, Electric Motors – Alignment of Standards and Best Practice Programmes with APEC, Final Report (2008).
  * [2] N. Van Oosterwyck, A. Ben yahya, A. Cuyt, S. Derammelaere, CAD Based Trajectory Optimization of PTP Motions using Chebyshev Polynomials, in: 2020 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM), IEEE, Boston, 2020, pp. 403–408.
  * [3] D. Richiedei, A. Trevisani, Analytical computation of the energy-efficient optimal planning in rest-to-rest motion of constant inertia systems, MECHATRONICS 39 (2016) 147–159.  doi:10.1016/j.mechatronics.2016.05.004.
  * [4] D. Ceulemans, N. van Oosterwyck, J. Demetsenaere, J. de Viaene, J. Steckel, S. Derammelaere, Time-optimal stepper motor motion profile through a novel load-angle-based step-command optimization, in: 2021 IEEE International Electric Machines & Drives Conference (IEMDC), IEEE, Hartford, 2021, pp. 1–7.  doi:10.1109/iemdc47953.2021.9449522.
  * [5] S. Derammelaere, M. Haemers, J. De Viaene, F. Verbelen, K. Stockman, A quantitative comparison between BLDC, PMSM, brushed DC and stepping motor technologies, 19th International Conference on Electrical Machines and Systems (2017).
  * [6] D. Lee, C. W. Ha, Optimization process for polynomial motion profiles to achieve fast movement with low vibration, IEEE Transactions on Control Systems Technology 28 (5) (2020) 1892–1901.  doi:10.1109/TCST.2020.2998094.
  * [7] G. Berselli, F. Balugani, M. Pellicciari, M. Gadaleta, Energy-optimal motions for Servo-Systems: A comparison of spline interpolants and performance indexes using a CAD-based approach, Robotics and Computer-Integrated Manufacturing 40 (2016) 55–65.  doi:10.1016/j.rcim.2016.01.003.
  * [8] J. S. Park, Motion profile planning of repetitive point-to-point control for maximum energy conversion efficiency under acceleration conditions, Mechatronics 6 (6) (1996) 649–663.  doi:10.1016/0957-4158(96)00012-8.
  * [9] G. Carabin, R. Vidoni, Energy-saving optimization method for point-to-point trajectories planned via standard primitives in 1-DoF mechatronic systems, International Journal of Advanced Manufacturing Technology (2021).  doi:10.1007/s00170-021-07277-y.
  * [10] C. Botan, M. Ratoi, F. Ostafi, V. Horga, Minimum energy control of servo drive systems with PMSM, in: SPEEDAM 2010 - International Symposium on Power Electronics, Electrical Drives, Automation and Motion, 2010, pp. 19–23.  doi:10.1109/SPEEDAM.2010.5542403.
  * [11] M. Pellicciari, G. Berselli, F. Balugani, On Designing Optimal Trajectories for Servo-Actuated Mechanisms: Detailed Virtual Prototyping and Experimental Evaluation, IEEE/ASME Transactions on Mechatronics 20 (5) (2015) 2039–2052.  doi:10.1109/TMECH.2014.2361759.
  * [12] Y. L. Hsu, M. S. Huang, R. F. Fung, Energy-saving trajectory planning for a toggle mechanism driven by a PMSM, Mechatronics 24 (1) (2014) 23–31.  doi:10.1016/j.mechatronics.2013.11.004.
  * [13] Y. L. Hsu, Adaptive Tracking Control of a PMSM-Toggle System with a Clamping Effect, International Journal of Mechanical Engineering and Applications 4 (1) (2016) 1.  doi:10.11648/j.ijmea.20160401.11.
  * [14] M. S. Huang, Y. L. Hsu, R. F. Fung, Minimum-energy point-to-point trajectory planning for a motor-toggle servomechanism, IEEE/ASME Transactions on Mechatronics 17 (2) (2012) 337–344.  doi:10.1109/TMECH.2010.2103366.
  * [15] J. L. Ha, R. F. Fung, K. Y. Chen, S. C. Hsien, Dynamic modeling and identification of a slider-crank mechanism, Journal of Sound and Vibration 289 (4-5) (2006) 1019–1044.  doi:10.1016/j.jsv.2005.03.011.
  * [16] K. Sollmann, M. Jouaneh, D. Lavender, Dynamic Modeling of a Two-Axis Parallel H-Frame-Type XY Positioning System, IEEE/ASME Transactions on Mechatronics 15 (2) (2010) 280–290.  doi:10.1109/TMECH.2009.2020823.
  * [17] F. Vanbecelaere, S. Derammelaere, N. Nevaranta, J. De, Online tracking of varying inertia using a SDFT approach, Mechatronics (2020).
  * [18] H. Walsch, I. Kolmanovsky, M. Steinbuch, L. del Re, Optimization and Optimal Control in Automotive Systems, Vol. 455 of Lecture Notes in Control and Information Sciences, Springer International Publishing, Cham, 2014.  doi:10.1007/978-3-319-05371-4.
  * [19] N. Oosterwyck, F. Vanbecelaere, M. Haemers, D. Ceulemans, K. Stockman, S. Derammelaere, CAD Enabled Trajectory optimization and Accurate Motion Control for Repetitive Tasks., in: IEEE International Conference on Control and Automation, ICCA, Vol. 2019-July, IEEE, Edinburgh, 2019, pp. 387–392.  doi:10.1109/ICCA.2019.8899728.
  * [20] Z. Shiller, Time-Energy Optimal Control of Articulated Systems With Geometric Path Constraints, Journal of Dynamic Systems, Measurement, and Control 118 (1) (1996) 139–143.  doi:10.1115/1.2801134.
  * [21] T. Chettibi, H. Lehtihet, M. Hadda, S. Hanchi, Minimum cost trajectory planning for industrial robots, European Journal of Mechanics - A/Solids 23 (4) (2004) 703–715.  doi:10.1016/j.euromechsol.2004.02.006.
  * [22] J. Huang, P. Hu, K. Wu, M. Zeng, Optimal time-jerk trajectory planning for industrial robots, Mechanism and Machine Theory 121 (2017) 530–544.  doi:10.1016/j.mechmachtheory.2017.11.006.
  * [23] G. Wenzhong, S. K. Porandla, Design optimization of a parallel hybrid electric powertrain, 2005 IEEE Vehicle Power and Propulsion Conference, VPPC 2005 (2005) 530–535.  doi:10.1109/VPPC.2005.1554609.
  * [24] J. T. Betts, A Survey of Numerical Methods for Trajectory Optimization, Journal of Guidance, Control, and Dynamics 21 (2) (1998) 193–207.  doi:10.2514/2.4231.
  * [25] A. Piazzi, A. Visioli, Global minimum-time trajectory planning of mechanical manipulators using interval analysis, International Journal of Control 71 (4) (1998) 631–652.  doi:10.1080/002071798221713.
  * [26] E. Kiel, Drive solutions, Springer, 2008.
  * [27] H. Dresig, F. Holzweißig, Dynamics of Machinery: Theory and Applications, Springer Berlin Heidelberg, 2010.
  * [28] G. Ellis, Nonlinear Behavior and Time Variation, in: Control System Design Guide, Butterworth-Heinemann, 2012, pp. 235–260.  doi:10.1016/b978-0-12-385920-4.00012-6.
  * [29] L. C. Westphal, Handbook of Control Systems Engineering, Springer US, 2001.  doi:10.1007/978-1-4615-1533-3.
  * [30] G. Rizzoni, J. Kearns, Principles and applications of electrical engineering, McGraw-Hill (London), 2003.
  * [31] M. Gadaleta, M. Pellicciari, G. Berselli, Optimization of the energy consumption of industrial robots for automatic code generation, Robotics and Computer-Integrated Manufacturing 57 (December 2018) (2019) 452–464.  doi:10.1016/j.rcim.2018.12.020.
URL https://doi.org/10.1016/j.rcim.2018.12.020

  * [32] W. J. Thompson, Chebyshev Polynomials: After the Spelling the Rest is Easy, Computers in Physics 8 (2) (2013) 161.  doi:10.1063/1.4823278.
  * [33] J. Nocedal, S. Wright, Numerical optimization, 2nd Edition, Springer Science & Business Media, 2006.
  * [34] J. H. Holland, Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control and Artificial Intelligence, MIT Press, Cambridge, MA, USA, 1992.
  * [35] H. Majidian, On the decay rate of Chebyshev coefficients, Applied Numerical Mathematics 113 (2017) 44–53.  doi:10.1016/j.apnum.2016.11.004.