# MAP-Elites based Hyper-Heuristic for the Resource Constrained Project
Scheduling Problem

Shelvin Chanda,∗, Kousik Rajeshb,∗, Rohitash Chandrac Commonwealth Scientific
and Industrial Research Organisation, Brisbane, Australia  Indian Institute of
Technology, Guwahati, India  University of New South Wales, Sydney, Australia

###### Abstract

The resource constrained project scheduling problem (RCPSP) is an NP-Hard
combinatorial optimization problem. The objective of RCPSP is to schedule a
set of activities without violating any activity precedence or resource
constraints. In recent years researchers have moved away from complex solution
methodologies, such as meta heuristics and exact mathematical approaches,
towards more simple intuitive solutions like priority rules. This often
involves using a genetic programming based hyper-heuristic (GPHH) to discover
new priority rules which can be applied to new unseen cases. A common problem
affecting GPHH is diversity in evolution which often leads to poor quality
output. In this paper, we present a MAP-Elites based hyper-heuristic (MEHH)
for the automated discovery of efficient priority rules for RCPSP. MAP-Elites
uses a quality diversity based approach which explicitly maintains an archive
of diverse solutions characterised along multiple feature dimensions. In order
to demonstrate the benefits of our proposed hyper-heuristic, we compare the
overall performance against a traditional GPHH and priority rules proposed by
human experts. Our results indicate strong improvements in both diversity and
performance. In particular we see major improvements for larger instances
which have been under-studied in the existing literature.

###### keywords:

Resource constrained project scheduling , MAP-Elites , Evolutionary algorithms
, Quality-Diversity , Hyper-Heuristics

##  1 Introduction

††* These authors contributed equally to the paper

The resource constrained project scheduling problem (RCPSP) is an NP-Hard
combinatorial optimization problem where the goal is to schedule a set of
activities without violating any of the activity precedence or resource
constraints [1]. Computational tools for solving RCPSP are well sought after
due to its wide applicability in areas such as manufacturing [2, 3],
construction and software development [4].

Recently there has been a surge in the use of genetic programming based hyper
heuristics (GPHH) [5, 6] for discovering efficient priority rules for RCPSP.
Compared to traditional metaheuristics [7, 8] and exact methods, these
priority rules [9], provide a much more simple, intuitive and computationally
inexpensive scheduling approach that is well suited for dynamic environments.
GPHH takes a set of problem characteristics/features as input and uses the
process of evolution and recombination to construct useful mathematical
expressions representing priority rules [5, 10]. These rules can then be used
to decide the order in which activities are processed or scheduled. Though
metaheuristics [11, 12, 13, 14] may generate better-performing schedules, the
computational costs often make them infeasible for quick generation of
schedules. The best performing metaheuristic [12] uses genetic algorithms with
neighborhood search. Each problem instance requires more than 1000 iterations
to arrive at a good solution. In contrast, hyper-heuristics provide a one-shot
solution. Once the operator tree of GPHH is evolved, it can be used to quickly
generate a well-performing schedule for any instance size. Thus, hyper-
heuristics provide a more scalable and generalizable approach to RCPSP.

Most GPHHs within literature follow a standard approach where they are
initialized with a set of randomly generated solutions and guided using an
elitist approach. The best performing solutions are carried forward across
generations and used as ’parents’ for generating new, hopefully better,
solutions. However, these hyper-heuristics often lose diversity within the
population quite early and suffer from pre-mature convergence and over-fitting
problems [5, 15]. Essentially, the elitist approach may often result in
insufficient exploration of the search space. When this is combined with the
bloat effect [16], a common problem affecting GPs, it often leads to poor
quality output in terms of the evolved rules.

One way to fix this problem is by using a quality-diversity approach such as
the Multidimensional Archive of Phenotypic Elites (MAP-Elites) algorithm [17].
In comparison to traditional elitist methods, MAP-Elites maintains a map or
grid of high performing solutions which are characterised along multiple
features of interest as defined by the decision maker [17]. While traditional
elitist methods output a single best solution, MAP-Elites outputs a set of
diverse, high-performing solutions [17]. This allows the decision maker to not
only have a wider range of choice but also to gain a better understanding of
the search space. But perhaps most importantly, in defining a feature map, the
decision maker is able to guide the search using domain specific knowledge. In
fact, prior research has shown that since MAP-Elites explores a larger part of
search space, it also tends to find better solutions than traditional elitist
methods [17].

In this paper, we study the automated evolution of priority rules for RCPSP
using a MAP-Elites based hyper-heuristic. In order to demonstrate the benefits
of our proposed hyper-heuristic, we compare the overall performance against a
traditional GPHH. We also compare the discovered rules against the state of
the art rules proposed by human experts. To the best of our knowledge, this is
first study looking into the use of a MAP-Elites based hyper-heuristic in
evolving rules for RCPSP. We also focus on larger instance sizes which have
been rarely studied in existing literature. More specifically, the key
intended contributions of this paper are:

  * 1.

Development and study of a quality-diversity hyper-heuristic framework to
evolve priority rules for the classical RCPSP.

  * 2.

In-depth analysis on quality-diversity elements such as MAP-Elites feature
set, grid size, etc.

  * 3.

Performance characterization of the evolved rules and benchmarking with the
existing rules in the literature on a large set of instances with as many as
300 activities.

The remainder of the paper is organized as follows. Section 2 provides a
discussion on existing literature. Section 3 gives details on the proposed
hyper-heuristic approach for evolving rules. Section 4 provides a discussion
on the results obtained while Section 5 summarizes the findings of this paper
and highlights some future research directions.

##  2 Background

###  2.1 Problem Definition

The RCPSP [18] involves a project with $M$ activities which need to be
scheduled while considering two types of constraints:

  * 1.

_Precedence constraints_ : These represent the interdependence between the
activities. If activity $j$ is a successor of activity $i$ then activity $i$
must be completed before activity $j$ can be started.

  * 2.

_Resource constraints_ : These represent limitations on resources such as
manpower, budget, etc. Each project is assigned a set of $K$ renewable
resources where each resource $k$ is available in $R_{k}$ units for the entire
duration of the project. Each activity may require one or more of these
resources to be completed. While scheduling the activities, the daily resource
usage for resource $k$ can not exceed $R_{k}$ units.

Each activity $j$ takes $d_{j}$ time units to complete. The overall goal of
the problem is to minimize the _makespan_ , i.e., the total duration from the
start of first to the end of last scheduled activity.

An example activity on node (AON) diagram for an RCPSP instance is given in
Fig 1. A possible schedule for this project that doesn’t violate the
precedence and resource constraints is shown in Fig 2.

Figure 1: An RCPSP instance visualised as a Directed Acyclic Graph (DAG) with
each node representing an activity. There is a single resource available (K =
1) and a node is represented as $i(d_{i},r_{i})$ where $d_{i}$ and $r_{i}$
represent the duration and resource requirement for activity $i$. The arrows
represent the precedence constraints between activities.

The RCPSP is considered to be NP-hard [19] and therefore the exact methods
[20] are only practical for smaller instances. For larger instances,
heuristic/meta-heuristic methods such as genetic algorithm [21], particle
swarm optimization [22], simulated annealing [23] and hybrid search [24], etc.
are utilized for finding ’good-enough’ or ’near-optimal’ solutions.

Figure 2: A possible schedule for the RCPSP instance shown in Fig 1.

###  2.2 Quality-Diversity based Hyper-Heuristics

The literature on quality-diversity based hyper-heuristics is rather scant.
There are a number of studies such as [25, 26, 27] which consider diversity in
one way or another to improve evolution. However, there are very few studies
that actually utilize a quality-diversity like framework. Nguyen et al. [28]
proposed the adaptive charting genetic programming (ACGP) hyper-heuristic (for
the dynamic flexible job shop scheduling problem) which incrementally maps the
explored areas of the search space and adaptively distributes its search
effort to efficiently explore the search space. Furthermore, growing neural
gas (GNG) and principle component analysis (PCA) utilize the phenotypic
characteristics of the evolved heuristics to efficiently generate and update
the map. Finally, a surrogate assisted model utilizes this map to determine
which heuristics are to be explored in the next generation.

While ACGP presented promising results, it also has its fair share of
challenges. Firstly, GNG and PCA introduce a number of extra parameters which
can be difficult to tune. Furthermore, the mapping strategy seems less
intuitive from the point of view of a decision maker since it is less
structured as opposed to the way MAP-Elites constructs its archive. This may
also lead to certain areas of the search space being left under-explored or
totally ignored. A variant of this algorithm has also been applied to the
resource constrained job scheduling problem [29]. The same authors [30] later
presented an extension of this algorithm with improved diversity control,
bloat control and human-driven interactive evolution.

###  2.3 GPHH for RCPSP

Priority rules are a much more preferred option, over meta-heuristics and
exact methods, when it comes to scheduling in dynamic environments. Meta-
heuristics and exact methods are usually time consuming to run and not to
mention, that tuning the parameters for such methods is often non-trivial. But
perhaps most importantly, these methods lack flexibility and generalization
capability needed for real world problems [5]. Priority rules are simple,
intuitive and easy to use. They also allow for quick decision making which is
especially important in dynamic environments where project elements are
constantly changing and evolving. Priority rules are also most often void of
any complex parameters that need tuning. Given their obvious advantages, the
use of priority rules for scheduling complex projects is on the rise. However,
designing such rules is often non-trivial. The literature contains a number of
rules which have been proposed over the years, by domain experts. Some of
these are outlined are in Table 1. A more comprehensive list of rules can be
found in [31].

Human experts are limited in their ability and scope when it comes to
proposing new and improved priority rules. As a result there has been a recent
push towards automating this process. Genetic programming based hyper-
heuristics (GPHH) are the most preferred option as they have had success in
evolving mathematical expressions [32].

To the best of our knowledge, the first study on evolving priority rules for
RCPSP was by Frankola et al. [33]. The authors proposed a GPHH which evolved
arithmetic priority rules constructed using binary operators and a limited set
of activity attributes which mostly focused on resource usage and network
structure. This would be seen more so an exploratory study since the authors
did not present any in-depth analysis of the results or any detailed
comparison with existing human designed priority rules. This was followed by a
more in-depth study by Šišejković in [34] in which he explored the use of a
much wider attribute and operator set. The author first used a greedy feature
selection approach to identify the best performing attributes. Later, these
attributes were combined to form arithmetic priority rules. The resulting
rules exhibited competitive performance but were unable to outperform the best
priority rules from literature.

