{'fraction_non_alphanumeric': 0.035203996669442135, 'fraction_numerical': 0.009625312239800166, 'mean_word_length': 4.8256693830034925, 'pattern_counts': {'":': 1, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning to Schedule Heuristics for the Simultaneous Stochastic Optimization of Mining Complexes

The simultaneous stochastic optimization of mining complexes (SSOMC) is a large-scale stochastic combinatorial optimization problem that simultaneously manages the extraction of materials from multiple mines and their processing using interconnected facilities to generate a set of final products, while taking into account material supply (geological) uncertainty to manage the associated risk. Although simulated annealing has been shown to outperform comparing methods for solving the SSOMC, early performance might dominate recent performance in that a combination of the heuristics' performance is used to determine which perturbations to apply. This work proposes a data-driven framework for heuristic scheduling in a fully self-managed hyper-heuristic to solve the SSOMC. The proposed learn-to-perturb (L2P) hyper-heuristic is a multi-neighborhood simulated annealing algorithm. The L2P selects the heuristic (perturbation) to be applied in a self-adaptive manner using reinforcement learning to efficiently explore which local search is best suited for a particular search point. Several state-of-the-art agents have been incorporated into L2P to better adapt the search and guide it towards better solutions. By learning from data describing the performance of the heuristics, a problem-specific ordering of heuristics that collectively finds better solutions faster is obtained. L2P is tested on several real-world mining complexes, with an emphasis on efficiency, robustness, and generalization capacity. Results show a reduction in the number of iterations by 30-50% and in the computational time by 30-45%.


Introduction
Solving a large-scale industrial mining complex problem is critical for decision-makers in the mining industry. An industrial mining complex is an engineering system that manages the extraction of materials from multiple mines, their processing in interconnected facilities and multiple waste sites to generate a set of products delivered to customers. These related operations require a strategic mine plan that maximizes the net present value (NPV) of mining operations while meeting various complex constraints. Mine scheduling optimization is critical as sub-optimal scheduling can adversely affect profit. Conventional mine planning methods break down the optimization into distinct steps to be solved sequentially. They typically optimize mine production scheduling for individual mines (defining the sequence of block extraction), then the downstream (destination) policies. This decomposition leads to suboptimal results because it does not account for the various interactions between processing, transportation, and destination policies contributing to the value of mined products. Conventional methods tend to ignore material supply uncertainty to reduce the model size to allow the use of commercial solvers, which is known to be the main cause of failure to meet production targets and forecasts [1,2]. Conventional methods also tend to aggregate blocks [3,4,5], leading to homogenization of spatial distribution of metal within an aggregate.
A new framework termed simultaneous stochastic optimization of mining complexes (SSOMC) has been proposed [6,7,8,9] to optimize the system by incorporating all components of the mining complex into a single mathematical model. The SSOMC model maximizes the NPV of the mining complex and includes geological uncertainty to manage the associated risk. This type of uncertainty is modeled using a set of geostatistically (stochastically) simulated scenarios of all relevant mineral deposit attributes [10,11,12,13]. Unlike conventional methods that consider the economic value of blocks, the SSOMC considers the simultaneous extraction of material from multiple mines and optimizes the economic value of the final products and all intermediate interactions. Since exact methods are impractical to use for SSOMC with more than a few thousand mining blocks, researchers have employed metaheuristics and hyper-heuristics to address the SSOMC (as well as mine production scheduling). Most notably, a simulated annealing approach proposed by Goodfellow and Dimitrakopoulos [8] outperforms all comparing methods. However, early performance may dominate recent performance in that the score function only uses a combination of past performance to sample the perturbations to be applied and solution time remains an issue.
Machine learning (ML) has shown promise for addressing small and large-scale operations research (OR) tasks [14,15,16,17,18]. Recent research even uses ML to guide OR algorithms to solve combinatorial optimization problems [19,20,21]. Given the complexity of these problems, state-ofthe-art algorithms rely on specific heuristics to make decisions otherwise computationally expensive or mathematically intractable. Therefore, ML emerges as a suitable predictive tool for making such heuristic scheduling decisions in a more reasoned and optimized manner [22]. In this manuscript, this line of work is followed by proposing, to the best of our knowledge, the first data-driven framework for heuristic scheduling in a fully self-managed hyper-heuristic-based solver. Indeed, in most of the data-driven AI for OR literature [23,22], the authors assume structural knowledge of the scheduling problem (e.g., optimal solutions or strong branching decisions). This assumption allows to apply ML to provide fast approximations of these decisions. Limited research has been conducted to solve large-scale stochastic combinatorial optimization problems such as the SSOMC where no structural knowledge pre-exists. In this case, ML models can be deployed to discover the structure of the problem. This is called policy learning and is an application of reinforcement learning (RL).
The proposed self-learning hyper-heuristic, L2P (learn-to-perturb), improves the search strategy by using past experience to better guide the exploration, thus accelerating the process. This hyperheuristic uses a multi-neighborhood simulated annealing algorithm in conjunction with RL. By defining a neighborhood structure, the RL agent learns to guide the search through the solution space landscape to find better solutions during the optimization process. This paper aims to contribute the following: (1) propose a self-managed hyper-heuristic for solving the SSOMC based on a multi-neighborhood simulated annealing algorithm with adaptive neighborhood search; (2) use RL to adapt the search and guide it towards better solutions; (3) test on several mining complexes with emphasis on efficiency, robustness, and generalization capacity. To test the proposed methodology, three policy-based agents are proposed: advantage actor-critic (A2C) [24], proximal policy optimization (PPO) [25], and soft actor-critic (SAC) [26]. Results demonstrate the effectiveness of L2P on a series of real-world mining complexes, where the proposed solution process takes only minutes, whereas CPLEX solves the linear relaxation in days or weeks. Using RL further reduces the computational time by 30-45%.
The rest of the paper is structured as follows: Section 2 reviews the literature on metaheuristics and hyper-heuristics in the context of the SSOMC. Section 3 presents an overview of the mathematical formulation of the SSOMC. Section 4 details the L2P solution approach. Then, Section 5 presents computational experiments on small-and large-scale instances of the SSOMC, along with comparisons for L2P variants (either using a standard score function or using various RL agents to guide the search), followed by conclusions in Section 6.

Mathematical programming for the simultaneous stochastic optimization of mining complexes
Strategic or long-term mine planning aims to optimize when to extract material from mineral deposits and where to send it under various technical constraints in order to maximize the economic value of extraction. Such planning uses simulated mineral deposits that are discretized into a threedimensional block model where each block is assigned values of spatially distributed attributes of interest (e.g., grade). Since this problem is of utmost importance to the mining industry, many papers address it [27,28,29,30,31,5]. Due to the large size of orebody models and the complexity of constraints within, the industry practice for many years has been to divide the strategic mine planning into several disconnected subproblems and solve them sequentially using approximate methods with strong assumptions, ultimately leading to sub-optimal solutions. Previous research highlights not only the importance of linking the multiple components of the mineral value chain using a single simultaneous optimization framework but also the incorporation of geological uncertainty in a stochastic formulation, which increases the net present value (NPV) by 15-30% [32,33,34,8], compared to assuming a deterministic input of grade distribution.
To overcome the above limitations, the SSOMC [6,7,8,9] accounts for geological uncertainty using two-stage stochastic programming models and links all components of a mining complex and can model non-linear interactions at each destination. Unlike previous work that considers the economic value of a block (due to dividing the problem into several sub-problems to be solved sequentially), the SSOMC considers revenues from final products and all intermediate interactions (blending, various costs along the value chain, any non-linear interactions such as in processors and stockpiles). As such, the SSOMC optimizes not only the block extraction sequence (production scheduling) but also the decision policies, downstream optimization (and potentially transportation), in a single model, capitalizing on the synergies and interactions between the different components.
Since the SSOMC integrates all components of the value chain into one model (from mines to products), the problem is too complex to be solved by exact methods. Thus, it becomes important to use heuristic methods, and mainly metaheuristics have been used to solve strategic mine planning in general, and the SSOMC in particular. Metaheuristics combine (problem-specific) heuristics into a more general framework. Thus, although they do not guarantee a global optimum in the presence of imperfect information and resource constraints, they may find a near-optimal solution.
Fundamental metaheuristics include particle swarm optimization [35] inspired by birds flocking, ant colony optimization [36] inspired by the interactive behavior of ants and simulated annealing [37] inspired by annealing process in metallurgy. Since these metaheuristics have been successfully applied to open-pit mine scheduling (i.e., production scheduling), such as Tabu search [32] and variable neighborhood descent algorithm [38], they were then extended to the SSOMC.
Two main approaches have been developed to address the SSOMC using metaheuristics. The first approach [6,7] proposes to simultaneously optimize the mining, processing, and transportation schedules by maximizing the NPV of the mining complex and minimizing the deviations from the ore and waste production targets. The proposed two-stage stochastic integer program (SIP) formulation is solved using an efficient meta-heuristic, which considers three sequentially executed perturbation strategies: block-based perturbations (modifying the extraction period and destination of each block), operating alternatives-based perturbations (allowing for better control of geometallurgical variables), and transpiration system perturbations (modifying the proportions of materials through the available logistics system to deliver materials from one destination to another).
In the second approach [8,9], the authors propose a formulation of a two-stage nonlinear SIP optimizing three decision variables. In addition to the scenario-independent block extraction sequence and scenario-dependent processing stream decisions (proportions sent from one location to another), the authors use k-means [39,40] to group blocks with similar characteristics (attributes) and send them to the same destination. This clustering allows destination policies to be dynamically changed during optimization and reduces the number of processing decisions that need to be made, in contrast to the blockwise destination SIP optimization in Montiel and Dimitrakopoulos [6] and Montiel and Dimitrakopoulos [7]. Goodfellow and Dimitrakopoulos [8,9] compare three metaheuristic algorithms, including simulated annealing, particle swarm optimization, and differential evolution, achieving an NPV up to 13% higher than standard deterministic industrial software. The framework is particularly flexible, allowing non-linear calculation of attributes (blending, grade recovery function, geometallurgical responses) at any stage of the mining complex.
A relatively recent development in metaheuristic optimization is that of hyper-heuristics [41]. Hyperheuristics have become particularly interesting since many existing approaches rely on metaheuristics, and no single approach works for various types of mining complexes. Hyper-heuristics typically combine low-level heuristics to perturb the current solution, so they operate at a higher level of abstraction than a metaheuristic [41]. Hyper-heuristics have been successfully employed for mine scheduling (but not yet for SSOMC). For example, Lamghari and Dimitrakopoulos [42] use 27 lowlevel heuristics to compare three different hyper-heuristics for mine scheduling. Lamghari et al. [43] combine machine learning (ML), exact algorithms, and heuristic methods to solve the SSOMC and capitalize on the synergies between artificial intelligence and optimization techniques. The authors show that using neural networks as a surrogate model for the downstream subproblem (instead of solving it to optimality) significantly reduces the computational time.

Machine learning for combinatorial optimization
Machine learning has shown promise for tackling small-scale operations research (OR) tasks, such as using graph convolutional neural networks to learn branch-and-bound variable selection policies [15], using hybrid architectures for efficient branching on CPU machines [17]. Moreover, ML has also been used to guide OR algorithms on large-scale tasks, such as using graph neural networks [44], neural networks [14,16,45], and structured convolutional kernel networks [18] to solve the airline crew scheduling problem [46] as well as the vehicle and crew scheduling problem in public transit and the vehicle routing problem with time windows [44,45]. When one has a structural understanding of the scheduling problem (i.e., a theoretical understanding of the decisions of the optimization algorithm), ML models can be used to provide fast approximations of these decisions and thus reduce the required computational time. This is referred to in the literature as "learning by imitation": samples of the expected behavior are available and used as demonstrations for learning in the ML model. As such, the ML model has the task of minimizing a cost function that compares the expert's decisions to its own. The application of ML, in this case, can be straightforward, but the assumption that enough demonstration samples are available for learning is a strong assumption that is not always true.
Recent research goes a step further by proposing data-driven methods that can guide OR algorithms to solve combinatorial optimization problems. To solve these problems, state-of-the-art algorithms use various heuristics to make decisions that would otherwise be computationally expensive or difficult to incorporate into the model. Therefore, ML emerges as a fitting predictive tool to make such heuristic scheduling decisions in a more principled and optimized manner [22]. As such, following Katz et al. [47], Lipets and Zadorojniy [21] use expert knowledge (also known as domainknowledge information) to generate a transition probability matrix [48] used to define a Markov decision process (MDP). The MDP was used as an orchestrator for scheduling heuristics to solve a complex airline crew scheduling problem (crew pairing and rostering). Concurrently, Zarpellon et al. [19] present a novel imitation learning framework and hypothesize that the state parameterization of the Branch-and-Bound search tree can help solve mixed-integer linear programming (MILP) problems. Chmiela et al. [20] propose the first data-driven framework for heuristic scheduling in an exact mixed-integer programming (MIP) solver, where they obtain problem-specific heuristic scheduling that collectively finds many solutions at minimum cost by learning from data describing the performance of primal heuristics.
This line of work is followed by proposing the first data-driven framework for heuristic scheduling in a hyper-heuristic solver that is fully self-managed. Since the SSOMC has no pre-established structural knowledge, ML models can be deployed to discover the problem structure using RL. The model can perform a set of actions in each state, with each action having a reward; the RL model then tries to maximize the expected sum of the rewards. A common challenge in implementing this approach is to find an appropriate reward function so that the learning algorithm does not fall into local minima or terminate without sufficiently exploring the search space. In this manuscript, the use of hyper-heuristics is justified by the extremely large size and complex constraints of the SSOMC and the demonstrated success of simulated annealing [49,37] and hyper-heuristics [1,50] in solving the SSOMC. Indeed, following the work of Goodfellow and Dimitrakopoulos [8,9] and Lamghari et al. [43], a self-learning hyper-heuristic called L2P (learn-to-perturb) is proposed.

Problem statement
This section presents an overview of the mathematical formulation [8,9,51] and the solution approach [8,42] that are used to flexibly model varying mining complexes and solve the SSOMC.
A mining complex is a mineral value chain, where materials extracted from the ground (multiple mines) flow through different processing streams and encounter multiple transformations into final sellable products that are delivered to the mineral market. The goal of the optimization of mining complexes is to provide the market with products, where the supply of materials extracted from the mines is characterized by the uncertainty and variability of the geological attributes.
This uncertainty is incorporated into the framework as a group of simulated orebody models with their respective relevant attributes, in the form of geostatistical simulations that generate multiple equiprobable scenarios of the attributes of interest (e.g., grades, densities, material types, etc.) As in Figure 1, the material extracted from mines flows through many different processes in the mining complex, propagating the uncertainty to all operations in the mining complex, which culminates in financial risk. In particular, a mining complex typically contains a set of mines (orebodies), stockpiles (to store incoming materials and meet blending constraints), and processors (to transform incoming products). As such, the SSOMC aims to optimize extraction sequences, destination policies, processing stream decisions, transportation alternatives altogether in one single mathematical model. This framework can be seen as a series of transfer functions that account for the transformations of in-situ material throughout the value chain, which can include various non-linear interactions (blending, non-linear recovery curves, stockpiling, etc.) In order to optimize the entire mining complex, the problem can be modeled using a two-stage SIP. The objective function maximizes the expected NPV of the mining complex (total expected discounted profit generated from processing the ore) and minimizes the expected recourse costs incurred when the stochastic (scenario-dependent) constraints are violated. Scenario-independent constraints include reserve constraints (ensuring that each block is mined at most once during the horizon), slope and precedence constraints (ensuring that each block can only be mined after all its predecessor blocks have been mined), and mining constraints (imposing a maximum amount of material to be mined per period). Scenario-independent constraints typically include processing constraints (imposing a minimum and maximum amount of material to be processed by each processor per period) and stockpiling constraints (balancing the flow at each stockpile and imposing a maximum amount of ore material to be stockpiled for each stockpile per period).
Due to the exceptionally large size that the model in question can reach and allow for the modeling of varying mining complexes, a generic modeling and solution approach was proposed by Goodfellow and Dimitrakopoulos [8,9]. To account for the flow of materials through the mining complex, block attributes are considered to be either primary or hereditary. Primary attributes are those where variables are sent from one location to another (e.g., metal tonnage). Note that these variables can be added directly. Hereditary attributes are variables that are not passed between locations in the mining complex (e.g., costs). These cost functions are assumed to be linear in order to simplify the solution process and to compare the results to the optimal solution (of the relaxed problem) when possible. In the general case, (non)linear formulations can be assigned to these hereditary attributes and are dynamically evaluated during the optimization process.
The flow of materials (and their respective attributes) through the mining complex is defined by three sets of decision variables (which the simultaneous stochastic optimizer can modify), namely mine production scheduling decisions, destination policy decisions, and processing stream decisions.
Production scheduling decisions define whether or not a certain mining block is extracted from a given mine in a certain production period. These variables are scenario-independent.
Once mining blocks are extracted, it is necessary to decide where to send each block. As a preprocessing step, multivariate distribution clustering divides blocks into groups based on multiple elements (e.g., grades) in each material type. K-means [39,40] is used to define the membership of a certain simulated orebody block for a certain scenario to a given cluster. These membership variables form a robust destination policy that decides where to send all blocks in a cluster (with similar attributes), rather than deciding the destination of individual blocks. Note that a block can be sent to different destinations in various scenarios, but the destination policy decisions are scenario-independent. The optimizer is responsible for deciding where each group (cluster) is sent per period. The use of these multivariate clusters makes the destination decisions more adept at creating destination policies for multi-element mining complexes.
Finally, the flow of materials through the rest of the mining complex is defined using processing stream decisions that define the proportion of a product sent from a location to a destination. These scenario-dependent variables act as recourse actions. For a detailed formulation of the model, see Appendix A. In short, it is a two-stage SIP, where first-stage decisions are mine production schedule and destination policies, and second-stage decisions are the processing stream decisions (and penalties for deviations from production targets). The objective function maximizes discounted cash flow and minimizes discounted risk-based deviation penalties, subject to the following constraints: (1) mine reserve and slope constraints, (2)

Solution methodology
This section presents the learn-to-perturb (L2P) solution approach that employs simulated annealing, Tabu search [52], and RL to solve the SSOMC. The proposed self-learning hyper-heuristic is constructed by hybridizing a multi-neighborhood simulated annealing algorithm with RL. A probabilistic approach aims at identifying a global optimum by identifying the neighboring states with lower energy or closer to an optimal solution in the solution space.

General structure of the learn-to-perturb (L2P) hyper-heuristic
Rather than directly solving an optimization problem, the proposed L2P hyper-heuristic attempts to identify the best heuristic for a given problem state and compensate for the weaknesses of different heuristic methods through simultaneous optimization [41]. Indeed, while a set of low-level heuristics (perturbations) work on finding optimality simultaneously, a high-level heuristic schedules their sequencing (ordering), ensuring that the best solution is found by combining the perturbations.
The proposed L2P hyper-heuristic is illustrated in Figure 2 (modified from [41,53,54]) and is divided into three components that are used simultaneously (rather than in sequential order): Simulated Annealing Acceptance Criterion Best solution x Collect problem independent information for domain barrier (e.g., the number of heuristics, the changes in evaluation function, a new solution or not, the distance between two solutions, etc.)

Simulated Annealing
Domain Barrier • Problem domain: includes the problem representation, the evaluation function, the initialization procedure, and the set of low-level heuristics. This is the lower block of Figure 2, and is based on Goodfellow and Dimitrakopoulos [9] and Goodfellow [51].
• Simulated annealing algorithm: consists of a stochastic heuristic selection mechanism, an acceptance criterion and a sampling function, modified from Lamghari and Dimitrakopoulos [42], Goodfellow and Dimitrakopoulos [9], and Goodfellow [51] • RL algorithm: predicts the future performance of low-level heuristics to guide the search towards better solutions. This self-managed component learns from past behavior to provide data-driven guidance for searching over the solution space landscape.
Following the work of Lamghari and Dimitrakopoulos [42], Goodfellow and Dimitrakopoulos [9], and Goodfellow [51], the pseudocode for the L2P hyper-heuristic is presented in Algorithm 1. The search is performed using a single candidate solution [42]. At each iteration, L2P improves the current solution x by iterating between heuristic selection and move acceptance (see boxes 1-3 of Figure 2). More precisely, once a heuristic is selected (see box 1), a candidate solution x is modified into a new solution using said heuristic (see box 2). The heuristic is chosen based on its past performance [42] and the expected performance of the RL agent, and simulated annealing [49,37,55] handles the heuristic selection mechanism and the acceptance criterion.
The algorithm starts with one solution, then moves to a neighboring solution in search of a better solution. This process is repeated until no significant improvement is possible. Heuristic selection is biased toward heuristics that have performed well in previous iterations, and it learns to direct solution sampling strategies to more promising regions of the heuristic space. Specifically, each time a heuristic h i (i = 1, . . . , n, where n is the number of low-level heuristics) is used, the algorithm returns the resulting change in the objective function value ∆f (h i ) and time required T (h i ). A heuristic is given more importance if it improves solutions quickly per unit time and it is considered less important if it deteriorates the solution and is more computationally expensive. Then, heuristic selection is performed based on the tabu status and performance of the heuristics.
To evaluate the tabu status, if the heuristic does not improve the current solution, it becomes tabu and is temporarily removed from the selection set for γ T iterations where γ T is randomly generated in the interval [Γ min , Γ max ], and the Tabu list is emptied when all heuristics are made tabu. This strategy allows for better exploration of the heuristic search space because it significantly reduces the probability of repeatedly choosing the same heuristics while leaving the possibility of not banning non-improving heuristics for too long by tuning Γ min and Γ max [32,42].
To evaluate the performance of (non-tabu) heuristics, L2P sampels towards the heuristics that are most likely to improve the objective function, by assigning a selection probability p i computed by normalizing the score function p i = In what follows, the computation of the score function S F is explained. Every ζ iterations, two measures π 1 (h i ) and π 2 (h i ) are updated for each heuristic h i using Eq. (1) and Eq. (2), respectively [42]. Then, S 1 is updated using Eq. (3), where η(h i ) is the number of times heuristic h i has been selected in the last ζ iterations, and α and β are two weight adjustment parameters in [0, 1] defining the importance given to recent performance and to improving heuristics, respectively. The value of β is initialized to 0.5 and is changed depending on whether a new best solution has been found in the last ζ iterations.
Input: if x is better than x best then x best ← x 6: Initialize the score function S 1

7:
Remove hj from the list L H , and set x ← x Stage 2: Selecting the heuristics based on their scores and tabu status 8: Initialize π 1 , π 2 , and η, and add all heuristics to a list L H (list of heuristics that are non-tabu) 9: while stopping criterion not met do 10: Generate a candidate solution x from x using (h i ) box 2

13:
Evaluate ∆f (h i ) and T (h i ) 14: if x is better than x best then Set x best ← x and newBest ← true Generate a random number γ T in [Γ min , Γ max ] and make h i tabu for γ T iterations 19: if if iter mod iter cooling = 0 then Temp ← Temp * k 23: if iter < ζ then iter ← iter +1 24: else 25: if newBest = true then β ← 1 26: Normalize score functions S 1 , π 1 , and π 2

32:
Normalize score function S F , revoke tabu status of all heuristics (reset L H ) box 6

33:
Set iter ← 1 and newBest ← false 34: return x best The computed score function S 1 is used along with S 2 (a score function provided by the RL agent) to calculate the final score function S F . To compute S 2 , all heuristic score measures for the last L w rounds are provided to the RL agent as input s t (in the form of a 3D matrix). As in box 5 of Figure 2, the agent outputs the action a t = π θ (s t ), where π θ is the agent's policy that is parameterized by the weights θ. The score function S 2 is obtained by normalizing a t . In order to guide the L2P hyper-heuristic to identify the best solution, the RL agent takes a particular action (a t ) in each state and transitions from one state to another, and based on the gain or loss occurring, the RL algorithm updates its policy weights in the direction that is likely to maximize the gain (reward) [22]. Once S 1 and S 2 are computed, a voting mechanism (box 6 in Figure 2) is used to compute the final score function S F based on S 1 and S 2 , as in Eq. (4), where λ RL is the RL agent's contribution to the final score function. While it is conceptually possible not to use voting, it is essential for the agent to learn good policy and find good solutions. The voting mechanism combined with sampling typically promotes exploration in RL to ensure the most viable outcome.
Note that the L2P variant with λ RL = 0 will be considered as a baseline in Section 5.
Using the computed S F to sample the heuristic to apply, the solution vector is updated. The reward r for the agent performing action S 2 is ∆f (h i ) (the change in the objective function value). The agent's performance is ultimately guided by the policy changes and experience gained by the agent over time. At each iteration t, before taking the action, the agent stores past experience in a replay buffer. Specifically, as in box 4, the tuple (S t−1 , a t−1 , r t , S t ) is saved in the buffer, which is used to update the agent's policy. And every few iterations (set by a hyperparameter), the agent updates its weights θ (see box 5). Note that as the algorithm progresses, the temperature of the simulated annealing selection mechanism is gradually reduced until only minor changes are accepted. The temperature is controlled by the initial temperature (Temp 0 ), and a cooling schedule, which is defined by a reduction factor, k ∈ [0, 1), and a number of iterations before the reduction factor is applied (see Algorithm 1). Note that the is repeated until a predefined stopping condition is met and only problem-independent information flow is allowed between the problem domain and other layers, as is standard in hyper-heuristics [41]. An pseudocode for L2P is given in Algorithm 1.

Reinforcement learning algorithm
Reinforcement learning (RL) is a goal-oriented algorithm that learns to achieve complex goals and maximize the objective towards a certain dimension, such as winning a game by accumulating the highest possible score in a certain number of moves. In this manuscript, the potential-based reward shaping [56,57] is used. The potential is obtained with a model-free RL agent (as opposed to a model-based agent) via policy iteration. Policy iteration is used as a form of optimal control to guide the search in the heuristic space, iterating between policy evaluation and improvement. Basically, the current policy is evaluated for each iteration, and the new cost is used to obtain a new improved Notation Description s ∈ S Set of states a ∈ A Set of actions r ∈ R Set of rewards s t , a t , r t Notation for state, action, and reward at a given time t γ Discount factor between 0 and 1, can be described as a penalty imposed due to the uncertainty of future rewards R t Return at time t, R t = ∞ i=0 γ i r t+i+1 prob(s , r|s, a) Probability of transitioning from current state s the next state s by taking the action a and receiving the reward r π θ (a|s) policy π parameterized with θ mapping each state with its optimal action V (s) State value function V w (.) State value function with parameter w V π (s) State value function while following the policy π for a state s; Action value function computing the expected return from (s, a) using policy π parameterized with parameter θ. It can also be referred to as Q(s, a) A(s, a)  Table 1.

The Markov Decision Process Formulation
A RL agent starts from an initial state and updates its weights in the direction of improving the policy, towards an optimal policy that yields the best possible performance. The agent's learning and state change mechanism are similar to those of humans in that it can be penalized for a series of incorrect decisions/actions and rewarded for correct decisions/actions. In what follows, the policy is denoted by π, states by s, actions by a, the discount factor by γ, and rewards by r. Figure 3, in a typical RL environment, an agent performs an action a according to policy π to transition from the current state s to the next state s . The probability prob(s , r|s, a) is the probability of transitioning from the current state s t to the next state s ∼ s t+1 at time t by taking action a for observation/tuple (s t , a t , r t , s t+1 ). The return R t at t is:

As in
On the one hand, for value-based methods, R t is computed to estimate π * that maximizes the expected return E π [R t ] 2 . This is achieved by exploiting value-based solutions, thus computing V π (s) using Eq. (5), or by computing Q π (s, a) in Eq. (6) by exploiting the temporal-difference learning to generate a series of greedy policies by choosing at each time-step the action: arg max a Q π (s, a ) [59]. Figure 3: An illustration of the interaction between the RL agent and the environment [58]. In the proposed L2P, the agent refers to the reinforcement learning algorithm (box 4-5 in Figure 2).
On the other hand, policy gradient methods typically model and update the policy to maximize the rewards. The quality of a policy is measured using the long-term cumulative reward function (formally known as average reward per time-step), as in Eq. (7), where ρ π (s) denotes the stationary distribution of the Markov chain for π θ . Thus, gradient descent is employed on Y (π θ ) to update the weight θ parameterizing the policy π. The agent outputs the normalized action (a t = π θ (s t )), where π is the agent's policy parameterized by the weight θ. The output action represents the predicted performance of each heuristic and is guided by updating the parameter θ of the stochastic policy π θ and the agent's exploration during training. The stochasticity comes from adding a Gaussian noise N (0, σ) to the (deterministic) output action a t to obtain the score function S 2 . This is used in continuous stochastic control problems (as in this case) to enrich the exploration behaviors. This process can be thought of as adding randomness to the output action in order to explore state and action spaces and thus escape local optima. Although σ can be learned during training, it is used as a hyperparameter in this case to avoid an overconfident agent (early convergence of the policy).

State and Action Encoding
After computing the score function S 2 , a voting mechanism is used to compute the final score function S F based on S 1 and S 2 , as formulated in Eq. (4), where λ RL denotes the contribution of the RL agent to the final score function that will be used for the next optimization iteration. In the next iteration, a heuristic h i is chosen using the computed score function S F , and a new (updated) solution vector is obtained, along with the resulting change in the value of the objective function ∆f (h i ) and the time required T (h i ). Next, the reward given to the agent is computed. Defining a reward is a rather difficult step in designing an RL framework. On the one hand, although defining a sparse reward is easy to encode (e.g., 1 when the problem is solved, 0 otherwise), this setting does not provide a very limited useful learning experience for the agent. On the other hand, when the reward is approximated (e.g., by a linear function approximator instead of solving a downstream process at each iteration), the learning is more likely to fall into local minima [22]. In this manuscript, the reward for action a t is defined as r t = ∆f (h i ). Note that to allow generalization for the learned policy across instances, the computed returns are standardized (transformed to have zero mean and unit standard deviation) before using them to update the agent's policy. In addition to ∆f (h i ), it is also conceptually possible to consider T (h i ) when computing the agent's reward r t . However, this may lead the agent to suggest perturbations that take very little time but provide no improvement in the value of the objective function, resulting in a sub-optimal policy.

Advantage Actor-Critic (A2C)
In policy-based methods, the policy is parameterized by θ (weights for action selection) and w (weights for state evaluation). These parameters are typically updated by gradient descent over E(R t ). For example, in the REINFORCE algorithm [60], one of the first policy-based methods, the parameters θ and w are updated in the direction of ∇ θ log π(a t |s t ; θ), which can be considered an unbiased estimate of ∇ θ E(R t ). The main limitation of policy-based methods is that this estimate can have a high variance, which hinders convergence to an optimal policy. To reduce the variance of the estimate while keeping it unbiased, a learned function b t (s t ) is subtracted from the reward estimate. The resulting equation that defines the gradient is formulated as follows: The learned function is referred to as the baseline (b t (s t ) ≈ V π (s t )). In the case where the approximate value function is used as the baseline function, the scaling quantity (R t − b t ) of the gradient can be described as A(a t , s t ) = Q(a t , s t ) − V (s t ) where A(a t , s t ) is referred to as the advantage of an action a t in a state s t , since R t is an estimate of Q(a t , s t ) and b t is an estimate of V (s t ).
The advantage evaluates the quality of an action with respect to the average action for a specific state. This method is called the actor-critic method since the policy π acts as an actor and the baseline b t acts as a critic [58,61]. The working of A2C is detailed in Appendix B. A2C maintains a policy π(a t , s t ; θ) and the state-value function V (s t ; θ s ). At each iteration, the algorithm updates the policy and state-value function using n-step returns, for example after reaching a terminal stage or after t max actions. The update can be thought of as ∇ θ log π(a t |s t ; θ )b(a t , s t ; θ, θ s ) where b can be described as an estimate of an advantage function described in Eq. (8), where the variable j varies from state to state and is bounded by t max .
b(a t , s t ; θ, θ s ) = The state value function parameter θ s and the policy parameter θ can be shared or completely separate (independent). In our implementation, these two networks share all hidden layers except the last hidden layer and the output layer, in order to speed up learning, while allowing different learning speeds for the actor and value networks. Furthermore, in order to avoid premature convergence, the policy entropy is added to the objective function. The resulting gradient ∆ θ is expressed in Eq. (9), where H denotes the value of the entropy and is a regularization parameter.

Soft Actor-Critic (SAC)
The soft actor-critic (SAC) [26] is an off-policy actor-critic algorithm that feeds the policy entropy measure H to the reward for robust state exploration while incorporating the maximum possible policy randomness. The training maximizes both H and expected return E π [R t ] 2 . By using the hyperparameter to control the entropy intensity, the reward can be defined as in Eq. (10). Using entropy maximization allows one to choose among optimal and equally good strategies to assign equal probability, and allows further exploration of near-optimal strategies.
In practice, SAC learns three functions: (i) the policy π with parameter θ (π θ ), (ii) the soft Q-value function Q parameterized by w (Q w ), and (iii) the soft state-value function V parameterized by ψ (V ψ ). Using the Bellman equation, the soft Q-value Q(s t , a t ) and soft state value V (s t ) can be defined as in Eq. (11), where ρ π (s) refers to the state and ρ π (s, a) refers to the state action pointers of the state distribution provided by the policy π(a|s).
First, to update the soft state function (Eq. (11)), the MSE (Eq. (12)) is minimized, resulting in the gradient in Eq. (13), where ∇ ψ Y V (ψ) is the gradient and Z denotes the reply buffer used.
Second, to update the soft-Q-function, the Bellman residual is minimized as in Eq. (14), where ∇ w Y Q (w) is the gradient and Vψ(s t+1 ) denotes the target value function. This estimator is updated periodically as an exponential moving average.
The update ensures that Q πnew (s t , a t ) ≥ Q π old (s t , a t ), where Π denotes the set of potential policies that can be used to model our policy. D π old (s t ) can be thought of as an intractable partition function that is used to normalize the distribution, but it never contributes to the calculation of the gradient.
The working of the SAC algorithm is detailed in Appendix B.

Proximal Policy Optimization (PPO)
The motivation behind incorporating proximal policy optimization into L2P comes from the fact that policy gradient algorithms can be unstable in continuous control problems, as is the case for the current RL environment. In order to improve the training process and stability and to avoid diverging from the optimal policy, the trust region policy optimization (TRPO) Schulman et al. [62], used the Kullback-Leibler divergence constraint [63] on the policy so that difference between the updated policy and the previous policy is relatively small.
The clipped advantage function becomes as in Eq. (17), where the clipping function CLIP(r(θ), 1 − , 1 + ) retains the ratio between (1 + ) and (1 − ). Note that the estimated advantage function W(.) is used here instead of true advantage function A(.) because true rewards are unknown in this case. The objective function to be maximized in order to update the weights θ becomes as in Eq. (18). The objective function employs an entropy term for sufficient state exploration and p 1 , p 2 are hyper-parameters controlling the exploration/exploitation tradeoff [25].
Following Hsu et al. [64], to avoid instability in continuous action spaces and (over)sensitivity to initialization (e.g., if the optimal local actions are close to initialization), the probability ratio is clipped to regularize the policy and a continuous Gaussian distribution based on the parameterized policy is used as in Section 4.2.2. The working of our PPO implementation is shown in Appendix B, where the surrogate loss is optimized using the Adamax optimizer [65] (or any similar optimizer).

Problem domain component
To produce new solutions, the hyper-heuristic described in the previous section uses 38 simple lowlevel perturbative heuristics from Goodfellow and Dimitrakopoulos [8,9], Goodfellow [51], referred to as h i (i = 1, . . . , n, where n is the number of low-level heuristics). Each heuristic aims to improve the current solution and examines a subset of one of the following four neighborhoods: 1. Block extraction sequence perturbations [8,9]: For each of this neighborhood's heuristics, a block is randomly selected. The block's extraction period is changed, according to the heuristic's rules, and any potentially violated slope constraints are fixed. Specifically, (1) the current block that is visited is set as the central block. Note that a predecessor of block i is a block that must be extracted to have access to i and block j is a successor of block i if and only if i is a predecessor of j. If the extraction period of a block i is advanced (resp. delayed), this process is equivalent to the predecessors (resp. successors) of i being moved to the same new extraction period as i. The set formed by a block i and its predecessors (resp. successors) mined in the same period is called the inverted cone (resp. cone), whose base (resp. apex) is i.

2.
Cluster destination policy perturbations [8]: A cluster is a group of blocks in a mine that can be grouped together based on certain characteristics. Clusters can be formed by grouping blocks with similar block attributes. The membership of a block in a cluster can change throughout simulations as block characteristics vary in different simulations. The overall optimal solution is essentially the sum of the optimal results from all simulations (scenarios). In the heuristics of this neighborhood, a cluster destination decision variable is randomly selected and sent to a different destination, if possible.
3. Destination policy perturbations [51]: A destination policy is often employed to group blocks into bins of ore and waste blocks. In this way, optimization is limited to looking at a collection of blocks instead of individual blocks. For example, in a time-varying grade policy, each block characteristic is assigned a value below which mining becomes infeasible. These heuristics will dynamically vary the cut-off range characteristics of clusters when clusterrelated perturbations occur, allowing for further optimization.

4.
Processing stream perturbations [9,51]: Open-pit mines can operate with one or multiple process streams, where materials are extracted in bulk, processed, and sent to specific destinations. To find these destination policies, the process stream variables are optimized to maximize the overall economic value of the operation. Process streams can have material quality constraints, although the current model uses the cut-off grade constraint. The transfer of materials between two locations is typically dependent on the process stream configuration and is, therefore, an integral part of the modeling and optimization process. A random normal number N (y i,j,t,s , 0.1) is used to change a randomly selected process stream variable. Keeping the distribution's variance small enough facilitates both global and local search.

Experimental setting, benchmark instances and parameter setting
To evaluate the performance of the proposed L2P hyper-heuristic, numerical experiments on five instances (of the simultaneous stochastic optimization of mining complexes) of various sizes and characteristics are reported, divided into three different stages of testing. The benchmark instances are briefly described below and detailed in Table 2 In the first stage of testing, three instances (mining complexes) I 1 -I 3 are used. This stage is intended to make an initial evaluation of the proposed L2P hyper-heuristic with emphasis on the generalization capability of the proposed method. The instances are small size with up to 15,000 blocks for one mine or up to 4,000 blocks per mine for two mines. Note that online learning (weight updating) is disabled in this stage, so the resulting performance is only obtained from the pretraining stage. Production, processing, and stockpile capacities are fixed for all periods/years. I 1 contains a copper mine (one type of metal), a processor, and a waste dump. I 2 contains two mines with two types of copper materials (oxides and sulfides), two processors (oxide leach pad and sulfide processor), and a waste dump. I 3 is similar to I 2 , except that it has ×2 times as many blocks and a sulfide stockpile. Each time, the method is initialized, pre-trained on one of the instances, and then tested on all three. I 3 is similar to I 2 , except that it has ×2 times as many blocks and a sulfide stock. To test, each L2P variant is trained using one instance, and then tested using all three instances. Since the L2P starts with a random initial solution and the heuristic selection mechanism is stochastic, all experiments are run 10 times and the results in this stage are reported as mean and standard deviation. To compare the different L2P variants, since finding the optimal value for any of the instances is not possible, the linear relaxation of the model is used, which leads to a weak bound on the objective value of the optimal solution. A time limit of four weeks is set to solve the linear relaxation for each instance using CPLEX 12.10.
In the second stage of testing, the agents are trained and tested on the large-scale instance of the SSOMC I 4 . The instance I 4 is similar to I 1 except that it has ×14 more blocks and imposes a lower processing limit, and this stage is intended to evaluate the performance of the L2P hyper-heuristic on a realistic case study. Each L2P variant is trained using I 4 , then tested on the same instance I 4 .
Since CPLEX was unable to solve the linear relaxation of I 4 within the time limit (four weeks), the objective value of the best solution found (after several trials) is therefore used as a comparison.
This point of reference is referred to as Z * . In this stage, all variants of the L2P hyper-heuristic will be compared based on the evolution of the objective function value (with respect to Z * in %) as a function of the number of iterations and the execution time (in minutes), in the form of confidence intervals for the estimated P10, P50, and P90 quantiles.
In the third stage of testing, the large-scale instance I 5 is used to evaluate the performance of the L2P hyper-heuristic on a realistic case study. The instance I 5 is similar to I 3 except that it has ×75 times more blocks. Each L2P variant is pre-trained on the instance I 4 , then the pre-trained weights are transferred and the variant is tested on the instance I 5 . Note that for both the second and third stages of testing, online learning is enabled as it would be for an end-user. The only difference with the second stage of testing is that the pre-training (warm-starting the policies) is performed here on a completely different instance.

First stage of testing
Using the value of the obtained linear relaxation, Tables 4-5 report the execution time (T gap=1% ) and the number of iterations (Iter gap=1% ), respectively, required by the variants to obtain a solution whose objective function value is 1% away from the linear relaxation value. The last row of Table 4 reports the time required by CPLEX to solve the linear relaxation. The best results obtained for each instance are indicated in bold.
The the heuristic choice and performance may vary slightly from instance to instance, which is a good  indicator of the generalization capacity of the proposed method, especially since the weights are not updated in real-time at this stage of testing (i.e., once the weights are pre-trained and transferred, no online learning is used). This difference in performance may only become smaller in an actual use case scenario, where online learning is enabled.
Results also show that although L2P-A2C outperforms the other variants, its performance is not as consistent across runs as L2P-PPO (the second-best performing variant). In particular, the standard deviations (both for Iter gap=1% and T gap=1% ) of L2P-A2C are larger than those of L2P-PPO. This may be explained by the fact that PPO takes much smaller steps to update its policy, and thus presents a more stable strategy than A2C (albeit less efficient). In addition, the standard deviation of L2P-Baseline is smaller than the standard deviation of all other variants. This is because, although the choice of heuristic is the result of a sampling process, the score function S 1 is biased towards the best performing heuristic, and S 2 might be more biased towards exploration, especially if the state given by the recent performance has not been seen before. Thus, using the score function S 1 solely may produce more "deterministic" (and less noisy) behavior.
Due to the large size of the instances in question, CPLEX takes much longer (150-250 times) to solve the linear relaxation than any of the L2P hyper-heuristic variants take to solve the instances.
Finding integer solutions using CPLEX would be impossible, let alone adding more components to the mining complexes presented in this manuscript. In summary, the best performing variant of the L2P hyper-heuristic is L2P-A2C, both for the instance on which it was trained and for the instances on which it was not trained. The RL-based variant with the smallest standard deviation is L2P-PPO, and the least-performing RL-based variant is L2P-SAC.

Second stage of testing
As shown in Figure  First, SAC uses entropy regularization to prevent premature convergence. This can make SAC take much longer to learn compared to A2C and PPO. Second, the trade-off between exploration and exploitation within SAC is controlled by the coefficient α H , which needs to be changed from one environment to another and requires careful tuning. This makes SAC difficult to include in an off-the-shelf solver with little to no pre-tuning and pre-training.
To compare the solutions obtained for the instance I 4 , Figure 5 shows the NPV and recovered copper for the solution obtained using L2P-A2C (the best performing variant of the L2P hyper-heuristic) and L2P-Baseline. The stopping time criterion used is equal to the execution time it took L2P-A2C to yield a solution value that is 1% away from Z * . The cash flow and recovered metal for the solution obtained using L2P-A2C are 2-3% higher than those obtained using L2P-Baseline. This is only observed when comparing L2P-Baseline and L2P-SAC and can be explained by the fact that the SAC policy is non-optimal and minimally differentiates states during optimization.
Although such non-optimality could be addressed by pre-training the L2P-SAC on I 5 instead of using pre-trained weights, it is crucial for users that the pre-trained solver be as generalizable as possible so that no tuning is required. Thus, using L2P-SAC in a real-life scenario is impractical.
In contrast, L2P-A2C and L2P-PPO emerge as better variants of the L2P hyper-heuristic.
To compare the solutions obtained for the instance I 5 , Figure 7 shows the NPV and recovered copper for the solution obtained using L2P-A2C (the best performing variant of the L2P hyper-heuristic) and L2P-Baseline. The stopping time criterion used is equal to the execution time it took L2P-A2C to yield a solution value that is 1% away from Z * . The cash flow and recovered metal (for the solution obtained using L2P-A2C) are 2−3% higher than those obtained using L2P-Baseline.

Conclusions
This paper presents a self-managed hyper-heuristic based on multi-neighborhood simulated annealing with adaptive neighborhood search, where state-of-the-art reinforcement learning agents are used to adapt the search and guide it towards better solutions. By defining a neighborhood structure, the reinforcement learning component learns to guide the search through the solution space landscape during the optimization process, while sampling towards the most promising heuristics.
We apply the proposed L2P hyper-heuristic to the simultaneous stochastic optimization of mining complexes (SSOMC), a large-scale stochastic combinatorial optimization problem of particular interest to the mining industry. Results comparing different variants of the L2P hyper-heuristic (with and without reinforcement learning) show its effectiveness on several case studies, reducing the number of iterations by 30-50% and the computational time by 30-45%, compared to the baseline using a standard score function. Since the SSOMC is similar to other scheduling problems, exploring this research direction to find more efficient methods to solve such problems becomes even more justified. Indeed, unlike problem-specific heuristics, the presented solution method can be generalized to other problems as in airlines (airline crew scheduling), bus (shift scheduling), trucking (routing problems), and rail industries. The only requirement would be to implement a set of low-level heuristics, none of which would need to be well-tuned.

Appendix A. Problem formulation of the two-stage stochastic integer program
In what follows, the mathematical formulation of the simultaneous stochastic optimization of mining complexes (SSOMC) [8,9] is outlined. All sets and indices used in the formulation are presented in Sets/Indices Description p ∈ P primary attributes, sent from one location in the mining complex to another (e.g., metal tonnages, or total tonnage) h ∈ H hereditary attributes, variables of interest at specific locations (nodes) in the mining complex, but not necessarily forwarded between locations (e.g., mining, stockpiling and processing costs) i ∈ S ∪ D ∪ M locations, within the set of mines M, stockpiles S, and other destinations D b ∈ B m block b within the set of blocks B m at mine m ∈ M s ∈ S a set of equally probable scenarios (simulations) describing the combined simulations from different sources of uncertainty (e.g., geological uncertainty) t ∈ T time period within the life-of-mine g ∈ G a group G is generated as a pre-processing step using clustering to reduce the number of policy decisions variables based on multiple elements within each material type.

O(b)
Block extraction precedence relationships. θ b,g,s Block group memberships O(i) and J (i), ∀i ∈ S ∪ D ∪ G A model of the mining complex, where O(i) is the set of locations that receive materials from location i and J (i) is the set of locations in the mining complex that send material to a location i. f h (p, i) A model of the hereditary attribute transformation functions p h,i,t Time-discounted price (or cost) per unit of attribute (typically a discount rate that is used to calculate the NPV) U h,i,t and L h,i,t Upper-and lower-bounds for each attribute (typically tonnage, metal production and product quality constraints).