# Solving large flexible job shop scheduling instances by generating a diverse
set of scheduling policies with deep reinforcement learning

Imanol Echeverria
TECNALIA, Basque Research and Technology Alliance (BRTA)
Mikeletegi Pasealeakua 7, 20009 Donostia-San Sebastián (Spain)
<EMAIL_ADDRESS>
&Maialen Murua
TECNALIA, Basque Research and Technology Alliance (BRTA)
Mikeletegi Pasealeakua 7, 20009 Donostia-San Sebastián (Spain)
<EMAIL_ADDRESS>
&Roberto Santana
Computer Science and Artificial Intelligence Department
University of the Basque Country
Lardizabal Pasealeakua 1, 20018 Donostia-San Sebastián (Spain).
<EMAIL_ADDRESS>
Corresponding author.

###### Abstract

The Flexible Job Shop Scheduling Problem (FJSSP) has been extensively studied
in the literature, and multiple approaches have been proposed within the
heuristic, exact, and metaheuristic methods. However, the industry’s demand to
be able to respond in real-time to disruptive events has generated the
necessity to be able to generate new schedules within a few seconds. Among
these methods, under this constraint, only dispatching rules (DRs) are capable
of generating schedules, even though their quality can be improved. To improve
the results, recent methods have been proposed for modeling the FJSSP as a
Markov Decision Process (MDP) and employing reinforcement learning to create a
policy that generates an optimal solution assigning operations to machines.
Nonetheless, there is still room for improvement, particularly in the larger
FJSSP instances which are common in real-world scenarios. Therefore, the
objective of this paper is to propose a method capable of robustly solving
large instances of the FJSSP. To achieve this, we propose a novel way of
modeling the FJSSP as an MDP using graph neural networks. We also present two
methods to make inference more robust: generating a diverse set of scheduling
policies that can be parallelized and limiting them using DRs. We have tested
our approach on synthetically generated instances and various public
benchmarks and found that our approach outperforms dispatching rules and
achieves better results than three other recent deep reinforcement learning
methods on larger FJSSP instances.

_Keywords_ Flexible job shop scheduling $\cdot$ Deep reinforcement learning
$\cdot$ Heterogeneous graph neural networks $\cdot$ Markov decision process

##  1 Introduction

Planning consists of assigning a series of operations to a set of resources
within a specified time period to optimize one or multiple objectives. It is a
widely studied problem in the literature due to its impact on the
manufacturing industry (Gupta and Sivakumar, 2006) and various other real-
world scenarios (Çaliş and Bulkan, 2015). One of the most well-known variants
is the job shop scheduling problem (JSSP) (Manne, 1960), in which a set of
jobs, composed of operations, must be assigned to a set of machines. Each one
of these operations can only be performed on a single given machine, and the
order of the operations within a job must be respected. The objective of the
problem is to create a resource allocation that optimizes one or multiple key
performance indicators, with makespan being the most commonly used in the
literature, defined as the time to complete all operations, as well as other
metrics related to job tardiness or production costs such as energy
consumption (Xiong et al., 2022). Due to its complexity, this problem is known
to be NP-Hard (Garey et al., 1976). A well-studied extension of the JSSP is
the Flexible Job Shop Scheduling Problem (FJSSP), where the difference lies in
the fact that each one of the operations can be processed on a set of machines
instead of a single one. This problem is also NP-Hard as it adds the
complexity of determining the optimal route for operations, resulting in a
significantly larger search space of possible solutions.

Multiple approaches have been proposed in the literature to find optimal
solutions for the FJSSP. Within the family of exact optimization methods,
dynamic programming-based methods have been proposed (Kim et al., 2022),
although the most common ones are based on linear programming or mixed linear
programming (Meng et al., 2020; Tutumlu and Saraç, 2023). However, due to the
NP-hardness of the FJSSP, these methods cannot handle medium-large instances
of the problem in a dynamic manner, that is, cannot generate optimal or
pseudo-optimal solutions in real-time, and the sizes of the instances
typically used for benchmarking are much smaller than those encountered in
real-world scenarios. For example, in Meng et al. (2020), the authors use the
first ten instances from the well-known Brandimarte benchmark (Brandimarte,
1993), with the largest instance consisting of only twenty jobs and fifteen
machines. This size is much smaller compared to real-world scenarios found in
the industry, which can involve hundreds or thousands of jobs and dozens of
machines.

On the other hand, metaheuristic algorithms, particularly genetic algorithms
(Zhang et al., 2020b; Chen et al., 2020), as well as approaches based on tabu
search (Hajibabaei and Behnamian, 2021) or simulated annealing (Defersha et
al., 2022), have been widely employed. Although metaheuristic methods are more
efficient than exact methods, both require a considerable amount of time to
obtain satisfactory solutions as the algorithms need to iterate to find these
solutions. Therefore, neither of these two approaches is suitable for
environments with disruptive events, such as machine breakdowns or new jobs
that need to be accommodated in real-time. Dispatching rules (DRs) (Calleja
and Pastor, 2014; Chen and Matis, 2013), methods framed within the family of
heuristics, are the most commonly used options in such situations. These rules
consist of two criteria: the sequencing order of operations and the assignment
of each operation to a machine. They also have the advantage of being capable
of generating solutions quickly for large instances of the FJSSP. However,
they present some disadvantages. Both criteria need to be manually designed by
an expert, which requires time and effort, and they are also static, making
them strategies that cannot adapt to diverse problem instances and cannot find
global optimal solutions (Luo et al., 2021).

Reinforcement learning (RL) and deep reinforcement learning (DRL) have emerged
as promising alternatives to tackle problems that can be modeled as a Markov
Decision Process (MDP), including other NP-Hard problems such as the Traveling
Salesman Problem (Cappart et al., 2021; Mazyavkina et al., 2021) or the
Vehicle Routing Problem (Joe and Lau, 2020). Scheduling problems such as the
FJSSP have received less attention in the literature but, in recent years,
various DRL approaches have been proposed to tackle the JSSP and the FJSSP,
with the former being more common due to be easier to model. To apply RL to
scheduling problems, the problem is first modeled as an MDP, and RL algorithms
are used to learn a policy that sequentially assigns operations to machines
while attempting to maximize the accumulated reward, which is usually designed
to minimize the makespan.

In general, the proposed DRL approaches outperform DRs, but there are still
limitations that need to be addressed to apply these methods in real-world
scenarios. One significant difference between current approaches and real-
world scenarios is the size of the benchmarks. For instance, in Song et al.
(2022) and Wang et al. (2023), two benchmark sets were used: the first ten
instances from Brandimarte (Brandimarte, 1993) and the sets of instances
proposed by Hurink (the rdata, vdata, and edata sets) (Hurink et al., 1994),
with a maximum of twenty machines and thirty jobs. These instance sets have a
limited number of elements compared to real-world scenarios, where there can
be thousands of jobs and dozens of machines. In this regard, Behnke and Geiger
(2012) proposed a benchmark of more realistic instances of larger sizes, with
up to one hundred jobs and sixty machines, which is utilized by Lei et al.
(2022) in their approach, in addition to the vdata set to evaluate the
performance on smaller instances. Although near-optimal results can be
obtained on small instances, there is more room for improvement on large
instances. Therefore, the objective of this paper is to propose a novel end-
to-end DRL framework that robustly solves FJSSP instances, outperforming DRs
and recent DRL-based methods. The aim is to narrow the gap between instances
typically used in the literature and problems that resemble real-world
scenarios. The contributions of this paper can be summarized as follows:

  * •

A novel representation of the state space is achieved by minimizing the
required information and appropriately distributing it among three types of
nodes (jobs, operations, and machines), along with a set of edges that
facilitate the transmission of information between nodes. As the instance
progresses, unnecessary edges and nodes are eliminated to enhance the
efficiency of the training and inference phases.

  * •

Limiting the action space using masks and employing DRs to achieve faster
model training and more robust inference behavior.

  * •

A method for generating a diverse set of scheduling policies through Bayesian
optimization and a method for selecting the best policies using K-Nearest
Neighbors, efficiently parallelizing models during inference to improve
performance.

##  2 Preliminaries

In this section, we will briefly describe how the FJSSP is formulated, as well
as introducing key concepts of RL and Graph Neural Networks (GNNs) that will
be used in our proposed method.

###  2.1 Problem formulation

An instance of the FJSSP problem is defined by a set of jobs
$\mathcal{J}=\\{j_{1},j_{2},\ldots,j_{n}\\}$, where each $j_{i}\in\mathcal{J}$
is composed of a set of operations
$\mathcal{O}_{j_{i}}=\\{o_{i1},o_{i2},\ldots o_{im}\\}$, and each operation
can be performed on one or more machines from the set
$\mathcal{M}=\\{m_{1},m_{2},\ldots,m_{p}\\}$. The processing time of operation
$o_{ij}$ on machine $m_{k}$ is defined as $p_{ijk}\in\mathbb{R}^{+}$. We
define $\mathcal{M}_{o_{ij}}\subseteq\mathcal{M}$ as the subset of machines on
which that operation $o_{ij}$ can be processed,
$\mathcal{O}_{j_{i}}\subseteq\mathcal{O}$ as the set of operations that belong
to the job $j_{i}$ where
$\bigcup\limits_{i=1}^{n}\mathcal{O}_{j_{i}}=\mathcal{O}$, and
$\mathcal{O}_{m_{k}}\subseteq\mathcal{O}$ as the set of operations that can be
performed on machine $m_{k}$. The execution of operations on machines must
satisfy a series of constraints:

  * •

All machines and jobs are available at time zero.

  * •

A machine can only execute one operation at a time, and the execution cannot
be interrupted.

  * •

An operation can only be performed on one machine at a time.

  * •

The execution order of the set of operations $\mathcal{O}_{j_{i}}$ for every
$j_{i}\in\mathcal{J}$ must be respected.

  * •

Job executions are independent of each other, meaning no operation from any
job precedes or has priority over the operation of another job.

In essence, the FJSSP combines two problems: a machine selection problem,
where the most suitable machine is chosen for each operation, a routing
problem, and a sequencing or scheduling problem, where the sequence of
operations on a machine needs to be determined. Given an assignment of
operations to machines, the completion time of a job, $j_{i}$, is defined as
$C_{j_{i}}$, and the makespan of a schedule is defined as
$C_{max}=\max\limits_{j_{i}\in\mathcal{J}}C_{j_{i}}$, which is the most common
objective to minimize.

###  2.2 Reinforcement Learning and Proximal Policy Optimization