Chand et al. [5] proposed a GPHH framework in which they explored the use of
different rule representations and a more structured selection of attributes.
The authors compared the classical arithmetic representation with a decision
tree like representation similar to the one proposed in [10]. They also
proposed an attribute set that focused on three key aspects of any project,
namely the resource requirements, the critical path and the precedence
network. The resulting rules outperformed all other existing human designed
rules including LST, WCS, ACS, etc. In terms of the study on representation
types, the authors found that the arithmetic representation performed the best
while the decision tree representation gave competitive results with greater
understandability.

The same authors [27] later proposed a multi-objective GPHH which tried to
optimize performance as well as rule complexity. All of the studies prior to
this had only focused on static instances. This study for the first time
focused on dynamic RCPSP in which resource availability varied with time. The
results indicated strong improvements over both, existing human designed rules
and also the results presented in [5].

To the best of our knowledge there have been no studies so far that deviate
from utilizing a traditional elitist approach to a more quality-diversity
based approach. Some of the studies mentioned above [5, 27] report cases of
pre-mature convergence and over-fitting which can be attributed to some extent
to the elitist nature of traditional GPs. Also, nearly all these studies only
focus on smaller instances (30 to 120 activities). However, the true benefit
of using a quality-diversity based approach over an elitist approach becomes
much more apparent when we work on large instances. In this study we aim to
explore the benefits of a quality-diversity based hyper heuristic which
utilizes MAP-Elites to search for high performing priority rules across a
structured feature map. We also specifically focus on evolving rules for
larger instances (300 activities) which are much more closer to the type of
instances encountered in real world project scheduling.

Figure 3: Solutions within a MAP-Elites grid. Table 1: Some of the prominent
existing priority rules.

Rule | Abbrv. | Extr. | Description
---|---|---|---
Early Start Time [35] | EST | Min | $ES_{j}$
Early Finish Time [35] | EFT | Min | $EF_{j}$
Late Start Time [35] | LST | Min | $LS_{j}$
Late Finish Time [35] | LFT | Min | $LF_{j}$
Most Total Successors [2] | MTS | Max | $|\bar{S_{j}}|$
First In First Out [35] | FIFO | Min | $ID_{j}$
Shortest Processing Time [35] | SPT | Min | $d_{j}$
Greatest Rank Position Weight[2] | GRPW | Max | $d_{j}+\sum_{i\in S_{j}}{d_{i}}$
Greatest Resource Demand [35] | GRD | Max | $d_{j}*\sum_{k=1}^{K}{r_{kj}}$

  * 1

$j\in\\{1...M\\}$

  * 2

$\bar{S_{j}}$ is the complete set of successor activities for activity $j$.

Table 2: Operator Set Add(a,b) | Mul(a,b) | Sub(a,b) | Div(a,b) | Max(a,b) | Min(a,b) | Neg1(a)
---|---|---|---|---|---|---
$a+b$ | $a*b$ | $a-b$ | $\begin{cases}a/b\quad if\ b>0\\\ 0\quad otherwise\end{cases}$ | $\begin{cases}a\quad if\ a>b\\\ b\quad otherwise\end{cases}$ | $\begin{cases}a\quad if\ a<b\\\ b\quad otherwise\end{cases}$ | $-1*a$

##  3 MAP-Elites for RCPSP

Our proposed framework is composed of a number of components which include the
MAP-Elites algorithm for evolving the rules, a set of operators (Table 2) and
attributes (Table 3) for constructing the rules and a schedule generation
scheme for generating complete schedules from the evolved rules. These are
individually outlined below.

###  3.1 MAP-Elites Hyper-Heuristic (MEHH)

We build our hyper-heuristic using the original version of MAP-Elites
(Algorithm 1) proposed by Mouret and Clune [17]. The features used in
constructing the feature map/grid are detailed in the next sub-section. The
algorithm starts by generating a set of random trees. These trees are
evaluated on the training set to determine their fitness. Next, their feature
values are calculated and these solutions are placed into the feature grid. If
two or more solutions fall within the same feature discretisation, then only
the solution with the lowest fitness value is retained.

Subsequently, crossover and mutation operators are applied to solutions within
the grid in order to generate new solutions. These solutions are similarly
evaluated and the grid is updated. This process is repeated until the
termination condition is reached. All solutions within the final grid are
evaluated on the validation set and the best performing solution is selected
to be evaluated on the test set.

The fitness of an individual solution is the average percentage deviation of
the calculated makespan ($ms_{i}$) from the instance lower bound ($lb_{i}$).
The lower bound used here is calculated using the critical path method, which
generates a precedence feasible schedule while ignoring the resource
constraints. The crossover and mutation operators are the same as the ones
used in our previous work [5, 36, 27].

Algorithm 1 Basic MAP-Elites Algorithm

1:for $gen=0$ to $g$ do

2: for $i=1$ to $\lambda$ do

3: if $gen=0$ then

4: $y\leftarrow random\\_individual()$

5: else

6: $y\leftarrow random\\_variation(P)$

7: $f\leftarrow calculate\\_features(y)$

8: $p\leftarrow calculate\\_performance(y)$

9: if $P(f)=\Phi$ or $F(P(f))>p$ then

10: $P(f)\leftarrow y$

###  3.2 MEHH Features

A given combination of feature values correspond to one unique individual in
the grid if the cell pointed to by the feature values is non empty. This is
illustrated in Figure 3. MAP-Elites was originally designed to operate on
phenotypic features, as the name suggests, however in this paper we have used
a combination of genotypic and phenotypic features. Our feature map/grid is
constructed using the following features:

  1. 1.

Number of nodes in an individual: This is the total number of nodes in the
tree representing the solution which include both internal (operations) and
leaf (activity attribute) nodes.

  2. 2.

Number of Resource related nodes: These are the leaf nodes which carry
information about resource usage of an activity; i.e
$RR_{j},AvgRReq_{j},MaxRReq_{j},MinRReq_{j}$.

  3. 3.

Slack: The slack provides a measure of how loosely packed the schedule is. The
slack value for a given activity within a generated schedule is taken to be
the maximum amount of time the activity can be delayed without causing a
violation of resource or precedence constraints. The total for a given
schedule is the sum of the individual activity slacks.

Note that the first two features are properties of the generated tree itself,
and independent of the performance on the instances in the training set. The
third feature is dependent on the training set and the resulting schedules
generated by the tree. In our experiments we limit the tree height. The number
of nodes in a tree and the number of resource related nodes are bounded due to
the fixed maximum height of generated trees; hence, these features are not
normalised. The slack value depends on the number of instances and the size of
the instances; hence, there is a need to normalise it. The total slack of an
individual $\zeta$ on the training set T is given by:

$\text{Total slack}=\frac{1}{|T|}\sum_{i\in T}\frac{slack_{i}}{n_{i}}$

where $T$ denotes the training set, $slack_{i}$ is the slack value for the
schedule generated on instance $i$ by the individual $\zeta$ and $n_{i}$ is
the number of activities in instance $i$.

The resulting slack was found experimentally which fell in the range [1.65,
2.00]; hence, the feature domain is chosen appropriately. Slack calculation is
outlined in Algorithm 2 and the variables involved are explained below:

  * 1.

slack stores the sum of slack over all activities

  * 2.

n is the total number of activities in the instance.

  * 3.

leastTime stores the least start time for the $i^{th}$ job over all its
successors. The $i^{th}$ job has to end before the leastTime.

  * 4.

startTimes is an array where the $i^{th}$ element represents the start time of
the $i^{th}$ job.

  * 5.

finishTimes is an array where the $i^{th}$ element represents the finish time
of the $i^{th}$ job.

  * 6.

jobRes is an array where the $i^{th}$ element represents the resource required
to complete the $i^{th}$ job.

  * 7.

totalRes is the amount of total resource available to the instance.

  * 8.

resCons represent an array for the amount of resource being consumed at time
$j$.

Algorithm 2 Calculating slack for a schedule

1:Set: $slack=0$

2:for $i=1$ to $n-1$ do

3: Set : leastTime= $\infty$

4: for $j$ in $successors[i]$ do

5: $leastTime=min(leastTime,startTimes[j])$

6: for $j=finishTimes[i]+1$ to $leastTime-1$ do

7: schedulable $\leftarrow$ True

8: for $k=1$ to $K$ do

9: if $\textit{jobRes}[i][k]+\textit{resCons}[j][k]>\textit{totalRes}[k]$ then

10: schedulable $\leftarrow$ False

11: break

12: if schedulable then

13: $slack\leftarrow slack+1$

14: else

15: break

###  3.3 Priority Rule Representation

In this study, we utilize the widely used arithmetic representation [10, 5]
for candidate solutions within our hyper-heuristic. The arithmetic
representation forms priority rules by combining arithmetic operators with
activity attributes. The main advantage of arithmetic representation is the
flexibility it provides in an evolutionary environment, allowing for easy
recombination without resulting in infeasible tree structures [36], as is the
case with other forms of representations (eg. decision tree representation
[5]). We utilize six binary operators and one unary operator which are listed
in Table 2.

###  3.4 Activity Attributes

As mentioned previously, a candidate solution within our hyper-heuristic is a
mathematical expression representing some combination of mathematical
operators and activity attributes. We consider 10 different activity
attributes in this research. The activity attributes that we consider are
listed below, and their calculations are detailed in Table 3:

  * 1.

TSC/TPC: The total count of the (immediate and non-immediate)
successors/predecessors of an activity.

  * 2.

ES/EF: The earliest start/finish time for an activity in the precedence
feasible schedule calculated by relaxing the resource constraints where each
activity is scheduled as early as possible.

  * 3.

LS/LF: The latest start/finish time for an activity in the precedence feasible
schedule calculated by relaxing the resource constraints where each activity
is scheduled as late as possible using backward recursion [37].

  * 4.

RR: The total number of resources required by an activity.

  * 5.

AvgRReq: The average resource requirement of an activity.

  * 6.

MinRReq: The minimum resource requirement of an activity.

  * 7.

MaxRReq: The maximum resource requirement of an activity.

