Adaptive Population-based Simulated Annealing for Uncertain Resource Constrained Job Scheduling

Transporting ore from mines to ports is of significant interest in mining supply chains. These operations are commonly associated with growing costs and a lack of resources. Large mining companies are interested in optimally allocating their resources to reduce operational costs. This problem has been previously investigated in the literature as resource constrained job scheduling (RCJS). While a number of optimisation methods have been proposed to tackle the deterministic problem, the uncertainty associated with resource availability, an inevitable challenge in mining operations, has received less attention. RCJS with uncertainty is a hard combinatorial optimisation problem that cannot be solved efficiently with existing optimisation methods. This study proposes an adaptive population-based simulated annealing algorithm that can overcome the limitations of existing methods for RCJS with uncertainty including the premature convergence, the excessive number of hyper-parameters, and the inefficiency in coping with different uncertainty levels. This new algorithm is designed to effectively balance exploration and exploitation, by using a population, modifying the cooling schedule in the Metropolis-Hastings algorithm, and using an adaptive mechanism to select perturbation operators. The results show that the proposed algorithm outperforms existing methods across a wide range of benchmark RCJS instances and uncertainty levels. Moreover, new best known solutions are discovered for all but one problem instance across all uncertainty levels.


I. INTRODUCTION
The mining supply chain in Australia is one of the largest in the world, where a significant proportion of the minerals mined are exported overseas.Mining is major contributor to the Australian economy, with the exporting minerals contributing around $AUD 200 billion per annum to the economy, employing 240,000 individuals directly, and around 850,000 individuals indirectly [1].Due to its role in the export of raw materials, Australia has developed a reputation in building mining technologies and is one of the largest suppliers of these technologies and software, worldwide. 1here are numerous problems within the mining supply chain.Of significant interest, is the transportation of ore from mines to ports.In particular, transporting ore from multiple mining sites to a port where resources (trains and trucks) must be shared, is a key problem.If this problem is solved effectively, it can lead to large increases in throughput, thereby leading to substantial gains in revenue for the mining companies.In an abstract form, this problem can be formulated as a resource constrained job scheduling (RCJS) problem.The transportation of ore in batches can be viewed as jobs, while different transport modes (trains or trucks), can be viewed as resources.In the real setting, certain batches are required to be loaded on ships before others, and this aspect can be modelled by precedences between jobs.The timely arrival of ore at the ports is very important, as batches of ore that unnecessarily wait will lead to demurrage costs.This aspect is incorporated as a weighted tardiness for each job (different batches of ore being delayed lead to different costs), and the aim is to ensure the total weighted tardiness across all jobs is kept at minimum.
In the literature, this problem is known as RCJS with a shared central resource [2], [3], [4].The complexity of RCJS makes it challenging for optimisation algorithms, and hence, several approaches have been proposed to tackle it.These include metaheuristics and matheuristics [4], [5], [6], [7], [8], [9], [10], [11], [12], [13].The problem was introduced by [4], who proposed a Lagrangian relaxation based heuristic, which was able to find good solutions and lower bounds.Singh and Weiskircher [2], [3] consider a similar problem to RCJS that is modelled as machine scheduling with fractional shared resources.They propose collaborative approaches via agent-based modelling, which prove to be effective in tackling these problems.Ernst and Singh [5] propose a Lagrangian relaxation based matheuristic that employs particle swarm optimisation, which showed improvements in solution quality across the RCJS benchmark dataset.Thiruvady et al. [6] consider a variant of RCJS with hard due dates, and they show that an ant colony optimisation (ACO) and constraint programming hybrid can efficiently find feasible and high quality solutions.The study by Thiruvady et al. [7] proposes a column generation approach for RCJS, which finds excellent solutions and lower bounds.An inherent issue is the time requirements in solving this problem, and hence, parallel implementations of metaheuristics and matheuristics have been attempted [8], [9].The best known solutions for RCJS have been obtained by Nguyen et al. [11] and Blum et al. [14].Nguyen et al. [11] develops a hybrid of differential evolution and column generation, and they show excellent results on small and medium sized problems.Blum et al. [14] devise a biased random key genetic algorithm, which produces excellent results, especially for large problem instances.In the interest of finding good solutions quickly, studies have shown that genetic programming can be effective for RCJS [10], [12], [15].Finally, Thiruvady et al. [13] propose Merge search and CMSA for RCJS.
RCJS under uncertainty (RCJSU) is a variant of RCJS that takes into consideration different practical random sources such as job release times, job processing times or the avail-ability of resources.Of vital importance in RCJS, and crucial in real-world settings, is the uncertain resource availabilities.That is, if trains and trucks breakdown or tracks and roads need maintenance leading to delays in the ore reaching ports, there can be substantial losses incurred.Therefore, planners need a solution method that can address disruptions due to random events to ensure that the solutions obtained can ensure smooth operations either by incorporating uncertainty into optimisation problems, reoptimisation, or reactive dispatching rules.Thiruvady et al. [16] developed the first RCJSU formulation, and proposed a surrogate-assisted ant colony system (SACS) to find robust solutions.The results show that SACS outperforms solution methods based on mixed integer programming (MIP), meta-heuristics and surrogate-assisted evolutionary algorithms.Further analyses show that SACS finds good solutions more quickly compared to standard ACO algorithms.
Although SACS has shown promising results, there are some key limitations.Firstly, surrogate models help SACS address the challenge of expensive solution evaluations, but the proposed algorithm tends to converge prematurely, especially for large instances.Therefore, SACS (as well as other existing algorithms in the study of Thiruvady et al. [16]) can waste a lot of time exploring the search space without making significant progress.Secondly, the proposed method has many pre-defined hyper-parameters such as the learning rate of ant colony systems (ACS) and pre-selection rate for the surrogate models.Because of this issue, the algorithm can have trouble in solving a wide range of instances without major modifications (e.g.reoptimise hyper-parameters).Finally, the performance of SACS is not consistent across different uncertainty levels and different problem sizes (i.e., number of machines).For medium problem sizes and large uncertainty levels, SACS does not show significantly superior performance compared to standard ACS.This paper addresses the three key limitations discussed above by developing a new adaptive population-based simulated annealing (APSA) algorithm for RCJSU.The main contributions of this paper are summarised as follows: (1) a new population-based simulated annealing algorithm with a balance between exploration and exploitation, (2) a new efficient Metropolis-Hastings algorithm extended from the simulated annealing algorithm proposed in [4], and (3) a new adaptive mechanism to update the probabilities for selecting perturbation operators to avoid premature convergence.APSA is similar to memetic algorithms, and population-based stochastic local search [17], [18] that combines evolutionary algorithms and local search heuristics; however, APSA adopts a much smaller population in which the population is improved mainly via local search and adaptive perturbations rather than crossover operators.The balance between exploration and exploitation within APSA allows the algorithms to perform efficiently against noisy and complex problems like RCJSU.Extensive experiments are conducted to validate the superior performance of the proposed algorithms for RCJSU and RCJS.
The paper is organised as follows.Section II provides a review of related work.In Section III, a formulation of the RCJSU problem is provided, including an integer programming model.Section IV discusses the adaptive population-based simulated annealing algorithm proposed in this study.The experimental evaluation and associated results are presented in Section V-A.Section VI concludes the paper and presents possibilities for future work.