MDPs provide a mathematical framework for modeling decision-making processes.
Formally, a MDP is defined as a tuple $(\mathcal{S},\mathcal{A},R,P)$, where
$\mathcal{S}$ is a finite set of elements called states, $\mathcal{A}$ is a
finite set of elements called actions,
$R:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ represents the reward
the agent receives when taking action $a$ in state $s$, and $P$ is a function
$P:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow[0,1]$ that
represents the probability of transitioning from state $s$ to $s^{\prime}$
when taking action $a$. The goal of RL is to learn a policy
$\pi_{\theta}:\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$, which is
parameterized by $\theta$, in order to maximize the expected cumulative reward
by mapping states to actions (Sutton and Barto, 2018). The future rewards that
the agent receives are usually multiplied by a discount factor
$\gamma\in[0,1]$.

There are multiple algorithms to find the parameters $\theta$ that generate an
optimal policy, and one of the most widely used is Proximal Policy
Optimization (PPO) (Schulman et al., 2017). PPO iteratively updates the policy
using a surrogate objective function based on the advantage function. The
advantage function, denoted by
$A:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$, measures how good an
action is compared to the average of all possible actions in a particular
state. The surrogate objective function encourages policy updates that
increase the probability of actions that maximize the expected cumulative
reward. At time step $t$, the policy will be updated to minimize the defined
loss function:

$L^{CLIP}(\theta)=\hat{E}_{t}\left[\min\left(r_{t}(\theta)\hat{A}_{t},\operatorname{clip}\left(r_{t}(\theta),1-\varepsilon,1+\varepsilon\right)\hat{A}_{t}\right)\right]$
(1)

where $\hat{E}_{t}$ denotes the empirical expectation over future time steps,
$\hat{A}_{t}$ represents the empirical advantage, $r_{t}$ is the ratio of
change between the current policy and the previous policy, the function
$\operatorname{clip}$ restricts the value of the first parameter within the
range defined by the second and third parameters it receives, and
$\varepsilon$ is a hyperparameter (usually 0.1 or 0.2) that is used in the
$\operatorname{clip}$ function and determines the extent to which the policy
is allowed to be updated.

###  2.3 Graph Neural Networks

GNNs (Scarselli et al., 2008) have received significant attention in recent
years (Zhang et al., 2021; Zhou et al., 2020) due to their ability to
efficiently extract information from structured data in the form of graphs,
where predictions are made based on the edges connecting the nodes. A directed
graph $\mathcal{G}$ is defined as a tuple of two sets
$(\mathcal{V},\mathcal{E})$ where $\mathcal{V}$ contains nodes,
$\mathcal{V}=\\{1,\ldots,n\\}$, and $\mathcal{E}$ edges,
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$, where $(i,j)\in E$
represents an edge from a node $i$ to a node $j$. It is assumed that every
node has an initial representation
$\boldsymbol{h_{i}}^{(0)}\in\mathbb{R}^{d_{0}}$, where $d_{0}$ is the number
of the initial features of the node, and that all edges are directed, so an
undirected graph can be represented with bidirectional edges. Multiple GNN
architectures have been proposed (Wu et al., 2020), but one of the most
commonly used ones is the Graph Attention Network (GAT) (Veličković et al.,
2017), where each node updates its representation by aggregating the
representations of its neighbors. Specifically, given a set of node
representations
$\mathcal{H}=\\{\boldsymbol{h}_{1},\boldsymbol{h}_{2}\ldots\boldsymbol{h}_{n}\\}$,
the update of the representation of a node $\boldsymbol{h}_{i}^{\prime}$ is
performed as follows:

$\boldsymbol{h_{i}}^{\prime}=f_{\theta}\left(\boldsymbol{h}_{i},\text{AGGREGATE}\left(\left\\{\boldsymbol{h}_{j}\mid
j\in\mathcal{N}_{i}\right\\}\right)\right)$ (2)

where $\mathcal{N}_{i}$ represents the set of nodes connected to node $i$,
while $f$ and AGGREGATE are functions that define the GNN architecture. Many
GNNs assign equal importance to all neighboring nodes, but in Veličković et
al. (2017), the authors proposed a weighted average of the representations of
$\mathcal{N}_{i}$ through a scoring function that learns the importance of
node $j$’s features to node $i$. The scoring function is defined as:

$e\left(\boldsymbol{h}_{i},\boldsymbol{h}_{j}\right)=\text{ LeakyReLU
}\left(\boldsymbol{a}^{\top}\cdot\left[\boldsymbol{W}\boldsymbol{h}_{i}\|\boldsymbol{W}\boldsymbol{h}_{j}\right]\right)$
(3)

where $\boldsymbol{a}\in\mathbb{R}^{2d^{\prime}}$ and
$\boldsymbol{W}\in\mathbb{R}^{d^{\prime}\times d_{0}}$ are learned vectors,
$d^{\prime}$ is the length of the vector $\boldsymbol{h_{i}}^{\prime}$ and
$||$ denotes vector concatenation. These scores are then normalized using the
softmax function across all neighbors of the node, resulting in:

$\alpha_{ij}=\operatorname{softmax}_{j}\left(e\left(\boldsymbol{h}_{i},\boldsymbol{h}_{j}\right)\right)=\frac{\exp\left(e\left(\boldsymbol{h}_{i},\boldsymbol{h}_{j}\right)\right)}{\sum_{j^{\prime}\in\mathcal{N}_{i}}\exp\left(e\left(\boldsymbol{h}_{i},\boldsymbol{h}_{j^{\prime}}\right)\right)}$
(4)

Finally, to obtain $\boldsymbol{h}_{i}^{\prime}$, a weighted sum of the
representation of neighboring nodes is computed using the attention
coefficients, followed by a non-linear activation function $\sigma$:

$\boldsymbol{h}_{i}^{\prime}=\sigma\left(\sum_{j\in\mathcal{N}_{i}}\alpha_{ij}\cdot\boldsymbol{W}\boldsymbol{h}_{j}\right)$
(5)

Brody et al. (2021) proposed GATv2, which introduces a change in the internal
order of these operations to achieve a more expressive attention mechanism.
The problem with the original order is that the vectors $\boldsymbol{W}$ and
$\boldsymbol{a}$ are applied consecutively, leading to them being collapsed
into a single linear layer. The proposed solution is to apply $\boldsymbol{a}$
after the non-linear function and $\boldsymbol{W}$ after the concatenation,
resulting in a more expressive score for each node-key pair. The updated
formula becomes:

$e\left(\boldsymbol{h}_{i},\boldsymbol{h}_{j}\right)=\boldsymbol{a}^{\top}\operatorname{LeakyReLU}\left(\boldsymbol{W}\cdot\left[\boldsymbol{h}_{i}\|\boldsymbol{h}_{j}\right]\right)$
(6)

and $\boldsymbol{W}$ can be divided into
$[\boldsymbol{W_{1}}\|\boldsymbol{W_{2}}]|\boldsymbol{W_{1}},\boldsymbol{W_{2}}\in\mathbb{R}^{d^{\prime}\times
d}$ where each matrix represents the left half and right half of the columns
of $\boldsymbol{W}$.

##  3 Related work

In this section, we will provide a brief analysis of the proposed methods for
addressing the JSSP and the FJSSP, focusing on how they represent states and
the algorithms they use to learn the scheduling policy.

Among the works addressing the JSSP, in Lin et al. (2019) a framework was
proposed for smart manufacturing based on an off-policy RL method, deep
Q-network (DQN), and edge computing is proposed to address the decisions of
multiple edge devices. They encapsulate the state information into a feature
vector that captures the relevant production state information and they employ
a multilayer perceptron (MLP) to determine the value of a state. Similarly, in
Luo (2020), a deep Q-network is employed to extract information from the
production status using seven features that are used to select from six
composite DRs. Instead of using vectors and off-policy methods, Liu et al.
(2020) used convolutional neural networks (CNNs) and the deep deterministic
policy gradient (DDPG) algorithm, which belongs to the actor-critic family,
was employed to tackle the problem. In this case, the input values are not a
vector but a matrix with three channels: the first channel contains the
process times of the operations, the second channel contains information about
whether a job has been assigned to a machine, and the third channel indicates
whether a job has been completed. This approach has the advantage of not
having to design complicated features that reflect the production state, which
can be time-consuming. One of the major drawbacks of these approaches is that
they are not size-agnostic, meaning that the number of jobs or machines must
be fixed since the inputs and outputs of an MLP cannot be variable.
Additionally, much of the information in the JSSP, especially in the FJSSP,
lies in its structure and the relationship between operations and machines,
which is not directly encapsulated in these state representations.

Building upon the work done in other domains such as vehicle routing (Tolstaya
et al., 2021), the information of a state can be represented using a graph,
where nodes possess features that reflect their state and share information
with connected nodes through edges. The state information can be extracted
using GNNs and attention mechanisms (Vaswani et al., 2017; Veličković et al.,
2017), which, unlike MLPs and CNNs, allow for capturing the structure of the
state and its relationships. Moreover, they are size-agnostic, meaning the
number of nodes in the problem is irrelevant. One way to represent JSSP and
FJSSP is through a disjunctive graph (Balas, 1969), where nodes represent
operations and two types of edges are employed. Directed edges indicate
precedence among operations belonging to the same job, while undirected edges
connect pairs of operations that can be performed on the same machine. With
this representation, a solution is defined as an acyclic orientation of the
undirected edges.

Using this representation, Zhang et al. (2020a) proposed a method to learn DRs
using an end-to-end DRL agent for the JSSP using the PPO algorithm to train
the policy. Song et al. (2022) extended this work to the FJSSP using a
Heterogeneous GNN to capture state information and the relationships between
operation and machine nodes. The authors in Lei et al. (2022) proposed to
model the FJSSP as a multiple MDP, and multi-pointer graph networks which is a
GNN that is employed to learn two sub-policies: one for operation selection
and another for machine assignment. Wang et al. (2023) used attention models
for deep feature extraction using a dual-attention network with operation
message attention blocks and machine message attention blocks. Lastly, Du et
al. (2022) modeled the multi-objective FJSSP with crane transportation and
setup times, and total processing time and total energy consumption are
optimized using a DQN model.

##  4 Proposed method

In this section, we will describe the proposed method to dynamically solve the
FJSSP. First, we will explain how we have formulated the problem as an MDP,
defining the state space, the action space, the reward function and state
transitions. Secondly, we will present how we have modeled the policy using a
GNN based on attention mechanisms and message passing between different nodes.
Thirdly, we will describe two methods to make these approaches more robust,
particularly in large instances, by applying masks based on DRs and using a
novel method for generating a diverse set of scheduling policies.

###  4.1 The FJSSP as a Markov Decision Process.

