Graph Neural Networks for the Offline Nanosatellite Task Scheduling
  Problem
This study investigates how to schedule nanosatellite tasks more efficiently
using Graph Neural Networks (GNNs). In the Offline Nanosatellite Task
Scheduling (ONTS) problem, the goal is to find the optimal schedule for tasks
to be carried out in orbit while taking into account Quality-of-Service (QoS)
considerations such as priority, minimum and maximum activation events,
execution time-frames, periods, and execution windows, as well as constraints
on the satellite's power resources and the complexity of energy harvesting and
management. The ONTS problem has been approached using conventional
mathematical formulations and exact methods, but their applicability to
challenging cases of the problem is limited. This study examines the use of
GNNs in this context, which has been effectively applied to optimization
problems such as the traveling salesman, scheduling, and facility placement
problems. More specifically, we investigate whether GNNs can learn the complex
structure of the ONTS problem with respect to feasibility and optimality of
candidate solutions. Furthermore, we evaluate using GNN-based heuristic
solutions to provide better solutions (w.r.t. the objective value) to the ONTS
problem and reduce the optimization cost. Our experiments show that GNNs are
not only able to learn feasibility and optimality for instances of the ONTS
problem, but they can generalize to harder instances than those seen during
training. Furthermore, the GNN-based heuristics improved the expected objective
value of the best solution found under the time limit in 45%, and reduced the
expected time to find a feasible solution in 35%, when compared to the SCIP
(Solving Constraint Integer Programs) solver in its off-the-shelf configuration
Bruno Machado Pacheco, Laio Oriel Seman, Cezar Antonio Rigo, Eduardo Camponogara, Eduardo Augusto Bezerra, Leandro dos Santos Coelho
# A Graph Neural Network Approach to Nanosatellite Task Scheduling:

###### Abstract

This study investigates how to schedule nanosatellite tasks more efficiently using Graph Neural Networks (GNN). In the Offline Nanosatellite Task Scheduling (ONTS) problem, the goal is to find the optimal schedule for tasks to be carried out in orbit while taking into account Quality-of-Service (QoS) considerations such as priority, minimum and maximum activation events, execution time-frames, periods, and execution windows, as well as constraints on the satellite's power resources and the complexity of energy harvesting and management. The ONTS problem has been approached using conventional mathematical formulations and precise methods, but their applicability to challenging cases of the problem is limited. This study examines the use of GNNs in this context, which has been effectively applied to many optimization problems, including traveling salesman problems, scheduling problems, and facility placement problems. Here, we fully represent MILP instances of the ONTS problem in bipartite graphs. We apply a feature aggregation and message-passing methodology allied to a ReLU activation function to learn using a classic deep learning model, obtaining an optimal set of parameters. Furthermore, we apply Explainable AI (XAI), another emerging field of research, to determine which features - nodes, constraints - had the most significant impact on learning performance, shedding light on the inner workings and decision process of such models. We also explored an early fixing approach, obtaining an accuracy above 80% both in predicting the feasibility of a solution and the probability of a decision variable value being in the optimal solution. Our results point to GNNs as a potentially effective method for scheduling nanosatellite tasks and shed light on the advantages of explainable machine learning models for challenging combinatorial optimization problems.

keywords: Scheduling, Graph Neural Network, Combinatorial Optimization, Nanosatellite, Quality of Service.

## 1 Introduction

Nanosatellites are gaining popularity for various applications, including Earth observation and scientific research. Due to its limited computational and energy resources, this spacecraft standard is associated with difficulties in mission planning despite the format's clear advantages, such as low cost and fast development time. Scheduling tasks is essential to mission planning, as it maximizes resource usage and increases data quality, cost savings, and mission success. The Offline Nanosatellite Task Scheduling (ONTS) problem is crucial in developing, deploying, and operating nanosatellites in orbit. It involves finding the best schedule for tasks execution in orbit, taking into account Quality-of-Service (QoS) factors such as priority, minimum and maximum activation events, execution time-frames, periods, and execution windows, as well as the limitations of the satellite's power resources and complexities of energy harvesting and management systems.

Traditional mathematical formulations and exact algorithms have been proposed to solve the ONTS problem, starting from Integer Programming (IP) [1] to Mixed Integer Programming (MILP) [2; 3] and Continuous-Time techniques [4]. More recently, given the difficulty in solving complex instances of the ONTS problem, Rigo et al. [5] proposed a Dantzig-Wolfe decomposition and a branch-and-price (B&P) methodology to build a unique
column-based formulation for producing feasible and optimal schedules. They also explored the Dynamic Programming (DP) technique to find optimal columns. Their computational experiments significantly improved overall solution time compared to a commercial MILP solver, with a 70% reduction in computation time. The evolution of these formulations and methodologies highlights the continuous efforts to find the most efficient and effective solutions to the problem.

Meanwhile, several recent investigations have considered machine learning tools to address combinatorial optimization problems [6; 7], such as the single machine problem [8], resource-constrained project scheduling [9], and knapsack problems [10]. Graph Neural Networks (GNN), in particular, has gained popularity in recent years to solve combinatorial optimization problems when the underlying structure may be represented as a graph [11]. The problem's graph structure is employed to transfer information between nodes, and several iterations of message passing update the node representations. The final representations of nodes can be utilized to generate predictions or solve optimization problems. GNNs are well-suited for handling combinatorial optimization problems because they can simulate complex problem structures and convey information across nodes. They have been successfully used for many optimization problems, such as traveling salesman problems [12], scheduling tasks, and facility location problems.

A popular approach has been the application of GNNs to learn variable selection for branching, directly [13], with the help of Markov Decision Process [14] or multi-layer perceptrons [15]. For instance, the authors of [16] were the first to suggest a novel graph convolutional neural network model that uses the inherent variable-constraint bipartite graph representation of MILPs, which is trained via imitation learning from a strong branching expert method. They show, using complex problems, that the technique can beat expert-designed branching rules applied in cutting-edge solvers, generating policies that outperform state-of-the-art machine-learning approaches for branching. They concluded that their model was a superior design choice before branching in MILP, and subsequent work demonstrated the feasibility of their strategy on a larger range of combinatorial problems tested with graph-based reinforcement learning techniques.

In the MILP framework, [15] presents a novel hybrid architecture for efficient branching on GPUs. For branching, they suggested an architecture to combine the capabilities of GNN with computationally cheap multi-layer perceptrons (MLP). The authors tested their technique on four classes of MILP problems and found that it reduced solver running time by up to 26% when compared to state-of-the-art solutions without a GPU, even when extrapolated to more complex problems than it was trained on. Further works explore GNNs in exact algorithms to select which cutting plane to add [17] or even learn a parallel Lagrangean decomposition by encoding the duals on a bipartite graph [18].

Several researchers have also explored GNNs to learn and improve on heuristic approaches for solving combinatorial optimization problems [19], such as the large neighborhood search [20; 21]. In [22], a Neural Improvement (NI) model for graph-based problems is presented that can efficiently guide various hill-climbing algorithms. The model leverages information stored in both nodes and edges and may be used to replace classic local-search algorithms while requiring less processing effort. The authors examine advanced models in order to avoid becoming locked in local optima, as well as to construct models for population-based metaheuristics. Experiments reveal that the NI model outperforms traditional variants for the preference ranking, traveling salesman, and graph partitioning problems.