II. RELATED WORK
A review of the literature related to RCJS is provided here, including methods and techniques that are effective at tackling problems similar to RCJS and RCJSU.

A. Optimisation techniques for tackling RCJS
Owing to the complexity of RCJS, mathematical programming methods on their own do not suffice in finding optimal or even feasible solutions to relatively small problems.Hence, metahuristics and matheuristics have been popular approaches for tackling the problem.Metaheuristics, including partical swarm optimisation [5], ACO [8], differential evolution [11], and genetic algorithms [14] have been known to find good solutions, albeit with large time overheads.RCJS can be modelled with graphs, and hence, ACO lends itself well to the problem.Nonetheless, past studies on RCJS have shown that ACO tends to converge slowly and gets trapped in local optima.Moreover, diversification strategies including local search (e.g., parallel implementations or niching techniques) help ACO, but the challenges still remain.Hybrid methods, especially those that can guarantee optimality, have also been explored [4], [7], [11], but the results show that only relatively small RCJS problem instances (up to 30 jobs) can be solved with these methods.

B. Scheduling on problems with uncertainty
Uncertainty is inherent in many real-world applications, and dealing with uncertainty is often challenging when devising scheduling algorithms.If uncertainty is not considered, the ensuing scheduling algorithms in the deterministic setting are often not effective at even finding feasible solutions (e.g. if a machine breaks down, a job will not be able to start at its pre-determined start time).Hence, a good understanding of the setting in which the schedules are applied is important to ensure success of the scheduling algorithms.Often, reactive scheduling strategies can be more effective in dealing with frequently occurring situations with uncertainties, since the latest information can be incorporated into the decisions.Nonetheless, if predicting the schedules is needed and knowledge of uncertainty is available, the preference is to develop robust schedules.The focus of this paper is on finding robust schedules for RCJSU, since predictable schedules are necessary for mining operations to coordinate supporting activities and to ensure that ore arrives at ports in a timely manner.Of particular interest is the timely arrival of ore at the ports, which relies greatly on the availability of resources.This is indeed the focus in the RCJSU problem and presented in the formulation in Equations ( 1) to (7) (see Section III).
When tackling scheduling problems with uncertainty, the aim is typically to devise methods that are capable of identifying good solutions under any uncertain scenario.These solutions are called robust solutions and the research area, robust proactive scheduling under uncertainty.There are numerous studies focusing on this area, which we summarise here.Beck and Wilson [19] investigate job shop scheduling, focusing on a variant which consists of stochastic processing times.Their results demonstrate that hybrids of tabu search and constraint programming with Monte Carlo simulation are effective in handling large-sized problems with uncertainty.The study by Song et al. [20], proposes branch-and-bound algorithms that generate robust schedules on a variant of resource-constrained project scheduling with processing times that are uncertain depending on time.Identifying robust solutions in uncertain problems introduces additional complexity, and for this purpose metaheuristcs such as swarm intelligence algorithms and evolutionary approaches are known to be effective.Bierwirth and Mattfield [21] study genetic algorithms for dynamic job arrivals in classic job shop scheduling.Their results show that their approaches are superior to priority dispatching rules, which are popular in dynamic production scheduling.Wang et al. [22] devise a surrogate-assisted multi-objective evolutionary method to tackle proactive scheduling where machines are considered to have stochastic breakdowns.To cope with the uncertainty, they propose a job compression and rescheduling strategy [23].Support vector regression in conjunction with surrogate models have been proposed to enhance the speed of time-intensive simulations.The study by Wang et al. [24] provides an in-depth analysis and discussion of scheduling with machine breakdowns.Systematic reviews of scheduling in problems with uncertainty have also been conducted [25], [26].
The least two decades have seen numerous studies that consider uncertainty in scheduling problems similar to that of RCJS.There are several similarities to project scheduling, and many studies have considered uncertainty within this setting [27], [28], [29], [30], [31], [32], [33].In these studies, uncertainty is typically considered in respect to the resources, though, a few consider processing time related uncertainty.Lambrechts et al. [29] study uncertainty with resources, and specifically modes of execution related to the activities.They develop the notion of "time buffering," that efficiently generates robust schedules.Masmoudi and Haït [30] consider resource levelling and resource constrained project scheduling, and they propose fuzzy modelling to tackle uncertainty.Their solution approaches are a genetic algorithm and a greedy method to identify robust schedules.Khodakarami et al. [31] investigate Bayesian networks to identify causal structures considering uncertainty in sources and parameters of a project.They find that Bayesian networks are indeed very effective at dealing with uncertainty, including aspects such as representation, elicitation and management.Uncertainty associated with processing times has also been investigated [32], [33].Moradi and Shadrokh [32] develop a problem-specific heuristic to deal with uncertainty, and the study by Chakrabortty et al. [33] propose several Branch & Cut heuristics for project scheduling. 2 2 Benchmark datasets obtained from the project scheduling problem library (PSPLIB) [34].