To define an MDP, we need to specify the space of states, the action space,
the reward and the state transitions. In broad terms, at each step, when an
operation is completed, the agent must decide which operation to perform next
and in which machine.

State space. Each state, which at step $t$ is denoted as $s_{t}$, is a
heterogeneous graph $\mathcal{G}_{t}=(\mathcal{V}_{t},\mathcal{E}_{t})$, where
the set $\mathcal{V}_{t}$ consists of three types of nodes: the set of
operations $\mathcal{O}_{t}$, the set of machines $\mathcal{M}_{t}$, and the
set of jobs $\mathcal{J}_{t}$. Each of these types of nodes has the following
features:

  * •

The operation nodes have a boolean feature indicating whether the operation
can be performed or not and the pending processing time. For an operation, it
is calculated as the sum of the average processing times of the remaining
operations.

  * •

The machine nodes have two features: the time at which the last assigned
operation on the machine is completed and the utilization percentage of that
machine.

  * •

The features of job nodes contain most of the information. Their features
include whether the job has been completed, the time at which its last
assigned operation is finished, the number of operations remaining for that
job and the pending processing time.

The set of edges $\mathcal{E}_{t}$ consists of the following types of edges:

  * •

Undirected edges between machines and operations, connecting operations that
can be executed in machines and indicating their processing time. An
undirected edge can be represented with two sets of directed edges, one
between the operations and the machines,
$\mathcal{OM}_{t}\subseteq\mathcal{O}_{t}\times\mathcal{M}_{t}$, and another
between the machines and the operations,
$\mathcal{MO}_{t}\subseteq\mathcal{M}_{t}\times\mathcal{O}_{t}$. Both have
features which are the processing time, the ratio of the processing time to
the maximum processing time of the operation and the ratio to the maximum of
all the options that the machine can process.

  * •

Directed edges between operations and jobs, indicating that an operation
belongs to a job,
$\mathcal{OJ}_{t}\subseteq\mathcal{O}_{t}\times\mathcal{J}_{t}$.

  * •

Directed edges between operations, indicating that an operation immediately
precedes the other,
$\mathcal{OO}_{t}\subseteq\mathcal{O}_{t}\times\mathcal{O}_{t}$.

  * •

Directed edges between machines and jobs, adding the processing time and the
gap that would be generated by executing that job (i.e., the first operation
that can be performed for that job) on that machine. It also has as features
the same ratios as the edges between the operations and machines. The set of
edges is defined as
$\mathcal{MJ}_{t}\subseteq\mathcal{M}_{t}\times\mathcal{J}_{t}$ and each edge
$(m_{k},j_{i})\in\mathcal{MJ}_{t}$ have features defined as $pg_{ki}$.

  * •

Jobs are connected to each other to enable information communication among
their states, $\mathcal{JJ}_{t}\subseteq\mathcal{J}_{t}\times\mathcal{J}_{t}$.

  * •

Machines are connected to each other to enable information communication among
their states, $\mathcal{MM}_{t}\subseteq\mathcal{M}_{t}\times\mathcal{M}_{t}$.

It should be noted that in our case, the number of features associated with
nodes is reduced, resulting in a lower computational cost as less information
needs to be calculated at each state. In comparison, in Song et al. (2022),
operations have six features and machines have three features and in Wang et
al. (2023) operations have ten features, machines eight and compatible
operation-machine pairs also eight features.

Furthermore, the inclusion of a job node offers two main advantages. Firstly,
it significantly reduces the action space, as the number of jobs is usually
much smaller than the number of operations. Secondly, it concentrates the
relevant information of the problem into a single type of node, making the
learning process more efficient and the inference faster.

Action space. At each step $t$, the set of actions $\mathcal{A}_{t}$ is
composed of all pairs of compatible jobs and machines. When choosing a job,
the first operation of the job that has not yet been scheduled is selected.
This set corresponds to the set of edges that connect jobs and machines,
$\mathcal{MJ}_{t}$. As explained in the next section, to improve learning
efficiency and stabilize inference, based on a DR, some of these actions will
be excluded from being eligible by the policy.

State Transitions. Since the MDP formulation of the FJSSP is deterministic,
performing an action $a_{t}$ in state $s_{t}$ transitions to a new state
$s_{t+1}$ where the operation has been executed. During each state transition,
the heterogeneous graph undergoes three transformations. Firstly, the
operation node that has been executed is removed, along with all the edges
connected to it. Secondly, the edges of the executed job are removed, and new
edges are created with the next operation that needs to be performed. Lastly,
the features of all nodes and edges are updated.

Figure 1: In the top figure, the original FJSSP is represented. In the bottom
figure, $o_{1}$ has been assigned to $m_{1}$ and some edges of the graph have
been removed.

In Figure 1, an example of an instance and state transition is shown. Blue
nodes represent jobs, red nodes represent operations, and green nodes
represent machines. The edges are directed, and the values on the edges
represent their processing times. In this instance, there are two jobs,
$j_{1}$ composed of $o_{1}$ and $o_{2}$, and $j_{2}$ composed of $o_{3}$.
$o_{1}$ can be processed on machine $m_{1}$ in five time units, $o_{2}$ on
machine $m_{2}$ in three time units, and $o_{3}$ on $m_{1}$ in eight units and
on $m_{2}$ in five units. The action taken is to select the edge
$(m_{1},j_{1})$, which corresponds to assigning operation $o_{1}$ to machine
$m_{1}$. After this action, several steps are performed. First, the node
$o_{1}$ and all its edges are removed. Then, the edges of node $j_{1}$ are
updated to correspond to its next operation, $o_{2}$. Finally, the edge of
$j_{2}$ connected to $m_{1}$ is updated because performing an operation on
this machine now requires waiting until $o_{1}$ is completed.

Reward function. The objective of the agent is to minimize the makespan. As
proposed in multiple papers (Zhang et al., 2020a; Lei et al., 2022), one way
to minimize the total makespan is by defining the reward function at step $t$
as $r(s_{t},a_{t},s_{t+1})=C(s_{t})-C(s_{t+1})$, where $C(s_{t})$ represents
the makespan at step $t$. If the discount factor, $\gamma$ is set to 1, the
accumulated reward becomes
$\sum_{t=0}^{|\mathcal{O}|}r\left(s_{t},a_{t},s_{t+1}\right)=C\left(s_{0}\right)-C\left(s_{|\mathcal{O}|}\right)$.
Assuming that the makespan is zero or a constant in the first step, minimizing
the accumulated reward is equivalent to minimizing $C_{s_{|\mathcal{O}|}}$,
which is the makespan in the final state.

###  4.2 Feature extraction

To extract the state information, which will be used to generate a policy, we
employed a GNN that utilizes the attention mechanisms of GATv2 explained in
Section 2.3. This approach differs from the one used by Song et al. (2022) and
Wang et al. (2023), where they employ a simpler attention mechanism. As
demonstrated by Brody et al. (2021), it is possible to obtain GNN
architectures with higher and more dynamic attention capacity without
significantly increasing the computational cost. This is crucial for a problem
such as the FJSSP. For each type of node, the initial representations are
denoted as follows: $\boldsymbol{h_{j_{i}}}\in\mathbb{R}^{d_{\mathcal{J}}}$
for jobs, $\boldsymbol{h_{o_{ij}}}\in\mathbb{R}^{d_{\mathcal{O}}}$ for
operations, and $\boldsymbol{h_{m_{i}}}\in\mathbb{R}^{d_{\mathcal{M}}}$ for
machines, where $d_{\mathcal{O}},d_{\mathcal{J}},d_{\mathcal{M}}$ are the
number of initial features of operations, jobs and machines respectively. The
following section describes how the embeddings of the three types of nodes
that compose the state are calculated.

Operation embeddings. Each operation receives information from its successor
and the machines on which it can be processed. Therefore, it is necessary to
calculate two attention weights. For each $o_{ij}\in{\mathcal{O}_{j_{i}}}_{t}$
with $j\leq|{\mathcal{O}_{j_{i}}}_{t}|-1$, the weight $e_{o_{ijj+1}}$ between
the preceding operations $o_{ij}$ and $o_{ij+1}$ is calculated as:

$e_{o_{ijj+1}}={\boldsymbol{a^{\mathcal{O}}}}^{\top}\text{ LeakyReLU
}(\boldsymbol{W_{1}^{\mathcal{O}}}\boldsymbol{h}_{o_{ij}}+\boldsymbol{W_{2}^{\mathcal{O}}}\boldsymbol{h}_{o_{ij+1}})$
(7)

where
${\boldsymbol{a^{\mathcal{O}}}}^{\top}\in\mathbb{R}^{2d_{\mathcal{O}}^{\prime}}$
and
$\boldsymbol{W_{1}^{\mathcal{O}}},\boldsymbol{W_{2}^{\mathcal{O}}}\in\mathcal{R}^{d_{\mathcal{O}}^{\prime}\times
d_{\mathcal{O}}}$ are learned linear transformations for all
$j\leq|{\mathcal{O}_{j_{i}}}_{t}|-1$ and $d_{\mathcal{O}}^{\prime}$ is the
size of the dimension of the embedding of the operation. We also calculate the
attention coefficient between a operation with itself, $e_{o_{ijj}}$. To
calculate the attention between the nodes of machines and operations,
considering that these edges also have features, we use $e_{om_{ijk}}$ to
represent the edge connecting operation $o_{ij}$ and machine $m_{k}$, and
$\boldsymbol{h_{p_{ijk}}}\in\mathbb{R}^{d_{\mathcal{OM}}}$ as the initial
representation of the features of that edge. The weight between nodes $o_{ij}$
and $m_{k}$ is calculated as follows:

$e_{om_{ijk}}={\boldsymbol{a^{\mathcal{OM}}}}^{\top}\text{ LeakyReLU
}(\boldsymbol{W_{1}^{\mathcal{OM}}}\boldsymbol{h}_{o_{ij}}+\boldsymbol{W_{2}^{\mathcal{OM}}}\boldsymbol{h}_{m_{k}}+\boldsymbol{W_{3}^{\mathcal{OM}}}\boldsymbol{h_{p_{ijk}}})$
(8)