A general framework for augmenting MILP solvers with data-driven insights by predicting variable biases using GNN topologies is proposed by [23]. The predicted biases are used to steer the solver, substituting heuristic components by storing the variable-constraint interactions as bipartite graphs. The framework is demonstrated to significantly enhance the solver's performance on two classes of difficult binary MILPs, and it is extendable to additional important solver components. The work of [24] explores learning to fix variables early in iterative approximation approaches applied to IP problems. The authors model the early fixing as a Markov decision process and train it through imitation learning. They undertake comprehensive experiments on three typical IP applications and demonstrate that their technique may dramatically accelerate prior approximation methods by up to ten times in most situations while generating comparable or even better solutions. The authors also analyze the use of their suggested learning-based early fixing approach and potential prospects for increasing its efficacy.

Beyond the pure application of GNN into optimization problems, this subfield of ML has been rapidly evolving as well, where unique new techniques have been proposed recently on how to explore the graph architecture in neural networks better [25], such as gated graphs [26], large graphs [27], and directional graphs [28]. In [29], a new pre-training strategy for graph datasets is introduced, named Graph Isomorphism Network with Edge Features (GINEConv). It involves training a GNN at both node and graph levels to learn good local and global representations simultaneously. The authors systematically studied pre-training on multiple graph classification datasets. Their proposed strategy significantly improves out-of-distribution generalization, achieving state-of-the-art performance. In contrast, [30] introduced the now widely used Graph Attention Network (GAT) technique
later improved in [31]. GAT uses masked self-attention layers to allow nodes to attend to their neighbor's features and implicitly assign weights to different nodes without costly matrix operations or prior knowledge of the graph structure.

Taking advantage of all this recent progress in GNN research and its successful application to optimization problems, this study proposes a novel solution methodology to the ONTs problem. By representing the problem as a bipartite graph, we leverage the robust representation learning capabilities of GNNs. The parameters of the MILP problem generate feature vectors fed into the model, allowing us to encode both the structure and parameters of each instance of the problem. GNNs can handle graphs of arbitrary size to handle optimization problems with varying numbers of variables and constraints. Our method employs two fully-connected, single-layer multilayer perceptron (MLP) networks with ReLU activations to encode the features of variables and constraints into hidden features that are updated using a two-step message-passing mechanism. The parameters can then be optimized similarly to conventional deep learning models. The proposed GNN model is then shown here to accelerate task scheduling by efficiently learning the relationships between tasks and resources and optimizing mission planning.

Furthermore, research in Artificial Intelligence (AI) as a whole has recently focused on Explainable Artificial Intelligence (XAI), which analyzes the factors that impact solution quality and their interconnections in AI methodologies [32; 33], such as those of traffic classification [34] of variable selection [35]. Similarly, Explainable Graph Neural Networks (XGNN) can offer insights into the underlying mechanisms of these black-box models and assist researchers and users in understanding better how these models are producing predictions [36; 37] through parameterized explanations [38], probabilistic explanations [39], or attribution evaluation [40]. A model-neutral method for explaining the predictions of any Graph Neural Network (GNN) on any graph-based machine learning is presented in [41]. The method uses the recursive neighborhood-aggregation methodology of GNNs to pinpoint significant graph paths and pertinent node characteristic data sent along each edge. The approach uses relational structures, including rich node-featured graphs, and offers an interface for analyzing GNN predictions, troubleshooting GNN models, and spotting systematic error patterns. In the study, GnnExplainer is defined as an optimization job that optimizes the mutual information between the prediction of a GNN and the distribution of potential subgraph structures.

In [42], the authors create graph analogs of three well-known explainability techniques for GNNs, including gradient-weighted CAM (Grad-CAM) and contrastive EB. These techniques include contrastive gradient-based saliency maps, class activation mapping, and excitation backpropagation (c-EB). Moreover, they examine the important sub-graphs derived from the explanations and note recurring trends. Our work also aims to explore explainability to understand better how a GNN and its attributes can contribute to the final quality of ONTS problem solutions, providing insights into this field of knowledge.

This paper is organized as follows: Section two describes the problem statement in detail, providing context and background information. The third section describes how the problem was approached and the methods employed for the investigation. The computational experiments are described in detail in the fourth section, along with a succinct summary of the findings. This study finishes with section five, summarizing and discussing the key findings.

## 2 Problem Statement

Given a set of jobs \(\mathcal{J}=\{0,...,J\}\) that represent a mission and a set of time units \(\mathcal{T}=\{0,...,T\}\) that represents the orbit period, the objective function (1) represents the goal of maximizing the mission quality of service (QoS) metric, which is represented as the sum of the priority values \(u_{j,t}\) for all the jobs \(j\) over all the periods \(t\).

\[QoS:\ \max_{x_{j,t}}\ \underbrace{\sum_{j=1}^{J}\sum_{t=1}^{T}u_{j,t}x_{j,t}}_{ \text{Quality of Service}} \tag{1}\]

Variable \(x_{jt}\) represents the binary decision of scheduling job \(j\) at time \(t\), which takes on value \(1\) if job \(j\) is scheduled to run at time \(t\) and \(0\) otherwise.

In constraints (2a) to (2d), the variable \(\phi_{j,t}\) is used to describe the period between task executions. In essence, these equations enforce the relationship between \(\phi_{j,t}\) and \(x_{j,t}\) such that \(\phi_{j,t}\) assumes value 1 only in the time step the task \(j\) started running and is later used to reflect the desired period between task executions in the scheduling problem.

\[\phi_{j,t}\geq x_{j,t}, \tag{2a}\]
\[\phi_{j,t}\leq x_{j,t}-x_{j,(t-1)}, \forall j\in\mathcal{J},\,\forall t\in\mathcal{T}:t>1 \tag{2b}\] \[\phi_{j,t}\leq x_{j,t}, \forall j\in\mathcal{J},\,\forall t\in\mathcal{T}\] (2c) \[\phi_{j,t}\leq 2-x_{j,t}-x_{j,(t-1)}, \forall j\in\mathcal{J},\,\forall t\in\mathcal{T}:t>1\] (2d) \[\sum_{t=1}^{w_{j}^{\min}}x_{j,t}=0, \forall j\in\mathcal{J}\] (2e) \[\sum_{t=w_{j}^{\max}+1}^{T}x_{j,t}=0, \forall j\in\mathcal{J}\] (2f) \[\sum_{l=t}^{t+t_{j}^{\min}-1}x_{j,l}\geq t_{j}^{\min}\phi_{j,t}, \forall t\in\{1,...,T-t_{j}^{\min}+1\}, \forall j\in\mathcal{J}\] (2g) \[\sum_{l=t}^{t+t_{j}^{\max}}x_{j,l}\leq t_{j}^{\max},\,\forall t \in\{1,...,T-t_{j}^{\max}\}, \forall j\in\mathcal{J}\] (2h) \[\sum_{l=t}^{T}x_{j,l}\geq(T-t+1)\phi_{j,t},\forall t\in\{T-t_{j}^ {\min}+2,...,T\}, \forall j\in\mathcal{J}\] (2i) \[\sum_{l=t}^{t+p_{j}^{\min}-1}\phi_{j,l}\leq 1,\forall t\in\{1,...,T-p_ {j}^{\min}+1\}, \forall j\in\mathcal{J}\] (2j) \[\sum_{l=t}^{t+p_{j}^{\max}-1}\phi_{j,l}\geq 1,\forall t\in\{1,...,T-p_ {j}^{\max}+1\}, \forall j\in\mathcal{J}\] (2k) \[\sum_{t=1}^{T}\phi_{j,t}\geq y_{j}^{\min}, \forall j\in\mathcal{J}\] (2l) \[\sum_{t=1}^{T}\phi_{j,t}\leq y_{j}^{\max}, \forall j\in\mathcal{J}\] (2m) \[\phi_{j,t}\in\{0,1\}, \forall j\in\mathcal{J},t\in\mathcal{T}\] (2n) \[x_{j,t}\in\{0,1\}, \forall j\in\mathcal{J},t\in\mathcal{T} \tag{2o}\]