C. Hybrid optimisation methods
Local search and population-based methods [35] have been extensively investigated in the literature to tackle challenging scheduling problems.While local search heuristics are best known for their ability to find local optima efficiently, population-based optimisation methods such as genetic algorithms are better known for their exploration ability.Empirical studies have demonstrated the need to maintain a balance between exploration and exploitation when dealing with difficult combinatorial optimisation problems such as job shop scheduling [36].Such hybrid methods are commonly referred to as memetic algorithms [17] or population-based stochastic local search algorithms [18] in the literature.While hybrid algorithms are effective in combinatorial optimisation, there are still a few shortcomings.Firstly, design choices to develop competitive hybrid algorithms are difficult to make and can be sensitive to different problem subsets.Secondly, the performance of hybrid algorithms depend strongly on an "ideal" setting of the algorithm parameters.
Hybrid algorithms have also been developed for RCJS.The study that introduced the problem [4] develop a Lagrangian relaxation heuristic, which combines MIP and simulated annealing.Following this, Thiruvady et al. [7] show that a column generation and ACO hybrid is very effective on this problem by finding good lower bounds.Parallel implementations of ACO [8] have also been effective, by finding highquality solutions quickly.[11] propose a hybrid algorithm combining differential evolution and local search with excellent results.The study by Thiruvady et al. [13] show that parallel matheuristics based on MIP and ACO can effectively tackle this problem.
Thiruvady et al. [16] has developed SACS, a sophisticated algorithm combining ACO, surrogate modelling, and local search to deal with RCJSU.The experiments show that SACS can find excellent solutions for small problems thereby outperforming other meta-heuristics.However, further analyses show that SACS tends to converge prematurely and is inconsistent when working with different problem subsets.Moreover, SACS has many hyper-parameters (for ACO, local search, surrogate models), for which the ideal settings are difficult to find.