where
${\boldsymbol{a^{\mathcal{OM}}}}^{\top}\in\mathbb{R}^{3d_{\mathcal{OM}}^{\prime}}$
and
$\boldsymbol{W_{1}^{\mathcal{OM}}}\in\mathcal{R}^{d_{\mathcal{OM}}^{\prime}\times
d_{\mathcal{O}}},\boldsymbol{W_{2}^{\mathcal{OM}}}\in\mathcal{R}^{d_{\mathcal{OM}}^{\prime}\times
d_{\mathcal{M}}},\boldsymbol{W_{3}^{\mathcal{OM}}}\in\mathcal{R}^{d_{\mathcal{OM}}^{\prime}\times
d_{\mathcal{OM}}}$. After normalizing the attention weights that the operation
has with its successor and itself, and the weight between the operation and
the machines where it can be processed using a softmax function, the
normalized attention coefficients $\alpha_{o_{ijj}}$, $\alpha_{o_{ijj+1}}$,
and $\alpha_{om_{ijk}}$ are obtained, and the new representation of the
operation $\boldsymbol{h_{o_{ij}}}^{\prime}$ is calculated as follows:

$\boldsymbol{h_{o_{ij}}}^{\prime}=\text{ELU}(\alpha_{o_{ijj}}\boldsymbol{W_{2}^{\mathcal{O}}}\boldsymbol{h_{o_{ij}}}+\alpha_{o_{ijj+1}}\boldsymbol{W_{2}^{\mathcal{O}}}\boldsymbol{h_{o_{ij+1}}}+\sum_{m_{k}\in\mathcal{M}_{o_{ij}}}\alpha_{om_{ijk}}(\boldsymbol{W_{2}^{\mathcal{OM}}}\boldsymbol{h_{m_{k}}}+\boldsymbol{W_{3}^{\mathcal{OM}}}\boldsymbol{h_{p_{ijk}}}))$
(9)

Machine embeddings. Machines have two connections from which they will receive
information, operations and other machines. The purpose of machines
communicating with each other is to be able to understand their state and
prioritize their tasks accordingly. The first weight is calculated similarly
to $e_{om_{ijk}}$, and it is denoted as $e_{mo_{ijk}}$. To calculate the
attention that a machine should assign to other machines, for all
$m_{i},m_{j}\in\mathcal{M}$, we compute:

$e_{m_{ij}}={\boldsymbol{a^{\mathcal{M}}}}^{\top}\text{ LeakyReLU
}(\boldsymbol{W_{1}^{\mathcal{M}}}\boldsymbol{h_{m_{i}}}+\boldsymbol{W_{2}^{\mathcal{M}}}\boldsymbol{h_{m_{j}}})$
(10)

where
${\boldsymbol{a^{\mathcal{M}}}}^{\top}\in\mathbb{R}^{2d_{\mathcal{M}}^{\prime}}$
and
$\boldsymbol{W_{1}^{\mathcal{M}}},\boldsymbol{W_{2}^{\mathcal{M}}}\in\mathcal{R}^{d_{\mathcal{M}}^{\prime}\times
d_{\mathcal{M}}}$. After normalizing the attention coefficients,
$\boldsymbol{h_{m_{i}}}^{\prime}$ is calculated as:

$\boldsymbol{h_{m_{i}}}^{\prime}=\text{ELU}(\sum_{m_{k}\in\mathcal{M}}\alpha_{m_{ik}}\boldsymbol{W_{2}^{\mathcal{M}}}\boldsymbol{h_{m_{k}}}+\sum_{o_{ij}\in\mathcal{O}_{m_{i}}}\alpha_{mo_{ijk}}(\boldsymbol{W_{2}^{\mathcal{MO}}}\boldsymbol{h_{o_{ij}}}+\boldsymbol{W_{3}^{\mathcal{MO}}}\boldsymbol{h_{p_{ijk}}}))$
(11)

where
$\boldsymbol{W_{2}^{\mathcal{MO}}}\in\mathcal{R}^{d_{\mathcal{O}}^{\prime}\times
d_{\mathcal{O}}},\boldsymbol{W_{3}^{\mathcal{MO}}}\in\mathcal{R}^{d_{\mathcal{MO}}^{\prime}\times
d_{\mathcal{MO}}}$.

Job embeddings. The calculation of the embeddings for jobs has three sources
of information: the information coming from the operations, the machines, and
the other jobs. We define ${\mathcal{MJ}_{j_{i}}}_{t}$ as the set of machines
that can be assigned the job $j_{i}$ at a step $t$ and
$\boldsymbol{h_{pg_{ki}}}\in\mathbb{R}^{d_{\mathcal{MJ}}}$ as the initial
representation of the features of that edge. In summary, by calculating the
attention coefficients in a similar manner to the previously described
emeddings, $\boldsymbol{h}_{j_{i}}^{\prime}$ is calculated as

$\boldsymbol{h_{j_{i}}}^{\prime}=\text{ELU}(\sum_{j_{j}\in\mathcal{J}}\alpha_{j_{ij}}\boldsymbol{W_{2}^{\mathcal{J}}}\boldsymbol{h_{j_{j}}}+\sum_{o_{ij}\in{\mathcal{O}_{j_{i}}}_{t}}\alpha_{jo_{ij}}\boldsymbol{W_{2}^{\mathcal{JO}}}\boldsymbol{h_{o_{ij}}}+\sum_{m_{k}\in{\mathcal{MJ}_{j_{i}}}_{t}}\alpha_{mj_{ki}}(\boldsymbol{W_{2}^{\mathcal{MJ}}}\boldsymbol{h_{m_{k}}}+\boldsymbol{W_{3}^{\mathcal{MJ}}}\boldsymbol{h_{pg_{ki}}}))$
(12)

where
$\boldsymbol{W_{2}^{\mathcal{J}}}\in\mathcal{R}^{d_{\mathcal{J}}^{\prime}\times
d_{\mathcal{J}}},\boldsymbol{W_{2}^{\mathcal{JO}}}\in\mathcal{R}^{d_{\mathcal{O}}^{\prime}\times
d_{\mathcal{O}}},\boldsymbol{W_{2}^{\mathcal{MJ}}}\in\mathcal{R}^{d_{\mathcal{M}}^{\prime}\times
d_{\mathcal{M}}},\boldsymbol{W_{3}^{\mathcal{MJ}}}\in\mathcal{R}^{d_{\mathcal{MJ}}^{\prime}\times
d_{\mathcal{MJ}}}$.

The process just described details how to perform the calculation for the
first layer of the GNN. However, this process can be repeated _L_ times,
generating the embeddings ${\boldsymbol{h_{j_{i}}}^{\prime}}^{(L)}$,
${\boldsymbol{h_{o_{ij}}}^{\prime}}^{(L)}$,
${\boldsymbol{h_{m_{k}}}^{\prime}}^{(L)}$, where each layer uses the
embeddings of jobs, operations, and machines from the previous layer, and all
layers make use of $\boldsymbol{h_{pg_{ki}}}$ and $\boldsymbol{h_{p_{ijk}}}$,
i.e., the original edge features.

###  4.3 Generating a policy

To generate a policy, we employed the well-known PPO algorithm (Schulman et
al., 2017). This involves generating a policy (actor)
$\pi_{\theta}(a_{t}|s_{t})$ and a value network (critic) $v_{\theta}(s_{t})$.
Both networks use the embeddings calculated by the GNN, although in different
ways. For the policy, we make predictions using the embeddings of the edges
$\mathcal{MJ}_{t}$. For each $(m_{k},j_{i})\in\mathcal{MJ}_{t}$, we
concatenate the embeddings of the machine
${\boldsymbol{h_{m_{k}}}^{\prime}}^{(L)}$, the job
${\boldsymbol{h_{j_{i}}}^{\prime}}^{(L)}$, and the edge feature
$\boldsymbol{h_{pg_{ki}}}$, and then apply $\text{MLP}_{1\theta}$, a mask,
that will be presented in the next section, and a softmax function to create a
probability distribution over all possible pairs of jobs and machines. For the
critic, we apply another $\text{MLP}_{2\theta}$ to all job embeddings, with
hyperbolic tangent as the activation function, and perform mean pooling to
obtain a single value, as described in the following formula:

$v_{\theta}(s_{t})=\frac{1}{|\mathcal{J}_{t}|}\sum_{j_{i}\in\mathcal{J}_{t}}\text{MLP}_{2\theta}({\boldsymbol{h_{j_{i}}}^{\prime}}^{(L)})$
(13)

To train both networks, we run the policy for $n_{eps}$ episodes, where in
each episode, the agent interacts with the environment by sampling from the
probability distribution $\pi_{\theta}(a_{t}|s_{t})$, and the experiences are
stored in a buffer. A schematic representation of the architecture of our
approach is presented in Figure 2. Every $n_{t}$ episodes, we use these
experiences to train the networks for $K$ epochs in batches of size $b$.
Additionally, every $n_{g}$ episodes, a new set of synthetic instances of
dimension $n_{ins}$ is generated for the agent to interact with.

Figure 2: The network architecture of our approach.

###  4.4 More Robust Inference

In some preliminary experiments, we observed that the results obtained from
the generated models were not very robust. Although two models may have
similar average results on a benchmark, there can be significant differences
between instances within the benchmark. To address this issue, we propose two
methods to make our approach more robust.

Dispatching Rules to Limit the Action Space. DRs have been widely used in
FJSSP, but their results can be improved, as shown in the next section.
However, they serve as a good baseline for the agent’s decision-making by
limiting the action space. In this research, we consider two simple heuristics
to limit the action space: allowing only the selection of the $k$ operations
that can start earlier (if there is a tie, all operations with the same value
are allowed) or allowing the $k$ operations that finish earlier. By
constraining the action space in this way, we only allow those actions that
typically lead to an optimal makespan, as it prevents selecting actions that
create large gaps in the scheduling. The value $k$ restricts the degree of
freedom the agent has in making decisions. We choose these simple heuristics
to avoid excessive computational costs that could increase inference time.
Moreover, limiting the action space accelerates the model training by reducing
the state space that the agent explores.