The constraints (2e) and (2f) are related to the execution of tasks in a given time window. The first, (2e), states that the sum of the binary variables \(x_{j,t}\) over the time interval \([1,w_{j}^{\min}]\) must be equal to zero, for all \(j\in\mathcal{J}\). Here, \(w_{j}^{\min}\) is the time when task \(j\) can start execution, meaning that it cannot run before this time. The second type of constraints, (2f), states that the sum of the binary variables \(x_{j,t}\) over the time interval \([w_{j}^{\max}+1,T]\) must also be equal to zero, for all \(j\in\mathcal{J}\). Here, \(w_{j}^{\max}\) is the maximum allowed time window for task \(j\), and \(T\) is the total number of time steps in the scheduling horizon. This means that if task \(j\) starts, it must finish before the time point \(w_{j}^{\max}\); otherwise the task cannot be executed. These constraints enforce that the tasks are executed only within the specified time windows, which can be used to ensure that a payload, for instance, runs only when passing above a certain territory.

Now, constraints (2g) ensure that if \(\phi_{j,t}\) is 1, meaning that task \(j\) started running at time \(t\), then at least \(t_{j}^{\min}\) units of \(x_{j,l}\) in the corresponding time window must also be \(1\). Similarly, (2h) ensures that the number of \(x_{j,l}\) values equal to \(1\) in the corresponding time window is limited by \(t_{j}^{\max}\). Complementary, (2i) ensures that if \(\phi_{j,t}\) is 1, then all \(x_{j,l}\) values from \(t\) to the end of the time horizon must also be 1 so that, if a task starts at the end of the orbit, then it executes until the final time step. Therefore, (2g) to (2i) ensure the task running time requirements are met. Constraints (2j) and (2k) states that the sum \(\phi_{j,t}\) over a window of size \(p_{j}^{\min}\) or \(p_{j}^{\min}\) must be equal to 1, ensuring that the period of execution of this task is respected.

Constraint (2l) specifies that the sum of all values of \(\phi_{j,t}\) for job \(j\) must be greater than or equal to a lower limit \(y_{j}^{\min}\). This means that the job must be performed a minimum number of times within the given time period. Constraint (2m) specifies that the sum of all values of \(\phi_{j,t}\) for job \(j\) must be less than or equal to an upper limit \(y_{j}^{\max}\). This means that the job must be performed at most the specified maximum number of times within the given time period.

Regarding the energy management formulations, equation (3b) calculates the balance energy at time step \(t\), \(b_{t}\)
by subtracting the total energy generated from the solar panels (\(q_{j}x_{j,t}\)) from the total energy demand (\(r_{t}\)). The second equation, (3c), calculates the energy required from or delivered to the battery (\(i_{t}\)) at time step \(t\).

\[\sum_{j=1}^{J}q_{j}x_{j,t}\leq r_{t}+\gamma\:V_{b}, \forall t\in\mathcal{T} \tag{3a}\] \[b_{t}=r_{t}-\sum_{j\in\mathcal{J}}q_{j}x_{j,t}, \forall t\in\mathcal{T}\] (3b) \[i_{t}=\frac{b_{t}}{V_{b}}, \forall t\in\mathcal{T}\] (3c) \[\text{SoC}_{t+1}=\text{SoC}_{t}+\frac{i_{t}\:e}{60\:Q}, \forall t\in\mathcal{T}\] (3d) \[\text{SoC}_{t}\leq 1, \forall t\in\mathcal{T}\] (3e) \[\text{SoC}_{t}\geq\rho, \forall t\in\mathcal{T} \tag{3f}\]

Equation (3d) establishes the state of charge (SoC) of the battery at every unit of time in the orbit. It is given as the sum of the state of charge at time \(t\) and the energy balance in this time step resulting from the current flowing in or out of the battery, expressed in terms of the battery capacity (\(Q\)). It also considers the battery charge and discharge efficiency (\(e\)). Constraints (3e) state that the State of Charge (SoC) at any time must be less than or equal to \(1\). This means that the battery can never be overcharged, Complementary constraints (3f) state that the SoC at any time must be greater than or equal to \(\rho\) -- it is a typical practice in such sensitive applications to impose large margins of safety. Finally, constraints (3a) ensure that the power demand does not exceed power availability. The battery can provide up to \(\gamma\cdot V_{b}\) Watts of power.

## 3 Methodology

The most traditional approach to use a deep learning model in a task in which the input is an optimization problem is to vectorize all problem parameters and feed them to a traditional neural network such as a multilayer perceptron (MLP) [43]. However, this might not even be possible, as different instances of a problem might have different constraints and/or variables and, thus, would yield vectors of varying sizes to the deep learning model.

Furthermore, even if the dimensions are fixed, traditional neural networks do not exploit the symmetries that exist in an optimization problem. There is no inherent ordering in the constraints nor in the variables of an optimization problem, i.e., an instance is not changed if we change the order of the rows (columns) of \(A\) along with the elements of \(b\) (\(c\)). By vectorizing these elements, we impose an order, upon which the output of the neural network will depend. Therefore, these symmetries are not embedded in the structure of a traditional deep learning model, and would need to be enforced during training, i.e., the model would need to learn these symmetries. This might slow down significantly the training and will not give any guarantees that the model will generalize with the given knowledge.

### Graph Neural Networks

Graph neural networks, or message-passing neural networks, are generalizations of convolutional neural networks from grid-structured data to graphs. GNNs work by propagating features between neighboring nodes recurrently. Let \(G=(V,E)\) be a graph and \(H^{(0)}\in\mathbb{R}^{n\times d}\) an initial feature matrix associated with the nodes, in which each row \(h_{v}^{(0)}\in\mathbb{R}^{d}\) is the feature vector of node \(v\in V\). At each layer \(l\in\{1,...,L\}\) of the GNN, and for each node \(v\in V\), we first compute the messages \(m_{u,v}^{(l)}\) propagated by its neighbors \(u\in\mathcal{N}(v)\) based on their features \(h_{u}^{(l-1)}\),

\[m_{u,v}^{(l)}=M_{l}(h_{u}^{(l-1)}),\:u\in\mathcal{N}(v), \tag{4}\]

where \(\mathcal{N}(v)\) represents the set of neighbors of \(v\), and \(M_{l}(\cdot)\) is the message function of layer \(l\). Then, the features of node \(v\) are updated with the information from these messages

\[h_{v}^{(l)}=U_{l}\left(h_{v}^{(l-1)},\texttt{Aggregation}\left(\{m_{u,v}^{(l) }:u\in\mathcal{N}(v)\}\right)\right), \tag{5}\]