III. PROBLEM FORMULATION
The formal definition of RCJSU as an integer program (IP) was provided by Thiruvady et al. [16].Past studies have considered uncertainty in different ways, including resource limits [29] or uncertainty in job/task processing times [33].However, in this study the focus is on uncertainty derived from the variability in resources, which in turn affects the processing times of jobs/tasks. 3.The problem is NP-hard as it is at least as hard as its deterministic counterpart, which is indeed shown to be NP-hard [4].This IP formulation extends the one of Singh and Ernst [4], proposed for RCJS.Uncertainty is introduced by considering several uncertain resource levels, that is, some proportion of the original resource limit is used for each uncertainty sample (details in Section V-A).
We are given several jobs J = {1, . . ., n}, which are associated with the following information.Each job i must execute on machine m i , (M = {1, . . ., l}), has a a release time r i , processing time p i , due time d i , weight w i and resource usage g i .Jobs on the same machine cannot execute concurrently, or in other words, only one job may execute at one time.A time horizon is given, consisting of discrete points T = {0, . . ., D}, and D is chosen to be sufficiently large to ensure a feasible solution will always be found.Precendences may exist between any two jobs i ∈ J and j ∈ J that execute on the same machine.This is represented as i → j, where job i completes before job j starts.
A central resource constraint exists across all machines, where the cumulative resource usage of jobs executing at any time point must be at most G.When considering uncertainty, G is allowed to vary.This leads to defining a sample, which has available a proportion of the resource G, and is an instance of the problem.As an input, the total number of samples u is provided, and the resource limit is G s applies to the s th sample U = {1, . . ., u S }.
The IP can be defined as follows.Let z sjt be a binary variable, which indicates in sample s, job j completes or has completed by time t, if z sjt = 1.By this definition, a job stays complete once completed.Compared to RCJS, the RCJSU IP formulation requires additional space complexity, which consists of O(n 3 ) variables, and each set of constraints with an additional dimension owing to the samples.The IP formulation is: Equation ( 1), the objective function, states that the average total weighted tardiness (TWT) across all samples must be minimised.Constraints (2) ensure that all jobs complete.Once a job completes it stays complete, which is enforced by Constraints (3).Constraints (4) ensure that all jobs satisfy release times and the precedences are imposed by Constraints (5).In this constraint, e j = r j + p j − 1. Constraints ( 6) require that only one job can execute on one machine at every time point.Constraints (7) ensure that the resource limits, considering each sample, are satisfied.For the RCJS, there are efficient scheduling heuristics that are known (see Section IV-D).These heuristics typically depend on an alternative formulation, where a sequence or permutation jobs are efficiently mapped to schedules.Previous studies have used this approach for RCJS (e.g.[7]), and it has also proved beneficial for the RCJSU [16].In this study, π represents a sequence of jobs while σ(π) represents an efficient mapping from π to a feasible schedule.Therefore, π is effectively a solution to the problem, where given a sequence, a deterministic schedule can be efficiently constructed.
An example RCJS problem instance is shown in Table I.This example shows how a solution is constructed and what the optimal solution is.The problem consists of three jobs, with the associated information in the tables.Moreover, the resource capacity is 10 units.Considering π, there are six potential sequences.Among the jobs, the resulting schedules from π will be efficient if Job 3 is scheduled early owing to an early due date and high weight.Additionally, Job 2 should be sequenced before Job 1 due to its higher weight.Given this information, the optimal solution is π = {3, 2, 1}.

IV. POPULATION-BASED SIMULATED ANNEALING
Several metaheuristics have been able to find excellent solutions to RCJS in a time efficient manner.One of these approaches was simulated annealing [4].However, to tackle RCJSU, the traditional simulated annealing method has a few shortcomings, and hence, we propose an adaptive populationbased simulated annealing (APSA).In particular, this method has three main differences from the traditional approach: (1) a population of solutions to track and exploit best known solutions, (2) a faster cooling schedule within the Metropolis-Hastings algorithm, and (3) adaptive probabilistic selection of the neighbourhood moves to perturb solutions in the population.Comparisons are made to ACS (ant colony system) and SACS (surrogate-assisted ant colony system) of [16], which have previously been effective on RCJSU.