Generating a Diverse Set of Scheduling Policies. We propose a method for
generating a diverse set of scheduling policies using Bayesian optimization
(BO) with the aim of achieving robust results in inference, as a single model
may not be robust. First, we define a search space for the BO algorithm,
denoted as $\mathcal{D}$, where we determine which variables to optimize.
These variables will be related to the structure of the GNN, the PPO
algorithm, the DRs used to limit the action space, or the parameters used to
generate the instances. Next, we generate a validation set $\mathcal{V}$ of
size $n_{\mathcal{V}}$ and obtain the optimal (or suboptimal) makespan of each
instance using exact methods, stored in the set $\mathcal{VM}$, though the
function
$\texttt{get\\_validation\\_makespans}(\mathcal{V})$.$n_{\mathcal{V}}$ should
not be very large to avoid excessive computational cost in the calculation of
makespans. For the validation and training of the models, small to medium-
sized instances are used, which can be solved without excessive computational
cost. We then sample the domain $\mathcal{D}$ using the function
$\texttt{sample}(\mathcal{D})$ to obtain hyperparameters $d$. Using $d$, we
train a policy using the method described in the previous section, implemented
in the function $\texttt{generate\\_policy}(d)$. The generated policy is
evaluated on $\mathcal{V}$, and the makespan it produces for each instance is
used to calculate the gap compared to the optimal makespan in $\mathcal{VM}$.
For a given instance $v\in\mathcal{V}$, the gap for policy $p$ is computed as
$g_{p_{v}}=\frac{(p_{v}-\mathcal{VM}_{v})}{\mathcal{VM}_{v}}$, where $p_{v}$
is the makespan achieved by policy $p$ on instance $v$ and $\mathcal{VM}_{v}$
is the result obtained using exact methods. These gap calculations are
implemented in the function
$\texttt{eval\\_policy}(\mathcal{V},\mathcal{VM},p)$. As the first policy is
generated, all its results are stored in $\mathcal{CG}$ as the current best
gaps achieved. Before starting the BO process, this policy is added to the set
of candidate policies to be selected, denoted as $\mathcal{C}$.

The BO process is executed for $n_{\mathcal{BO}}$ iterations. In each
iteration $j\in\\{1,2,\ldots n_{\mathcal{BO}}\\}$, a set of hyperparameters
$d_{j}$ is chosen, and a new policy $p_{j}$ is generated and evaluated to
obtain the gaps $\mathcal{G}_{p_{j}}$. The maximum difference between the
current gaps and those generated by $p_{j}$ is calculated. This difference
serves as the objective function for the BO algorithm, encouraging the
generation of policies that achieve good results on validation instances with
more room for improvement. If this difference is positive, it means that the
new policy has obtained a better result on at least one instance. Therefore,
this new policy is added to the set of candidate policies, and the current
gaps are updated.

Once a set of diverse scheduling policies is generated, a subset of these
policies is selected to run in parallel during inference. The number of
policies to parallelize will be determined based on the available computing
capacity. To accomplish this, the K-Nearest Neighbors algorithm is applied to
the validation set results for each policy to group in $n_{\mathcal{P}}$
clusters the policies based on similar outcomes in the validation set. The
centroid of each cluster is calculated, and the closest policy, using the
euclidian distance, within each cluster is selected as a final candidate,
using the function $\texttt{cluster\\_center}(\mathcal{K},\mathcal{C})$, and
generating the set of policies $\mathcal{P}$ that will be parallelized during
inference. This approach ensures that similar policies are not executed in
inference. The algorithm for creating a diverse set of scheduling policies
(DSP) is shown in Algorithm 1.

Algorithm 1 Algorithm to generate a set of diverse scheduling policies.

Input: BO domain $\mathcal{D}$, BO iterations $n_{BO}$, number of policies
$n_{\mathcal{P}}$, validation set size $n_{\mathcal{V}}$
Output: Policies $\mathcal{P}$

1:function DSP($\mathcal{D}$,$n_{\mathcal{P}}$, $n_{B}$, $n_{\mathcal{V}}$)

2: Generate a validation set $\mathcal{V}$

3: $\mathcal{VM}$ $\leftarrow$
$\texttt{get\\_validation\\_makespans}(\mathcal{V})$

4: $d\leftarrow$ $\texttt{sample}(\mathcal{D})$

5: $p\leftarrow$ $\texttt{generate\\_policy}(d)$

6: $\mathcal{CG}$ $\leftarrow$ $\texttt{eval\\_policy}(\mathcal{V}$,
$\mathcal{VM}$, p)

7: $\mathcal{C}$ $\leftarrow$ $\\{p\\}$

8: for $j=1$ to $n_{BO}$ do

9: $d_{j}$ $\leftarrow$ sample($\mathcal{D}$)

10: $p_{j}$ $\leftarrow$ generate_policy($d_{j}$)

11: $\mathcal{G}_{p_{j}}$ $\leftarrow$ eval_policy($\mathcal{V}$,
$\mathcal{VM}$, $p_{j}$ )

12: $g_{j}\leftarrow\text{max}(\mathcal{CG}-\mathcal{G}_{p_{j}})$

13: if $g_{j}>0$ then

14: $\mathcal{C}$ $\leftarrow$ $\mathcal{P}\cup\\{p_{j}\\}$

15: $\mathcal{CG}_{v}$ $\leftarrow$
$\min\limits_{v\in\mathcal{V}}(\mathcal{CG}_{v},\mathcal{G}_{p_{j_{v}}})$
Update surrogate model with $(d_{j},g_{j})$

16: $\mathcal{K}\leftarrow\texttt{KNN}(\mathcal{C},n_{\mathcal{P}})$

17: $\mathcal{P}\leftarrow\texttt{cluster\\_center}(\mathcal{K},\mathcal{C})$
return $\mathcal{P}$

Figure 3: The gap of two instances of the validation set of all the generated
candidate policies is shown, colored by the cluster they belong to. All the
policies have been divided into 9 clusters. It can be observed that the yellow
and purple clusters achieve the best results for these two instances.

##  5 Experimental results

In this section, we present the experiments conducted to evaluate the validity
of the proposed methodology, which we refer to as Diverse Scheduling Policies
(DSP), we study its capability to outperform DRs and its competitiveness with
other state-of-the-art approaches. For this purpose, comparisons have been
made with synthetically generated instances and publicly available FJSSP
benchmarks 111We obtained the benchmarks from https://openhsu.ub.hsu-
hh.de/handle/10.24405/436., measuring both the quality of the solutions and
the time required. The method was implemented using Python 3.10, PyTorch
Geometric (Fey and Lenssen, 2019) was utilized for GNNs, Optuna (Akiba et al.,
2019) was used for hyperparameter optimization and OR-Tools for the exact
method 222 We obtained the implementation from https://github.com/google/or-
tools/blob/stable/examples/python/flexible_job_shop_sat.py.. The hardware is a
machine with AMD Ryzen 5 5600X processor and an NVIDIA GeForce RTX 3070 Ti.
The code and the synthetically generated test instances will be made publicly
available upon acceptance of the paper.

###  5.1 Generation of Synthetic Instances

To train the models, synthetic instances must be generated, and for this
purpose, the instance generation method proposed in Brandimarte (1993) and
Song et al. (2022) has been adapted. The generation of synthetic instances is
governed by various parameters that determine their size and structure. These
parameters include the minimum and maximum number of jobs ($j_{\text{min}}$
and $j_{\text{max}}$), the minimum and maximum number of machines
($m_{\text{min}}$ and $m_{\text{max}}$), the minimum and maximum number of
operations per job ($o_{\text{min}}$ and $o_{\text{max}}$), the maximum number
of options each operation can have ($op_{\text{max}}$, with a minimum fixed at
four), the maximum processing time for an operation $\overline{p}$ (with a
minimum defined as one), and the deviation from the mean $d$. To generate an
instance, the uniform distribution within the defined ranges is sampled. The
processing time on machine $m_{k}$ is calculated by sampling from the uniform
distribution $\mathcal{U}(\overline{p_{ij}}(1-d),\overline{p_{ij}}(1+d))$.
Training on large instances can be time-consuming, so smaller instances are
used for that purpose, while larger instances are used for testing.

###  5.2 Configuration

For the experiments, the validation set size, $n_{\mathcal{V}}$, has been set
to 100, the number of BO iterations $n_{\text{BO}}$ has been set to 200 and
the synthetic instances have been generated by randomly selecting from the
number of jobs, machines and operation ranges defined below. Tree-structured
Parzen Estimator has been used as the sampler. The size of the policy set
$n_{\mathcal{P}}$ has been set to 6, to enable efficient parallelization of
the models. Regarding the domain of BO, the hyperparameters can be divided
into categories related to: the structure of the GNN, the training mode of the
PPO algorithm, the size of the generated instances, and the heuristic and
degrees of freedom for decision-making. The following ranges have been
considered:

  * •

Related to the GNN: the number of layers $L\in[1,3]$ and the number of hidden
channels $d^{\prime}\in\\{32,64,128\\}$. For simplicity, the dimension of the
node embeddings has been considered the same for all nodes. $MLP_{1\theta}$
and $MLP_{2\theta}$ of the actor and the critic have a single layer of the
same dimension.

  * •

The defined ranges for generating instance sets in each iteration of BO were:
$j_{\text{min}}\in[5,9]$, $j_{\text{max}}\in[10,15]$,
$m_{\text{min}}\in[4,7]$, $m_{\text{max}}\in[9,13]$,
$op_{\text{max}}\in[5,9]$, $p\in[8,25]$. $d$ has been set to $0.2$.

  * •

Related to network training: the maximum number of episodes during training
$n_{eps}\in[10000,20000]$, learning rate $lr\in[0.00001,0.001]$, batch size
$b\in\\{64,128,256\\}$, model training frequency $n_{t}\in[10,100]$, and
frequency of creating new instances $n_{ins}\in[1,10]$.

  * •

The heuristic used by the environment, as well as the possible options
$k\in[1,3]$.

Allowing the BO algorithm to have many hyperparameter options aims to obtain
diverse models that can be effective on different subsets of instances. For
the PPO algorithm, the discount factor $\gamma$ was fixed to 1, the clipping
parameter $\epsilon$ to 0.2, and the coefficients for policy, value, and
entropy in the loss function were set to 1, 0.5, and 0.01, respectively. The
number of epochs per episode was set to 3, and Adam optimizer (Kingma and Ba,
2014) was used.

###  5.3 Results in synthetically generated instances

First, we have validated our approach by comparing it with various well-known
DRs. Additionally, optimal or near-optimal solutions have been generated using
OR-Tools, running it for 1800 seconds, using the CP-SAT solver. OR-Tools has
been used in multiple papers for this purpose (Song et al., 2022; Wang et al.,
2023; Tassel et al., 2021), because it is one of the best open-source tools
for solving scheduling problems (Da Col and Teppan, 2019). The objective is to
study if our approach outperforms the DRs, as their execution times are
similar to our method, and to measure the gap between the solution obtained by
our approach and the optimal solution. For this purpose, four sets of
instances larger than the instances used during the training phase have been
generated. The metrics used to measure the results have been the gap between
the method and the solution proposed by OR-Tools (the lower the value, the
better, and a negative value indicates that it is better than OR-Tools), and
the execution time.

Regarding the heuristic rules that have been used, they have been shown to be
effective in the literature (Jun et al., 2019) and have recently been used to
compare them with other deep learning-based methods (Lei et al., 2022; Chen
and Matis, 2013). To define a heuristic for solving the FJSSP, two rules need
to be defined: one for job selection and sequencing, and another for operation
assignment to machines. Regarding the former, the following rules have been
considered:

  * •