where \(U_{l}(\cdot)\) is the update function of layer \(l\), and \(\texttt{Aggregation}\) is a function that receives multiple message vectors and returns a single vector. The most usual choice for \(\texttt{Aggregation}\) is the sum, but many are the possibilities,
such as

\[\texttt{Aggregation}=\begin{cases}\frac{1}{|\mathcal{N}(v)|}\sum\limits_{u\in \mathcal{N}(v)}m_{u,v}^{(l)},&\text{if mean}\\ \max\limits_{u\in\mathcal{N}(v)}m_{u,v}^{(l)},&\text{if max}\\ \sum\limits_{u\in\mathcal{N}(v)}m_{u,v}^{(l)},&\text{if sum}\\ \sum\limits_{u\in\mathcal{N}(v)}\alpha_{u,v}m_{u,v}^{(l)},&\text{if attention}\end{cases} \tag{6}\]

where \(\alpha_{u,v}\) is the attention weight for node \(u\) given node \(v\). The attention weights can be learned using a neural network or other techniques. Furthermore, the message function can easily be extended to consider the edge weight (or even edge features) along with the feature vectors of the neighbors.

A common approach is to define the message functions \(M_{l},l=1,\ldots,L\) as linear operators over the hidden features of the neighbors, aggregate these messages by summing, and use the ReLU activation function with a bias as the update functions \(U_{l},l=1,\ldots,L\). We use the approach of [44] as a reference point and write

\[\begin{split} m_{u,v}^{(l)}&=\frac{1}{c_{vu}}W^{(l) }h_{u}^{(l-1)},\,u\in\mathcal{N}(v)\\ h_{v}^{(l)}&=\text{ReLU}\left(b^{(l)}+\sum\limits_{u\in \mathcal{N}(v)}m_{u,v}^{(l)}\right)\end{split} \tag{7}\]

where \(c_{vu}=\sqrt{|\mathcal{N}(u)|}\sqrt{|\mathcal{N}(v)|}\) with \(|\mathcal{N}(v)|\) denoting the number of neighbors, and \(W^{(l)}\in\mathbb{R}^{d\times d},b^{(l)}\in\mathbb{R}^{d}\) are (learnable) parameters.

A more recent method was proposed by [27] and named SAGE (SAmple and aGgrEgate). The authors propose to directly aggregate the features of the neighbors, i.e., to use the identity as the message function and apply a linear operator with a nonlinear activation as the update function. Putting it into terms,

\[\begin{split} m_{u,v}^{(l)}&=h_{u}^{(l-1)},\,u\in \mathcal{N}(v)\\ h_{v}^{(l)}&=\text{ReLU}\left(b^{(l)}+W_{1}^{(l)}h_ {v}^{(l-1)}+W_{2}^{(l)}\texttt{Aggregation}(m_{u,v}^{(l)},\,u\in\mathcal{N}( v))\right)\end{split} \tag{8}\]

where \(W_{1}^{(l)},W_{2}^{(l)}\in\mathbb{R}^{d\times d},b^{(l)}\in\mathbb{R}^{d}\) are the parameters. The authors suggest using more complex aggregation operators, such as an LSTM and a fully-connected single-layer neural network followed by a pooling operation (element-wise maximum).

After recurrent message passing operations through the \(L\) layers of a GNN, \(H^{(L)}\) can be further aggregated to generate a single feature vector of the entire graph. The GNN can be trained end-to-end by minimizing a prediction loss based on its outputs, optimizing its parameters (_e.g._, \(W^{(l)}\) and \(b^{(l)}\) of (7)) in the same way as a traditional deep learning model.

### GNNs for Combinatorics

Given a linear problem, we can build a graph \(G=(V,E)\) in which we add one node for each variable of the problem, one node for each constraint, and connect each variable node to constraint nodes whenever the coefficient of the respective variable is not null in the respective constraint. More precisely, given the problem of the form

\[\begin{split}\max&\quad c^{T}x\\ \text{s.t.:}&\quad Ax\geq b\end{split} \tag{9}\]

where \(x\in\mathbb{R}^{n}\) and \(b\in\mathbb{R}^{m}\), we can build a graph \(G=(V_{\text{var}}\cup V_{\text{con}},E)\), in which \(|V_{\text{var}}|=n\), \(|V_{\text{con}}|=m\), and \(E=\{(v_{var,i},v_{con,j}):A_{i,j}\neq 0\}\). Intuitively, the graph represents the structure of the problem at hand, the relationship between variables and constraints. Note that this approach yields a bipartite graph, that is, a graph in which the nodes are separated into two disjoint sets, \(V_{\text{var}}\) and \(V_{\text{con}}\), with edges connecting only nodes from different sets.

For illustration purposes, consider an optimization problem like Eq. 9 with

\[c=\begin{bmatrix}1\\ 2\\ 3\end{bmatrix};\;A=\begin{bmatrix}1&2&0\\ 0&1&-1\\ 3&0&1\end{bmatrix};\;b=\begin{bmatrix}2\\ 1\\ 4\end{bmatrix} \tag{10}\]
and \(x=[x_{1},x_{2},x_{3}]^{T}\). The bipartite graph can be represented as in Figure 1.

By representing an optimization problem as a graph, we can feed it to a GNN. The parameters of the optimization problem can be used to generate the feature vectors fed to the model, enabling us to codify not only the structure but also the parameters of any given instance of the problem. Because of the convolutional nature of the message-passing iterations of the GNNs, the model can deal with arbitrary-sized graphs, which enables us to handle optimization problems with varying variables and constraints with the same GNN. Furthermore, the message function is invariant to the ordering of the neighboring nodes (see Eq. (27)), which are precisely the symmetries of the optimization problem (order of variables and constraints).

### SatGNN

We name _SatGNN_ the network that serves as a basis for the experiments reported in the next section. To encode the features associated with the variables and the constraints into the hidden features of the first layer \(H^{(0)}\in\mathcal{R}^{(n+m)\times d}\), we use two fully-connected, single-layer MLPs, \(\text{NN}_{\text{var}}\) and \(\text{NN}_{\text{con}}\) with ReLU activations. We can write

\[h_{v}^{(0)}=\begin{cases}\text{NN}_{\text{var}}(f_{v}),&v\in V_{\text{var}}\\ \text{NN}_{\text{con}}(f_{v}),&v\in V_{\text{con}}\end{cases},\]

where \(f_{v}\) is the vector of features associated with each constraint or variable node. In our experiments, for a node \(v_{var,i}\in V_{\text{var}}\) associated with variable \(x_{i}\), the feature vector is \(f_{v_{var,i}}=(\hat{x}_{i},c_{i})\), where \(\hat{x}\) is a candidate solution and \(c_{i}\) is the weight of variable \(x_{i}\) in the objective function. Likewise, for a constraint node \(v_{con,i}\in V_{\text{con}}\) associated with the \(i\)-th constraint, \(f_{v_{\text{var}},i}=(b_{i},s_{i})\), where \(s_{i}\in\{=,\geq,\leq\}\) models the constraint type.

At the core of the model are \(MP\) operators, which perform the update of the hidden features of the nodes through message-passing, as described in Section 3.1,

\[H_{\mathcal{N}(v)}^{(l)}=\text{MP}_{l}(h_{v}^{(l-1)},H_{\mathcal{N}(v)}^{(l-1) },w_{\mathcal{N}(v)}),\]