A. The APSA Algorithm
Algorithm 1 shows the high-level APSA procedure for RCJSU.As input, the algorithm takes the following: (a) T 0 : an initial temperature, (b) t lim : time limit for the execution of the algorithm, (c) Iter: number of Metropolis-Hastings iterations, (d) s: population size, (e) γ: a parameter to control the cooling schedule and (f) ρ: a parameter that specifies how to adapt the probabilities (ρ = 0.9 for this study).In Line 2, the population Π is initialised, where each solution consists of a random list of Normalise(p b , p j , p r ) 20: end for 21: end while 22: π f := SwapAllJobPairs(π bs ) 23: return π f all the tasks.The best known solution, π bs , is initialised to be empty.In Line 3, the probability distributions associated with three neighbourhood moves (discussed in the following) are initialised.The main procedure takes place between Lines 4-21, and the algorithm completes execution when the time limit is reached.
Within the main loop, the following is applied to each solution in the population.First, in Line 6, the Metropolis-Hastings algorithm is applied (Algorithm 2) using π as the best solution and for Iter iterations.The resulting solution is stored in π t .π bs is updated to π, if π is an improvement, i.e., f (π) < f (π bs ).Following this, between Lines 9-19, π is perturbed using a neighbourhood move (details follow in Section IV-C).As part of these steps, the probabilities associated with choosing one of the perturbation moves are adapted or updated.That is, if a move is selected, its corresponding probability reduced by an adaptation factor ρ = 0.9.
Every solution obtained (by applying any neighbourhood move, etc.) requires the computation of its objective value.As mentioned earlier, a solution is represented by a sequence of the tasks.For each sequence of tasks π, a feasible schedule (σ(π)) is obtained by constructing u s schedules corresponding to each uncertainty sample.Note, each schedule is resource feasible, ensuring that at most G s resources are used at each time point for sample s.The result is that each sample has its own TWT.The objective value of π is computed as the average of the TWTs across all samples.

B. The Metropolis-Hastings Algorithm
The Metropolis-Hastings method relies on the cooling schedule to ensure a gradual convergence to high-quality Algorithm 2 MetropolisHastings end if 14: end if 16: end for 17: return π ib regions of the search space.However, due to the overheads introduced by uncertainty, the traditional gradual reduction in temperature will not suffice, and hence we propose a modified version of the method for RCJSU.
Algorithm 2 shows the Metropolis-Hastings procedure.It takes as input: (a) π: current best solution, (b) T 0 : initial temperature, (c) Iter: number of iterations that the algorithm will execute for and (d) γ: a parameter to determine the cooling schedule.C represents the best known objective value, and is set to a large number.The temperature T c is set to the initial temperature T 0 , and the best solution π ib is initialised to π.The algorithm executes for Iter iterations between Lines 3-16.Using π ib , the neighbourhood move of swapping a pair of jobs is applied.That is, two jobs within the sequence are selected at random and swapped (see Figure 2).If the new solution found π b is an improvement, C is updated to the objective value of π b .Moreover, π ib is updated to π b and the iteration count is reset.Otherwise, a move to the non-improving solution π b is still chosen with probability e −(f (π ib )−C)/T c .The cooling schedule takes place in Line 14.

C. Adaptive Perturbation
Key to the success of APSA is perturbing the best known solution, via three neighbourhood moves.In particular, when applied carefully at different stages of the search, the overall APSA algorithm is able to achieve excellent trade-offs between intensification and diversification.
The first neighbourhood move is β-sampling, which was first studied in the context of project scheduling [37].The basic idea underlying the method is that consecutive jobs or tasks with a sequence are inter-dependent, and moving whole sub-sequences (maintaining certain dependencies) can lead to significant improvements in solutions.For the purposes of the current study, starting at a randomly selected position in the sequence, a sub-sequence of fixed length (5 jobs) is moved to the back of the permutation.
Figure 1 provides an illustration of this method.π consists of 11 jobs.A position p ∈ 1, . . ., n − 5 is chosen randomly, Fig. 1: An example that demonstrates β-sampling.From the sequence π, a subset of jobs is selected (7,9,10,3,11) and moved to the end of the sequence.The second neighbourhood move is Job swapping.This consists of swapping pairs of jobs in π. Figure 2 shows an example of swapping one pair of jobs.Two positions or indices are chosen at random in the sequence π, and the jobs in these positions are swapped.In the example, Jobs 2 and 8 are swapped.
The third neighbourhood move is to re-initialise the current solution to a randomly chosen sequence.Not surprisingly, the resulting solution can belong to vastly different areas of the search space with a very different objective function value.Hence, this neighbourhood move is expected to produce the largest amount of diversification.