First In First Out (FIFO), which processes the first operation to arrive in
the queue of a machine.

  * •

Most Work Remaining (MWKR), which processes a job with the highest remaining
workload.

  * •

Least Work Remaining (LWKR), which processes a job with the lowest remaining
workload.

  * •

Most Operation Number Remaining (MOPNR), which processes a job with the
highest number of remaining operations.

For the machine selection criteria, the following rules have been considered:

  * •

Shortest Processing Time (SPT), which assigns a machine with the shortest
processing time for an operation.

  * •

Earliest Ending Time (EET), which selects the machine where an operation can
start earliest.

To compare our method, we generated four sets of instances, each consisting of
50 instances, which are larger than the ones used to train our model. The
creation of these sets followed an incremental approach, where they became
progressively more complex. The parameters used for each set are detailed in
Table 1. For the first test set, $\mathcal{S}_{1}$, we doubled the number of
jobs and operations compared to the ones used in the training set. In
$\mathcal{S}_{2}$, we increased the number of jobs and the processing times,
allowing for more extensive and unseen ranges during the training phase. In
$\mathcal{S}_{3}$, we further increased the number of jobs and, finally, in
$\mathcal{S}_{4}$, we increased the number of possible machines.

Table 1: Summary of the size of the test instances Name | Range jobs | Range machines | Ops per job | Processing time
---|---|---|---|---
$\mathcal{S}_{1}$ | $\mathcal{U}(20,30)$ | $\mathcal{U}(8,16)$ | $\mathcal{U}(4,8)$ | $\mathcal{U}(1,7)$
$\mathcal{S}_{2}$ | $\mathcal{U}(40,50)$ | $\mathcal{U}(8,16)$ | $\mathcal{U}(4,8)$ | $\mathcal{U}(5,20)$
$\mathcal{S}_{3}$ | $\mathcal{U}(50,60)$ | $\mathcal{U}(8,16)$ | $\mathcal{U}(4,8)$ | $\mathcal{U}(5,20)$
$\mathcal{S}_{4}$ | $\mathcal{U}(50,60)$ | $\mathcal{U}(16,32)$ | $\mathcal{U}(8,16)$ | $\mathcal{U}(5,20)$

In Table 2, the average gaps and execution times obtained by the eight DRs and
our approach are shown for the synthetically generated test sets. As can be
observed, although our approach takes slightly more time, it consistently
outperforms the eight DRs. Furthermore, the performance of our proposed method
is more robust than that of the DRs as the problem scales increases. More
interestingly, the results obtained for large instances are better than those
obtained by OR-Tools (as the gap achieved is negative), although the execution
time of our methods is much less. This demonstrates that our approach has been
able to generalize the knowledge learned from small instances to larger
instances with more machines, jobs, operations, and average processing times.

Table 2: Summary of the result in test instances of all DRs and our approach
(DSP) with respect to OR-Tools.

Name | Obj. | FIFO | LWKR | MWKR | MOPNR | DSP | OR-Tools
---|---|---|---|---|---|---|---
EET | SPT | EET | SPT | EET | SPT | EET | SPT
$\mathcal{S}_{1}$ | Gap(%) | 27.28 | 61.30 | 92.56 | 86.03 | 26.05 | 52.80 | 28.08 | 53.36 | 10.84 | -
| Time(s) | 0.11 | 0.13 | 0.13 | 0.13 | 0.13 | 0.14 | 0.14 | 0.13 | 8.54 | 1591.63
$\mathcal{S}_{2}$ | Gap(%) | 13.40 | 33.25 | 58.24 | 48.07 | 12.07 | 30.62 | 13.25 | 32.29 | 6.96 | -
| Time(s) | 0.26 | 0.30 | 0.29 | 0.28 | 0.27 | 0.29 | 0.29 | 0.34 | 32.90 | 1800.38
$\mathcal{S}_{3}$ | Gap(%) | 10.36 | 27.56 | 50.31 | 36.93 | 9.95 | 26.35 | 10.77 | 26.41 | 6.05 | -
| Time(s) | 0.36 | 0.40 | 0.36 | 0.37 | 0.37 | 0.39 | 0.39 | 0.43 | 57.61 | 1800.45
$\mathcal{S}_{4}$ | Gap(%) | 7.43 | 29.99 | 95.92 | 69.29 | 4.43 | 26.82 | 7.85 | 26.11 | -3.06 | -
| Time(s) | 0.94 | 1.08 | 0.93 | 0.98 | 0.97 | 0.95 | 1.04 | 1.08 | 166.56 | 1801.52

###  5.4 Results in publicly available FJSSP benchmarks

To study the generalization capability of the proposed method on real-world
instances, 45 large instances from the Behnke benchmark (Behnke and Geiger,
2012) have been used and, although our main focus is large instances, we have
also used the instances proposed by Hurink et al. (1994) (40 small and medium-
sized instances from the vdata set) to assess that there still is a
competitive alternative for small ones. We have chosen these subsets from
these benchmarks in order to compare our results with the approach presented
by Lei et al. (2022). In addition to comparing it with this approach, we
compare our method with the methods proposed by Wang et al. (2023) and Song et
al. (2022), both also recent DRL methods, using their publicly available open-
source code and using their default settings to train and test the models. To
perform the inference, the sampling strategy has been used, as described in
Song et al. (2022). This strategy involves running the instance 100 times by
sampling from the distribution $\pi_{\theta}$ instead of selecting the maximum
probability, and choosing the smallest makespan value among all runs. This
strategy has been chosen because it yields better results than selecting the
maximum probability from $\pi_{\theta}$, which is called the greedy approach.
In each benchmark, we have also compared the DRL methods with the best DR
among the proposed ones.

Results in vdata benchmark. Table 3 shows the results obtained by the four
methods grouped according to the size of the instances. Since the optimal or
sub-optimal solutions for these instances are available, as published by
Behnke and Geiger (2012), we will use them instead of OR-Tools in this case,
and also because the other approaches also use them. Nevertheless, the results
obtained with OR-Tools are very similar to those obtained from Behnke and
Geiger (2012).

Table 3: Comparison of the mean gap grouped by size obtained by DRs, the methods proposed by Song et al., Wang et al. and Lei et al. and our approach (DSP) in vdata instances. Size | Obj. | DR | Song | Wang | Lei | DSP
---|---|---|---|---|---|---
10x5 | Gap(%) | 41.40 | 6.89 | 3.09 | 8.00 | 4.14
| Time (s.) | 0.03 | 4.79 | 0.72 | - | 2.67
10x10 | Gap(%) | 49.38 | 0.28 | 0.27 | 0.32 | 1.24
| Time (s.) | 0.06 | 9.10 | 1.75 | - | 2.77
15x5 | Gap(%) | 30.34 | 2.99 | 1.13 | 3.24 | 2.10
| Time (s.) | 0.04 | 6.38 | 1.05 | - | 2.38
15x10 | Gap(%) | 63.06 | 7.35 | 5.13 | 11.44 | 4.07
| Time (s.) | 0.10 | 13.19 | 3.47 | - | 4.17
15x15 | Gap(%) | 57.26 | 0.88 | 0.88 | 3.25 | 0.91
| Time (s.) | 0.16 | 24.11 | 7.85 | - | 7.52
20x5 | Gap(%) | 24.32 | 2.97 | 0.94 | 1.85 | 0.88
| Time (s.) | 0.07 | 7.90 | 1.58 | - | 2.92
20x10 | Gap(%) | 50.01 | 5.02 | 1.79 | 4.08 | 1.50
| Time (s.) | 0.14 | 18.14 | 5.29 | - | 6.11
30x10 | Gap(%) | 33.08 | 3.25 | 1.16 | 1.91 | 0.57
| Time (s.) | 0.27 | 30.36 | 11.48 | - | 12.10

As expected, the results of the best DR, which is FIFO for job selection and
EET for machine selection, yield the worst results. Regarding the DRL-based
methods, all three of them achieve a similar average gap, with the best mean
gap achieved by Wang et al. (2023), $1.81\%$. Our DSP method follows closely
with the second-best results, featuring an average gap of $1.94\%$ nearly
identical to the previous. This is followed by the method of Song et al.
(2022) at $3.72\%$ and Lei et al. (2022) at $4.27\%$. We emphasize that, as
seen in Table 3, the average results for larger instances are better with our
method, suggesting that our approach is more robust with this kind of
instances. We speculate that the better outcomes in smaller instances by Wang
et al. (2023) may be attributed to their state representation featuring more
features and their GNN utilizing a dual attention heads architecture. As for
the average execution times, which are shown in Table 3, we compared the
proposed methods except for the one by Lei et al. (2022) as they did not
provide their exact times. As can be seen in Table 3, the DRL-based methods
achieve similar execution times.

Results in Behnke benchmark. Table 4 shows the results grouped according to
the size of the instances from the benchmark obtained by the best DR, the four
DRL methods, and OR-Tools, which was used to obtain pseudo-optimal results by
running it for 1800 seconds. In A, the results for each of the instances are
displayed. It is important to show the results obtained by the solver on each
instance because, for example, in comparison, Lei et al. (2022) achieve their
results by running the Gurobi solver for 3600 seconds and obtain considerably
worse results.

In this case, our approach outperforms Wang et al. (2023) in 40 out of 45
instances, Lei et al. (2022) in 34 and Song et al. (2022) and the DRs in all
instances. Furthermore, our approach performs better on more challenging
instances. In terms of the average gap, notable differences exist among the
DRL methods in this benchmark. The method proposed by Song et al. (2022)
achieves the worst results among the DRL methods, $38.48\%$, close to the DRs,
$55.65\%$, indicating instability issues for large instances. The other two
DRL approaches are more robust, although our method, $2.91\%$, is clearly
superior to those proposed by Lee et al. (2021), $7.27\%$, and Wang et al.
(2023), $6.42\%$. Regarding execution time, as shown in Table 4, it can be
observed that as the instance dimension increases, the processing time
significantly especially increases for especially Wang et al. (2023), which is
slower than our approach. This could be because they calculate more features,
particularly within the operation-type nodes, which are more numerous.