where \(H_{\mathcal{N}(v)}^{(l-1)}=\{h_{u}^{(l-1)}:u\in\mathcal{N}(v)\}\) is the set of hidden features of the neighbors of the target node and \(w_{\mathcal{N}(v)}=\{w_{u,v}:u\in\mathcal{N}(v)\}\) is the set of edge weights. However, we generalize this operator to apply multiple convolutions at each node instead of applying a single convolution, considering the same context. The application of an MP operator with \(K\) convolutions is illustrated in Algorithm 1. This allows for more complex features to be extracted at each model layer.

```
Data: Target node \(v\), node features \(h_{v}\), neighbors' features \(H_{\mathcal{N}(v)}\), edge weights \(w_{u,v},\forall u\in\mathcal{N}(v)\). Result: Updated node features \(h_{v}^{*}\). \(h^{(0)}\gets h_{v}\) for\(k\gets 1\)to\(K\)do \(h^{(k)}\gets U^{(k)}\left(h^{(k-1)},\texttt{Aggregation}^{(k)}\left(M^{(k)}(h_{u},w_{u,v}):u\in \mathcal{N}(v)\right)\right)\)  end for\(h_{v}^{*}\gets h^{(K)}\)
```

**Algorithm 1**Application of an \(MP\) operator with \(K\) convolutions to update the features of a node \(v\) through message-passing.

The message-passing is split into two steps, one for each set of nodes, similar to the approach of [16]. At each layer, the messages are propagated first from the variable nodes to the constraint nodes and then from the constraint nodes to the variable nodes, exploiting the bipartite nature of the graph. Algorithm 2 describes this process with further detail.

Figure 1: Bipartition graph representation of Eq. (10)
Finally, the output of model is generated from the last hidden feature vectors of the variable nodes, generating an output in the same shape as the problem's variable vector. The hidden features are fed to an MLP with two hidden layers and ReLU activations \(\text{NN}_{\text{out}}\), which maps each \(d\)-dimensional vector into a single output, _i.e._,

\[\hat{y}_{v}=\text{NN}_{\text{out}}(h_{v}^{(L)}),\forall v\in V_{\text{var}}.\]

Figure 2 shows an overview of the architecture.

### XAI for GNN

Explainable Artificial Intelligence (XAI) refers to the development of AI models that can explain their predictions and human-comprehensible decisions. As AI models are increasingly utilized in high-stakes domains, such as healthcare and finance, where understanding why a model makes a specific prediction is critical, XAI is gaining importance.

XAI can be performed in the context of GNNs by employing interpretable models or interpretability approaches that reveal the inner workings of the GNN model. For example, feature importance analysis can be used to determine which graph features are most important for the prediction made by the GNN. The intermediate representations learned by the GNN can also be visualized using the activations visualization technique. These illustrations help show how the GNN analyzes data and generates predictions.

For instance, the authors in [45] introduce GNNExplainer, a new approach for explaining Graph Neural Network (GNN) predictions in graph-based machine learning tasks. GNNs are powerful, but explaining their predictions is challenging due to their complexity. GNNExplainer identifies a compact subgraph, and a subset of node features important for GNN predictions, allowing for consistent and concise explanations across instances. The approach maximizes the mutual information between GNN predictions and subgraph structures. Experiments show that GNNExplainer outperforms other methods and can identify important graph structures and node features, providing interpretability and insights into faulty GNNs.

Figure 2: Overview of the components of _SatGNN_ and the operations it performs given an optimization problem (embedded as a bipartite graph \(G\)). \(H\) variables represent the sets of hidden features of the nodes. The connection of \(G\) and both \(MP\) operators represents both the weights of the edges as well as the neighborhood information.
### Hyperparameter Optimization

Hyperparameter optimization is an essential task in machine learning that involves adjusting the hyperparameters of a model so that its performance is maximized. Hyperparameters are not learned during the training and must be defined before the training begins, e.g., learning rate, number of layers, number of learnable parameters, and regularization factor. Adjusting the hyperparameters can be done manually, through trial-and-error iterations, using an expert's intuition, or automatically treating the relationship between the hyperparameters and the model performance as a black-box function.

Optuna is a Python library that provides a flexible and efficient platform for hyperparameter optimization using a variety of algorithms, including TPE (Tree-structured Parzen Estimator). TPE is a Bayesian optimization variant that models the hyperparameter distribution using a tree-structured Parzen estimator. Let \(\mathbf{x}=(x_{1},x_{2},...,x_{d})\) denote a set of \(d\) hyperparameters, and let \(f(\mathbf{x})\) be the cost function to be optimized.

The TPE algorithm divides the hyperparameter space into two regions: a region containing the hyperparameters that have been observed to result in a good performance, denoted by \(\mathbf{u}\), and a region containing the hyperparameters that have not been observed to result in a good performance, denoted by \(\mathbf{v}\). TPE constructs the acquisition function \(a(\mathbf{x})\) as follows:

\[a(\mathbf{x})=\frac{p_{\mathbf{u}}(\mathbf{x})}{p_{\mathbf{v}}(\mathbf{x})}, \tag{11}\]

where \(p_{\mathbf{u}}(\mathbf{x})\) and \(p_{\mathbf{v}}(\mathbf{x})\) are the density functions estimated by the tree-structured Parzen estimator for the regions \(\mathbf{u}\) and \(\mathbf{v}\), respectively. Specifically, \(p_{\mathbf{u}}(\mathbf{x})\) and \(p_{\mathbf{v}}(\mathbf{x})\) are estimated as follows:

\[p_{\mathbf{u}}(\mathbf{x})=\frac{1}{|S_{\mathbf{u}}|}\sum_{\mathbf{y}\in S_{\mathbf{u}}}K\left( \frac{\mathbf{x}-\mathbf{y}}{\sigma_{\mathbf{u}}}\right), \tag{12}\]

and

\[p_{\mathbf{v}}(\mathbf{x})=\frac{1}{|S_{\mathbf{v}}|}\sum_{\mathbf{y}\in S_{\mathbf{v}}}K\left( \frac{\mathbf{x}-\mathbf{y}}{\sigma_{\mathbf{v}}}\right), \tag{13}\]

where \(S_{\mathbf{u}}\) and \(S_{\mathbf{v}}\) are the sets of hyperparameters observed in regions \(\mathbf{u}\) and \(\mathbf{v}\), respectively, \(K(\cdot)\) is a kernel function, and \(\sigma_{\mathbf{u}}\) and \(\sigma_{\mathbf{v}}\) are bandwidth parameters. The TPE algorithm selects the next set of hyperparameters to evaluate by maximizing the acquisition function:

\[\mathbf{x}_{\text{next}}=\arg\max_{\mathbf{x}}a(\mathbf{x}). \tag{14}\]

This process is repeated until the optimal set of hyperparameters is found, or a stopping criterion is met.

### Data

To supply the models with data suitable for learning the tasks of interest, we generate new instances of the OTS problem on demand, including the energy input and task QoS parameters. These methodologies are briefly presented in the following topics and were previously published in [46].

#### 3.6.1 Power Input Vector