1) Convergence of Probabilities:
Assuming the probability values of p b , p j and p r at time step t are p t b , p t j and p t r , the expected probabilities at time step t + 1 can be calculated as Hence the governing differential equations for evolving the probabilities p b , p j and p r can be derived: Theorem 1.For any initial condition of p 0 b , p 0 j , p 0 r ∈ [0, 1] and p 0 b + p 0 j + p 0 r = 1, the probability values converge to p b = p j = p r = 1/3.
Proof.At equilibrium of the differential equations ( 11)- (13), Solving this system of equations with p b + p j + p r = 1, we can obtain the following seven equilibria from the differential equations: However, the first six equilibria are unstable, in the sense that a small disturbance will perturb the system away from those equilibria.The last equilibrium (p b = 1/3, p j = 1/3, p r = 1/3) is the only stable equilibrium.Hence, after a sufficient number of time steps, the probability values will converge to this stable equilibrium.
Figure 3 shows a simulation of how the probabilities associated with each neighbourhood move (Algorithm 1) evolve across 200 simulations.We see that three probabilities converge to almost equal levels by 125 simulations.

D. Scheduling Heuristic
The aim of the scheduling heuristic is to generate a resource feasible schedule σ(π) from the sequence π.For this purpose, the serial scheduling scheme is employed [38], which is adapted from previous studies on RCJS [4], [6], [7], [8].Note, Thiruvady et al. [16] showed that this scheduling scheme can indeed be effective for RCJSU.Fig. 4: An example of a problem with three machines and 15 jobs.From π, jobs are selected in order (left to right) and scheduled on machines that they belong to (e.g.orange jobs on m 2 ).Those jobs that have predecessors that are not yet scheduled will be put on the waiting list π; (a) Job 8 requires Job 2 to be scheduled before it can be scheduled.(b) Job 8 is scheduled once Job 2 has been scheduled, and this will be done as early as possible considering resources.
The scheduling heuristic works as follows.Jobs in π are selected in the order in which they appear.Each job (j) is tested to see if any of its preceding jobs have not yet been scheduled, and if so, j is placed on the waiting list.Otherwise, j is scheduled as early as possible to ensure that the resource limits are not violated.Once j is scheduled, any job that depends on j (and only j) on the waiting list is scheduled immediately.All such jobs on the waiting list are scheduled in order.Note, for RCJSU, schedules using the procedure above are created for each sample.
Figure 4 shows how the scheduling heuristic works.The problem consists of 15 jobs to be scheduled on three machines.The height of a machine represents the amount of resources available, while a job's height represents its resource requirement.A colour implies a job's machine, i.e., jobs of one colour must go on one machine (e.g."orange" jobs should be scheduled on machine m 2 ).
Figure 4a shows a permutation π, and a solution that is partially completed with a few jobs allocated to each machine.The waiting list also consists of some waiting jobs π, where for example, Job 8 waits for Job 2 to be scheduled, considering Job 2 precedes Job 8. Figure 4b shows one step in the method, where the next job in π, Job 2, is chosen and scheduled in the first available resource feasible time slot on m 1 .Once Job 2 is scheduled, the waiting list can be examined, and it is found that Job 8 can be released.Job 8 is scheduled as early as possible considering the resource availability.