Table 4: Comparison of the mean gap grouped by size obtained by DRs, the methods proposed by Song et al., Wang et al. and Lei et al. and our approach (DSP) in Behnke instances. Size | Obj. | DR | Song | Wang | Lei | DSP | OR-Tools
---|---|---|---|---|---|---|---
20x20 | Gap(%) | 74.07 | 42.65 | 17.55 | 11.93 | 10.53 | -
| Time (s.) | 0.08 | 10.29 | 2.78 | - | 3.60 | 1481.77
20x40 | Gap(%) | 69.09 | 50.57 | 13.86 | 16.19 | 12.53 | -
| Time (s.) | 0.09 | 11.39 | 6.43 | - | 6.67 | 607.57
20x60 | Gap(%) | 67.38 | 50.22 | 14.36 | 20.23 | 12.64 | -
| Time (s.) | 0.10 | 14.38 | 12.35 | - | 12.06 | 539.37
50x20 | Gap(%) | 62.98 | 34.44 | 8.11 | 0.01 | 1.94 | -
| Time (s.) | 0.33 | 29.16 | 14.63 | - | 15.26 | 1800.30
50x40 | Gap(%) | 54.14 | 40.59 | 7.15 | 11.88 | 3.27 | -
| Time (s.) | 0.36 | 30.67 | 37.505 | - | 26.07 | 1800.57
50x60 | Gap(%) | 53.19 | 44.09 | 6.13 | 12.27 | 3.78 | -
| Time (s.) | 0.39 | 37.82 | 65.824 | - | 41.53 | 1800.83
100x20 | Gap(%) | 44.03 | 23.70 | -4.73 | -11.02 | -9.36 | -
| Time (s.) | 1.12 | 68.88 | 54.794 | - | 61.05 | 1800.62
100x40 | Gap(%) | 38.63 | 28.21 | -2.69 | 2.03 | -4.95 | -
| Time (s.) | 1.27 | 82.27 | 134.12 | - | 93.44 | 1801.15
100x60 | Gap(%) | 37.37 | 31.85 | -1.92 | 1.92 | -4.21 | -
| Time (s.) | 1.28 | 102.07 | 235.10 | - | 130.15 | 1801.69

##  6 Conclusions and future work

This paper has introduced a novel method based on DRL and GNNs to dynamically
and robustly solve the FJSSP. To achieve this, we first proposed a new
representation of states, action space, and transition function, where nodes
and edges are progressively removed to make our method more efficient.
Secondly, we designed a message passing scheme between nodes and a more
expressive attention mechanism, not employed by other methods in the
literature. Moreover, we presented two new approaches to enhance inference,
making it more robust. This was achieved by creating a set of diverse policies
that can be executed in parallel and by constraining the action space using
DRs, resulting in faster model training and more robust inference.

To validate our approach, we generated new instances larger than those used in
the training phase, demonstrating the generalization capability of our
architecture. The results revealed that our method outperforms eight DRs on
large instances. Furthermore, we compared our method with public benchmarks
and recent DRL-based approaches. In small to medium-sized instances, our
method outperformed two of the three methods and was slightly worse than the
third. In large instances, it significantly outperformed all three methods.
Notably, our method even outperformed OR-Tools on instances with a higher
number of jobs and machines, all while maintaining significantly lower
execution times.

For future work, we plan to extend this method to handle additional
complexities in the FJSSP, such as new job insertions or setup times.
Moreover, we aim to apply our approach to the low-carbon FJSSP, which
considers both the makespan and energy consumption of machines. Additionally,
it will be interesting to study the possibility to extend this approach to
other combinatorial optimization and operations research problems, such as hub
location problems or knapsack problems.

## Acknowledgements

This work was partially financed by the Basque Government through their
Elkartek program (SONETO project, ref. KK-2023/00038). R. Santana thanks the
Misiones Euskampus 2.0 programme for the financial help received through
Euskampus Fundazioa. Partial support by the Research Groups 2022-2024
(IT1504-22) and the Elkartek Program (KK-2020/00049, KK-2022/00106, SIGZE,
K-2021/00065) from the Basque Government, and the PID2019-104966GB-I00 and
PID2022-137442NB-I00 research projects from the Spanish Ministry of Science is
acknowledged.

## References

  * Akiba et al. (2019) Akiba, T., Sano, S., Yanase, T., Ohta, T., Koyama, M., 2019.  Optuna: A next-generation hyperparameter optimization framework, in: Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 2623–2631.
  * Balas (1969) Balas, E., 1969.  Machine sequencing via disjunctive graphs: an implicit enumeration algorithm.  Operations Research 17, 941–957.
  * Behnke and Geiger (2012) Behnke, D., Geiger, M.J., 2012.  Test instances for the flexible job shop scheduling problem with work centers.  Arbeitspapier/Research Paper/Helmut-Schmidt-Universität, Lehrstuhl für Betriebswirtschaftslehre, insbes. Logistik-Management .
  * Brandimarte (1993) Brandimarte, P., 1993.  Routing and scheduling in a flexible job shop by tabu search.  Annals of Operations Research 41, 157–183.
  * Brody et al. (2021) Brody, S., Alon, U., Yahav, E., 2021.  How attentive are graph attention networks?  arXiv preprint arXiv:2105.14491 .
  * Çaliş and Bulkan (2015) Çaliş, B., Bulkan, S., 2015.  A research survey: review of AI solution strategies of job shop scheduling problem.  Journal of Intelligent Manufacturing 26, 961–973.
  * Calleja and Pastor (2014) Calleja, G., Pastor, R., 2014.  A dispatching algorithm for flexible job-shop scheduling with transfer batches: an industrial application.  Production Planning & Control 25, 93–109.
  * Cappart et al. (2021) Cappart, Q., Moisan, T., Rousseau, L.M., Prémont-Schwarz, I., Cire, A.A., 2021.  Combining reinforcement learning and constraint programming for combinatorial optimization, in: Proceedings of the AAAI Conference on Artificial Intelligence, pp. 3677–3687.
  * Chen and Matis (2013) Chen, B., Matis, T.I., 2013.  A flexible dispatching rule for minimizing tardiness in job shop scheduling.  International Journal of Production Economics 141, 360–365.
  * Chen et al. (2020) Chen, R., Yang, B., Li, S., Wang, S., 2020.  A self-learning genetic algorithm based on reinforcement learning for flexible job-shop scheduling problem.  Computers & Industrial Engineering 149, 106778.
  * Da Col and Teppan (2019) Da Col, G., Teppan, E.C., 2019.  Industrial size job shop scheduling tackled by present day cp solvers, in: Principles and Practice of Constraint Programming: 25th International Conference, CP 2019, Stamford, CT, USA, September 30–October 4, 2019, Proceedings 25, Springer. pp. 144–160.
  * Defersha et al. (2022) Defersha, F.M., Obimuyiwa, D., Yimer, A.D., 2022.  Mathematical model and simulated annealing algorithm for setup operator constrained flexible job shop scheduling problem.  Computers & Industrial Engineering 171, 108487.
  * Du et al. (2022) Du, Y., Li, J., Li, C., Duan, P., 2022.  A reinforcement learning approach for flexible job shop scheduling problem with crane transportation and setup times.  IEEE Transactions on Neural Networks and Learning Systems .
  * Fey and Lenssen (2019) Fey, M., Lenssen, J.E., 2019.  Fast graph representation learning with PyTorch Geometric, in: ICLR Workshop on Representation Learning on Graphs and Manifolds.
  * Garey et al. (1976) Garey, M.R., Johnson, D.S., Sethi, R., 1976.  The complexity of flowshop and jobshop scheduling.  Mathematics of Operations Research 1, 117–129.
  * Gupta and Sivakumar (2006) Gupta, A.K., Sivakumar, A.I., 2006.  Job shop scheduling techniques in semiconductor manufacturing.  The International Journal of Advanced Manufacturing Technology 27, 1163–1169.
  * Hajibabaei and Behnamian (2021) Hajibabaei, M., Behnamian, J., 2021.  Flexible job-shop scheduling problem with unrelated parallel machines and resources-dependent processing times: A tabu search algorithm.  International Journal of Management Science and Engineering Management 16, 242–253.
  * Hurink et al. (1994) Hurink, J., Jurisch, B., Thole, M., 1994.  Tabu search for the job-shop scheduling problem with multi-purpose machines.  Operations-Research-Spektrum 15, 205–215.
  * Joe and Lau (2020) Joe, W., Lau, H.C., 2020.  Deep reinforcement learning approach to solve dynamic vehicle routing problem with stochastic customers, in: Proceedings of the International Conference on Automated Planning and Scheduling, pp. 394–402.
  * Jun et al. (2019) Jun, S., Lee, S., Chun, H., 2019.  Learning dispatching rules using random forest in flexible job shop scheduling problems.  International Journal of Production Research 57, 3290–3310.
  * Kim et al. (2022) Kim, M.S., Oh, S.C., Chang, E.H., Lee, S., Wells, J.W., Arinez, J., Jang, Y.J., 2022.  A dynamic programming-based heuristic algorithm for a flexible job shop scheduling problem of a matrix system in automotive industry, in: 2022 IEEE 18th International Conference on Automation Science and Engineering (CASE), IEEE. pp. 777–782.
  * Kingma and Ba (2014) Kingma, D.P., Ba, J., 2014.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980 .
  * Lee et al. (2021) Lee, K., Laskin, M., Srinivas, A., Abbeel, P., 2021.  Sunrise: A simple unified framework for ensemble learning in deep reinforcement learning, in: International Conference on Machine Learning, PMLR. pp. 6131–6141.
  * Lei et al. (2022) Lei, K., Guo, P., Zhao, W., Wang, Y., Qian, L., Meng, X., Tang, L., 2022.  A multi-action deep reinforcement learning framework for flexible job-shop scheduling problem.  Expert Systems with Applications 205, 117796.
  * Lin et al. (2019) Lin, C.C., Deng, D.J., Chih, Y.L., Chiu, H.T., 2019.  Smart manufacturing scheduling with edge computing using multiclass deep q network.  IEEE Transactions on Industrial Informatics 15, 4276–4284.
  * Liu et al. (2020) Liu, C.L., Chang, C.C., Tseng, C.J., 2020.  Actor-critic deep reinforcement learning for solving job shop scheduling problems.  IEEE Access 8, 71752–71762.
  * Luo (2020) Luo, S., 2020.  Dynamic scheduling for flexible job shop with new job insertions by deep reinforcement learning.  Applied Soft Computing 91, 106208.
  * Luo et al. (2021) Luo, S., Zhang, L., Fan, Y., 2021.  Dynamic multi-objective scheduling for flexible job shop by deep reinforcement learning.  Computers & Industrial Engineering 159, 107489.
  * Manne (1960) Manne, A.S., 1960.  On the job-shop scheduling problem.  Operations Research 8, 219–223.
  * Mazyavkina et al. (2021) Mazyavkina, N., Sviridov, S., Ivanov, S., Burnaev, E., 2021.  Reinforcement learning for combinatorial optimization: A survey.  Computers & Operations Research 134, 105400.
  * Meng et al. (2020) Meng, L., Zhang, C., Ren, Y., Zhang, B., Lv, C., 2020.  Mixed-integer linear programming and constraint programming formulations for solving distributed flexible job shop scheduling problem.  Computers & Industrial Engineering 142, 106347.
  * Scarselli et al. (2008) Scarselli, F., Gori, M., Tsoi, A.C., Hagenbuchner, M., Monfardini, G., 2008.  The graph neural network model.  IEEE Transactions on Neural Networks 20, 61–80.
  * Schulman et al. (2017) Schulman, J., Wolski, F., Dhariwal, P., Radford, A., Klimov, O., 2017.  Proximal policy optimization algorithms.  arXiv preprint arXiv:1707.06347 .
  * Song et al. (2022) Song, W., Chen, X., Li, Q., Cao, Z., 2022.  Flexible job-shop scheduling via graph neural network and deep reinforcement learning.  IEEE Transactions on Industrial Informatics 19, 1600–1610.
  * Sutton and Barto (2018) Sutton, R.S., Barto, A.G., 2018.  Reinforcement Learning: An introduction.  MIT press.
  * Tassel et al. (2021) Tassel, P., Gebser, M., Schekotihin, K., 2021.  A reinforcement learning environment for job-shop scheduling.  arXiv preprint arXiv:2104.03760 .
  * Tolstaya et al. (2021) Tolstaya, E., Paulos, J., Kumar, V., Ribeiro, A., 2021.  Multi-robot coverage and exploration using spatial graph neural networks, in: 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), IEEE. pp. 8944–8950.
  * Tutumlu and Saraç (2023) Tutumlu, B., Saraç, T., 2023.  A mip model and a hybrid genetic algorithm for flexible job-shop scheduling problem with job-splitting.  Computers & Operations Research 155, 106222.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I., 2017.  Attention is all you need.  Advances in Neural Information Processing Systems 30.
  * Veličković et al. (2017) Veličković, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., Bengio, Y., 2017.  Graph attention networks.  arXiv preprint arXiv:1710.10903 .
  * Wang et al. (2023) Wang, R., Wang, G., Sun, J., Deng, F., Chen, J., 2023.  Flexible job shop scheduling via dual attention network based reinforcement learning.  arXiv preprint arXiv:2305.05119 .
  * Wu et al. (2020) Wu, Z., Pan, S., Chen, F., Long, G., Zhang, C., Philip, S.Y., 2020.  A comprehensive survey on graph neural networks.  IEEE Transactions on Neural Networks and Learning Systems 32, 4–24.
  * Xiong et al. (2022) Xiong, H., Shi, S., Ren, D., Hu, J., 2022.  A survey of job shop scheduling problem: The types and models.  Computers & Operations Research 142, 105731.
  * Zhang et al. (2020a) Zhang, C., Song, W., Cao, Z., Zhang, J., Tan, P.S., Chi, X., 2020a.  Learning to dispatch for job shop scheduling via deep reinforcement learning.  Advances in Neural Information Processing Systems 33, 1621–1632.
  * Zhang et al. (2020b) Zhang, G., Hu, Y., Sun, J., Zhang, W., 2020b.  An improved genetic algorithm for the flexible job shop scheduling problem with multiple time constraints.  Swarm and Evolutionary Computation 54, 100664.
  * Zhang et al. (2021) Zhang, X.M., Liang, L., Liu, L., Tang, M.J., 2021.  Graph neural networks and their current applications in bioinformatics.  Frontiers in Genetics 12, 690049.
  * Zhou et al. (2020) Zhou, J., Cui, G., Hu, S., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C., Sun, M., 2020.  Graph neural networks: A review of methods and applications.  AI Open 1, 57–81.