An analytical model has been used to determine the power input vector of each instance of the OTS problem considered in this study. Once orbits are stable and solar flux constant - \(1360W/m^{2}\) - one can calculate this vector by knowing the spacecraft orbit, attitude - its kinematics - and size. We have taken the FloripaSat-I mission as a parameter for orbital data, which has an altitude of 628 kilometers and an orbital period of 97.2 minutes [47]. The attitude considered here is the Nadir, in which the satellite turns at the same rate around the Earth, so one side (or axis) always faces the Earth's surface. This analytical model then utilizes a rotation matrix to simulate the satellite's dynamics and can be adapted for larger or different geometries by adjusting the normal vectors representing the body. For this study, we considered a 3U nanosatellite size. The power generated by photovoltaic panels on each of the CubeSat's six sides depends on the efficiency of the cells, the area of the cells, the view factor of the surface to the Sun, and a step function that accounts for the satellite's location concerning Earth's shadow, as:

\[P_{k}=\eta A_{pv_{k}}I_{Sun}F_{k\to Sun}\Psi, \tag{15}\]

where \(\eta\) and \(A_{pv}\) are the solar cell efficiency and area, respectively; \(I_{Sun}\) is the solar flux; \(F_{k\to Sun}\) is the cell projection to the Sun; and \(\Psi\) is a step function that, when the spacecraft is in the shade of the Earth, takes a value of zero. More details about the equations used can be found in [48].
#### 3.6.2 Tasks Parameters

For any particular mission size and orbital length, the main objective is to generate a realistic ONTS case using random data. The number of tasks or time units can be increased to make the instances distinct. Algorithm 3 presents an instance generator technique that has been used here to accomplish this. It requires two inputs: the number of tasks (\(J\)) and the time units (\(T\)) to be taken into account. The process produces nine variables for each task, including \(u_{j}\), \(q_{j}\), \(y_{j}^{\min}\), \(y_{j}^{\max}\), \(t_{j}^{\min}\), \(t_{j}^{\max}\), \(p_{j}^{\min}\), \(p_{j}^{\max}\), \(w_{j}^{\min}\), and \(w_{j}^{\max}\).max. These parameters completely describe an instance of the ONTS regarding the QoS aspects.

```
Input : Number of jobs \(J\), number of time periods \(T\) Output : Initial values for \(u_{j},q_{j},y_{j}^{\min},y_{j}^{\min},t_{j}^{\min},t_{j}^{\max},y_{j}^{\min},p_{ j}^{\max},w_{j}^{\min},w_{j}^{\max}\) for\(j\gets 1\)to\(J\)do \(u_{j}\leftarrow\text{U}(1,J)\); \(q_{j}\leftarrow\text{U}(3,2,5)\); \(y_{j}^{\min}\leftarrow\text{U}[1,[T/45]]\); \(y_{j}^{\max}\leftarrow\text{U}[y_{j}^{\min},[T/15]]\); \(t_{j}^{\min}\leftarrow\text{U}[1,[T/10]]\); \(t_{j}^{\max}\leftarrow\text{U}[t_{j}^{\min},[T/4]]\); \(p_{j}^{\min}\leftarrow\text{U}[t_{j}^{\min},[T/4]]\); \(p_{j}^{\max}\leftarrow\text{U}[p_{j}^{\min},T]\); \(w_{j}^{\min}\leftarrow\text{U}[0,[T/5]]\); \(w_{j}^{\max}\leftarrow\text{U}[1-[T/5]],T]\);  end while
```

**Algorithm 3**Instance Generator Algorithm

## 4 Computational Experiments

The experiments for this paper were conducted in Python, using PyTorch and the DGL libraries, and the Gurobi solver, on a server with an Intel i7-12700 16-Core (12 cores, 20 threads), 16 GB of RAM, and Ubuntu 22.04.1 LTS 64 bits. An NVIDIA RTX A4000 was used to speed up the DGL library calculations.

In the following sections, we present three experiments that utilize GNN for different optimization problems regarding ONTS instances. In the first experiment, we propose a GNN-based approach to classify the feasibility of candidate solutions given problem instances. In the second experiment, we take a step further and propose a GNN to predict the consistency of each variable from a candidate solution with the maximization of the objective. Finally, in the last experiment, we used the GNN from experiment 2 to generate candidates suitable for early fixing binary variables of the problem.

### Experiment 1 - Feasibility Classification

In the first experiment, we aim to predict the feasibility of a candidate solution using GNNs. First, we tackle the ONTS problem for a single job, without constraints (3a) to (3f), which turns it into an example of a task scheduling problem. Then, we generalize this approach to the complete ONTS problem.

#### 4.1.1 Single Job Scheduling

For an instance \(I\) of the ONTS problem, we train SatGNN on data from multiple jobs \(j\in\mathcal{J}\), to learn the problem's underlying structure and to generalize to unseen jobs. As detailed in Section 3.3, we represent the task scheduling problem for a given job in a CubeSat as a bipartite graph. More specifically, we focus on an instance of the ONTS problem with 97-time steps and 9 jobs; therefore, each of the nine job scheduling problems has 194 variables and several constraints ranging from X to Y.

The bipartite graph for a given job of the instance and the feature vectors are fed to the SatGNN. We implement the model with a single regular convolution (see Eq. 7) in the message-passing operators (\(K=1\)), weight sharing between the operators (\(MP_{\text{con}}=MP_{\text{var}}\)), and a single message-passing iteration (\(L=1\)). Furthermore, as the task at hand requires a classification of the entire candidate solution, we aggregate the outputs into a single value

\[\hat{y}=\sigma\left(\frac{1}{n}\sum_{v\in V_{\text{var}}}\hat{y}_{v}\right),\]

where \(\sigma:\mathbb{R}\rightarrow[0,1]\) is the sigmoid function.
We first build a dataset of random candidate solutions alongside their feasibility to train the model. More precisely, we build a dataset \(\mathcal{D}\) composed of tuples \((\hat{x},j,y)\in\mathbb{Z}^{n}\times\mathcal{J}\times\{0,1\}\) in which \(\hat{x}\) is a candidate solution, and \(y\) takes the value \(1\) whenever \(\hat{x}\) is feasible for the task scheduling problem defined by job \(j\), being \(0\) otherwise. For each of the 9 jobs of the selected instance of the OTS problem, we generate 1000 pseudo-random candidate solutions, half of which are feasible. The feasible candidate solutions are generated by solving the optimization problem with Gurobi and retrieving a sample of solutions near the optima. The infeasible candidate solutions are generated by perturbing the decision variables (and updating the non-decision variables accordingly) of the feasible solutions until they violate some constraint. Therefore, \(\mathcal{D}\) is a balanced dataset with 9000 elements.

For training, we first split the dataset into training, validation, and test sets \(\mathcal{D}_{\text{train}}\cup\mathcal{D}_{\text{val}}\cup\mathcal{D}_{\text{ test}}=\mathcal{D}\) in such a way that \(\mathcal{D}_{\text{test}}\) contains all elements of \(\mathcal{D}\) associated to one of the jobs, \(\mathcal{D}_{\text{val}}\) contains all elements associated to a different job, and \(\mathcal{D}_{\text{train}}\) contains the elements associated to the remaining 7 jobs, i.e., no job is present in both datasets. We then optimize the parameters of the GNN to minimize the binary cross-entropy between the predicted feasibility and the actual feasibility of the candidate solutions

\[\sum_{(\hat{x},j,y)\in\mathcal{D}_{\text{train}}}-y\log(\text{GNN}(\hat{x},j) )-(1-y)\log(1-\text{GNN}(\hat{x},j)),\]