Table 3: Activity Attribute Set

Attribute | Details
---|---
Early Start ($ES_{j}$) | $\frac{1}{max(ES)}ES_{j}$
Early Finish ($EF_{j}$) | $\frac{1}{max(EF)}EF_{j}$
Late Start ($LS_{j}$) | $\frac{1}{max(LS)}LS_{j}$
Late Finish ($LF_{j}$) | $\frac{1}{max(LF)}LF_{j}$
Total Predecessor Count ($TPC_{j}$) | $\frac{1}{M-1}|\bar{P_{j}}|$
Total Successor Count ($TSC_{j}$) | $\frac{1}{M-1}|\bar{S_{j}}|$
Resources Required ($RR_{j}$) | $\frac{1}{K}\sum_{k=1}^{K}\begin{cases}1\quad if\ r_{jk}>0\\\ 0\quad otherwise\end{cases}$
Avg. Resource Requirement ($AvgRReq_{j}$) | $\frac{1}{K}\sum_{k=1}^{K}\\{\frac{1}{R_{k}}r_{jk}\\}$
Max Resource Requirement ($MaxRReq_{j}$) | $Min(\frac{1}{R_{k}}r_{jk},\ k\in 1...K)$
Min Resource Requirement ($MinRReq_{j}$) | $Max(\frac{1}{R_{k}}r_{jk},\ k\in 1...K)$

  * 1

$\bar{P_{j}}$ is the complete set of predecessor activities for activity $j$.

These attributes have been used in our previous works [5, 36] and have shown
to be useful in constructing good priority rules. All attributes are also
scaled between 0 and 1 to avoid any undesirable outcomes resulting from
attributes with particularly high or low values.

###  3.5 Schedule Generation Scheme

There are two different schedule generation schemes detailed in [38] which are
widely used within the RCPSP literature. These are the parallel schedule
generation scheme (PSGS) and the serial schedule generation scheme (SSGS)
[18]. It has been shown in previous research [2, 5] that PSGS performs better
when scheduling with priority rules. PSGS is also based on time increment
which makes it suitable for dynamic online scheduling [5]. Hence, for these
reasons, we use PSGS within our hyper-heuristic framework. We use the
implementation detailed in [18].

When generating a schedule with a priority rule, we always pick the activity
with minimum value of the priority function to be scheduled next. If two
activities have the same priority value, the tie is broken using their
activity ID, that is, the one with the lower ID value is picked first.

##  4 Experiment and Results

###  4.1 Experiment setup

Our key algorithmic comparison is between a standard GPHH as detailed in [5]
and our proposed MEHH with varying archive sizes. Table 4 gives a summary of
the algorithm parameters. These are standard values which have been used in
our previous research [5]. The feature domains have been determined using
preliminary runs. Both methods were implemented using the DEAP
111https://deap.readthedocs.io/en/master/api/algo.html framework. All data and
code from our research is available on github 222https://github.com/sydney-
machine-learning/MEHH_RCPSP.

For all our experiments we use a training, validation and test set. For the
training set we use the entire J30 instance set from the Project Scheduling
Library (PSPLib) [39]. J30 set consists of 480 instances with 30 activities
each. The instances are constructed using three complexity parameters, namely
network complexity (NC), resource factor (RF) and resources strength (RS). The
instances were constructed using 48 different parameter combinations and each
combination consists of 10 instances. For the validation and test, we use the
RG300 instances. The RG300 [8] set was constructed using three complexity
parameters, namely order strength (OS), resource utilization (RU), resource
constrainedness (RC). The instances were constructed using 48 different
parameter combinations and each combination consists of 10 instances. For
validation, we pick the first instance from each parameter combination,
resulting in 48(10%) instances in total. The remaining 90% of the instances
are used for testing. Having an equal number of instances from each parameter
combination ensures that instances of different hardness levels are included
in training and validation. Using a validation set helps us reduce the effect
of over-fitting and also lets us take advantage of the diversity in the
archive. Training, validation and test performance is based on average
percentage deviation from the lower bounds, with a lower value considered to
be better.