V. EXPERIMENTAL SETTING
We investigate APSA on RCJSU, and compare its performance to current state-of-the-art approaches, including ACS and SACS.Additionally, we investigate APSA against standard simulated annealing (SA) with the same settings, and the SA of Singh and Ernst [4].Finally, we examine APSA's performance on the deterministic problem (RCJS), and we compare APSA to ACS, SACS and a column generation and ACO hybrid by Thiruvady et al. [7].The experiments are conducted on problem instances from the literature. 4hile this data was originally created for the RCJS, each problem instance is modified to generate the data for RCJSU.Specifically, there are 10 samples generated for each problem instance, where each one has resource limits in the range [U min , U l × G].U min is chosen to be the largest resource required by any job, i.e. max i∈J g i , which ensures a resource feasible schedule can always be found.Additionally, the multiplier U l ∈ {0.5, 0.6, 0.7, 0.8, 0.9, 1.0} is used to select a proportion of the resources.Typically, a smaller value (e.g.0.5 or 0.6) means that the problem instance is tightly constrained, where the schedules will need larger completion times, leading to a higher overall TWTs.
The parameter settings for APSA were selected through tuning by hand and by using past studies as a guide.We refer to Algorithms 1 and 2. The time limit t l im is set to 600 seconds.The initial temperature T 0 was set to 1500 [4].Due to uncertainty and the overheads required to construct schedules with a number of samples, the parameters within the Metropolis-Hastings algorithm have been set to allow for much quicker cooling schedule.This is achieved by allowing only 1000 iterations Iter = 1000 and γ = 0.5.The population size was found through tuning by hand s = 10.The subsequence size used within β-sampling of length 5 jobs, which is obtained from [4].Moreover, in Line 14 of Algorithm 1, there are m = rand() * n applications of job swapping, where rand() is a random number between 0 and 1.The parameter settings for ACS and SACS were obtained from [16], to which we refer the reader for complete details.
We carry out experiments on all problem instances considering all uncertainty levels.There are 25 repetitions of for each algorithm, with each one given 10 minutes of wall clock time.The memory limit was 2 GB, which proved to be sufficient for all algorithms.The algorithms were implemented in C++ with GCC-5.4.0.All experiments were carried out on Monash University's Campus Cluster -MonARCH. 5he IP was implemented in Gurobi 9.0.1 [39].It was run for 10 minutes per problem instance and allowed a larger amount of memory of 10GB.Across all the runs, only one solution was found for one problem instance at the uncertainly of 0.8.This is not surprising given the size of the model, which is unsuitable in practical settings.