where \(\text{GNN}(\hat{x},j)\) is the model's predicted probability of \(y=1\), and we evaluate the models based on their accuracy on data unseen during training.

The training is performed with the Adam optimizer [49] with a budget of 100 epochs. We observed that the model's performance was highly dependent on the initialization. Therefore, multiple models were trained, with random Glorot uniform weight initialization, as described in [44]. We then select the best model on the validation set and evaluate it on the test set. The performance of the model on \(\mathcal{D}_{\text{test}}\) can be seen in Table 1.

#### 4.1.2 Ots Problem

Seeing that the feasibility classifier GNN could learn and perform on unseen samples of the task scheduling problem, we generalize the approach described above to the complete OTS problem. The graph is constructed in the same way, but now considering all jobs and all constraints described in Section 2. Therefore, the problem now has 9 times more integer variables and several others necessary for the coupling constraints (3a) to (3f).

The architecture is also identical, with the only exception being the output aggregation. Once the continuous variables can be determined entirely from the values of the binary variables, the prediction is made solely with features from the nodes that correspond to the binary variables. In other words, let \(V_{\text{var}}=V_{\text{int}}\cup V_{\text{cont}}\) where \(V_{\text{int}}\) is the set of nodes associated with the integer variables and \(V_{\text{cont}}\) the set of nodes associated with the continuous variables. Then, the output is computed as

\[\hat{y}=\sigma\left(\frac{1}{n}\sum_{v\in V_{\text{in}}}\hat{y}_{v}\right)\]

The dataset is built with 21 different instances of the OTS problem. For each instance, we generate 1000 pseudo-random candidate solutions following the same approach as previously. Once again, we split the data into training, validation, and test sets, with the samples from two instances on the test set, the samples from two other instances on the validation set, and the samples from the remaining 17 instances on the training set. As a similar initialization impact was observed, we followed the same procedure for model selection based on validation performance. The performance of the final model can be seen in Table 1.

For an additional analisys, as can be seen in Figure 3, the more relevant variables for the decision process were regarding \(\phi\), which determine the exactly startup time of a task; since \(x\) is heavily dependent on \(\phi\), this aids the model to maker better overall choices.

#### 4.1.3 XAI on Feasibility Classification

In this sub-experiment, we investigate which inputs had the greatest impact on classifying a solution as feasible or not feasible for both the single job problem and the complete OTS problem. First, GNNExplainer was used

\begin{table}
\begin{tabular}{l c c} \hline \hline Task & Accuracy & F1 score \\ \hline Single Job Feasibility Classification & 83.0\% & 0.8172 \\ Full ONTS Feasibility Classification & 75.2\% & 0.7956 \\ \hline \hline \end{tabular}
\end{table}
Table 1: SatGNN’s test set performance on the feasibility classification task.
on the trained SatGNN model to determine which variable nodes, constraint nodes or even features in the bipartite graph had the most influence on the feasibility categorization of a solution.

By identifying the most pertinent nodes and restrictions, we better understood the GNN model's decision-making process, which could improve the model's interpretability and lead to the design of more efficient solutions for the task scheduling problem. Considering the graph explanation, more specifically regarding edges types "var2con" and "con2var", it was possible to notice that no edge type had significant more importance than the other; also, no linear correlation was observed between them, as shown in Figure 5. On the other hand, when trying to summarize the most important nodes and edges involved in the output decision as a sub-graph, the number of edges can be reduced from \(7160\) to \(5722\) in the resulting sub-graph, and the number of variables from \(194\) to \(139\); which once again emphasizes how most of the variables and connections are important in the final decision.

Figure 4: Edge mask considering graph explanation for “var2con” and “con2var”, considering the correlation between the two edges types in a mean of 50 instances.

Figure 3: Feature importance regarding the model decisions variables; it is possible to observe that \(\phi\) is more easily distinguished by the GNN, which directly influence the feasible region of \(x\).
### Experiment 2 - Optimality Classification

In our second experiment, we aim to predict the probability of each integer variable in a given candidate solution being consistent with the maximization of the result for the ONTS problem. In other words, given a candidate solution \(\hat{x}\in\{0,1\}^{n}\), the ideal output would be a vector \(y\in[0,1]^{n}\) in which \(y_{i}=\mathbf{1}(x_{i}^{*}=\hat{x}_{i})\), where \(x^{*}\) is the optimal solution.

We represent the problem as a bipartite graph and apply SatGNN to classify the variable. Differently from the experiments above, we do not need to aggregate the output of SatGNN. Instead, we apply the sigmoid function directly to the output of each node associated with the integer variables. We build a dataset with the same 21 instances of the ONTS problem as in the previous experiment. We generate pseudo-random candidate solutions and their label for each instance, which is computed given the optimal solution. The optimal solution for each of the 21 instances was found using the Gurobi solver.

The dataset was divided into training, validation, and test sets, with validation and test having the data from two instances each and the training set with the data from the remaining 17 instances. The instances in the test set will be referred to as instances A and B. Multiple models with different random initializations were trained on the validation set, which was also used to perform hyperparameter tuning, as described in sec. 3.5. The hyperparameters selected for tuning alongside the best configuration found can be seen in table 2. The importance of each hyperparameter was assessed by training a random forest on the task of predicting the performance measure based on the values of the hyperparameters, upon which the importance is taken as the Gini importance of each hyperparameter [50]. The hyperparameter importance can be seen in Figure 6, while the coordinate plot of the hyperparameters interconnections is presented in Figure 7.

On the test set, the best model could correctly predict, on average, 85.4% of the variables (with a standard deviation of 2.0 p.p.) and achieved an average F1 score of 0.8531 (standard deviation of 0.02). A summary of the best model's performance can be seen in Table 3. We also analyze the output of the model for each variable. Overall, the histogram of the model output of Figure 7(a) indicates that the prediction is usually close to the interval's limits, i.e., it is approximately binary. Furthermore, we evaluate the accuracy for each of the 1746 variables on all candidate solutions of instances A and B. These accuracies can be seen in Figure 7(b). The model could correctly predict the optimality of most variables in both instances of the test set.

Furthermore, we evaluate the out-of-distribution generalization capacity of the SatGNN model in the optimality classification task by feeding it with larger instances of the problem. More specifically, we generate two new instances: instance C has the same time horizon but 11 jobs instead of 9; instance D has the same number of jobs but requires scheduling over 120-time steps instead of 97. Then, we generate new random candidate solutions in the same way as previously described. The performance of SatGNN on the larger instances can be seen in Table 3. Not only was the model able to handle the new instances without any modifications, but it also achieved an average accuracy of 80.8% and an average F1 score of 0.8037 on instance C (more jobs), and an average accuracy of

Figure 5: When trying to explain the output with a resulting subgraph, the number of edges can be reduced from \(7160\) to \(5722\) in the resulting sub-graph, and the number of variables from \(194\) to \(139\), which once again emphasizes how most of the variables and connections are important in the final decision.
73.6% and an average F1 score of 0.7405 on instance D (more time steps). More details on the output of SatGNN model for the larger instances can be seen in Figure 9.

### Experiment 3 - Early Fixing