Table 4: Parameters for GPHH and MEHH. Parameter | Value
---|---
Population size | 1024
Stopping condition | 25 Generations
Fitness function | Average % Deviation
Selection method | Tournament Selection
Selection pool size | 7
Gen height | [2,5]
Height limit | 7
Crossover probability | 0.8
Mutation probability | 0.2
Number of runs | 31
Max items per bin (MEHH) | 1
Fitness domain (MEHH) | [0,1]
$Feature_{1}$ (# of nodes) domain (MEHH) | [4,127]
$Feature_{2}$ (# of res. nodes) domain (MEHH) | [0, 30]
$Feature_{3}$ (Slack) domain (MEHH) | [1.65, 2.00]

###  4.2 Comparison Against GP

We start by comparing our proposed MEHH with a standard GPHH [5]. For MEHH we
consider 4 different variants in which everything is kept constant except for
archive size. These values are chosen based on desired feature discretization.
We consider 4 different feature discretization values which include 5, 10, 15
and 20. Since we have three different features, this results in a total
archive size of 125 (5*5*5), 1000, 3375 and 8000. From here on we denote all
our MEHH variants as $MEHH_{<archive\\_size>}$.

Both methods are evaluated using 31 independent runs. For each run we evaluate
the final archive or population (in the case of GP) on the validation set and
pick the best performing solution to be the representative for that run. This
solution is then evaluated on the completely unseen test set. The final
statistics given in Table 5 are calculated based on the test performance of
the 31 representative rules from each of the runs.

The results clearly indicate the benefit of using MEHH over GPHH in finding
high performing priority rules for larger instances. For all archive sizes,
our proposed MEHH presents better mean and median performance. $MEHH_{8000}$
performs the best with a 2.53% improvement over GPHH. It is also more reliable
in comparison to a GPHH due to low standard deviation in the results across
multiple runs. In fact, all our MEHH variants have lower standard deviation in
comparison to GPHH. We also see a correlation between increase in archive size
and improvement in overall performance. This would make sense, since a larger
archive would allow MEHH to maintain a larger set of useful genetic material
that could greatly assist the evolutionary search.

GPHH loses diversity over generations. Even if the solutions in GPHH are
different in terms of the tree structure, they tend to be variants of one
another resulting in the same performance. Essentially the bloat effect takes
over resulting in longer trees with little to no improvement in terms of
objective value. This is further illustrated in Figure 4 which shows the
percentage of unique solutions in the archives or the population for the
respective median runs (based on test performance). As the run progresses,
GPHH continues to lose diversity which greatly hinders its ability to find
better solutions. MEHH on the hand continues to improve in terms of diversity
as more and more cells in the archive get occupied. Table 6 sheds more light
on this issue. It tabulates the number of unique individuals across all runs.
In our experiments GPHH was run with a population size of 1024 yet on average
it only has about 37 unique individuals at the end of a run. Even $MEHH_{125}$
manages to outperform GPHH in this aspect with more than twice as many unique
individuals in the end. This again highlights the advantage and need for
diversity preserving hyper-heuristics as early convergence, overfitting and
bloat effect are all real issues affecting the quality of solutions in a
standard GPHH.

Table 5: Comparison of test performance (average percentage deviation ) for GPHH and MEHH over all 31 runs.  Measure | GPHH | $MEHH_{125}$ | $MEHH_{1000}$ | $MEHH_{3375}$ | $MEHH_{8000}$
---|---|---|---|---|---
Mean | 1005.93 | 1004.00 | 1003.61 | 1003.41 | 1003.40
Median | 1004.92 | 1003.42 | 1003.39 | 1003.22 | 1003.41
Best | 1002.93 | 1002.36 | 1001.72 | 1001.62 | 1001.47
Worst | 1010.15 | 1009.25 | 1006.34 | 1007.34 | 1004.64
St. Dev. | 2.01 | 1.50 | 0.87 | 1.25 | 0.80
Figure 4: Change in diversity over generations for median runs of GPHH and MEHH. Table 6: Unique individuals at termination of training. | GPHH | $MEHH_{125}$ | $MEHH_{1000}$ | $MEHH_{3375}$ | $MEHH_{8000}$
---|---|---|---|---|---
Mean | 37 | 83 | 469 | 1124 | 1905
Median | 23 | 83 | 471 | 1124 | 1891
Min | 2 | 77 | 436 | 1045 | 1817
Max | 133 | 91 | 509 | 1190 | 2069
St. Dev. | 36 | 3 | 17 | 43 | 59

###  4.3 MEHH Archive Coverage

An important aspect to consider is the coverage of the MEHH grid. Coverage is
calculated as the percentage of cells in the grid which contain a solution.
Since it isn’t necessary that a solution exists for each feature value, some
cells will inevitably be empty in the grid.

Table 7: Coverage at termination for median run of MEHH for different grid sizes. Grid size | Coverage (%)
---|---
$MEHH_{125}$ | 66.4
$MEHH_{1000}$ | 47.1
$MEHH_{3375}$ | 33.3
$MEHH_{8000}$ | 23.6

A greater coverage is beneficial since it provides the decision maker with a
greater number of choices. In our case, it allows us to have a better chance
of avoiding over-fitted solutions. Table 7 gives a summary of archive coverage
by size for a typical run. We observe that archive coverage actually decreases
with an increase in archive size. This may be due to a number of reasons.

First contributing factor is the correlation between two of the selected
features. Number of nodes and number of resource nodes are positively
correlated (Figure 5). In fact one is a subset of the other and hence this
means that certain feature combinations are realistically not possible.
Furthermore, the same two features can only have integer values. This can
create gaps in the feature space due to the non-continuous nature and once
again lead to combinations or feature ranges which realistically cannot have
any corresponding solutions. This problem can get worse with an increase in
archive size. However, as noted earlier, we see a positive correlation between
increase in archive size and improvement in results. Hence, even though the
coverage for larger archive sizes is less, the quality is sufficient to ensure
good results.

Figure 5 shows pair-wise feature-fitness plots. We can see that smaller trees
generally have better training performance. This is a trend that has been
reported in our previous work [27] as well. Next we see that solutions with
lower slack values also tend to perform better. This is again expected as less
slack corresponds to more compact schedules with lower makespans.

(a) $Feature_{1}$ vs $Feature_{2}$

(b) $Feature_{2}$ vs $Feature_{3}$

(c) $Feature_{2}$ vs $Feature_{3}$

Figure 5: Feature wise plots of the performance grid. The fitness of
individuals corresponding to the third feature (in each of the plots) is
averaged along the z axis

###  4.4 Comparison with existing rules

Now that we have established that our MEHH outperforms a classical GPHH, we
turn our attention to the quality of individual rules. It is important to
analyze the performance of the evolved rules against existing human designed
priority rules. Table 8 lists the performance of existing priority rules from
literature on the test set. The best performing rule is MTS closely followed
by LFT, FIFO and EST.

Table 8: Performance of human designed priority rules on test set. Rule | Percentage Deviation | Makespan
---|---|---
EST | 1020.19 | 313714
EFT | 1028.43 | 316192
LST | 1020.43 | 313813
LFT | 1014.39 | 311633
SPT | 1077.04 | 330034
FIFO | 1016.30 | 312323
MTS | 1013.72 | 311336
RAND | 1071.38 | 328539
GRPW | 1044.91 | 323228
GRD | 1081.61 | 331584
IRSM | 1079.58 | 330329
WCS | 1023.52 | 314897
ACS | 1050.69 | 322928

As mentioned earlier, for each run we pick the best performing rule on the
validation set to be the representative for that run and to be subsequently
evaluated on the test set. From these 31 rules we select a sub-set of rules
for further analysis and comparison with MTS. We pick 3 rules from GPHH and
each of the MEHH variants. For each algorithmic variant we select the best,
worst and median rules from the set of 31 representative rules selected across
runs. These are of course selected based on test performance. Each of the MEHH
rules are denoted as $MEHH_{<grid\ size>-<metric>}$ where metric can be either
B (Best), W(worst) or M(median) while the grid_size can be either 125, 1000,
3375, 8000. For GPHH, the rules are denoted as $GPHH_{<metric>}$. The
performance of the selected rules is given in Table 9.

Based on the results, we can clearly see that all of our selected rules
outperform the existing the human designed priority rules. Even in the worst
case, our evolved rules are still better than MTS. Rule $MEHH_{8000-B}$
performs best in terms of average percentage deviation while $MEHH_{1000-B}$
performs the best in terms of cumulative makespan. In all cases, the rules
selected from MEHH outperform the GPHH evolved rules. This indicates the
strength of our proposed MEHH in discovering priority rules which perform
exceptionally well on new unseen cases.

Table 9: Performance of selected evolved rules on the test set. Rule | Deviation | Makespan
---|---|---
$GPHH_{B}$ | 1002.93 | 308460
$MEHH_{125-B}$ | 1002.36 | 308510
$MEHH_{1000-B}$ | 1001.72 | 308141
$MEHH_{3375-B}$ | 1001.62 | 308199
$MEHH_{8000-B}$ | 1001.47 | 308199
$GPHH_{M}$ | 1004.92 | 309125
$MEHH_{125-M}$ | 1003.42 | 309072
$MEHH_{1000-M}$ | 1003.39 | 309000
$MEHH_{3375-M}$ | 1003.22 | 308922
$MEHH_{8000-M}$ | 1003.41 | 308958
$GPHH_{W}$ | 1010.15 | 310597
$MEHH_{125-W}$ | 1009.25 | 310290
$MEHH_{1000-W}$ | 1006.34 | 309714
$MEHH_{3375-W}$ | 1007.34 | 310190
$MEHH_{8000-W}$ | 1004.64 | 309636

We continue our analysis by comparing the performance of the same set of rules
on only hard instances from the test set. We characterise hard instances as
those having RC $\geq$ 0.6 and RU $\geq$ 3\. This essentially refers to
instances where the resource requirements are high and the availability is
relatively limited. The results for this comparison are given in Tables 10 and
11. We can see that once again MTS outperforms all other human designed
priority rules. Similarly, $MEHH_{8000-B}$ performs best in terms of average
percentage deviation while $MEHH_{3375-B}$ performs the best in terms of
cumulative makespan. However, the gap in performance between the evolved rules
and MTS is now much larger. This further emphasizes the usefulness of our
evolved rules which perform equally well on new unseen difficult instances.
For a potential user or decision maker, this information is handy as it gives
an indication of how well the rules perform under difficult circumstances.

Table 10: Performance of selected evolved rules on hard instances. Hard instances were classified as having RC $\geq$ 0.6 and RU $\geq$ 3\.  Rule | Deviation | Makespan
---|---|---
$GPHH_{B}$ | 2219.69 | 161659
$MEHH_{125-B}$ | 2218.42 | 161551
$MEHH_{1000-B}$ | 2218.64 | 161576
$MEHH_{3375-B}$ | 2217.96 | 161502
$MEHH_{8000-B}$ | 2217.70 | 161523
$GPHH_{M}$ | 2224.90 | 161998
$MEHH_{125-M}$ | 2219.55 | 161684
$MEHH_{1000-M}$ | 2220.31 | 161656
$MEHH_{3375-M}$ | 2220.16 | 161656
$MEHH_{8000-M}$ | 2220.33 | 161722
$GPHH_{W}$ | 2234.72 | 162631
$MEHH_{125-W}$ | 2230.40 | 162246
$MEHH_{1000-W}$ | 2227.37 | 162196
$MEHH_{3375-W}$ | 2227.69 | 162220
$MEHH_{8000-W}$ | 2224.58 | 161935

Finally, we compare the performance of the same set of rules on the remaining
PSPLib instances. As mentioned earlier, we used the entire J30 instance set
for training. Here we compute the performance of these rules on the J60, J90
and J120 set. J60 and J90 both consist of 480 instances which are constructed
using the same parameter combinations as J30. J120 consists of 600 instances
which are constructed using a slighter harder set of parameter combinations
[39]. The results for this comparison are given in Tables 12 and 13.
$MEHH_{8000-M}$ is the best performing rule for J60, while $MEHH_{1000-B}$ and
$MEHH_{3375-M}$ perform the best on J90 and J120, respectively. Majority of
our evolved rules out-perform all of the human designed rules. The difference
in performance is a lot smaller in this case but still significant if we
consider the number of instances involved.

Table 11: Performance of human designed priority rules on hard instances. Rule | Deviation | Makespan
---|---|---
EST | 2242.00 | 163232
EFT | 2251.00 | 163826
LST | 2249.78 | 163724
LFT | 2236.70 | 162772
SPT | 2297.57 | 166841
FIFO | 2237.51 | 162870
MTS | 2234.50 | 162596
RAND | 2284.17 | 166059
GRPW | 2265.79 | 165272
GRD | 2301.08 | 167284
IRSM | 2297.44 | 166848
ACS | 2255.45 | 164369
WCS | 2250.37 | 163797
Table 12: Performance of selected evolved rules on smaller instance sets.  Rule | J60 | J90 | J120
---|---|---|---
$GPHH_{B}$ | 17.06 | 15.54 | 43.04
$MEHH_{125-B}$ | 16.89 | 15.53 | 42.93
$MEHH_{1000-B}$ | 16.97 | 15.26 | 43.21
$MEHH_{3375-B}$ | 17.30 | 15.56 | 43.00
$MEHH_{8000-B}$ | 17.09 | 15.62 | 43.16
$GPHH_{M}$ | 16.94 | 15.46 | 42.99
$MEHH_{125-M}$ | 17.08 | 15.49 | 43.08
$MEHH_{1000-M}$ | 17.02 | 15.49 | 43.42
$MEHH_{3375-M}$ | 16.87 | 15.42 | 42.83
$MEHH_{8000-M}$ | 16.82 | 15.52 | 42.88
$GPHH_{W}$ | 17.32 | 15.63 | 43.33
$MEHH_{125-W}$ | 17.74 | 15.54 | 43.31
$MEHH_{1000-W}$ | 17.01 | 15.46 | 43.80
$MEHH_{3375-W}$ | 16.88 | 15.50 | 43.39
$MEHH_{8000-W}$ | 17.32 | 15.62 | 42.93
Table 13: Performance of human designed priority rules on smaller instance sets. Rule | J60 | J90 | J120
---|---|---|---
EST | 21.68 | 21.45 | 55.78
EFT | 22.46 | 21.82 | 55.77
LST | 17.12 | 15.80 | 44.04
LFT | 17.46 | 15.90 | 43.86
SPT | 23.77 | 23.60 | 60.33
FIFO | 20.38 | 19.47 | 51.57
MTS | 17.98 | 16.70 | 46.03
RAND | 23.53 | 22.81 | 59.08
GRPW | 22.32 | 21.83 | 57.87
GRD | 24.11 | 23.18 | 61.30
IRSM | 20.73 | 19.14 | 53.62
ACS | 20.41 | 19.09 | 51.98
WCS | 18.48 | 16.92 | 46.43

###  4.5 Discussion

The results presented above highlight the strength of our MEHH. MEHH produces
rules which on average generalize significantly better to an unknown test set.
GPHH starts with a diverse population which quickly erodes away as the
evolutionary selection pressure starts filling up the population with copies
of few ’good’ solutions and their variants. MEHH on the other hand maintains
an archive that has been partitioned based on selected features of interest. A
new solution is only added into the archive if it exhibits features which the
archive has not yet encountered or if it outperforms an existing solution with
the same characteristics. This approach ensures and encourages a diverse
population which in turn provides a wider range of options when picking the
final representative rule.

We observed the same pattern when comparing individually selected rules. Rules
selected from MEHH always outperformed rules selected from GPHH in the best,
median and worst case scenarios. An illustration of this is given in Figure 6.
The gap in performance increases with an increase in instance size. As we
tackle larger more difficult instances, the benefits of MEHH evolved rules
become more apparent. Our evolved rules also outperform existing human
designed rules by a significant margin. Prior to this, majority of the studies
on hyper-heuristic focused solely on smaller instances (30 to 120 activities).
For the first time, we presented an in-depth study with larger and more
difficult instances.

The rules evolved in this study will be highly desirable in dynamical
scheduling environments where sudden and unexpected changes require quick
scheduling responses. In such situations the goal is to get a ’good enough’
solution in a reasonable amount of time as opposed to focusing on optimality.
Furthermore the evolved rules can also be used to seed or augment existing
meta-heuristic [40] and exact methods to improve their performance in cases
where optimality and not speed is the desired goal. Some of the rules
discovered in this study are simplified and presented in Table 14. A further
comparison in terms of performance and rule complexity is given in Figure 7.
Rule $MEHH_{125-B}$ offers the best trade-off between performance and
complexity. In fact, except for $MEHH_{1000-B}$, all other evolved rules
exhibit a reasonable level of complexity and can be easily utilized by
decision makers.

Furthermore, MEHH allows the decision maker to customize the search process to
their liking by choosing features that are more relevant to their scheduling
environment. For instance, in some real world applications incorporating slack
is an important aspect of the scheduling process. This is done in order to
’cushion’ the schedule from unplanned disruptions. The decision maker can thus
refer to the archive and pick a heuristic that appropriately reflects their
desired quality indicators(eg. a decent trade-off between performance and
slack).

Figure 6: Improvement in percentage deviation given by MEHH rules over $GPHH_{M}$. Figure 7: Comparison of evolved and existing rules. Table 14: Selected evolved rules and their performance on the test set.  Rule | Description | Deviation
---|---|---
$GPHH_{B}$ |  10cm-AvgRReq - EF - ES - LF - 2*LS - Max(AvgRReq, TSC) - Min(AvgRReq, -TSC) + Min(EF, LS) | 1002.93
$MEHH_{125-B}$ |  10cm-2*AvgRReq + EF - LF - Max(EF, TSC) + Max(LS, RR/TPC - Max(MaxRReq + MinRReq, MaxRReq + TSC)) | 1002.36
$MEHH_{1000-B}$ |  10cm-AvgRReq + LF*LS - TSC + (EF*MaxRReq - MaxRReq*TPC)*(LS*MaxRReq + MaxRReq/TSC)*(-MaxRReq*RR - 1/(LS*MaxRReq))*(-MaxRReq*Min(AvgRReq, EF) + Min(ES, RR) + Min(Min(1, AvgRReq*TSC)/(EF*TPC + 2*MinRReq), 1/(-AvgRReq - Min(AvgRReq, TPC))) + Min(1/(MinRReq + RR), Max(TPC, TSC)))*Max(ES*LF, Min(LS, RR))*Min(ES*MinRReq, Max(EF, MinRReq))*Min(MinRReq + TSC, Max(MinRReq, RR))*Min(AvgRReq, TSC, Max(EF, LS)) | 1001.72
$MEHH_{3375-B}$ |  +10cmMax(LS, MinRReq) + Min(-TSC*(AvgRReq*TSC + AvgRReq + ES + 1 + 1/(-EF - LS)), ES**2*LS*MinRReq*Max(AvgRReq, EF, AvgRReq*ES) + Min(-MinRReq, MaxRReq*RR*(AvgRReq + 1/MaxRReq), LS - RR) - 1/(EF*ES*(EF + MaxRReq))) | 1001.62
$MEHH_{8000-B}$ |  10cmLF*LS - Max(-LF, Min(AvgRReq, MaxRReq)) - Min(-AvgRReq, AvgRReq)*Min(LF, -MaxRReq)*Min(-LS, ES + MaxRReq) + Min(-TSC, 2*LF*(-RR - TPC + Min(MinRReq, RR))*Min(AvgRReq, EF)) | 1001.47

##  5 Conclusion

In this study we explored the use of MAP-Elites based Hyper Heuristics (MEHH)
for evolving priority rules for the classical RCPSP. A systematic
investigation was presented comparing MEHH with a standard GPHH. MEHH
outperformed GPHH by maintaining a diverse repertoire of solutions which
allowed for better generalization on new unseen instances. In particular, this
study focused on larger instances which have often been ignored/under-studied
in the existing literature. The rules discovered through MEHH were extensively
compared with human designed priority rules. The comparisons revealed major
performance gaps and even cases where the worst case outputs from MEHH
performed better than the best human designed rules.

In future research we will consider improving the proposed MEHH by
investigating the impact of each of the individual features and also
considering newer more advanced features. We will also consider applying our
proposed MEHH to dynamic problem instances and evaluating its efficacy.

## References

  * [1] R. Kolisch, Serial and parallel resource-constrained project scheduling methods revisited: Theory and computation, European Journal of Operational Research 90 (2) (1996) 320 – 333.
  * [2] R. Kolisch, Efficient priority rules for the resource-constrained project scheduling problem, Journal of Operations Management 14 (3) (1996) 179 – 192\.
  * [3] J. Li, P. Duan, H. Sang, S. Wang, Z. Liu, P. Duan, An efficient optimization algorithm for resource-constrained steelmaking scheduling problems, IEEE Access 6 (2018) 33883–33894.
  * [4] P. Brucker, A. Drexl, R. Möhring, K. Neumann, E. Pesch, Resource-constrained project scheduling: Notation, classification, models, and methods, European journal of operational research 112 (1) (1999) 3–41.
  * [5] S. Chand, Q. Huynh, H. Singh, T. Ray, M. Wagner, On the use of genetic programming to evolve priority rules for resource constrained project scheduling problems, Information Sciences 432 (2018) 146 – 163.
  * [6] M. Dumic, D. Sisejkovic, R. Coric, D. Jakobovic, Evolving priority rules for resource constrained project scheduling problem with genetic programming, Future Generation Computer Systems 86 (2018) 211 – 221.
  * [7] R. Zamani, A competitive magnet-based genetic algorithm for solving the resource-constrained project scheduling problem, European Journal of Operational Research 229 (2) (2013) 552 – 559.
  * [8] D. Debels, M. Vanhoucke, A decomposition-based genetic algorithm for the resource-constrained project-scheduling problem, Operations Research 55 (3) (2007) 457–469.
  * [9] A. A. B. Pritsker, L. J. Waiters, P. M. Wolfe, Multiproject scheduling with limited resources: A zero-one programming approach, Management Science 16 (1) (1969) 93–108.
  * [10] S. Nguyen, M. Zhang, M. Johnston, K. C. Tan, A computational study of representations in genetic programming to evolve dispatching rules for the job shop scheduling problem, Trans. Evol. Comp 17 (5) (2013) 621–639.
  * [11] R. Pellerin, N. Perrier, F. Berthaut, A survey of hybrid metaheuristics for the resource-constrained project scheduling problem, European Journal of Operational Research 280 (2) (2020) 395–416.
  * [12] S. Proon, M. Jin, A genetic algorithm with neighborhood search for the resource-constrained project scheduling problem, Naval Research Logistics (NRL) 58 (03 2011).  doi:10.1002/nav.20439.
  * [13] E. Goncharov, V. Leonov, Genetic algorithm for the resource-constrained project scheduling problem, Automation and Remote Control 78 (2017) 1101–1114.  doi:10.1134/S0005117917060108.
  * [14] D. Debels, M. Vanhoucke, A bi-population based genetic algorithm for the resource-constrained project scheduling problem, in: O. Gervasi, M. L. Gavrilova, V. Kumar, A. Laganá, H. P. Lee, Y. Mun, D. Taniar, C. J. K. Tan (Eds.), Computational Science and Its Applications – ICCSA 2005, Springer Berlin Heidelberg, Berlin, Heidelberg, 2005, pp. 378–387.
  * [15] J. Branke, S. Nguyen, C. W. Pickardt, M. Zhang, Automated design of production scheduling heuristics: A review, IEEE Transactions on Evolutionary Computation 20 (1) (2016) 110–124.
  * [16] M. Affenzeller, S. M. Winkler, G. Kronberger, M. Kommenda, B. Burlacu, S. Wagner, Gaining deeper insights in symbolic regression, in: Genetic Programming Theory and Practice XI, Springer, 2014, pp. 175–190.
  * [17] J.-B. Mouret, J. Clune, Illuminating search spaces by mapping elites (2015).  arXiv:1504.04909.
  * [18] R. Kolisch, S. Hartmann, Heuristic Algorithms for the Resource-Constrained Project Scheduling Problem: Classification and Computational Analysis, Springer US, Boston, MA, 1999, pp. 147–178.
  * [19] J. Blazewicz, J. Lenstra, A. Kan, Scheduling subject to resource constraints: classification and complexity, Discrete Applied Mathematics 5 (1) (1983) 11–24.
  * [20] J. A. Carruthers, A. Battersby, Advances in critical path methods, Journal of the Operational Research Society 17 (4) (1966) 359–380.
  * [21] R. Zamani, A competitive magnet-based genetic algorithm for solving the resource-constrained project scheduling problem, European Journal of Operational Research 229 (2) (2013) 552–559.
  * [22] G. Koulinas, L. Kotsikas, K. Anagnostopoulos, A particle swarm optimization based hyper-heuristic algorithm for the classic resource constrained project scheduling problem, Information Sciences 277 (2014) 680–693.
  * [23] V. Valls, F. Ballestin, S. Quintanilla, Justification and rcpsp: A technique that pays, European Journal of Operational Research 165 (2) (2005) 375–386, project Management and Scheduling.
  * [24] S. Asta, D. Karapetyan, A. Kheiri, E. Özcan, A. J. Parkes, Combining monte-carlo and hyper-heuristic methods for the multi-mode resource-constrained multi-project scheduling problem, Information Sciences 373 (2016) 476–498.
  * [25] S. Wang, Y. Mei, M. Zhang, Novel ensemble genetic programming hyper-heuristics for uncertain capacitated arc routing problem, in: Proceedings of the Genetic and Evolutionary Computation Conference, GECCO ’19, Association for Computing Machinery, New York, NY, USA, 2019, p. 1093–1101.
  * [26] S. Nguyen, M. Zhang, K. C. Tan, Adaptive charting genetic programming for dynamic flexible job shop scheduling, in: Proceedings of the Genetic and Evolutionary Computation Conference, GECCO ’18, Association for Computing Machinery, New York, NY, USA, 2018, p. 1159–1166.
  * [27] S. Chand, H. Singh, T. Ray, Evolving heuristics for the resource constrained project scheduling problem with dynamic resource disruptions, Swarm and Evolutionary Computation 44 (2019) 897–912.
  * [28] S. Nguyen, M. Zhang, K. C. Tan, Adaptive charting genetic programming for dynamic flexible job shop scheduling, in: Proceedings of the Genetic and Evolutionary Computation Conference, GECCO ’18, Association for Computing Machinery, New York, NY, USA, 2018, p. 1159–1166.
  * [29] S. Nguyen, D. Thiruvady, Evolving large reusable multi-pass heuristics for resource constrained job scheduling, in: 2020 IEEE Congress on Evolutionary Computation (CEC), 2020, pp. 1–8.
  * [30] S. Nguyen, M. Zhang, D. Alahakoon, K. C. Tan, People-centric evolutionary system for dynamic production scheduling, IEEE Transactions on Cybernetics 51 (3) (2021) 1403–1416.
  * [31] R. Klein, Bidirectional planning: improving priority rule-based heuristics for scheduling resource-constrained projects, European Journal of Operational Research 127 (3) (2000) 619–638.
  * [32] Q. N. Huynh, S. Chand, H. K. Singh, T. Ray, Genetic programming with mixed-integer linear programming-based library search, IEEE Transactions on Evolutionary Computation 22 (5) (2018) 733–747.
  * [33] T. Frankola, M. Golub, D. Jakobovic, Evolutionary algorithms for the resource constrained scheduling problem, in: ITI 2008 - 30th International Conference on Information Technology Interfaces, 2008, pp. 715–722.
  * [34] D. Šišejković, Evolution of scheduling heuristics for the resource constrained scheduling problem, Master’s thesis, Fakultet elektrotehnike i računarstva, Sveučilište u Zagrebu (2016).
  * [35] E. W. Davis, J. H. Patterson, A comparison of heuristic and optimum solutions in resource-constrained project scheduling, Management Science 21 (8) (1975) 944–955.
  * [36] S. Chand, H. Singh, T. Ray, Evolving rollout-justification based heuristics for resource constrained project scheduling problems, Swarm and Evolutionary Computation 50 (2019) 100556.
  * [37] S. E. Elmaghraby, Activity networks: Project planning and control by network models, John Wiley & Sons, 1977.
  * [38] R. Kolisch, Serial and parallel resource-constrained project scheduling methods revisited: Theory and computation, European Journal of Operational Research 90 (2) (1996) 320 – 333.
  * [39] R. Kolisch, A. Sprecher, Psplib - a project scheduling problem library: Or software - orsep operations research software exchange program, European Journal of Operational Research 96 (1) (1997) 205 – 216.
  * [40] C. Fang, L. Wang, An effective shuffled frog-leaping algorithm for resource-constrained project scheduling problem, Computers & Operations Research 39 (5) (2012) 890–901.


COMPARING FINGERS AND GESTURES FOR BCI CONTROL USING AN OPTIMIZED CLASSICAL
MACHINE LEARNING DECODER
D. Keller1, M. J. Vansteensel1, S. Mehrkanoon2, M. P. Branco1
1UMC Utrecht Brain Center, Department of Neurology and Neurosurgery, Utrecht,
The Netherlands
2Department of Information and Computing Sciences, Utrecht University,
Utrecht, The Netherlands
E-mail: {d.keller, m.pedrosobranco<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>
ABSTRACT: Severe impairment of the central motor network can result in loss of
motor function, clinically recognized as Locked-in Syndrome. Advances in
Brain-Computer Interfaces offer a promising avenue for partially restoring
compromised communicative abilities by decoding different types of hand
movements from the sensorimotor cortex. In this study, we collected ECoG
recordings from 8 epilepsy patients and compared the decodability of
individual finger flexion and hand gestures with the resting state, as a proxy
for a one-dimensional brain-click. The results show that all individual finger
flexion and hand gestures are equally decodable across multiple models and
subjects (>98.0%). In particular, hand movements, involving index finger
flexion, emerged as promising candidates for brain-clicks. When decoding among
multiple hand movements, finger flexion appears to outperform hand gestures
(96.2% and 92.5% respectively) and exhibit greater robustness against
misclassification errors when all hand movements are included. These findings
highlight that optimized classical machine learning models with feature
engineering are viable decoder designs for communication-assistive systems.

INTRODUCTION

Dysfunction of the neuromotor system may precipitate transient or, in severe
cases, enduring global loss of motor control. Global dysfunction may be
referred to as Locked-In Syndrome (LIS) [1], often characterized by
quadriplegia and aphonia. In recent decades, efforts to replace dysfunctional
motor control have seen pioneering developments in Brain-Computer Interfaces
(BCIs) [2]. BCIs extract information directly from cortical activity to
control mechanical or digital effectors without relying on neuromuscular
activation, essentially bypassing the muscular output. Restoration of effector
control can serve several purposes, ranging from object manipulation [3, 4,
5], locomotion and mobility [6] and speech production [7, 8, 9]. However, for
individuals with severe impairment, the restoration of the communicative
agency has been identified as one of the most urgent needs [10].

A simple approach towards communication BCI is automatic letter selection on a
digital keyboard [11]. An attractive signal recording modality for
communication BCIs is electrocorticography (ECoG) due to its high
spatiotemporal precision, good signal-to-noise ratio, and reliable signal
stability over extended periods [11, 12]. Several studies have demonstrated
that hand movement recognition from ECoG recordings can be performed with high
accuracy. Consequently, a variety of hand movements have been explored for
this purpose, including but not limited to finger flexion [13, 14, 15, 16,
17], reaching and grasping [3, 4, 5, 18, 19], and wrist flexion and extension
[15], more complex hand gestures [20, 21, 22, 15] and handwriting [23]. To
identify a reliable motor signature for a unidimensional BCI control signal
(i.e., ’brain-click’) many studies have examined different types of hand
movements in isolation [20, 21, 22, 13, 14], but few have compared different
hand movements against each other within a unified framework.

Figure 1: The four gestures executed by subject S1 - S5.

The goal of this work is to contribute to a deeper understanding of the
decodability of individual finger flexion and hand gestures against the
resting state. Specifically, we explored which hand strategy is the most
promising for a reliable brain-click and which is more transferable across
subjects. In addition, we aim to extend our analysis to the prospect of
multidimensional control with four and eight degrees of freedom (DoF) to
investigate which type of hand movement intrinsically yields a better within
and across-category discriminability. To overcome the notorious data sparsity
in this domain, we employed an optimized feature selection decoder with
different classification models and assessed which (offline) machine learning
approach yields the best performance on individual hand movements (2-DoF),
within hand movement types (4-DoF), and within all hand movements (8-DoF).

MATERIALS AND METHODS

Data description: This study is based on two ECoG data sets, consisting of
finger flexion, including the index finger, little finger, and thumb, and hand
gestures associated with the American Sign Language letters D, F, V, and Y
[24, 21] (see Figure 1).

The data were collected from 8 subjects (NG = 5, NF = 6) in an epilepsy
monitoring unit of the University Medical Center Utrecht (see Table 1). These
subjects had 32, 64 or 128 high-density subdural ECoG electrodes with an
inter-electrode distance of 3 or 4mm and an exposed diameter of 1 to 1.3mm
(AdTech, Racine, USA; or PMT Corporation, Chanhassen, MN, USA) implemented
over the hand-knob region of the sensorimotor cortex. The ECoG data were
recorded using a 128-channel Micromed LTM system (subjects 1 - 5; Treviso,
Italy; 22 bits, hardware bandpass filter 0.15–134.4Hz; sampling frequency
512Hz) and a Blackrock system (subjects 6 - 8; Microsystems LLC, Salt Lake
City, USA, digital bandpass filter 0.3 - 500Hz; sampling frequency 2000Hz).
Data were converted to the BIDS standard format [25].

Table 1: Subject Details

Sub- | Task | Trials | Age | Sex | Hand | Hand- | Hemi- | Grid
---|---|---|---|---|---|---|---|---
ject | (C / T) | edness | sphere | (incl.)
S1 | G | 37 / 74 | 19 | F | Right | Right | Left | 4x8 (32)
F | 90 / 181
S2 | G | 68 / 138 | 45 | F | Left | Left | Right | 8x8 (59)
S3 | G | 34 / 69 | 29 | M | Right | Right | Left | 4x8 (29)
S4 | G | 32 / 67 | 19 | M | Right | Right | Left | 4x8 (31)
F | 90 / 181
S5 | G | 34 / 69 | 42 | M | Right | Right | Left | 4x8 (32)
F | 88 / 177
S6 | F | 89 / 179 | 30 | F | Left | Right | Right | 16x8 (123)
S7 | F | 85 / 171 | 20 | F | Right | Right | Left | 8x8 (64)
S8 | F | 84 / 169 | 36 | F | Right | Right | Left | 16x8 (128)

  * •

Note.  Trials are presented as the ratio of hand movement trials per condition
(C) out of all trials (T; including the rest trials). In Grid, (incl.)
indicates the number of channels included. Abbreviations: Gesture, G; Finger,
F, Male, M; Female, F.

Experimental Design: Subjects were instructed to initiate movements based on
visual cues that were presented in a randomized, event-driven design. For the
gestures, the subjects imitated the depicted gesture after stimulus onset and
maintained the posture until the end of the trial before returning to a
resting position. Each subject performed 10 trials with an intertrial interval
of 4.4s and a run duration of 6.7m. Rest trials were implicitly calculated
from a small time interval before the onset of the next movement. For finger
flexion, the subjects performed two finger flexions immediately after cue
onset and then returned to a resting position afterward. In contrast to the
gestures, each movement was interleaved with an explicit resting trial. The
design consisted of 30 trials with an intertrial interval of 7s and a run
duration of 8.2m. In both experiments, each subject performed the tasks with
the hand contralateral to the grid location, and subject 2 performed the task
twice. In addition, a data glove (5DT, Irvine CA, USA, 20 ms sampling time)
was used during both experiments to record motor activity.

Preprocessing: Data preprocessing included the removal of bad trials and
channels (identified by [24, 21] based on data glove data and raw signal
inspection), followed by common average referencing, notch, and bandpass
filtering (56 Hz - 130 Hz) to remove artifacts. Finger flexion data, sampled
at 2000 Hz, were downsampled to 512 Hz for consistency across subjects. The
data were then subsequently aligned with movement onset markers obtained from
data glove recordings and segmented accordingly. For decoding individual hand
movements and within hand movement types a segmentation window of WF = [-0.5,
1.5s] and WG = [-0.5s, 2.5s] was used. For decoding all hand movements, the
two 4-DoF settings for subject 5 were combined with a common segmentation
window of WFG = [-0.5s, 2s]. In all three settings t = 0 represents the
motion-aligned stimulus presentation.

Features were extracted using a continuous Morlet wavelet transformation,
which produced spectral power features for the high-frequency band (60 Hz -
126 Hz) in 2 Hz frequency bins. To reduce the feature space the power was
averaged and the time dimension was decimated to TG = 154 and TF = 102 time
points per channel for the fingers and gestures, respectively. The resulting
feature vectors were used for subsequent model training. Preprocessing was
conducted in Python (v3.9) using the MNE library (v1.16).

Decoder: The architecture of the decoder, depicted in Figure 2, revolves
around an optimized data-driven feature engineering approach for conventional
classical machine learning classifiers. The decoder encompasses four modules:
Normalization, Incremental Feature Selection (IFS), Feature Reduction (FR),
and Classification.

Figure 2: General decoder architecture. The pipeline is composed of four
modules: Normalization, Incremental Feature Selection (IFS), Feature Reduction
(FR), and a Classification Model. A Box-Cox (BC) transformation and a Standard
Scalar normalize the signal. The incremental selection procedure can recruit
different combinations of Variance Thresholding (VT), Mutual Information
Criterion (MIC), and Recursive Feature Elimination (RFE) to select the
relevant features from the spatio-temporal representation. RFE iteratively
removes the set of S least important features from the pool N until it reaches
KRFE features. After feature selection, the feature space can be further
compressed with various Feature Reduction (FR) methods. The resulting vector
is forwarded to one of four classifiers. A Bayesian optimization algorithm
orchestrates the order and method of feature selection (OIFS, MIFS), other
pipeline configuration ($\lambda$BC, TrVT, KMIC, KRFE, CRFE, CFR), and various
model-specific hyperparameters. Dotted lines represent meta-routing processes;
single and double-lined boxes represent dynamic and predefined processes,
respectively. Abbreviations: LBC, Lower Confidence Bound; EL, Negative
Expected Improvement; PI Negative Probability Improvement

The initial step of the decoder normalises the spectral power of the spatio-
temporal feature vector (N) via a Box-Cox transformation [26], followed by
mean centering and unit variance scaling to ensure data normality and variance
stabilisation. Feature selection employs Variance Thresholding (VT), Mutual
Information Criterion (MIC) [27, 28], and Recursive Feature Elimination (RFE).
RFE iteratively removes a set of features (S = 1e-3 * N) corresponding to the
least important coefficients, and akin to MIC, retains a subset of the best K
temporal features across all channels. Each selection method generates a
binary mask indicating the retained features. The individual modules are
applied incrementally, with the specific combination and its parameter
configuration being delegated to a Bayesian optimization algorithm, which
avoids manual tuning and efficiently navigates through the high-dimensional
parameter space. In particular, this approach aims to balance the advantages
of filter and wrapper methods [29] to remove noisy and redundant features
while prioritising discriminative ones. The Feature Reduction (FR) method can
be applied in isolation or in conjunction with IFS processing. In the
classification phase, four algorithms were selected to compete with each
other: Logistic Regression (LR), Linear Discriminant Analysis (LDA), Boosted
Decision Trees (BDT), and Support Vector Machines (SVM). A Majority Class
Predictor, which predicts the most frequent hand movement, was evaluated on
the data to establish a ’chance’ baseline. Performance is assessed using the
F1 scoring metric, adjusted for label imbalance with the inversely weighted
class distribution.

For this ’black-box’ optimization problem, the Bayesian algorithm [30]
approximates an expensive non-smooth objective function by inference,
essentially guiding the search process based on prior results. To find an
optimal decoder candidate in the large parameter space, a Gaussian Process
model [31] with a hedging portfolio strategy [32]is used, where hedging
probabilistically choose the best acquisition function from three candidates:
Lower Confidence Bound, Negative Expected Improvement or Negative Probability
Improvement. The search is restricted to a maximum of 256 candidates,
exploring a hyperparameter space, ranging from 9 (LDA) to 20 (BDT)
configurations for different algorithms, of which up to 60% are conditional
hyperparameters; the number-of-components hyperperameter for the FR step was
shared among all three methods. Model performance is evaluated using
stratified 10-fold crossvalidation, with the best candidate further assessed
through leave-one-out crossvalidation. The decoder pipeline adheres to the
scikit-learn architecture, ensuring compatibility with the scikit-learn
library and its derivatives. The implementation is in Python 3.9, using
scikit-learn (v1.4.0) and xgboost (v2.0.3).

Statistical Analysis: The analysis relies on a Friedman ANOVA to identify a
general effect and Dunn’s test with Benjamini-Hochberg’s false discoveries
rate correction the post hoc analysis and pairwise comparison. The statistical
analysis was performed in Python 3.9, using scipy (v1.13.0) and scikit-
posthocs (v0.9.0).

RESULTS

Individual Hand Movements: In the context of individual hand movements, all
classification models exhibited a high F1 performance (averaged across
subjects), exceeding 98.0% for each finger flexion and gesture (details
summarised in Table 2 and Figure 3 A). Notably, the Index finger (99.59%),
Gesture V (99.11%), and Gesture F (98.09%) were the most promising candidates
for brain-click BCI control. Interestingly, they share a commonality in index
finger flexion. Within the hand gestures, no gesture significantly
outperformed the others (Friedman ANOVA and Dunn’s test; ns). Similarly,
within the finger flexions, the overall effect was significant, $\chi^{2}$(3)
= 15.32, p = $4.71e^{-4}$), but no finger flexion was significantly different
from the others (Dunn’s test). Moreover, decoding performance remained
remarkably stable for each subject across all hand movements and models, with
consistent trends in variability observed for each subject (i.e., subject 1
consistently had the lowest and highest scores for the fingers and the
gestures, respectively).

Table 2: Mean performance across all subjects of the Optimize Feature
Selection Decoder for Individual Movements (2-DoF), Within Types (4 DoF) and
Within all Hand Movements (8-DoF).

Models | Fingers (vs. rest) | Gestures (vs. rest) | Multi-DoF (incl. rest)
---|---|---|---
Index | Little | Thumb | Gesture D | Gesture F | Gesture V | Gesture Y | 4-Finger | 4-Gesture* | 8-Hand**
| A | F1 | A | F1 | A | F1 | A | F1 | A | F1 | A | F1 | A | F1 | A | F1 | A | F1 | A | F1
Chance | 75.3 | 64.6 | 75.5 | 64.9 | 74.9 | 64.2 | 83.3 | 75.7 | 79.8 | 70.8 | 79.4 | 70.3 | 79.8 | 70.9 | 50.3 | 33.7 | 56.6 | 40.9 | 50.4 | 33.8
LR | 98.2 | 98.3 | 98.3 | 98.4 | 95.9 | 96.0 | 97.4 | 97.5 | 98.2 | 98.2 | 98.7 | 98.7 | 96.6 | 96.8 | 95.7 | 95.7 | 91.3 | 91.5 | 87.4 | 87.7
BDT | 99.6 | 99.6 | 98.1 | 98.1 | 98.1 | 98.0 | 97.6 | 97.7 | 98.7 | 98.7 | 98.4 | 98.5 | 98.2 | 98.2 | 95.8 | 95.8 | 92.0 | 92.0 | 88.6 | 88.7
LDA | 99.3 | 99.3 | 98.7 | 98.8 | 97.5 | 97.6 | 98.1 | 98.1 | 98.7 | 98.7 | 98.0 | 98.0 | 96.4 | 96.6 | 96.1 | 96.1 | 91.1 | 91.1 | 93.1 | 92.9
SVM | 99.5 | 99.4 | 98.8 | 98.8 | 97.7 | 97.6 | 97.9 | 98.0 | 98.8 | 98.89 | 99.1 | 99.1 | 98.4 | 98.4 | 96.7 | 96.6 | 92.5 | 92.5 | 92.7 | 92.6

  * •

Note.  Values are in %. *The 4-DoF gesture decoding includes Gesture F, Y, and
V. ** The 8-DoF decoding of all hand movements was only obtained from subject
5. Abbreviations: Accuracy, A; Versus, vs; Inclusive, incl.

Figure 3: Box plots for different machine learning models of the decoder for
(A) all individual fingers and (B) the three multi-DoF comparisons. Each point
represents one subject, with the upper and lower error bars representing an
interquartile range of 25 and 75, respectively, and where omitted when
performance for one or more subjects exceeded this range. The 4-DoF gesture
decoding includes Gesture F, Y, and V. The red line represents the highest
chance level among all subjects.

Hand Movements Types: For the hand movement types, a different trend emerged.
For 4-DoF classification, finger flexion (95.7% - 96.6%) outperformed the hand
gestures (91.1% - 92.5%; average across subjects), which could be
statistically verified (Dunn’s test), p = 0.031; we excluded the worst
decodable gesture (Gesture D) to ensure a similar task complexity. An
extension to 8-DoF classification preserves a high F1 score for Subject 5
(88.7% - 92.9%), with remarkably minimal confusion between gestures and
fingers. On visual inspection, fingers exhibit more confusion with the resting
state, while gestures are more often confused among themselves (as depicted in
Figure 4).

Figure 4: Confusion matrices of the LR model for subjects 1, 4, and 5. In each
confusion matrix, the horizontal axis represents the predicted hand movement
(or rest), and the vertical axis represents the ground truth hand movement (or
rest). Henceforth, the diagonal elements represent the total correct values
predicted per hand movement. The colour coding represents the proportion of
absolute counts for each hand movement ranging from black (0%) to white
(100%).

Classification Models: When training classification models within an optimised
feature engineering framework, no model emerges as significantly superior to
the others (Friedman ANOVA), $\chi^{2}$(4) = 60.48, p = $2.3e^{-12}$, (Dunn’s
test; ns), although, all perform significantly above chance level (Dunn’s
Test), p < $1.28e^{-8}$. In general, Boosted Decision Trees (BDT) and Support
Vector Machines (SVM) demonstrate the highest classification performance
across all conditions, except for Little Finger and Gesture D.

DISCUSSION

The current work demonstrates that optimised spatio-temporal feature
engineering of finger flexion and hand gestures, recorded from high-density
ECoG, enables reliable decoding for one-dimensional brain-click, 4- and 8-DoF
decoding tasks, even with very small data volumes. Notably, within each
category, no single hand movement emerged as superior decodable. However, upon
qualitative inspection, the index finger and Gestures V and F, all sharing
index finger flexion, appeared as the most promising candidates. Moreover, the
fingers exhibited a better performance in more complex 4- and 8-DoF decoding;
finger flexion may possess more discriminative properties for multi-DoF tasks.
Finally, no classical machine learning model outperformed the others, but BDT
and SVM may have a small advantage.

Hand gestures entail a more complex interplay of motor components than simple
finger flexion, including wrist and finger flexion, and lateral extension,
among others. However, our results revealed that individual finger flexion
alone yields a near perfect neuroelectrical signature. Consequently, the
addition of supplementary motor components may not increase decodability but
rather confront the decoder with a motor signature that inherently has a
higher variability in its signal. As we evaluated the decoder on small data
volumes, an effect of additional components on decoder performance may emerge
with larger sample sizes. Extending this rationale, gestures might possess a
more intricate spatial and temporal pattern, recruiting various neuronal
populations from a larger effector field in the sensorimotor cortex. The
absence of advanced feature transformation techniques in our classical machine
learning approach might hinder the decoder from exploiting the full potential
of the gesture’s electrical signature. In addition, differences in sample
size, with fewer than 30 trials for fingers and only 10 trials for gestures,
may impact the comparison, while the overall small sample size may not provide
sufficient power to detect a potential statistical difference.

In line with the evident discernible difference between the 4-DoF types,
misclassification errors were more pronounced for gestures, with frequent
confusion among different fingers and the resting state. Notably, substantial
confusion between hand movements and the resting state suggests a potential
contamination of the rest periods with unintended movement. This may be
attributed to the design of the experimental design of the gestures, which
lacked separate explicit rest trials. Incorporating a threshold based on data
glove recordings could be beneficial. However, defining true rest may not be
practically feasible for real BCI applications, as it necessitates extensive
subject training to suppress such activity [11, 12, 19], and a more
naturalistic approach would be to build a decoder that can successfully
discriminate between meaningful and non-task-related sporadic motor activity.
Furthermore, some features along the temporal dimension may not reflect actual
motor activity, but ’resting’ activity, especially before the movement onset
and towards the end of the segmentation window.

The proposed decoder design surpassed prior approaches evaluated on gesture
data for four out of five subjects [20, 21], albeit a direct comparison is
difficult due to differences in task complexity. For communication assistance
systems, the proposed decoder design can offer a viable alternative to deep
learning approaches for one-dimensional brain-click tasks [19] and even larger
DoF applications [14, 17], where data acquisition is challenging. Although
trained offline, the decoder can process individual segments of preprocessed
data in as little as 2 - 10ms. Consequently, although theoretically deployable
in an online setting, regular offline retraining on new data is necessary to
address concept drifts for ensuring long-term stability, in particular for
individuals with neurodegenerative diseases. Importantly, the experiments
involved movement execution by epileptic individuals, further validation in
attempted movement is imperative to extend applicability to online BCI
settings for individuals with LIS.

CONCLUSION

Electrocorticography data provides a high-resolution spatiotemporal feature
representation, forming a suitable foundation to tailor an optimised classical
machine learning decoder with automatic feature engineering to the large
feature space. We demonstrate that within this framework, both finger flexion
and hand gestures enable reliable decoding across multiple subjects, and when
extended to a multiple degrees of freedom, maintain high discriminability
between hand movements. For click-based letter selection in communication-
assistive BCI systems, the index finger flexion emerges as an optimal
candidate. Moreover, all tested models consistently exhibit high
classification performance across multiple subjects - a comparable performance
to deep learning approaches.

REFERENCES

## References

  * [1] Fred Plum and Jerome B Posner  “The diagnosis of stupor and coma”  Oxford University Press, USA, 2000
  * [2] Nitish V Thakor  “Translating the brain-machine interface”  In _Science translational medicine_ 5.210 American Association for the Advancement of Science, 2013, pp. 210ps17–210ps17
  * [3] Leigh R Hochberg et al.  “Reach and grasp by people with tetraplegia using a neurally controlled robotic arm”  In _Nature_ 485.7398 Nature Publishing Group UK London, 2012, pp. 372–375
  * [4] Jennifer L Collinger et al.  “High-performance neuroprosthetic control by an individual with tetraplegia”  In _The Lancet_ 381.9866 Elsevier, 2013, pp. 557–564
  * [5] Chad E Bouton et al.  “Restoring cortical control of functional movement in a human with quadriplegia”  In _Nature_ 533.7602 Nature Publishing Group UK London, 2016, pp. 247–250
  * [6] Henri Lorach et al.  “Walking naturally after spinal cord injury using a brain–spine interface”  In _Nature_ Nature Publishing Group UK London, 2023, pp. 1–8
  * [7] David A Moses et al.  “Neuroprosthesis for decoding speech in a paralyzed person with anarthria”  In _New England Journal of Medicine_ 385.3 Mass Medical Soc, 2021, pp. 217–227
  * [8] Francis R Willett et al.  “A high-performance speech neuroprosthesis”  In _Nature_ 620.7976 Nature Publishing Group UK London, 2023, pp. 1031–1036
  * [9] Sean L Metzger et al.  “A high-performance neuroprosthesis for speech decoding and avatar control”  In _Nature_ 620.7976 Nature Publishing Group UK London, 2023, pp. 1037–1046
  * [10] Yu Kageyama et al.  “Nationwide survey of 780 Japanese patients with amyotrophic lateral sclerosis: their status and expectations from brain–machine interfaces”  In _Journal of Neurology_ 267 Springer, 2020, pp. 2932–2940
  * [11] Mariska J Vansteensel et al.  “Fully implanted brain–computer interface in a locked-in patient with ALS”  In _New England Journal of Medicine_ 375.21 Mass Medical Soc, 2016, pp. 2060–2066
  * [12] Elmar GM Pels et al.  “Stability of a chronic implanted brain-computer interface in late-stage amyotrophic lateral sclerosis”  In _Clinical Neurophysiology_ 130.10 Elsevier, 2019, pp. 1798–1803
  * [13] Tao Xie et al.  “Classifying multiple types of hand motions using electrocorticography during intraoperative awake craniotomy and seizure monitoring processes—case studies”  In _Frontiers in neuroscience_ 9 Frontiers Media SA, 2015, pp. 353
  * [14] Ziqian Xie, Odelia Schwartz and Abhishek Prasad  “Decoding of finger trajectory from ECoG using deep learning”  In _Journal of neural engineering_ 15.3 IOP Publishing, 2018, pp. 036009
  * [15] Tessy M Thomas et al.  “Simultaneous classification of bilateral hand gestures using bilateral microelectrode recordings in a tetraplegic patient”  In _MedRxiv_ Cold Spring Harbor Laboratory Press, 2020, pp. 2020–06
  * [16] Jathurshan Pradeepkumar et al.  “Decoding of Hand Gestures from Electrocorticography with LSTM Based Deep Neural Network”  In _2021 43rd Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_ , 2021, pp. 420–423
  * [17] Lin Yao, Bingzhao Zhu and Mahsa Shoaran  “Fast and accurate decoding of finger movements from ECoG through Riemannian features and modern machine learning techniques”  In _Journal of Neural Engineering_ 19.1 IOP Publishing, 2022, pp. 016037
  * [18] Tianxiao Jiang et al.  “Characterization and decoding the spatial patterns of hand extension/flexion using high-density ECoG”  In _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ 25.4 IEEE, 2017, pp. 370–379
  * [19] Nathan Crone et al.  “A click-based electrocorticographic brain-computer interface enables long-term high-performance switch-scan spelling”  In _Research Square_ American Journal Experts, 2023
  * [20] Martin G Bleichner et al.  “Give me a sign: decoding four complex hand gestures based on high-density ECoG”  In _Brain Structure and Function_ 221 Springer, 2016, pp. 203–216
  * [21] Mariana P Branco et al.  “Decoding hand gestures from primary somatosensory cortex using high-density ECoG”  In _NeuroImage_ 147 Elsevier, 2017, pp. 130–142
  * [22] Yue Li et al.  “Gesture decoding using ECoG signals from human sensorimotor cortex: a pilot study”  In _Behavioural neurology_ 2017 Hindawi, 2017
  * [23] Mingyue Xu et al.  “Swin-TCNet: Swin-based temporal-channel cascade network for motor imagery iEEG signal recognition”  In _Biomedical Signal Processing and Control_ 85 Elsevier, 2023, pp. 104885
  * [24] Jeroen CW Siero et al.  “BOLD matches neuronal activity at the mm scale: a combined 7 T fMRI and ECoG study in human sensorimotor cortex”  In _Neuroimage_ 101 Elsevier, 2014, pp. 177–184
  * [25] Cyril R Pernet et al.  “EEG-BIDS, an extension to the brain imaging data structure for electroencephalography”  In _Scientific data_ 6.1 Nature Publishing Group UK London, 2019, pp. 103
  * [26] George EP Box and David R Cox  “An analysis of transformations”  In _Journal of the Royal Statistical Society Series B: Statistical Methodology_ 26.2 Oxford University Press, 1964, pp. 211–243
  * [27] Alexander Kraskov, Harald Stögbauer and Peter Grassberger  “Estimating mutual information”  In _Physical review E_ 69.6 APS, 2004, pp. 066138
  * [28] Brian C Ross  “Mutual information between discrete and continuous data sets”  In _PloS one_ 9.2 Public Library of Science San Francisco, USA, 2014, pp. e87357
  * [29] Ioannis Tsamardinos and Constantin F Aliferis  “Towards principled feature selection: Relevancy, filters and wrappers”  In _International Workshop on Artificial Intelligence and Statistics_ , 2003, pp. 300–307  PMLR
  * [30] Jonas Mockus and Jonas Mockus  “The Bayesian approach to local optimization”  Springer, 1989
  * [31] Carl Edward Rasmussen and Christopher KI Williams  “Gaussian processes for machine learning”  Springer, 2006
  * [32] Matthew Hoffman, Eric Brochu and Nando De Freitas  “Portfolio Allocation for Bayesian Optimization.”  In _UAI_ , 2011, pp. 327–336