A. Results
The results that follow are reported as a percentage difference from the average solution quality over 25 runs for each algorithm to the best solution found by any algorithm.Consider a problem instance and let B * be the best solution obtained by any algorithm for that problem instance.Let A k be the best solution found by APSA for instance k.Then the percentage difference for APSA is A k −B * A k .Tables II and III show the results of ACS, SACS and APSA for low and high uncertainty levels U L, respectively.The tables show the best result obtained by any algorithm (which is highlighted in bold if it is found by APSA) for each uncertainty level, followed by the percentage differences of each algorithm to the corresponding best solution.The best result found by any algorithm is highlighted in bold, and those results that are statistically significant (determined using a Wilcoxon ranked-sum test) are marked with a "*".To summarise the performance of the algorithms across the problem instances, the last row of the tables report the number of times (# best) the algorithm finds the best solution.
The results in both tables show that APSA finds best known solutions for all problem instances and uncertainty levels (except instance 15-2 at uncertainty level 0.9), and the best average results across all problem instances for all uncertainty levels.Looking closely, we see that the differences between the algorithms are smallest for the smallest problems (3 machines) and largest problems (20 machines).Everywhere else, APSA has significant advantages.For the problem instances with 3 machines, this is not surprising, because it is expected that all algorithms find solutions that are close to optimal.However, APSA has the advantage that, for every run, it finds a solution close to the best.On investigating the individual runs on the large problem instances, we find that there are relatively few iterations conducted (< 10).This leads to APSA not converging sufficiently well, thereby not finding as large improvements are seen on other smaller problem instances.
The results in Tables II and III are summarised in Figure 6.The figure shows the performance of each algorithm (ACS, SACS and APSA) by uncertainty levels.The bars represent the average of percentage difference to best across all instances in an uncertainty level.As expected, APSA clearly outperforms the other two methods, with a slight advantage with increasing uncertainty levels.

B. Convergence Characteristics of APSA
APSA outperforms SACS and ACS across all problem instances for all uncertainty levels, due to its ability to diversify well.In particular, APSA's ability to stop getting stuck in local optima leads to substantial gains.
In the following, we examine the convergence of APSA, and compare it to ACS and SACS.We choose four problem instances, which are sufficiently diverse in size to allow making somewhat general inferences about APSA.The results are shown in Figure 6, where each sub-figure corresponds to one problem instance.The plots show time in seconds on the x-axis and TWT (logarithmic scale) on the y-axis.Again, APSA is clearly the best performing method across all four problem instances.There are two main advantages of the algorithm, especially related to the population and diversification.Firstly, the starting solution, irrespective of the problem instance, is better than ACS or SACS.This can be attributed to the population, where at the beginning of the search process, there is a substantial amount of diversity in the starting solutions.Secondly, and more importantly, there is continuous improvement in APSA over-time.This can be clearly seen for medium to large problem instances (6-28, 9-47 and 12-14).The increased diversity allowed by the reinitialisation of the solutions in the population ensures that the population does not get stuck in local optima.

C. Comparisons of APSA to the original SA ([4])
The original SA (OSA) algorithm [4] can be adapted to tackle the RCJSU problem. 6The main difference is that every time a new solution is created, a new objective needs to be computed, which involves creating multiple schedules corresponding to the samples.Otherwise, all the settings for the algorithm remain the same.
Table IV shows the results for APSA and OSA considering the uncertainty levels 0.6 and 0.9.Other than the smallest problem instances (3 machines), APSA clearly is superior, finding better average solutions (> 85%) across both uncertainty levels.

D. Comparison of APSA with SA
The novel components of APSA including the population of solutions, adaptive probabilistic selection and quicker cooling schedule clearly lead to significant gains over traditional OSA.In this section, we aim to understand the impacts of the population on APSA.This essentially involves an implementation of APSA with a population size of 1.
The results can be seen in Table V, where comparisons are made for the uncertainty levels 0.6 and 0.9.We see that there is a clear advantage of the population, and the advantage increases with the higher uncertainty levels.

E. Comparisons on the Deterministic Problem
APSA's advantages obtained from the population, a quicker cooling schedule and adaptive probabilistic selection for the perturbation of solutions, make for an effective algorithm in the uncertain setting.Its performance on the deterministic problem or RCJS is worth investigating, to see if any of the novelties proposed in the algorithm can be beneficial for the deterministic problem.The results are shown in table VI, where APSA to ACS, SACS, and a column generation and ACO hybrid (CGACO) are compared.
CGACO is also run 25 times on each problem instance.A single sample was used for ACS, SACS and APSA, leading to the original resource limit (effectively RCJS).We see in Table VI that beyond small problem instances (> 4 machines), APSA has a clear advantage.It generally performs best, but without the same large gains seen on RCJSU.

VI. CONCLUSION
This study investigates an optimisation problem that originates in the mining industry, specifically transporting ore from mines to ports.Past studies have considered the problem as deterministic, however, in real settings, the resources needed to transport the ore are subject to uncertainty.The problem is also called resource constrained job scheduling with uncertainty, where for this study, we consider multiple uncertain scenarios.To tackle this problem, we propose an adaptive populationbased simulated approach, which has significant advantages over previously proposed ACO based approaches.We find that using a population, increasing the rate of the cooling schedule within the Metropolis-Hastings algorithm and adaptively selecting intensification and diversification via the neighbourhood moves, proves advantageous.These modifications, compared to traditional simulated annealing, lead to finding robust solutions, and hence, best results on a benchmark dataset for RCJSU.
Despite the success of APSA on RCJSU, investigating how the approach performs on other optimisation problems with uncertainty will be of interest in future work.In particular, it is expected an adaptation of APSA is likely to be effective for problems that consist of highly constrained disjoint feasible regions, such as RCJSU.
The MIP model in Section III is clearly too complex and large for existing solvers, but decomposition approaches such as column generation or Benders' decomposition may prove effective.The advantage of such methods is that they can provide guarantees of solution quality (via lower bounds to RCJSU), which can help to measure the quality of the solutions found by APSA.

Fig. 2 :
Fig. 2: An example that demonstrates job swapping.Two indices are selected randomly, and the corresponding jobs are swapped.In this example, indices 3 and 9 are selected, and the jobs within those indices, 2 and 9, are swapped.and this sub-sequence (p to p + 5) moves to the end (bottom part of the figure).Here, p = 4, which starts at Job 7, and Jobs 7, 9, 10, 3 and 11 are chosen to move to the end of the sequence.Subsequent jobs (Jobs 8,4, and 5) move up in the sequence to position 4. β-sampling is applied in the above manner in Line 10 of Algorithm 1.

TABLE I :
An example with 3 jobs and 10 units of resource.Proc. is the processing time.

TABLE II :
Results of ACS, SACS and APSA shown as the percentage difference to the best solution for uncertainty levels 0.5, 0.6 and 0.7; Results in column "Best" are highlighted in bold if APSA finds the best solution; Results for each algorithm are highlighted in bold when they achieve best average results, while statistically significant results (pair-wise Wilcoxon ranked-sum test) at a 95% confidence interval are marked with a *.

TABLE III :
Results of ACS, SACS and APSA shown as the percentage difference to the best solution for uncertainty levels 0.8, 0.9 and 1.0; Results in column "Best" are highlighted in bold if APSA finds the best solution; Results for each algorithm are highlighted in bold when they achieve best average results, while statistically significant results (pair-wise Wilcoxon ranked-sum test) at a 95% confidence interval are marked with a *.

TABLE IV :
[4]omparison of PSA and the original SA[4](SAO) for the uncertainty levels of 0.6 and 0.9.

TABLE V :
A comparison of PSA and SA for the uncertainty levels of 0.6 and 0.9.