Using the outcomes of the previous experiment, we tackle the task of early fixing variables of the ONTS problem. For this, we point out that it is possible to recover the optimal solution to the problem given any \(\hat{x}\) candidate solution along with its associated \(y\) label. Therefore, we get a predicted optimal solution by using the predicted optimality of the candidate solution (output of the SatGNN model from experiment 2). Additionally, as the predicted optimal solution can be generated from any random candidate solution, we use a set of random candidates and average the predicted optimals. Figure 10 illustrates how to use SatGNN for early fixing.

Therefore, given a set \(\hat{X}\) of random candidate solutions for a given problem instance, we compute

\[\hat{x}^{*}=\frac{1}{|\hat{X}|}\sum_{\hat{x}\in\hat{X}}\hat{x}\odot\hat{y}(\hat {x})+(1-\hat{x})\odot(1-\hat{y}(\hat{x})),\]

\begin{table}
\begin{tabular}{l l l} \hline Hyperparameter & Value range & Final value \\ \hline \# of MP operations & 1..3 & 3 \\ Convolutions’ type & [Regular, SAGE] & Regular\(\rightarrow\) SAGE\(\rightarrow\) Regular \\ SAGE aggregation & [lstm, pool] & pool \\ SAGE feature drop & 0.0..0.5 & 0.09 \\ \# of hidden features & 2..20 & 19 \\ Share weights between \(MP_{\text{var}}\) and \(MP_{\text{con}}\) & yes/no & no \\ \# layers & 1..3 & 1 \\ \# random samples per instance & \(2^{6}..2^{10}\) & \(2^{9}\) \\ batch size & \(2^{2}..2^{7}\) & \(2^{2}\) \\ \hline \end{tabular}
\end{table}
Table 2: Hyperparameters of SatGNN selected for hyperparameter optimization using Optuna, along with the values of the best model found. For “Type of each MP operation”, _Regular_ represents the convolution using a linear operator to combine the neighbor features and a ReLU activation function, as in eq. 7; _SAGE_ represents the convolution operator used in SageGNN, as described in sec. 3.1.

\begin{table}
\begin{tabular}{l c c c c} \hline \hline  & \multicolumn{4}{c}{Accuracy} \\ Task & Instance A & Instance B & Instance C & Instance D \\ \hline Optimality Classification & 83.4\% & 87.3\% & 80.8\% & 73.6\% \\ Early Fixing & 82.5\% & 87.9\% & 66.9\% & 70.9\% \\ \hline \hline \end{tabular}
\end{table}
Table 3: SatGNN’s performance on new instances (unseen during training) of the ONTS problem. For “Optimality Classification”, average accuracy over all samples is reported.

Figure 6: Importance of SatGNN hyperparameters tuned for the optimality classification. The score is calculated from the Gini importance of a random forest fitted on predicting the model’s performance from the hyperparameters’ values. The values were normalized to sum 1.
where \(\odot\) is the element-wise product and \(\hat{y}(\hat{x})\) is the predicted optimality of candidate solution \(\hat{x}\in\hat{X}\) generated using the model from the previous experiment. In the results reported

Furthermore, we can say that the closer a given predicted optimal variable \(\hat{x}_{i}^{*}\) is to 1 (resp. 0), the more certain the model is that that variable should be fixed at 1 (resp. 0). Therefore, we use the model's certainty to select the

Figure 8: (a) Histogram of the outputs given the samples from the test set. Each dimension (variable node) of the output is treated as a different occurrence; (b) Accuracy for each binary variable of the problem over all candidate solutions on the test set. Orange and blue (with transparency) were used to distinguish between the two instances in the test set. Over 99% accuracy was observed for the same 1208 variables on both instances A and B (darker rows), while less than 1% accuracy was observed on both instances for 52 variables (white rows), over 1746 variables in total.

Figure 7: Hyperparameters interconnections during the trials.
variables to be fixed; that is, if we want to fix 50 binary variables, we will choose the 50 variables that the model is most certain of. We evaluate the accuracy of the SatGNN for early fixing as a function of the number of fixed variables on the two instances of the test set. These results can be seen in figure 11. As expected, the accuracy decreases as we include variables for which the model is less certain, to the limit of 82.5% and 87.9% accuracy on the two instances, which is the accuracy of the predicted optimal solution over the 1745 variables. A summary of the model's performance when fixing all variables can be seen in Table 3.

In face of these results, we evaluate how early fixing using SatGNN impacts the optimization performance both in terms of runtime and maximum objective value. Specifically, we solve the two instances of the ONTS problem on the test set using Gurobi under an increasing number of fixed variables. The results can be seen in figure 12.

As expected, correctly fixing the variables positively impacts the optimization, while wrongly fixing variables may decrease the runtime but often impacts the objective negatively. However, we see that, at the limit, a substantial runtime reduction is achieved (90% and 28%, for instances 6 and 9, resp.) with a negligible objective cost (1.3% and 0.3%, resp.). Beyond that, fixing more than 500 variables for instance 6 and more than 200 for instance

Figure 10: Early fixing with SatGNN.

Figure 9: SatGNN optimality classification output for larger instances. Over 99% accuracy was observed for 1556 and 1230 variables (solid rows), and less than 1% accuracy was observed for 227 and 178 variables (white rows) for instances C and D, respectively. Instance C has 2134 binary variables, while instance D has 2160.
9 deemed the problems infeasible within a 5 minutes budget.

Given that the SatGNN model could generalize the optimality classification for larger instances of the problem, we also evaluate the impact of early fixing based on our model for the same two larger instances used in the previous experiment. The performance on the larger instances can be seen in Figure 13 and in Table 3. Even though these instances' sizes were not seen during training (not even during validation), the model was still able to handle them and provide sensible early fixing candidates. The performance drop is significant in terms of accuracy and optimization performance. In most configurations, however, the model was still able to reduce the runtime with little to no objective value reduction.

## 5 Conclusion

This work has proposed a novel approach to tackle the ONTS problem using graph neural networks. Our experiments showed that our proposed architecture, SatGNN, is successful in classifying the feasibility and the optimality of candidate solutions to varied instances of the ONTS problem. Not only the model was able to generalize to unseen instances, but it also showed promising results on out-of-distribution instances, which were larger than the ones seen during training. This shows how the inherent symmetries of graph neural networks make them suitable for dealing with the structures of optimization problems.

By leveraging on the optimality classification results, we used the SatGNN to generate candidate solutions to the binary variables of the problem. Through these candidate solutions, we were able to fix the variables, reducing the size of the problem and, consequently, the time the solver takes to converge. This approach for early fixing outperformed using the Gurobi solver alone, even when no tolerance is permitted for the optimal. Furthermore, generalization to larger instances is still a challenge for early fixing, even though promising results were observed.

In summary, this work has shown how graph neural networks can be used to improve nanosatellite task scheduling. As we propose a supervised learning method, our approach is still limited to the availability of labeled data, which can be costly to obtain. Nonetheless, our results suggest that using graph neural networks for combinatorial optimization problems holds excellent promise and opens up new avenues for future research.

Figure 11: Early fixing accuracy for the two instances of the ONTS problem in the test set.

Figure 12: Optimization results of the two ONTS instances with SatGNN-based early fixing. The objective is plotted with respect to the maximum of the original problem (without any fixed variables). Accuracy is measured with respect to the optimal value of the fixed variables.
## Acknowledgments

The authors acknowledge support from CNPq (Conselho Nacional de Desenvolvimento Cientifico e Tecnologico) under grant number 150281/2022-6, 404576/2021-4 as well as FAPESC under grant number 2021TR001851.