##  Appendix A Behnke results

In the following table, the results obtained by the different methods and OR-
Tools are shown. The gap relative to OR-Tools is presented in parentheses as a
percentage.

Table 5: Comparison of the best makespan obtained by DRs, the methods proposed
by Song et al., Wang et al. and Lei et al., our approach (DSP) and OR-Tools in
instances.

Index | Size | DR | Song | Wang | Lei | DSP | OR
---|---|---|---|---|---|---|---
06 | 20x20 | 176 (64.84) | 179 (39.84) | 148 (15.62) | 143 (11.72) |  141 (10.16) | 128
07 | 20x20 | 174 (65.89) | 176 (36.43) | 152 (17.83) | 142 (10.08) |  138 (6.98) | 129
08 | 20x20 | 200 (83.33) | 185 (46.83) | 148 (17.46) |  139 (10.32) | 142 (12.7) | 126
09 | 20x20 | 169 (83.2) | 180 (44.0) | 147 (17.6) | 144 (15.2) |  142 (13.6) | 125
10 | 20x20 | 181 (73.08) | 190 (46.15) | 155 (19.23) | 146 (12.31) |  142 (9.23) | 130
11 | 50x20 | 321 (59.53) | 337 (31.13) | 276 (7.39) |  250 (-2.72) | 251 (-2.33) | 257
12 | 50x20 | 282 (57.37) | 336 (33.86) | 267 (6.37) |  247 (-1.59) | 250 (-0.4) | 251
13 | 50x20 | 315 (67.34) | 343 (38.31) | 271 (9.27) |  249 (0.4) | 265 (6.85) | 248
14 | 50x20 | 352 (66.14) | 342 (36.25) | 275 (9.56) |  257 (2.39) | 258 (2.79) | 251
15 | 50x20 | 324 (64.54) | 333 (32.67) | 271 (7.97) |  255 (1.59) | 258 (2.79) | 251
16 | 100x20 | 547 (45.06) | 598 (23.05) | 460 (-5.35) |  437 (-10.08) | 444 (-8.64) | 486
17 | 100x20 | 562 (50.32) | 597 (25.68) | 449 (-5.47) |  430 (-9.47) | 439 (-7.58) | 475
18 | 100x20 | 524 (42.68) | 600 (25.52) | 466 (-2.51) |  428 (-10.46) | 440 (-7.95) | 478
19 | 100x20 | 539 (41.19) | 594 (21.72) | 468 (-4.1) |  423 (-13.32) | 436 (-10.66) | 488
20 | 100x20 | 524 (40.91) | 593 (22.52) | 454 (-6.2) | 427 (-11.78) |  426 (-11.98) | 484
26 | 20x40 | 167 (74.34) | 173 (53.1) | 132 (16.81) | 129 (14.16) |  126 (11.5) | 113
27 | 20x40 | 176 (64.23) | 180 (46.34) |  137 (11.38) | 140 (13.82) | 140 (13.82) | 123
28 | 20x40 | 152 (71.93) | 182 (59.65) | 132 (15.79) | 133 (16.67) |  121 (6.14) | 114
29 | 20x40 | 145 (64.96) | 177 (51.28) |  131 (11.97) | 138 (17.95) | 134 (14.53) | 117
30 | 20x40 | 185 (70.0) | 171 (42.5) |  136 (13.33) | 142 (18.33) | 140 (16.67) | 120
31 | 50x40 | 306 (55.08) | 326 (38.14) | 253 (7.2) | 271 (14.83) |  242 (2.54) | 236
32 | 50x40 | 287 (56.17) | 326 (38.72) | 254 (8.09) | 267 (13.62) |  241 (2.55) | 235
33 | 50x40 | 274 (51.28) | 314 (34.19) | 248 (5.98) | 261 (11.54) |  247 (5.56) | 234
34 | 50x40 | 280 (53.02) | 338 (45.69) | 250 (7.76) | 250 (7.76) |  239 (3.02) | 232
35 | 50x40 | 290 (55.16) | 326 (46.19) | 238 (6.73) | 249 (11.66) |  229 (2.69) | 223
36 | 100x40 | 500 (44.21) | 575 (33.1) | 429 (-0.69) | 442 (2.31) |  414 (-4.17) | 432
37 | 100x40 | 504 (37.19) | 567 (26.28) | 428 (-4.68) | 444 (-1.11) |  423 (-5.79) | 449
38 | 100x40 | 490 (35.51) | 570 (28.09) | 427 (-4.04) | 454 (2.02) |  422 (-5.17) | 445
39 | 100x40 | 498 (39.6) | 574 (28.41) | 432 (-3.36) | 471 (5.37) |  412 (-7.83) | 447
40 | 100x40 | 518 (36.64) | 567 (25.17) | 450 (-0.66) | 460 (1.55) |  445 (-1.77) | 453
46 | 20x60 | 176 (74.56) | 172 (50.88) |  130 (14.04) | 145 (27.19) | 131 (14.91) | 114
47 | 20x60 | 152 (61.54) | 177 (51.28) | 131 (11.97) | 144 (23.08) |  127 (8.55) | 117
48 | 20x60 | 165 (64.0) | 182 (45.6) | 143 (14.4) |  141 (12.8) | 142 (13.6) | 125
49 | 20x60 | 157 (72.57) | 171 (51.33) | 131 (15.93) | 134 (18.58) |  126 (11.5) | 113
50 | 20x60 | 170 (64.23) | 187 (52.03) | 142 (15.45) | 147 (19.51) |  141 (14.63) | 123
51 | 50x60 | 280 (50.44) | 328 (45.13) |  238 (5.31) | 253 (11.95) | 240 (6.19) | 226
52 | 50x60 | 287 (52.65) | 323 (42.92) | 236 (4.42) | 242 (7.08) |  233 (3.1) | 226
53 | 50x60 | 280 (55.41) | 323 (45.5) | 236 (6.31) | 256 (15.32) |  231 (4.05) | 222
54 | 50x60 | 302 (55.7) | 343 (44.73) | 251 (5.91) | 268 (13.08) |  243 (2.53) | 237
55 | 50x60 | 297 (51.74) | 327 (42.17) | 250 (8.7) | 262 (13.91) |  237 (3.04) | 230
56 | 100x60 | 499 (36.64) | 569 (31.11) | 425 (-2.07) | 439 (1.15) |  411 (-5.3) | 434
57 | 100x60 | 526 (43.29) | 565 (30.79) | 427 (-1.16) | 442 (2.31) |  412 (-4.63) | 432
58 | 100x60 | 475 (34.09) | 592 (34.55) | 428 (-2.73) | 442 (0.45) |  424 (-3.64) | 440
59 | 100x60 | 458 (36.26) | 567 (30.95) | 426 (-1.62) | 443 (2.31) |  424 (-2.08) | 433
60 | 100x60 | 493 (36.57) | 584 (31.83) | 434 (-2.03) | 458 (3.39) |  419 (-5.42) | 443
Mean Gap (%) | 55.65 | 38.48 | 6.42 | 7.27 | 2.91 |