Capacity planning in supply chains of mineral resources

Accepted Manuscript Capacity Planning in Supply Chains of Mineral Resources Joey Fung, Gaurav Singh, Yakov Zinder PII: DOI: Reference:

S0020-0255(14)01085-8 http://dx.doi.org/10.1016/j.ins.2014.11.015 INS 11251

To appear in:

Information Sciences

Received Date: Revised Date: Accepted Date:

29 November 2013 18 September 2014 6 November 2014

Please cite this article as: J. Fung, G. Singh, Y. Zinder, Capacity Planning in Supply Chains of Mineral Resources, Information Sciences (2014), doi: http://dx.doi.org/10.1016/j.ins.2014.11.015

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Capacity Planning in Supply Chains of Mineral Resources Joey Funga , Gaurav Singhb , Yakov Zindera a

School of Mathematical Sciences, University of Technology, Sydney, PO Box 123, Broadway NSW 2007, Australia b CSIRO Computational Informatics, Private Bag 33, South Clayton VIC 3169, Australia

Abstract The paper addresses the existing gap in the literature on the optimisation in capacity planning for mineral supply chains. The presented optimisation procedure aims at minimising the cost of infrastructure expansion for any given scenario of future demand. The optimisation procedure is designed as a matheuristic – a hybridisation of a mixed integer linear program (MILP), and a simulated annealing based scheduler. The optimisation procedure is iterative in nature and has the following distinctive features. Each iteration starts with generating a MILP model and finding a minimal cost infrastructure expansion for this model. Then, the scheduler analyses the MILP solution by constructing a schedule. In constructing this schedule, the scheduler reduces its search space using the MILP solution. The scheduler identifies bottlenecks in the infrastructure, which are used for generating a new MILP model at the next iteration. The MILP and the scheduler use different levels of data aggregation and their interaction mechanism is designed as a solution process of a bi-criteria optimisation problem. The computational experiments on data, originating from the world’s largest coal exporter, shows the ability of the developed matheuristic to solve industrial-scaled instances of the problem. Keywords: Capacity Planning, Matheuristic, Hybrid Heuristic, Mixed Integer Linear Program, Simulated Annealing

Email addresses: [email protected] (Joey Fung), [email protected] (Gaurav Singh), [email protected] (Yakov Zinder) Preprint submitted to Information Sciences

November 21, 2014

1. Introduction Coal, iron and other mineral resources have a vital role in the world economy. For example, in 2012 coal provided 29.9% of the global primary energy needs and generated 41% of the world’s electricity1 . One of the difficult managerial problems in supply chains of mineral resources is capacity planning, which aims at finding a cost effective expansion of the infrastructure in order to satisfy a forecasted increase in demand. The difficulty of the problem is due to the complexity of this type of supply chains which involve mines that produce the commodity; trains which transport the commodity to the sea terminals through a railway network; train-loading facilities at the mines; building stockpiles on the limited strips of land at the sea terminals, with each stockpile having a pre-determined composition of shipments from different mines; as well as groups of equipment at the sea terminals used to unload trains, to build stockpiles and to load the ships. Capacity planning is crucial for supply chains of mineral resources due to the multi-million dollar cost associated with capacity expansion, the long lead-time for building the additional infrastructure and the difficulty of withdrawing from building projects and purchasing commitments. In the considered class of supply chains, the demand is generated by prespecified arrival and departure times of ships and their cargo requirements. The cargo requirement of a ship is the set of stockpiles to be loaded onto the ship. One of the approaches to capacity planning is to base the decision making process on a set of scenarios of ships’ arrivals, departures, and their cargoes. According to the scenario-based approach, for each scenario, it is necessary to minimise the cost of expansion of the supply chain infrastructure in order to meet the scenario requirements. In this paper, we present an optimisation procedure which, for a given scenario, optimises the cost of infrastructure expansion. In practice, the capacity planning is an extremely complex process that involves different methods and analysis from different viewpoints. The cost minimisation problem considered in this paper is the core of this planning process. The contribution of this paper is two-fold. Firstly, the paper fills the gap in the literature on optimisation in capacity planning for supply chains of mineral resources, since very little has been published on this very challenging problem. Secondly, the paper contributes a new optimisation procedure with 1

http : //www.worldcoal.org/bin/pdf /original pdf f ile/coal f acts 2013(11 09 2013).pdf

2

the design motivated by the necessity to consider the entire supply chain and a long planning horizon, and at the same time to account for the impact of queues and congestions. The optimisation procedure, presented in this paper, is based on the hybridisation of two optimisation engines: mixed integer linear programming (MILP) and a metaheuristic based scheduler. The hybridisation of a mathematical program and a metaheuristic is commonly referred to as a matheuristic [22, 5]. As far as the field of hybrid optimisation algorithms is concerned, the paper contributes an optimisation procedure with the following distinct features: • the engines use different levels of data aggregation and their interaction is designed as an iterative optimisation of the constraints in the MILP; • the MILP, which determines the infrastructure expansion at each iteration, is specifically designed for this type of supply chains and, to the authors’ best knowledge, has not been previously reported in the literature; • the simulated annealing based scheduler uses the flows of commodity specified by the MILP in order to reduce the search space in the process of analysing the infrastructure expansion; • the scheduler identifies bottlenecks in the supply chain, and therefore what constraints in the MILP to be modified at the next iteration in order to make the MILP model more accurate. The above features of our optimisation procedure can be viewed as a framework which leads to different algorithms depending on the choice of MILP model and/or scheduler. Various implementations of this framework and their comparison is one of the directions of future research. The presented implementation of this framework has been justified by means of computational experiments with real world data provided by the Hunter Valley Coal Chain2 , which is the world’s largest exporter of coal. To the authors’ knowledge, this is the only attempt of modelling with such level of accuracy whilst enabling the solution of industrial-scaled instances of the problem. Also, even though the optimisation procedure was tested using data from 2

http : //www.hvccc.com.au

3

coal supply chains, our procedure is equally applicable to supply chains of other mineral resources, like iron ore. The remainder of the paper is organised as follows. Section 2 provides a motivation and a comparison with the existing literature. Section 3 provides more details on the supply chains of the considered type. Section 4 gives an overview of the optimisation procedure. Section 5 presents the mixed integer linear program and Section 6 describes the scheduler. Section 7 contains the results of computational experiments. The conclusions can be found in Section 8. 2. Motivation and Comparison with the Existing Literature Two factors have motivated this research. The first is the lack of publications on the very challenging optimisation aspect of the capacity planning in supply chains of mineral resources, which partly can be attributed to the complexity of these supply chains. For example, in regards to the Hunter Valley Coal Chain, [4] acknowledges that “To the best of our knowledge, there are only a few papers in the literature on the use of optimization models in bulk goods supply chains of this complexity”. The complexity of the cost minimisation problem, considered in our paper, stems from the following requirements. The capacity planning as a long-term planning problem needs data aggregation (see for example [30, 7]). However, this aggregation should not diminish the significance of the multi-stage nature of supply chains of mineral resources. It is well known (see for example [32]) that neglecting the scheduling can distort capacity estimates so far as to make the results of capacity planning invalid. This leads to a difficult trade-off – aggregation versus the necessity to take care of queues and congestions. The gap in the literature which has motivated our research is not compensated by the vast body of knowledge on the related topics such as capacity planning in manufacturing or supply chain network design. The majority of these publications are not relevant to the optimisation of capacity planning for supply chains of mineral resources since they do not address to the required extent the necessity of considering the long-term strategic decisions in conjunction with the impact of possible queues and congestions. Thus, an obvious related area is manufacturing with its multi-stage production processes. However, these publications either do not take into account enough details of transportation which, in supply chains of mineral resources, involves trains and a railway network, or are concerned with cyclical production (see for 4

example the survey [10]). Also, the common practice in capacity planning in manufacturing settings is to replace a sequence of operations by a set of operations or even to aggregate sequences of operations into a single operation. This approach ignores the order of operations and hence congestion (see for example [12, 6, 33]). Another closely related area is capacity investment and supply chain network design. However, the rich literature in this area focuses on highly aggregated models, which ignores the impact of queues and congestions, hence rendering these publications inapplicable for the problem considered in our paper. For example, in [28, 13], the models do not go below the plant level. If transportation is considered at all, then it is considered either as traveling salesman type models (see for example [14]) or as very aggregated flows of commodity. For instance, [11] aggregates transportation between pairs of locations as flows and [26] aggregates all operations such as loading and unloading and all details of transportation under the umbrella of the term “trip”. The second factor motivating our paper was the necessity to design a new approach to the considered cost minimisation problem. To the authors’ knowledge, the only publication which is devoted to optimisation of capacity planning for supply chains of mineral resources is [31]. This paper attempts to determine the desired levels of resources whilst accounting for the transportation and multistage material handling. The method suggested in [31] is solely MILP, which resulted in a very large mixed integer linear program and caused difficulties in obtaining even a feasible solution. In comparison with [31], the MILP model in our paper reflects the real-world situation more accurately and permits to avoid some restrictions which may not occur in the real-world. For instance, the MILP model in [31] discretises the planning horizon into days and assumes that each delivery starts and completes on the same day. The deficiency of linear programming in modelling systems with queues and congestions is well known, and unfortunately also affects the accuracy of the modelling in [31]. The above mentioned difficulties have motivated our research on an alternative approach to this important cost minimisation problem. As far as the second motivating factor is concerned – the necessity to develop a suitable approach to the very challenging cost minimisation problem, one of the closest fields of research is the field which is commonly referred to in the literature as simulation optimisation (see for example [15, 24, 9, 8]). The publications on simulation optimisation are normally cast into one of two 5

categories. The first category is comprised of simulation algorithms where optimisation is used as a subroutine, predominantly for determining the simulation rules in the course of simulation (see for example [29]). Given that decisions on capacity planning affect the entire planning horizon and that the nature of discrete event simulation is a step-by-step application of certain dispatching rules, the first category of publications on simulation optimisation is not suitable for the purpose of our paper. The second category uses simulation as a subroutine for calculating the objective value in the process of optimisation search (see for example [25]). In the publications comprising this category, the simulation does not have a direct impact on the optimisation search besides providing an estimate of the objective value. The only publication known to the authors which acknowledges the possibility of changing the optimisation model based on the results of simulation is [1], although this publication leaves any discussion on the interaction mechanism between optimisation and simulation outside the scope of the paper. Instead, [1] outlines a simple example where simulation is used for replacing the initial estimates of waiting times by the estimates obtained through simulation. Furthermore, in publications comprising the second category, no optimisation occurs inside the simulation process. These make the the second category not suitable for the purpose of this paper either. Saying the above, it should be stressed that simulation models have been used to assess the merits of different capacity configurations in supply chains of mineral resources under the assumption that the capacity configuration is given [36]. In other words, such publications do not consider the optimisation aspect of the problem, which is the main focus of our paper. According to [36], even the assessment of the given capacity configuration requires significant efforts itself. The necessity to make decisions on capacity expansion by analysing the entire system and the impact of these decisions over the entire planning horizon suggests the use of MILP. On the other hand, as has been mentioned above, the capacity planning in supply chains of mineral resources should take into account the impact of possible queues and congestions, which suggests some scheduling optimisation procedure. Given the scale and complexity of the problem on hand, it is reasonable to resort to metaheuristic based scheduling, which is commonly used in practical applications [27]. These observations have motivated the design of our optimisation procedure as a hybridisation of MILP and simulated annealing. The development of hybrid optimisation algorithms is an active area of 6

research [3, 34]. If a mathematical program is hybridised with a metaheuristic, the resultant algorithm is commonly referred to as a matheuristic [22, 5]. According to the classification in [34], the matheuristic presented in our paper is a so-called high–level relay hybridisation. In contrast to the low–level hybridisation, where one component is a sub-routine of another (see for example [2, 20, 35]), the high–level relay hybridisation is where both components are self-contained optimisation engines which are executed sequentially – in our matheuristic even iteratively. The hybridisation implemented in our paper differs from the commonly used approach in the high–level hybridisation where an optimisation engine can affect the other components only by providing a starting point in the optimisation search (see for example [14, 12, 21, 23]). In contrast, in our optimisation procedure, the solution of the MILP reduces the search space of the simulated annealing based scheduler, whereas the bottlenecks in the supply chain infrastructure, identified by the scheduler, result in the change of constraints in the MILP. Furthermore, although it is not unusual in the hybridised algorithms found in the literature that the optimisation engines may focus on different parts of the overall problem, for example by freezing the values of some decision variables, they are effectively dealing with the same problem (see for example [14]). On the other hand, the optimisation engines in our paper not only use different objective functions, but also different levels of data aggregation. Another differentiating factor from the common hybridisation approaches is the method of interaction between the hybridised optimisation engines, which in our procedure is designed as an optimisation problem itself. 3. Problem Description Even though the problem considered in this paper is motivated by the Hunter Valley Coal Chain of Australia, the optimisation procedure presented is not limited to the operations of Hunter Valley Coal Chain or to coal supply chains in general. Our procedure is equally applicable to other commodities, such as iron ore. Thus, in this section we provide a general description of the problem and the commodity transported in the supply chain will be referred to interchangeably as the mineral or commodity. Figure 1 shows a simple schematic of the mineral supply chains and various resources involved in transporting the mineral. The supply chain is comprised of several production facilities, referred to as mines. Each mine 7

Figure 1: Various processes within mineral supply chains

supplies minerals using a train loading station, referred to as a loading point. A mine can use only one loading point, but several mines can use the same loading point. A loading point can load one train at a time. A railway network transports the mineral from the loading points to one or more terminals. Each ship arrives at its designated terminal according to a predetermined schedule to pick up minerals. At each terminal, the ships normally berth in non-decreasing order of their arrivals [31]. The mineral for each ship is stored at its designated terminal in one or several stockpiles. The mineral delivered by each train is intended for only one stockpile. A stockpile is built using minerals from one or several mines, and the amount required from each mine is predetermined. Each stockpile occupies a part of a strip of land called a pad, and cannot be split between several pads or between several disjoint locations on the same pad. A pad can contain several stockpiles. The sizes of pads are an important characteristic of the considered supply chains. In order to minimise disruptions, it is necessary that the train trips associated with all stockpiles of a ship must be completed before the ship berths at the terminal [31]. Furthermore, prior to departure of the first train trip associated with a stockpile, the entire space needed for this stockpile is reserved. There are several types of wagons and a train is formed by wagons of 8

the same type. A train type is defined by the amount of mineral that it can transport, and the type and number of wagons that form it. As mentioned, a railway network is used to transport mineral from loading points to terminals. Each node of the railway network is either the loading point of a mine, the unloading point of a terminal, or a junction. A junction is a point where different routes merge and/or split. Each loading point prefers trains of a certain type due to technical restrictions such as track length. For any loading point and any terminal, there is a unique route (sequence of nodes in the railway network) which connects them. When a train unloads mineral at a terminal, two types of machines, dumper (or in-loader) and stacker, are used to unload this train. A dumper permits a train to unload its cargo on to the conveyor belt of a stacker, which in-turn builds a stockpile. A dumper and stacker combination can serve only one train at a time. At any point in time, a dumper can be paired with only one stacker, and a stacker can be paired with only one dumper. Each stacker has access to one or several pads in the terminal. For loading, a ship needs a berth. Each terminal has several berths. A berth can only be used by one ship at a time. During loading of a ship, a special device, reclaimer (out-loader), takes mineral from a stockpile and sends it to another device, shiploader, using a conveyor belt. A reclaimer can work on only one stockpile at a time. Stockpiles for the same ship are loaded in a pre-determined order. Each reclaimer can transfer mineral to several shiploaders, but is connected to only one shiploader at a time. Each reclaimer has access to one or several pads. There are also machines (stackerreclaimers) that can be operated as either a stacker or reclaimer, but not both at the same time. The data which is used for validation of the presented optimisation procedure was provided by the Hunter Valley Coal Chain, which is the world’s largest exporter of coal at more than 100 million tonnes of coal per year. [4] describes this coal supply chain as “. . . a large, multifaceted, and dynamic system consisting of several interacting subsystems that are in themselves highly complex.” The reader is referred to Section 7 for further information on the scale of the problem. 4. Optimisation Procedure This section presents the overall design of the optimisation procedure, describing the input and output of both optimisation engines, and the hy9

bridisation mechanism – the method of interaction between the MILP and the metaheuristic based scheduler. Figure 2 complements the text below by presenting the procedure graphically. As has been mentioned above, the aim of the entire procedure is to minimise the cost of capacity expansion required to meet the demand, which is specified by the information of ships’ arrivals, the due dates of their departures, and their cargo in terms of stockpiles and the content of each stockpile. Of course, if a scenario is infeasible, i.e. there exist no suitable expansion of infrastructure, the developed optimisation procedure terminates without specifying the capacity expansion. The optimisation procedure is iterative in nature and, at each iteration, sequentially executes both optimisation engines – first, the MILP, and then the metaheuristic based scheduler.

4.1. Groups of Resources and Units of Expansion The developed optimisation procedure determines the capacity expansion for various parts of the supply chain infrastructure, where each part is characterised by its type and location and is referred to as a group of resources. Thus, for example, a group of resources could be all wagons of the same type, all equipment at a loading point, or all stackers of a terminal. Expansion of capacity, at a cost, is allowed for each group of resources, and once added, the expansion is assumed to be available from the beginning of the planning horizon. Different groups of resources have different measurement of their capacity and different forms of expansion. For example, the expansion can be an additional piece of equipment, or additional working hours, or a replacement of the existing equipment by a more efficient one. Therefore, the presented approach utilises the notion of a unit of expansion. The contribution of a unit of expansion for each group of resources is specified in the input by the end–user, and can be measured, for example, in wagons, hours, tonnes, etc. This flexibility that allows the end–user to define the contribution of a unit of expansion makes our optimisation procedure applicable to a variety of mineral supply chains and alternatives for capacity expansion. 4.2. Mixed Integer Linear Programming Engine The MILP model in the MILP engine is specifically designed for the considered type of supply chains and is based on the partitioning of the 10

Figure 2: The structure of the hybrid optimisation procedure.

11

planning horizon into disjoint time periods which may have different lengths and depend on the ships’ arrivals and departures. The detailed description of the partitioning method is presented in Section 5. The MILP optimisation engine, with its global perspective of the supply chain, determines the cost effective capacity expansion under the requirement that any ship’s departure is not delayed. The input of the MILP engine is data on existing infrastructure, data specifying the demand, and the output of the hybridisation mechanism. The data on demand is ship arrivals and departures and their cargo. For each group of resources, the data on the infrastructure provides the existing capacity and the contribution of one unit of expansion. Using the input, the MILP engine partitions the planning horizon; computes for each group of resources g and each period p a parameter Fgp – the contribution made in period p by a unit of expansion; and generates the MILP model and solves it. In this model, the additional capacity in period p is expressed as eg × Fgp , where eg is a variable specifying the number of units of expansion for the group g. The MILP model is designed in such a way that its output is not only the cost effective expansion, but also the corresponding flows of commodity in the form of train trips and their loads. The detailed description of the MILP model is given in Section 5. For each period p and each group of resources g, the MILP model contains a constraint which restricts the total utilisation of this group to be no greater than the available capacity in the considered period. For example, for each group of dumpers and each period, there is a constraint which limits the total tonnage of trainloads by the available capacity of the group in the considered period. In this constraint, the available capacity is measured as the sum of the original capacity plus eg × Fgp . As another example, the total utilisation can be the number of trains which are assigned to pass a junction in the considered period. In this constraint, the availability of the resource is the originally given maximum number of trains that can pass the junction, plus the additional number due to the expansion. 4.3. Scheduling Engine The deficiency of MILP in modeling sequences of operations is well known [27]. Thus, the above mentioned constraints, which are standard in the literature on capacity planning with MILP, lead to the situation where a MILP solution, which is perfectly feasible for the corresponding MILP model, may not withstand more precise scrutiny. For example, suppose that according 12

to the MILP solution, several trains should pass some junctions during a certain period of time. Being part of the solution, this flow of commodity is feasible, which in terms of the MILP, means that the number of trains which pass each junction does not exceed the maximum number for this junction during the considered period of time. However, in reality, the trains may form a queue at a particular junction. In other words, the simple summation of the trains overestimates the available capacity since it does not take into account the unavoidable idle times. Another simple example that demonstrates the overestimation of the available capacity by ignoring the unavoidable idle times is the assignment of three equal trainloads to a group, comprising two dumpers, for processing during a time period of duration T . If the duration of processing for each trainload is 23 T , then the three trainloads require 2T of processing time. Since the group is comprised of two dumpers, this is exactly the dumper time available in the considered time period, which has length T . Therefore, from the viewpoint of the MILP model, this assignment of trainloads to the dumpers is feasible. However, since two dumpers cannot process the same trainload, there is no feasible schedule which permits to process these trainloads in T units of time. Since a MILP solution may prove to be inadequate in reality, the optimisation procedure uses the output of the MILP as the input to the metaheuristic based scheduler, which performs more detailed analysis. In contrast to the MILP, the scheduler permits violation of the ships’ departure due dates and aims at minimising the total tardiness with respect to these due dates. The scheduler output is the total tardiness and the detected bottlenecks, where a bottleneck is a group of resources and a time period in which the scheduler detected insufficient capacity for this group. The presented implementation of the scheduler is based on simulated annealing and its detailed description is given in Section 6. 4.4. Hybridisation of Two Engines The presented hybridisation of the two engines – MILP and scheduler, has two aspects: the impact of the MILP on the scheduler, and the impact of the scheduler on the MILP. The MILP impacts the scheduler by the fact that its output is used as input for the scheduler (see Subsections 6.1 and 6.2). As has been mentioned above, the MILP output is comprised of the information on the expansion of the existing infrastructure as well as the corresponding flows of commodity (train trips). The scheduler uses these 13

flows in two ways: in constructing an initial schedule, and in reducing the search space. A description of the methods used for constructing an initial schedule and the reduction of the search space requires detailed information on the scheduling engine, and therefore is postponed until Section 6. The impact of the scheduler on the MILP engine is the modification of constraints in the MILP model, based on bottlenecks detected by the scheduler. The part of the optimisation procedure which is responsible for modifying the constraints will be referred to as the constraints modification module. The input for this module is the bottlenecks produced by the scheduler. The output of the constraints modification module is a set of parameters which are used by the MILP engine in generating the MILP model. Each parameter, which will be referred to as a feedback parameter, is associated with a particular group of resources and a particular time period, and specifies how the upper limit on the total utilisation of this group in the considered time period is calculated. The remainder of this section describes the constraints modification in detail. 4.5. Feedback Matrix As has been mentioned, each bottleneck, detected by the scheduler, specifies a group g with insufficient capacity, and a period p where this insufficiency was detected. Hence, a bottleneck can be associated with a pair (g, p). The bottleneck (g, p) causes a modification of the constraint which limits the total utilisation of group g in period p. In the presented optimisation procedure, the modification of this constraint is the replacement of the available capacity by its product with 0.9. In the course of optimisation, the same constraint can be modified several times. Therefore, the constraint modification module keeps, for each constraint that limits the capacity of a group g in a time period p, a feedback parameter δgp = (0.9)k , where k is the number of times the constraint has been modified. In other words, at each iteration, the MILP model uses (Cgp + eg Fgp ) × δgp as the available capacity of group g in period p, where Cgp is the original capacity of group g in period p; eg is the expansion measured in units; and Fgp is the contribution of one unit of expansion to the period p. The formula (Cgp + eg Fgp ) × δgp serves as the available capacity only in the MILP model, whereas the scheduler uses the actual capacity Cgp + eg Fgp . The parameters δgp form a two-dimensional matrix. Since every element of this matrix is a feedback parameter, the entire matrix will be referred to as a feedback matrix. At the start of the first iteration, all feedback parameters 14

are set to 1, signifying that no modifications have been made. Observe that this initial setting is in accord with the formula δgp = (0.9)k , which for k = 0, i.e. no changes have been made, gives δgp = (0.9)0 = 1. 4.6. Hybridisation as a Bi-criteria Optimisation Problem As has been mentioned in the Introduction, one of the distinct features of the presented optimisation procedure is the treatment of the hybridisation as an optimisation problem in its own right. More specifically, the hybridisation is viewed as a bi-criteria problem on the set of feedback matrices where the first objective function is the cost of capacity expansion determined by the MILP engine, and the second objective function is the total tardiness determined by the scheduler for this capacity expansion. Hence, the goal is to find a feedback matrix which corresponds to a point on the Pareto curve with acceptable expansion cost and acceptable total tardiness. Of course, the meaning of the word ‘acceptable’ may change from implementation to implementation. The difficulty of this problem arises, in particular, from the fact that computing both objective values requires the execution of both optimisation engines. In this paper, the hybridisation is implemented as follows. The heuristic is a search procedure, where the search is comprised of several blocks. Each block is a sequence of consecutive iterations, where each iteration is associated with a corresponding feedback matrix. Within a block, each feedback matrix but the first is generated by using the previous feedback matrix and the bottlenecks corresponding to this previous feedback matrix. Thus, each iteration in the block is comprised of the generation of a new MILP model, using the current feedback matrix; solution of this model; the subsequent execution of the scheduler; and the modification of the feedback matrix using the bottlenecks detected by the scheduler. For the first iteration of the first block, the current feedback matrix is comprised of all its entries as 1. The choice of feedback matrices for the first iteration of remaining blocks is described later. The scheduler may identify more than one bottlenecks. In the presented implementation of the hybridisation mechanism, only one bottleneck is selected to modify the feedback matrix, since the selection of one bottleneck may fix the others. For the chosen bottleneck (g, p), it may be reasonable to multiply by 0.9 not only δgp , but also some other feedback parameters δgp , where p is a time period in some neighbourhood of time period p. The programming implementation of our optimisation procedure allows to change 15

the definition of this neighbourhood which, at the extreme, can be reduced to p. The choice of the bottleneck is made at random. The resultant new matrix is the current matrix at the next iteration provided that the block is not terminated. The termination of blocks and the choice of a current matrix for the first iteration in each block is based on the idea of dominance. This idea emanates from the bi-criteria nature of the hybridisation mechanism. A feedback matrix dominates another if its associated expansion cost and total tardiness are not greater than those of the second, and one of these two values of the first is strictly less than the corresponding value of the other. Our heuristic stores all non-dominated feedback matrices generated so far. Each non-dominated feedback matrix is stored together with the associated bottlenecks which have not been used previously in the modification of this matrix. If a pre-determined number of consecutive iterations fails to produce a new non-dominated matrix, the current block is terminated. If the block is terminated, one of the stored non-dominated matrices is chosen at random and modified using one of the associated bottlenecks. The resultant modified matrix is used as the current matrix for the first iteration of the new block. The list of the stored non-dominated matrices and the associated sets of bottlenecks are updated in the course of optimisation. Thus, if a new matrix dominates the previously stored non-dominated matrices, the latter is eliminated from the list. The matrix is also eliminated from the list if all bottlenecks associated with this matrix have been already used for its modification. Each time when a bottleneck is used for the modification of the corresponding previously stored matrix, this bottleneck is removed from the associated set of bottlenecks. The optimisation procedure in our paper allows different stopping criteria, such as maximum possible number of iterations, maximum possible computer time or total tardiness within a pre-specified tolerance. The procedure also terminates if the list of non-dominated matrices is empty at the step that requires a non-dominated feedback matrix to be chosen. In the computational experiments presented in this paper, a limit was imposed on the duration of the optimisation procedure, namely that the constraints modification at an iteration could commence only if the time lapsed from the beginning of computations is less than 36 hours.

16

5. Mixed Integer Linear Programming Model The MILP model used by the MILP engine is presented in this section. This model aims to determine, via MILP’s global view of both the planning horizon and supply chain, a cost effective capacity expansion, under the requirement that all ship departures are not delayed. The presented MILP model partitions the planning horizon into non-overlapping time periods, with the partition designed to reduce the number of variables and constraints used by the model. No artificial restrictions are placed on the model. In particular, the model permits a train trip to lapse over more than one time period. The formulation of this MILP is based on the following observation. Points in time when ships arrive and depart split the planning horizon into non-overlapping time periods. Within such a time period the availability of ships for loading and the list of stockpiles which can be built remain unchanged. The main idea is to associate each operation (and therefore the corresponding variable) with a period or adjacent periods where the operation occurs. Since some periods can be very short, a special procedure, described below, modifies (aggregates) the partitioning of the planning horizon. 5.1. Partition of the planning horizon To avoid the necessity of considering a big number of short periods and therefore a big number of variables and constraints associated with these periods, a short period is merged with a following period. An aggregated period, in turn, can be merged again with the next period if required. A time period is deemed short if its length is less than some chosen value ε. In our paper, ε is the duration of the longest round trip of a train, including its loading time. Thus, each round trip will occur either within a period or over two adjacent periods. For each period p, the aggregation procedure determines the earliest arrival or departure time which exceeds tstart + ε, where tstart is the starting time of period p. This earliest arrival p p or departure time determines the end of the current period, and the start of the next period. The aggregation starts with the first period, whose starting point is the start of the planning horizon, and terminates when the end of the planning horizon is reached. If the last period has duration shorter than ε, this period is merged with the preceding one.

17

Set L J W D ST SR PA RC SH SL

– – – – – – – – – –

SLg

–

Qg P Ts P Ls

– – –

Description the set of all loading points the set of all junctions the set of all groups of wagons the set of all groups of dumpers the set of all groups of stackers the set of all groups of stacker-reclaimers the set of all groups of pads the set of all groups of reclaimers the set of all groups of shiploaders the set of all pairs (s, l), where stockpile s contains commodity from loading point l is a subset of SL, comprised of all pairs where the delivery of commodity from loading point l to stockpile s requires the group of resources g the set of stockpiles which utilise the group of resources g the set of all periods in which building of stockpile s is permitted the set of all periods in which loading of stockpile s is permitted Table 1: Sets used in the MILP model

Although each end point of the new aggregated periods are points of arrival or departure of one or more ships, some ships may arrive or depart within a new aggregated period. For a ship which arrives within an aggregated period, this period is the last in which transportation is allowed, and the first period where loading of the ship is permitted. For a ship which departs within an aggregated period, this is the last period in which loading is allowed. 5.2. Capacity Expansion and the Objective Function Capacity expansion is determined for various parts of the supply chain infrastructure. As has been mentioned above, each part is characterised by its type and location and is referred to as a group of resources. The sets formed by the groups of resources are shown in Table 1. Any expansion of capacity comes into effect from the beginning of the planning horizon and is measured in units, which can have different meaning for different 18

groups (see Subsection 4.1). As has been mentioned in Subsection 4.2 the MILP engine calculates the contribution Fgp of one unit to the capacity of group g in period p, based on what the unit is and the duration of this period. Hence, as in Subsection 4.5, the available capacity is computed as (Cgp + eg Fgp ) × δgp , where Cgp is the original capacity of group g in period p; eg is the expansion measured in units; and δgp is a feedback parameter provided by the constraints modification module. The cost of one unit of expansion also depends on what the unit is. For example, the acquisition of a new dumper may result in quite different cost in comparison to the introduction of a number of extra working hours per week. Let τg denote the cost of one unit of expansion for the group of resources g. Then, the objective function can be written as      min τg eg + τg eg + τg eg + τg eg + τg eg g∈L



+

g∈J

g∈W

τg eg +

g∈SR



τg eg +

g∈P A

g∈D



g∈ST

τg eg +

g∈RC



τg eg ,

(1)

g∈SH

or in other words, to minimise the total cost of capacity expansion over all groups of resources (see Table 1). 5.3. Constraints The constraints of the MILP model, described below, uses the notations in Table 1 and Table 2.  p∈P Ts



aslp +

bslp = Msl

∀(s, l) ∈ SL

(2)

p∈P Ts −{pft s}

Constraints (2) ensure that for each stockpile s the required amount of commodity Msl from loading point l is delivered during the set of periods where transportation is permitted for stockpile s, denoted P Ts (see Table 1). The first period pfts where transportation to stockpile s is permitted is excluded from the second summation because train trips which span over two periods cannot finish in this first period. 

usp = Ms

∀s

p∈P Ls

19

(3)

Variable aslp – bslp

–

αslp βslp eg usp zgp xstk gp

– – – – – –

xrcl gp

–

Definition tonnes delivered from loading point l to stockpile s with transportation started and ended in period p tonnes delivered from loading point l to stockpile s with transportation started in period p − 1 and ended in period p the number of train trips used to transport aslp the number of train trips used to transport bslp the number of units of expansion for the group of resources g tonnes loaded onto ship from stockpile s during period p tonnes stored on the group of pads g at the end of period p tonnes stacked by the stacker-reclaimers at the same terminal as the terminal of the group of resources g during period p tonnes reclaimed by the stacker-reclaimers at the same terminal as the terminal of the group of resources g during period p Table 2: Variables in the MILP model

In constraints (3), Ms is the size of the stockpile s in tonnes. These constraints ensure that the entire stockpile is loaded onto the ship during the set of periods where loading is permitted for the stockpile, denoted P Ls . Sslmin αslp ≤ aslp ≤ Sslmax αslp Sslmin βslp

≤ bslp ≤

Sslmax βslp

∀(s, l) ∈ SL,∀p ∈ P Ts ∀(s, l) ∈ SL,∀p ∈ (P Ts −

(4) {pfts })

(5)

Constraints (4) and (5) determine the number of train trips used to transport aslp and bslp . The parameters Sslmin and Sslmax refer to the minimum and maximum tonnes of commodity which can be carried by the type of train that is preferred for transporting commodity from loading point l to stockpile s. In order to simplify the presentation and facilitate the reading, the description of constraints below assumes that each sum contains only variables which are defined for the considered indices. More specifically, for each stockpile s • aslp and αslp are defined only where (s, l) ∈ SL and p ∈ P Ts , i.e. where the stockpile requires commodity from loading point l and transportation is permitted for the stockpile in period p. 20

• bslp and βslp are defined only where (s, l) ∈ SL and p ∈ (P Ts − {pfts }), where pfts is the first period in which transportation is permitted for stockpile s. As has been mentioned above, the elimination of pfts is based on the fact that bslp and βslp correspond to transportation over two adjacent periods and hence cannot complete in the first period pfts . • usp are defined only where p ∈ P Ls , i.e. where loading is permitted for the stockpile in period p. Therefore, for example in (6), the sum in the left-hand side contains, for each p, only variables aslp which are defined for the considered combinations of indices. In other words, for each p, the left-hand side of (6) contains only the variables which are associated with the stockpiles s that are permitted to be built in period p. 

aslp ≤ (Cgp + eg Fgp )δgp

∀g ∈ L, ∀p

(6)

aslp ≤ (Cgp + eg Fgp )δgp

∀g ∈ D, ∀p

(7)

aslp − xstk gp ≤ (Cgp + eg Fgp )δgp

∀g ∈ ST, ∀p

(8)

αslp ≤ (Cgp + eg Fgp )δgp

∀g ∈ J, ∀p

(9)

aslp × Rsl ≤ (Cgp + eg Fgp )δgp Sg

∀g ∈ W, ∀p

(10)

(s,l)∈SLg





(s,l)∈SLg

(s,l)∈SLg





(s,l)∈SLg

(s,l)∈SLg

Constraints (6) to (10) restrict the total capacity in each period for the loading points, dumpers, stackers, junctions and wagons. In each constraint, the left-hand side is the total utilisation of the group of resources in the considered period. The parameter Cgp is the original capacity of the corresponding group of resources in the period p; Fgp is the contribution of a unit of expansion to the capacity in period p; and δgp is the feedback parameter. Hence, the right-hand side is the available capacity in period p (see Subsections 4.2 and 4.5 for more details). According to the definition in Table 1, SLg is a specific subset of SL. For example, in (6), SLg is the set of all pairs (s, l) such that the loading point 21

l is g (each loading point is considered as a separate group of resources). As another example, SLg in (10) is the set of all pairs (s, l) for which the transportation of commodity from loading point l to stockpile s requires wagons from the group of wagons g. The units of measurement in (6) to (8) are in tonnes, while in (9) it is the number of trains during the considered period. The left- and right-hand side of (10) are measured in wagon-hours. The wagon-hours are computed based on the length of the considered period. As has been mentioned in Subsection 4.1, all stackers at a terminal form a distinct group of resources. All stacker-reclaimers of a terminal form another distinct group of resources. Observe that in (8), the first term in the lefthand side is the amount of commodity which should be stacked in the period p by stackers and stacker-reclaimers together at the terminal corresponding to the group of stackers g. The second term is the part of this amount which is assigned to stacker-reclaimers at this terminal rather than to the stackers. In the left hand side of (10), i.e. in the expression measuring the number of wagon-hours utilised, the parameter Sg denotes the tonnes per wagon of wagon-group g and Rsl denotes the shortest possible duration of a round trip between loading point l and stockpile s, including the lower bound estimate of the time to load trains. In what follows, we assume that the periods of the planning horizon are numbered starting from the first period, i.e. the first period is numbered as 1. Therefore, for example, aslp−1 is the total amount of commodity transported from loading point l to stockpile s in the period which immediately precedes period p. 

(aslp−1 + aslp + bslp ) ≤

p 

(Cgh + eg Fgh )δgh ∀g ∈ L, ∀p ≥ 2 (11)

(s,l)∈SLg

h=p−1



p 

(aslp−1 + aslp + bslp ) ≤

(s,l)∈SLg



(Cgh + eg Fgh )δgh ∀g ∈ D, ∀p ≥ 2 (12)

h=p−1

stk (aslp−1 + aslp + bslp ) − (xstk gp−1 + xgp ) ≤

p 

(Cgh + eg Fgh )δgh

h=p−1

(s,l)∈SLg

∀g ∈ ST, ∀p ≥ 2 22

(13)



(αslp−1 + αslp + βslp ) ≤

(s,l)∈SLg

(Cgh + eg Fgh )δgh ∀g ∈ J, ∀p ≥ 2 (14)

h=p−1

(s,l)∈SLg



p 

p  aslp−1 + aslp + bslp × Rsl ≤ (Cgh + eg Fgh )δgh Sg h=p−1

∀g ∈ W, ∀p ≥ 2

(15)

Constraints (11) to (15) are similar to constraints (6) to (10) and differ only in that they consider pairs of adjacent periods. Thus, these constraints restrict the total capacity, in each pair of adjacent periods, for the loading points, dumpers, stackers, junctions and wagons. Each right-hand side is the total available capacity of group g in periods p − 1 and p. Each lefthand sides is the total utilisation of the group over the periods p − 1 and p. Since these constraints consider pairs of adjacent periods, the left-hand sides of the constraints contain the variables bslp and βslp , which correspond to transportation starting in period p − 1 and completing in p (see Table 2). Otherwise, the constraints (11) to (15) are similar to constraints (6) to (10). For any three time periods, p−1, p and p+1, the above constraints restrict the utilisation of a group of resources in the adjacent periods p − 1 and p, and in the adjacent periods p and p + 1. Therefore, the available capacity in period p is considered twice, which may lead to the overestimation of available capacity. In order to reduce the occurrence of overestimation of the available capacity, the MILP model contains also constraints (16) – (20). p  

p  

aslh +

(s,l)∈SLg h=1

bslh ≤

(s,l)∈SLg h=2

p 

(Cgh + eg Fgh )δgh ∀g ∈ L, ∀p ≥ 3

h=1

(16) p  

p  

aslh +

(s,l)∈SLg h=1

bslh ≤

(s,l)∈SLg h=2

p 

(Cgh + eg Fgh )δgh ∀g ∈ D, ∀p ≥ 3

h=1

(17) p   (s,l)∈SLg h=1

aslh +

p  

bslh −

p 

(s,l)∈SLg h=2

h=1

xstk gh ≤

p 

(Cgh + eg Fgh )δgh

h=1

∀g ∈ ST, ∀p ≥ 3 23

(18)

p  

p  

αslh +

(s,l)∈SLg h=1

βslh ≤

(s,l)∈SLg h=2

p 

(Cgh + eg Fgh )δgh ∀g ∈ J, ∀p ≥ 3

h=1

(19) p   aslh (s,l)∈SLg h=1

Sg

× Rsl +

p   bslh (s,l)∈SLg h=2

Sg

× Rsl ≤

p 

(Cgh + eg Fgh )δgh ∀g ∈ W, ∀p ≥ 3

h=1

(20) Constraints (16) – (20) are similar to (6) – (10) and (11) – (15). In contrast to (6) – (10), which restrict capacity in single time periods, and in contrast to (11) – (15), which restrict the capacity in pairs of adjacent time periods, the constraints (16) – (20) restrict the capacity utilisation from the beginning of the planning horizon till each period p in the planning horizon. In other words, the capacity utilisation is restricted over the periods {1, 2, 3}, {1, 2, 3, 4}, {1, 2, 3, 4, 5}, . . .. Observe that the capacity restriction for period {1} and periods {1, 2} has already occurred in the previously considered constraints. 

bslp ≤ 2Rgmax Ag

∀g ∈ L, ∀p ≥ 2

(21)

(s,l)∈SLg

Constraints (21) are based on the following observation. If a train trip, transporting commodity from loading point g, spans over two adjacent time periods p and p + 1, then this trip should start no earlier than tend − Rgmax , p max where Rg is the duration of the longest round trip associated with loading point g and tend is the point in time when period p ends. Similarly, this train p max . Hence, all such train trips should belong trip is completed before tend p +Rg to the time interval of length 2Rgmax . If Ag is the loading rate at loading point g, then the right-hand side of (21) is an upper bound on the total amount of commodity which can be transported by such train trips. zgp ≤ (Cgp + eg Fgp )δgp

∀g ∈ P A, ∀p

(22)

Constraints (22) limit the utilisation of each group of pads g in each time period p. Observe that for each g ∈ P A, Cgp are identical for all periods p. Furthermore, Fgp are also identical for all period p. 24



zgp = zgp−1 +

aslp +

(s,l)∈SLg



bslp −

(s,l)∈SLg



usp

∀g ∈ P A, ∀p

(23)

s∈Qg

Constraints (23) are balance equations for pads. For each group of pads g and each period p, the left-hand side is the total amount of commodity stored, at the end of period p, on pads constituting group g. The right-hand side is the amount of commodity stored at the end of the previous period, plus total delivery during the period p, minus total amount loaded onto the ships during the period p. Observe that for p = 1, the variable zgp−1 does not exist and therefore, according to the assumption stated earlier, this variable does not contribute to the right-hand side. As per Table 1, Qg is the set of stockpiles that should be built on the pads constituting g.  

usp ≤ (Cgp + eg Fgp )δgp

∀g ∈ SH, ∀p

(24)

usp − xrcl gp ≤ (Cgp + eg Fgp )δgp

∀g ∈ RC, ∀p

(25)

∀g ∈ SR, ∀p

(26)

s∈Qg

s∈Qg

xstk gp s δ  Cgp gp

+

xrcl gp ≤ (Cgp + eg Fgp )δgp , r δ  Cgp g p

Constraints (24), (25) and (26) impose capacity restrictions on equipment used to load the commodity onto the ships. This equipment is comprised of reclaimers, stacker-reclaimers and shiploaders. Constraints (24) restrict the amount of commodity which is loaded onto ships (left-hand side) to the available capacity of the considered group of shiploaders g (right-hand side). In constraints (24), the unit of measurement is tonnes. All commodity which is loaded should be reclaimed by reclaimers, stackerreclaimers or both. The left-hand side of (25) is the amount reclaimed by reclaimers, i.e. is the amount loaded onto ships minus the amount reclaimed by stacker-reclaimers. The right-hand side is the available capacity of the considered group of reclaimers g in the period p. The units of measurement for these constraints are also tonnes. Constraints (26) restrict the capacity of the stacker-reclaimers. In contrast to (24) and (25) where units of measurement are tonnes, the units of s measurement in (26) are number of stacker-reclaimers. The parameter Cgp 25

in the left hand side is the amount of commodity which one stacker-reclaimer xstk can stack during period p when operating as a stacker only. Hence, Cgps is gp

the number of stacker-reclaimers needed to stack xstk gp tonnes of commodity r in period p. Similarly, Cgp is the amount of commodity which one stackerreclaimer can reclaim during period p when operating as a reclaimer only, xrcl and Cgp is the number of stacker-reclaimers needed to reclaim xrcl r gp tonnes of gp commodity in period p. The feedback parameter δgp in the right-hand side of (26) remains equal to 1 at all iterations of the optimisation procedure, and the detection of bottlenecks impacts constraints (26) through another two feedback parameters, δg p and δg p , in the left hand side. In these two feedback parameters, g  are all stackers and g  are all reclaimers at the same terminal as the considered group of stacker-reclaimers g. If, for example, a bottleneck in stacking is detected and δg p is reduced, the contribution of one stacker-reclaimer to stacking in the time period p is also reduced. Therefore, the same amount of commodity xstk gp would require a larger number of stacker-reclaimers. 6. Scheduling Engine The scheduling engine analyses the capacity expansions provided by the MILP engine and aims to identify bottlenecks in the supply chain which remained undetected by the MILP engine. The underlying scheduling problem solved by the scheduling engine is described in Subsection 6.1. The simulated annealing scheduling algorithm is based on the idea of specifying task priorities by means of two lists. This idea is presented in Subsection 6.2. The simulated-annealing based scheduling algorithm is described in Subsection 6.3. The form in which input data is provided significantly affects the actual implementation of the algorithm described in Subsections 6.2 and 6.3. Some data-specific technical details on the implementation of the presented scheduling algorithm are provided in Subsection 6.4. As has been mentioned in Section 3, the building of a stockpile can commence only if pad space for the entire stockpile is reserved. Furthermore, according to the same section, the loading of a ship can commence only after building all stockpiles for this ship. These two requirements can lead to a situation where neither ship loading nor new deliveries of commodity to the terminal are possible. Such a situation is referred to as a deadlock. The following simplified example describes a deadlock. 26

Suppose that there are two ships, each requiring three stockpiles. Assume that each stockpile occupies exactly a quarter of the total pad space. If the pad space is occupied by two stockpiles of each ship, then the pad space is completely occupied and no more deliveries can be made. On the other hand, none of the ships has all its stockpiles built, and therefore the ship loading cannot commence, and hence pad space cannot be released for new stockpiles. The deadlock prevention mechanism, described in Subsection 6.5, allows to avoid deadlocks. Subsection 6.6 presents the identification of bottlenecks. 6.1. Scheduling Problem Solved by the Scheduling Engine The underlying scheduling problem solved by the scheduling engine depends on several factors, including the MILP model in the MILP engine and the available data. Even the specification of all factors does not exclude the possibility of alternative formulations of the underlying scheduling problem. The solution method also can be chosen differently even for the same formulation. One of the important factors is the choice of units of expansion (variables eg in the MILP model). The choice of the unit of expansion in turn depends on the available data and the needs of the planners. In this paper, the unit of expansion for groups of dumpers, stackers, stacker-reclaimers, reclaimers and shiploaders is the number of additional machines for that group. The unit of expansion for loading points is the additional tonnes per day. For pads, the expansion specifies additional size in metres. For junctions, the expansion is specified in additional trains which can pass the junction per day. The expansion of wagons is specified in additional number of wagons. In our paper, the underlying scheduling problem is based on the output of the MILP engine. Each train trip specified in this output is translated into a linearly-ordered set of tasks (chain of tasks). Each of these tasks represents either the passage of a junction; or loading at the corresponding loading point; or stacking at the terminal using a pair comprised of a dumper and a stacker (or stacker-reclaimer). The tasks in a chain are ordered according to the sequence of events in the course of the train trip. The first task in the chain has a release time – the given time when building of the stockpile can begin. The tasks in each chain can be processed only if the required number of wagons of the required type is allocated to the chain together with the pad space for the entire stockpile (see Section 3 for the rationale behind this 27

Figure 3: Simplified example of a railway network.

assumption). The wagons are released (and can be again allocated) at the start of the last task in the chain, i.e. at the beginning of stacking, since it is assumed that unloading of a train is instantaneous. The release of pad space will be described below. Each stockpile also corresponds to a task – the task of loading the stockpile onto the ship using a pair comprised of a shiploader and a reclaimer (or stacker-reclaimer). The order in which stockpiles should be loaded onto the ship is given and therefore the tasks of loading the stockpiles form another chain. All chains generated by the train trips for the same ship precede the chain of loading tasks, i.e. the tasks of loading the stockpiles onto this ship. The first task in this chain also has a release time which is specified by the ship arrival, and the last task of the chain has the due date – the desired ship departure. For each stockpile, the pad space is released at the end of loading this stockpile. The objective is to minimise the total tardiness with respect to the due dates specified for each ship (in practice, this may mean the total tardiness within the adopted tolerance). The following simplified example illustrates the above text. Suppose that there is only one terminal, three junctions and two loading points. The structure of the railway network for the simplified example is depicted in Figure 3. Suppose also that there is only one ship which will carry two stockpiles, SP 1 and SP 2, and that according to the output of the MILP engine, each stockpile is formed by a single trainload – say, SP 1 from loading point 1 and SP 2 from loading point 2. Then, the output of the MILP engine results in the partially-ordered set of tasks, in the form of an in-tree, depicted 28

Figure 4: The partially ordered set of tasks for two train trips and two stockpiles. The tasks labeled “j1” and “J1” is, respectively, the passage of the first junction when traveling to and returning from loading points (see Figure 3). The task labeled “Stk” is stacking of the delivered commodity.

in Figure 4. 6.2. Scheduling Using Two Lists In contrast to the MILP, which considers infrastructure as groups of resources, the scheduling engine distinguishes the individual members of each group, e.g. individual dumpers, individual stackers, individual shiploaders, etc. For each output of the MILP engine, the minimisation of total tardiness is an iterative procedure (this paper uses simulated annealing). At each iteration of this iterative procedure, a so-called “list schedule” is constructed. As the term list schedule indicates, the construction of this schedule uses the technique commonly referred to as list scheduling [27]. In list scheduling, the construction of the schedule starts at the beginning of the planning horizon and the resources are allocated to the jobs which are ready at this point in time. If the resources are scarce, then the conflict for resources is resolved by allocating the resources to the jobs with higher priority, where the priority is given by the position of the job in some list. After allocating the resources at a point in time, the next point in time is selected where the current situation changes (e.g. new jobs become available, or some previously allocated resources are released) and the allocation of resources repeats. More specifically, in constructing a schedule, for the purpose of assigning wagons and pad space, all train trips are arranged in a list which specifies their priorities. The priority of a train trip is given by its position in the list, where a train trip closer to the beginning of the list has a higher priority. Observe that, in the terminology of Subsection 6.1, wagons and pad space 29

are assigned not to the individual tasks, but to the chains of individual tasks. If several trips (chains of tasks) compete for the same wagons or pad space, the resource is assigned to the chain of highest priority. The output of the MILP engine determines the relative position of the train trips in the list, but leaves some freedom in their ordering. By restricting the train trip positions in the list, the output of the MILP engine is used for the reduction of the search space. More specifically, as has been mentioned in Section 5, the planning horizon is partitioned into time periods, and in the output of the MILP engine, each train trip is associated either with a time period, or with a pair of adjacent time periods. The association with the earlier time period puts a train trip closer to the beginning of the list. Any list used in constructing a schedule for the purpose of allocating wagons and pad space satisfies the following conditions. (R1) A train trip associated with period p − 1 precedes in the list any train trip associated with period p. (R2) A train trip associated with period p − 1 precedes in the list any train trip associated with adjacent periods p and p + 1. (R3) A train trip associated with adjacent periods p − 1 and p precedes in the list any train trip associated with period p. (R4) A train trip associated with adjacent periods p − 1 and p precedes in the list any train trip associated with adjacent periods p and p + 1. The construction of a schedule requires the resolution of conflicts not only related to the allocation of wagons and pad space, but also to other types of resources, i.e junctions, loading points, dumpers, stackers, and stackerreclaimers. In order to diversify the search for a desired schedule, these conflicts are resolved with a separate list. This separate list also specifies the order of the train trips and satisfies the same conditions (R1) – (R4). In the conflict situation, i.e. when several tasks compete for the same resource, such as a junction or a dumper, the resource is allocated to the task that corresponds to the higher priority train trip, i.e. the trip with the left-most position in the list amongst all train trips associated with the competing tasks. In the construction of a schedule, the conflicts for shiploaders, reclaimers and stacker-reclaimers for reclaiming are resolved based on the ships’ due dates, i.e. a smaller due date gives a higher priority. 30

6.3. Simulated Annealing Scheduling Algorithm The description of the simulated annealing scheduling algorithm requires the notions of a feasible pair of lists and the definition of the neighbourhood of a feasible pair. A list is feasible if it satisfies conditions (R1) – (R4). A pair of lists is feasible if both lists are feasible. The neighbourhood of a feasible pair is the set of all feasible pairs such that each is obtained from the original pair by a single interchange of two train trips in one of the original lists. As has been mentioned above, the simulated annealing scheduling algorithm is an iterative procedure. Each iteration uses a feasible pair of lists, which will be referred to as the current pair. At the first iteration, the current pair contains identical lists generated at the initialisation phase which precedes the first iteration. The initialisation phase of the scheduling algorithm generates these identical lists by randomly choosing the order of train trips which is not determined by conditions (R1) – (R4). The initialisation phase also finds the total tardiness of the schedule for the initial current pair. If this schedule has a desired total tardiness, then the scheduling algorithm terminates. The implemented scheduling algorithm follows the standard simulated annealing scheme, which is governed by a parameter commonly referred to as a temperature, which will be denoted by ψ. The value of ψ is updated during the course of optimisation using a certain rule which is commonly referred to as the cooling schedule. Thus, each iteration consists of the following major steps. Step 1 Randomly choose a new pair from the neighbourhood of the current pair and find the total tardiness of the new pair. Step 2 If the schedule for the new pair has the desired total tardiness, terminate the scheduling algorithm. Step 3 If the total tardiness of the schedule for the new pair is less than that for the current pair, choose the new pair as the current pair, and go to Step 7. Step 4 If the number of consecutive iterations without improvement of the total tardiness has reached the threshold, then terminate the scheduling algorithm. Step 5 With some probability which depends on ψ, choose the new pair as the current pair for the next iteration and go to Step 7. 31

Step 6 Choose the current pair as the current pair for the next iteration. Step 7 Update ψ. In Step 1, a new pair is chosen in several steps. First, a list is selected randomly from two lists of the current pair. Second, a pair of train trips is selected randomly in such a way that their interchange in the selected list also results in a feasible list. If possible, these two train trips are chosen from the trips which competed for a resource in the schedule for the current pair. Hence, the new pair is formed by the chosen list with interchanged train trips, and the remaining list of the current pair. In Step 5, the probability of choosing the new pair as the current pair in Tn −Tc the next iteration is computed as e ψ , where Tn and Tc , respectively, is the total tardiness of the schedule for the new pair and the current pair. The cooling schedule (Step 7), implemented in this paper, decreases ψ if one of the following two conditions is satisfied: the number of new pairs accepted at the current value of ψ has hit a specified threshold; or if the number of iterations using the current value of ψ has hit another specified threshold. 6.4. Resource Allocation The actual computations in the course of constructing a schedule depends on the form of the input data. For example, the data provided by the Hunter Valley Coal Chain specifies, for each junction, the total number of trains that can pass through this junction per day (throughput capacity). In addition, the data provided by the Hunter Valley Coal Chain also specifies the time required to pass through each junction. The structure of the provided data dictated the following method of calculating the corresponding processing times. In the presented implementation, after a train enters a junction, this junction is blocked for the duration calculated as a ratio of the duration of a day and the throughput capacity. Furthermore, in the presented implementation, once a train enters a junction, it starts to compete for the next junction after the time required for this train to pass the current junction. Another example is loading points. The data provided by the Hunter Valley Coal Chain specifies, for each loading point, two parameters: the loading rate, in tonnes per hour, and the daily capacity, in tonnes per day. Correspondingly, the duration of loading is calculated using the first of the above-mentioned parameters, whereas the total load permitted, and therefore the number of trains which can be loaded, is determined using also the second parameter. 32

Observe that the number of wagons required for the train trip is determined by the information provided by the MILP. More specifically, the number of wagons is calculated using the train load, which in turn is calculated as a ratio of the variables aslp /αslp or bslp /βslp . 6.5. Deadlock Prevention Mechanism This subsection describes the deadlock prevention mechanism adopted in the presented implementation of the scheduling engine. Computational experiments with data provided by the Hunter Valley Coal Chain showed that most of the time, the deadlock prevention mechanism is not required. However, if the simulated annealing based scheduling algorithm is completely unable to produce a resultant schedule, then the deadlock prevention mechanism is activated. If a deadlock occurs for the new pair of lists in Step 1 of the simulated annealing scheduling algorithm (see Subsection 6.3), the deadlock prevention mechanism identifies the earliest arriving ship, say ship i, amongst all ships for which commodity cannot be delivered. It also identifies all ships which arrive later and use the same terminal. The set of these ships will be denoted V . If several ships at the same terminal have the same arrival time, the tie is broken arbitrarily, and this ordering of the ships is maintained through the entire optimisation process. The schedule is constructed again, with the additional restriction that any attempt to reserve pad space for stockpiles of any ship in V immediately triggers an attempt to reserve pad space for all stockpiles for ship i. Hence, the allocation of the pad space to any stockpiles of ships in V is possible only if the ship i is allocated space sufficient for all its stockpiles. The deadlock prevention mechanism is illustrated by considering again the example at the start of this section. Suppose that Ship 1 arrives before Ship 2. If the deadlock, described in this example, occurs, then the additional restriction is introduced that the reservation of pad space for the stockpiles of Ship 2 immediately triggers the reservation of pad space for all stockpiles of Ship 1. Thus, as soon as the first stockpile of Ship 2 is scheduled for delivery, the space for the entire three stockpiles required by Ship 1 is reserved. Therefore, all three stockpiles for Ship 1 will be delivered, and when loaded, will release pad space for the remaining stockpiles of Ship 2. If a new deadlock occurs, similar deadlock prevention restrictions are added to the existing ones. At each iteration of the simulated annealing

33

scheduling algorithm, once the total tardiness is found for the new pair of lists, Step 1 ends and all deadlock prevention arrangements are discarded. 6.6. Identifying Bottlenecks Each bottleneck is specified as a group of resources g and a time period p. In the presented implementation, the bottlenecks are detected by analysing a predetermined number of the best schedules, obtained in the course of the simulated annealing search. The analysis based on several schedules (ten, in the presented implementation) instead of only the best one makes the identification of bottlenecks more reliable. In detecting bottlenecks, the average total waiting time over all analysed schedules is computed for each group of resources. In this computation, only waiting times related to the ships which are late are considered. In other words, if a train waits at a junction, but the corresponding ship nevertheless is not late, then the waiting time of the train does not contribute to the waiting time at this junction. Furthermore, the waiting times are calculated only after a certain point in time, which is referred to the as the end of the warm up period. The warm up period is introduced in order to consider only the typical behaviour of the supply chain which is not influenced by the initial conditions. Several groups with the largest average waiting times are selected for the identification of bottlenecks. This determines the groups of resources for the bottlenecks. The calculation of waiting times caused by insufficient capacity of shiploaders, pads and insufficient number of wagons require additional explanation. For shiploaders, the waiting time of a ship is calculated only from the completion of all its stockpiles or the ship arrival, whichever occurs later. According to the scheduling algorithm, a large number of train trips (chains of tasks, in the terminology of Subsection 6.1) can wait simultaneously for the allocation of pad space. This can occur if at some point in time there are several stockpiles ready for building. Although the simultaneous departure of any number of trains is feasible from the viewpoint of the scheduling algorithm, this never can happen in the real world and will immediately create queues at junctions, loading points, and other groups of resources. So, if all train trips, waiting for pad space simultaneously, contribute to the total waiting time, then the necessity to expand the pad space may be exaggerated. Furthermore, after the departure of the first train trip for some stockpile, all other train trips for this stockpile cease waiting for the pad space. Thus, their contribution to the total waiting time can also 34

exaggerate the necessity of expansion. Therefore, the train waiting time is calculated only for the train trip of highest priority. The point in time when the reservation of pad space is made is the starting point for calculating the waiting time of the next highest priority train amongst those waiting for this pad space. The train trip may also be delayed by the absence of the required wagons. For the same reason as above, the waiting time caused by the absence of wagons is calculated only for the trip of highest priority. For the trips which are the first delivery to a stockpile, the waiting time caused by the absence of wagons is computed only from the point in time at which pad space is allocated. The contribution of the waiting times to each period p is calculated as follows: • For loading points, junctions, dumpers, stackers, pads, and wagons, the waiting time incurred by a train trip associated with period p contributes to the waiting time of period p. A half of the waiting time incurred by a train trip associated with the adjacent periods p − 1 and p contributes to period p − 1 and the other half of this waiting time is contributes to period p. • If a stockpile waits for a reclaimer, its waiting time is assigned proportionally to all periods in which the loading of this stockpile occurred in the solution of the MILP. • For shiploaders, waiting time of each ship is distributed proportionally to all periods in which the loading of this vessel occurred in the solution of the MILP. 7. Computational Experiments The optimisation procedure presented in this paper was tested using data provided by the Hunter Valley Coal Chain. The results of the experiments show that the procedure is capable of solving industrial size problems. Table 3 gives information on initial capacities. The numbers in this section are accurate to two decimal places or three significant places, or is otherwise exact if there are no decimal digits.

35

Dumpers: 8 Stackers: 13 Stack-Reclaim: 3 Reclaimers: 9 Shiploaders: 7

Total Pad Capacity (tonnes): 4,000,000 Loading points: 33 Loading point avg tonnes per day: 40,909.09

Junctions: 15 Avg trains daily per junction: 62.34 Type 1 wagons: 2046 Type 2 wagons: 248

Table 3: Initial capacity of the supply chain

7.1. Setup of Computational Experiments The computational experiments were conducted for four scenarios, each with 100 ships. Each scenario was transformed into a benchmark scenario as follows. For each scenario, one iteration of the optimisation procedure, described in Section 4, was executed without the phase of bottleneck detection and constraint modification. This results in a schedule S in the output of the scheduler. After that, the due date of each ship departure in the scenario was increased by the tardiness of this ship in the schedule S. This transformation changed the planning horizon of the original scenario and led to a different partition of the planning horizon by the MILP engine, therefore resulting in a new MILP model. Hence, the transformation produces benchmark scenarios with the following important properties: 1. each scenario has a known upper bound on the optimal cost of capacity expansion; 2. for each scenario, there exists a schedule which has a zero total tardiness at the upper bound cost of expansion, i.e. the schedule used to construct the benchmark scenario. As has been described in Section 6, the scheduling engine in the current implementation of the optimisation procedure is based on simulated annealing. In order to address the stochastic nature of simulated annealing, the optimisation procedure was applied to each scenario forty times. Each application of the optimisation procedure results in one or more non-dominated pairs of expansion cost and associated total tardiness. The discussion on the bi-criteria nature of the optimisation procedure can be found in Section 4. The output of the optimisation procedure was compared with the known benchmark results, i.e. the upper bound on the cost of expansion and zero total tardiness. The computational experiments were conducted using a random number generator commonly referred to as Mrg32k3a and originally published in 36

Scenario 1 2 3 4

Number of Stockpiles 158 150 160 152

Horizon Length (days) 43.94 52.33 70.73 53.40

Benchmark Cost ($) 9.26E+06 2.67E+07 2.09E+07 6.00E+06

Table 4: Characteristics of the Benchmark Scenario

[16]. One of the biggest benefits of this generator is its ability to produce, using only a single initial seed, multiple sequences of random numbers each of length up to 2121 [19]. The Mrg32k3a generator is widely used in research and in industrial applications and can be found, for example, in such leading software as Matlab, Arena and R [17]. The choice of the Mrg32k3a generator for the computational experiments presented in this paper permitted to avoid distortion in the results caused by possible correlation in the sequences of random numbers. The implementation of the generator utilised in the presented experiments was from the “Stochastic Simulation in Java” package [18]. The optimisation procedure was set to terminate when it required a nondominated feedback matrix and yet the list of non-dominated feedback matrices has been emptied. In other words, the procedure terminates when all non-dominated feedback matrices have been completely explored. The discussion on exploring the non-dominated feedback matrices can be found in Section 4. Additionally, a limit was imposed on the duration of the optimisation procedure, namely that the constraints modification at an iteration could commence only if the time lapsed from the beginning of computations is less than 36 hours. The computational experiments showed that out of 160 runs of the optimisation procedure only 15 were terminated due to time limits. The testing of the optimisation procedure was conducted on 64-bit 28x dual 3.2 GHz Xeon machines with 8GB of virtual memory. The optimisation procedure was implemented in Java 1.7.0, and the MILP solved using CPLEX Concert Technology 12.5. CPLEX was set to use the deterministic branch and bound mode, with an optimality gap tolerance of 5%, and limited to use at most 4 CPU threads. 37

7.2. Results As has been discussed above, each scenario is associated with a benchmark cost. Thus, in the text which follows, any cost of capacity expansion produced by a run of the optimisation procedure can be represented as the difference from the corresponding benchmark cost, and expressed as a percentage of - benchmarkCost said benchmark cost, i.e. producedCost × 100. This representation benchmarkCost has the benefit of giving a direct comparison against the upper bound cost of capacity expansion. For convenience, such a percentage will be referred to as a Percentage Indicator. As mentioned in Subsection 7.1, each run of the optimisation procedure produces one or more non-dominated pairs, where each pair is an expansion cost and the total tardiness associated with the expansion. Depending on the requirements of a decision maker, the definition of acceptable total tardiness may differ, and so a useful capacity expansion might not be restricted to the non-dominated pair associated with zero total tardiness, and instead might be another non-dominated pair whose associated total tardiness is within the acceptable threshold. In the text which follows, let the least-cost pair of a run of the optimisation procedure refer to the non-dominated pair with the smallest cost of capacity expansion amongst the non-dominated pairs (produced by this run) with associated total tardiness within the acceptable threshold. Of course, if a run of the optimisation procedure fails to produce any non-dominated pairs with associated total tardiness within the acceptable threshold, then such a least-cost pair does not exist for this run at the specified threshold. Statistics such as the mean and median over all runs of a scenario will be used to describe the following results. Thus, Table 5 provides, for each scenario, a summary of the costs of capacity expansion (expressed as a percentage indicator) corresponding to the least-cost pairs that were produced by the runs for the scenario, for different thresholds of total tardiness (first column of the table). The second column of Table 5 shows the number of runs, out of the 40 for each scenario, that were able to produce a least-cost pair for the specified threshold of total tardiness. The statistics for each threshold were calculated only over the runs which produced a least-cost pair for the particular threshold. Thus, for example, the mean for Scenario 1 at a threshold of zero hours of total tardiness (first row) is calculated over the 26 runs which produced a least-cost pair at zero hours of total tardiness.

38

Total Tardiness =0 ≤1 ≤2 ≤3 ≤4 ≤5 ≤ 10

hrs hrs hrs hrs hrs hrs hrs

Number of Runs

Median Cost

Mean Cost

Standard Deviation

26 33 36 36 36 36 38

1.26% 0.70% -0.46% -0.58% -0.58% -0.58% -0.58%

26.60% 13.16% 10.05% 9.49% 9.47% 6.33% 3.41%

52.00% 29.62% 26.45% 26.52% 26.51% 19.76% 12.87%

Total Tardiness =0 ≤1 ≤2 ≤3 ≤4 ≤5 ≤ 10

hrs hrs hrs hrs hrs hrs hrs

(a) Scenario 1 Total Tardiness =0 ≤1 ≤2 ≤3 ≤4 ≤5 ≤ 10

hrs hrs hrs hrs hrs hrs hrs

Number of Runs

Median Cost

Mean Cost

Standard Deviation

3 12 18 22 29 33 39

11.49% 6.21% 3.11% 2.45% 0.65% 0.59% 0.34%

9.49% 8.33% 8.69% 8.63% 7.34% 6.53% 4.74%

8.34% 7.74% 12.58% 12.01% 11.55% 10.20% 7.92%

(b) Scenario 2

Number of Runs

Median Cost

Mean Cost

Standard Deviation

37 40 40 40 40 40 40

5.64% 0.16% -0.07% -0.23% -0.61% -1.06% -1.17%

11.26% 5.03% 4.95% 4.90% 4.77% 4.27% -0.48%

19.14% 14.61% 14.62% 14.64% 14.62% 14.63% 2.00%

(c) Scenario 3

Total Tardiness =0 ≤1 ≤2 ≤3 ≤4 ≤5 ≤ 10

hrs hrs hrs hrs hrs hrs hrs

Number of Runs

Median Cost

Mean Cost

Standard Deviation

40 40 40 40 40 40 40

1.43% 1.43% 1.43% 1.43% 1.43% 1.43% 1.43%

1.43% 1.43% 1.43% 1.43% 1.43% 1.43% 1.43%

0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00%

(d) Scenario 4

Table 5: A summary of the smallest expansion costs (expressed as a percentage indicator) found by the runs of each scenario with associated total tardiness below different thresholds.

39

Considering the first scenario, Table 5 indicates that 26 of the 40 runs for the scenario were capable of producing non-dominated pairs which permitted zero hours of total tardiness. Of these 26 runs, 13 achieved the zero hours of tardiness without exceeding the benchmark cost by any more than 1.26%, as indicated by the median cost column. If even 2 hours of total tardiness is permitted, then the number of runs which found non-dominated pairs with acceptable total tardiness increases to 36 of the 40 runs (90% of all runs for the scenario). Furthermore, half of these 36 runs actually produced an expansion cost that is less than that of the benchmark cost by 0.46%, as indicated again by the median cost column. Taking into account the length of the planning horizon which, according to Table 4, was 43.94 days for this scenario, and the 100 ships in each scenario, 2 hours, or even 10 hours, of total tardiness can be viewed as negligible, equating to at most six minutes of delay per ship on average. For Scenarios 1 and 2, 38 runs and 39 runs, out of the total of 40 runs for each scenario, permitted a total tardiness of 10 hours or less. For Scenarios 3 and 4, all but 3 runs amongst the 40 runs conducted for each scenario permitted zero hours of total tardiness, and every run permitted 1 hour of total tardiness or less. Observe that Table 5 was compiled by extracting the relevant data from the results of computational experiments, rather than running the procedure for each threshold. Observe also that the least-cost pair produced by a same run, but for different thresholds of total tardiness, may be different, since each run of the optimisation procedure can produce a number of non-dominated pairs. The number of runs, aggregated over all scenarios, that produced nondominated solutions within different thresholds for total tardiness are shown in Table 6. According to this table, even if the strict threshold of zero total tardiness is enforced, more than 66% of the runs of the procedure were capable of producing a non-dominated pair satisfying this strict restriction. If the threshold is increased to 5 hours, which is only, on average, a three minute increase in the tardiness per ship, then according to Table 6, 149 runs out of 160 (93.13% of total runs) were capable of producing non-dominated pairs satisfying this threshold. Another increase, on average, of three minutes of tardiness per ship, i.e. the threshold of 10 hours, results in an increase of the percentage of 93.13% to 98.13% of all runs. It is important to stress that column 3 of Table 5 shows that the total tardiness, which is an indicator of the feasibility of a capacity expansion, is accompanied by the median costs 40

Total Tardiness

= 0 hrs

≤ 1 hrs

≤ 2 hrs

≤ 3 hrs

≤ 4 hrs

≤ 5 hrs

≤ 10 hrs

Number of Runs Percent of Runs

106 66.25%

125 78.13%

134 83.75%

138 86.25%

145 90.63%

149 93.13%

157 98.13%

Table 6: The number and percentage of runs which found expansions of capacity that permitted total tardiness within different thresholds.

Percentage Indicator

Number of Runs

Median T.T.

Mean T.T.

Standard Deviation

Percentage Indicator

Number of Runs

Median T.T.

Mean T.T.

Standard Deviation

≤ 0% ≤ 5% ≤ 10% ≤ 20% ≤ 30%

40 40 40 40 40

1.88 0.21 0.21 0.10 0

14.68 7.50 7.50 3.57 3.37

23.97 13.87 13.87 8.90 8.89

≤ 0% ≤ 5% ≤ 10% ≤ 20% ≤ 30%

0 40 40 40 40

4.16 4.01 2.43 2.43

16.40 12.12 3.98 3.65

34.16 26.23 5.60 5.40

(a) Scenario 1

(b) Scenario 2

Percentage Indicator

Number of Runs

Median T.T.

Mean T.T.

Standard Deviation

Percentage Indicator

Number of Runs

Median T.T.

Mean T.T.

Standard Deviation

≤ 0% ≤ 5% ≤ 10% ≤ 20% ≤ 30%

40 40 40 40 40

1.43 0.06 0 0 0

4.54 2.02 0.72 0.44 0.42

6.56 3.89 2.05 1.78 1.78

≤ 0% ≤ 5% ≤ 10% ≤ 20% ≤ 30%

0 40 40 40 40

0 0 0 0

0 0 0 0

0 0 0 0

(c) Scenario 3

(d) Scenario 4

Table 7: A summary of the smallest total tardiness (T.T.) found by the runs of each scenario with the expansion cost, expressed as a percentage indicator, not exceeding different thresholds (specified in each row).

which are either very close to the benchmark cost or, in most cases, even better than the benchmark cost. Table 7 also reflects the interplay between the expansion cost and the total tardiness, but in contrast to Table 5, the thresholds are introduced for the expansion cost (expressed as a percentage indicator). Thus, the table provides information on the smallest total tardiness associated with capacity expansions with the cost within the specified threshold. In contrast to Tables 5 and 7, Table 8 specifies a threshold on both the total tardiness (different threshold tardiness per row) and on the expansion cost (different threshold cost per column, expressed as a percentage indicator). The table shows the number of runs (out of the forty for each scenario) which found one or more non-dominated pair that satisfied both thresholds, 41

and Table 9 shows these results graphically. Thus, according to Table 8 and considering Scenario 1 for example: • the first row and column indicates that 7 of the 40 runs found solutions that required no additional cost compared to the benchmark solution and had an associated total tardiness of zero; and • the second row and column indicates that, by permitting to exceed the benchmark cost by no more than 1%, and the total tardiness by no more than 1 hour, the number of runs which found non-dominated pairs for these thresholds increased to 19 runs (almost half). Naturally, these 19 runs include the 7 runs mentioned in the point above. The remainder of this section presents information on the size and durations of the optimisation procedure and the optimisation engines. Although the constraints of the MILP model are constantly modified in the course of optimisation, the number of variables and the number of constraints remain the same for all iterations of the optimisation procedure. The size of the MILP in terms of the number of variables and constraints is given in Table 10. As has been discussed in Subsection 6.1, the scheduling problem solved by the scheduling engine depends on the output of the MILP engine. In particular, the number of train trips in these scheduling problems varies from one iteration of the optimisation procedure to another. Table 10 gives information on the smallest (column 5) and largest (column 6) number of train trips in the formulations of the scheduling problem solved by the simulated annealing scheduling algorithm in the course of optimisation. The average iterations per run of the optimisation procedure, and the average duration of each execution of the MILP engine and scheduling engine are given in Table 11. The average durations are calculated by considering all of the iterations executed by all runs of the scenario. The percentage of time spent by the optimisation procedure on the MILP engine and the scheduling engine are calculated by taking the second and third row from the table, respectively, and dividing by the fourth row. Observe that the time required by the constraints modification module was negligible. The durations in Table 11 were obtained with the CPLEX software, used to solve the MILP model, set to not keep any previous solutions. This setting aimed to facilitate the generation of different MILP solutions and therefore different flows of commodity. As can be seen in the durations in the same table, the majority of the time spent by the optimisation procedure is devoted 42

0 1 2 3 4 5 10

hrs hrs hrs hrs hrs hrs hrs

0%

1%

2%

3%

4%

5%

10%

15%

20%

25%

30%

7 15 20 23 23 24 27

11 19 24 27 27 28 31

15 23 27 28 28 29 32

16 23 27 28 28 29 32

16 24 27 28 28 29 32

16 24 27 28 28 29 32

16 24 27 28 28 29 32

17 27 31 32 32 33 36

20 28 32 32 32 33 37

20 28 32 32 33 34 37

21 29 33 33 33 34 37

(a) Scenario 1

0 1 2 3 4 5 10

hrs hrs hrs hrs hrs hrs hrs

0%

1%

2%

3%

4%

5%

10%

15%

20%

25%

30%

0 0 0 0 0 0 0

1 3 6 9 15 17 21

1 4 8 11 17 19 24

1 4 8 11 18 21 27

1 5 10 12 19 22 27

1 5 10 12 19 23 28

1 7 12 14 20 24 31

2 9 14 16 22 26 36

3 11 17 21 27 31 37

3 12 17 21 27 31 38

3 12 17 21 28 32 38

(b) Scenario 2

0 1 2 3 4 5 10

hrs hrs hrs hrs hrs hrs hrs

0%

1%

2%

3%

4%

5%

10%

15%

20%

25%

30%

11 20 21 22 23 25 35

14 24 24 24 25 27 37

16 26 26 26 26 29 38

16 27 27 27 27 30 38

18 30 31 31 31 32 38

18 30 31 31 31 32 38

25 36 36 36 36 37 40

30 37 37 37 37 37 40

32 38 38 38 38 38 40

32 38 38 38 38 38 40

33 38 38 38 38 38 40

(c) Scenario 3

0 hrs

0%

1%

2%

0

0

40

(d) Scenario 4 Table 8: The number of runs which produced non-dominated pairs that were within the different thresholds for total tardiness (each row) and expansion cost (each column, expressed as a percentage indicator). A plot of these numbers can be found in Table 9.

43

(a) Scenario 1

(b) Scenario 2

(c) Scenario 3

(d) Scenario 4

Table 9: A plot of the data from Table 8. Each type of symbol represents a different threshold on total tardiness, whilst the horizontal axis represents the different thresholds for expansion cost (as a percentage indicator). The vertical axis shows the number of runs satisfying the different pairs of thresholds.

Scenario

Number of Periods

Number of Variables

Number of Constraints

Minimum Number of Train trips

Maximum Number of Train Trips

1 2 3 4

61 63 71 60

19932 17390 19194 16500

33746 31043 34125 29334

1703 1670 1689 1601

2078 2022 2010 2038

Table 10: Size of the MILP and minimum and maximum number of train trips over all runs for each scenario.

44

Scenario Average Iterations Per Run Average MILP Duration Average Scheduler Duration Average Iteration Duration % of Average Time for MILP % of Average Time for Scheduler

1

2

3

4

66.63 3.46 8.71 12.23 28.30% 71.23%

117.05 3.18 9.45 12.68 25.08% 74.55%

99.75 3.16 4.60 7.82 40.37% 58.85%

2.13 2.36 2.21 4.63 51.07% 47.69%

Table 11: Iterations and durations (minutes) for each scenario.

to the scheduling engine, which is currently implemented as a single-threaded simulated annealing algorithm. Thus, a direction of future research could be the design and implementation of multi-threaded algorithms for the scheduling engine. 8. Conclusion The optimisation procedure presented in this paper addresses the lack of literature on the optimisation aspect of capacity planning in supply chains of mineral resources. The procedure is a hybridisation of a mixed integer linear programming and simulated annealing. The use of the MILP for determining the capacity expansions permits optimisation with a global view of the supply chain and planning horizon. On the other hand, the simulated annealing based scheduling engine analyses the effect of queuing and congestions, which is difficult to model adequately with MILP. The novelty of the presented optimisation procedure is the method of modelling the considered type of supply chains by MILP and the mechanism of hybridisation. This mechanism utilises the identification of bottlenecks and the modification of MILP constraints, complemented by the reduction of search space using the information on the flows of commodity in the MILP solution. The optimisation procedure was tested on data originating from the Hunter Valley Coal Chain. The testing demonstrated the ability of the procedure to solve problems of sizes faced by the industry. Further research may comprise the development and comparison of different mathematical programming models and different metaheuristic based schedulers under the presented hybridisation framework.

45

9. Literature References [1] Bachelet, B., Yon, L., 2007. Model enhancement: improving theoretical optimization with simulation. Simulation Modelling Practice and Theory 15, 703–715. doi:10.1016/j.simpat.2007.02.003. [2] Bai, R., Kendall, G., Qu, R., Atkin, J.A., 2012. Tabu assisted guided local search approaches for freight service network design. Information Sciences 189, 266–281. doi:10.1016/j.ins.2011.11.028. [3] Blum, C., Puchinger, J., Raidl, G.R., Roli, A., 2011. Hybrid metaheuristics in combinatorial optimization: A survey. Applied Soft Computing 11, 4135–4151. doi:10.1016/j.asoc.2011.02.032. [4] Boland, N.L., Savelsbergh, M.W., 2012. Optimizing the hunter valley coal chain, in: Gurnani, H., Mehrotra, A., Ray, S. (Eds.), Supply Chain Disruptions. Springer London, pp. 275–302. doi:10.1007/978-0-85729778-5 10. [5] Boschetti, M.A., Maniezzo, V., Roffilli, M., R¨ohler, A.B., 2009. Matheuristics: Optimization, simulation and control, in: Blesa, M., Blum, C., Di Gaspero, L., Roli, A., Sampels, M., Schaerf, A. (Eds.), Hybrid Metaheuristics. Springer Berlin Heidelberg. volume 5818 of Lecture Notes in Computer Science, pp. 171–177. doi:10.1007/978-3-642-049187 13. [6] Chen, T.L., Chen, Y.Y., Lu, H.C., 2013. A capacity allocation and expansion model for tft-lcd multi-site manufacturing. Journal of Intelligent Manufacturing 24, 847–872. doi:10.1007/s10845-012-0634-9. [7] Crainic, T.G., Laporte, G., 1997. Planning models for freight transportation. European Journal of Operational Research 97, 409–438. doi:10.1016/S0377-2217(96)00298-6. [8] Fu, M.C., 2002. Optimization for simulation: Theory vs. practice. INFORMS Journal on Computing 14, 192–215. doi:10.1287/ijoc.14.3.192.113.

46

[9] Fu, M.C., Glover, F.W., April, J., 2005. Simulation optimization: a review, new developments, and applications, in: Proceedings of the 37th conference on Winter simulation, Winter Simulation Conference. pp. 83– 95. doi:10.1109/WSC.2005.1574242. [10] Geng, N., Jiang, Z., 2009. A review on strategic capacity planning for the semiconductor manufacturing industry. International Journal of Production Research 47, 3639–3655. doi:10.1080/00207540701871051. [11] Georgiadis, M.C., Longinidis, P., 2011. Optimal design and operation of supply chain networks under demand uncertainty, in: Minis, I., Zeimpekis, V., Dounias, G., Ampazis, N. (Eds.), Supply Chain Optimization, Design, and Management: Advances and Intelligent Methods. IGI Global, pp. 73–108. doi:10.4018/9781615206339.ch004. [12] Ghezavati, V., Sadjadi, S., Dehghan Nayeri, M., 2011. Integrating strategic and tactical decisions to robust designing of cellular manufacturing under uncertainty: Fixed suppliers in supply chain. International Journal of Computational Intelligence Systems 4, 837–854. doi:10.1080/18756891.2011.9727835. [13] Hsu, C.I., Li, H.C., 2009. An integrated plant capacity and production planning model for high-tech manufacturing firms with economies of scale. International Journal of Production Economics 118, 486–500. doi:10.1016/j.ijpe.2008.09.015. [14] Javid, A.A., Azad, N., 2010. Incorporating location, routing and inventory decisions in supply chain network design. Transportation Research Part E: Logistics and Transportation Review 46, 582–597. doi:10.1016/j.tre.2009.06.005. [15] Law, A., McComas, M., 2000. Simulation-based optimization, in: Simulation Conference, 2000. Proceedings. Winter, pp. 46–49. doi:10.1109/WSC.2000.899696. [16] L’Ecuyer, P., 1999. Good parameters and implementations for combined multiple recursive random number generators. Operations Research 47, 159–164. doi:10.1287/opre.47.1.159.

47

[17] L’Ecuyer, P., 2012. Random number generation, in: Gentle, J.E., H¨ardle, W.K., Mori, Y. (Eds.), Handbook of Computational Statistics. Springer Berlin Heidelberg. Springer Handbooks of Computational Statistics, pp. 35–71. doi:10.1007/978-3-642-21551-3 3. [18] L’Ecuyer, P., Meliani, L., Vaucher, J., 2002a. SSJ: a framework for stochastic simulation in Java, in: Simulation Conference, 2002. Proceedings of the Winter, pp. 234–242 vol.1. doi:10.1109/WSC.2002.1172890. [19] L’Ecuyer, P., Simard, R., Chen, E.J., Kelton, W.D., 2002b. An object-oriented random-number package with many long streams and substreams. Operations Research 50, 1073–1075. doi:10.1287/opre.50.6.1073.358. [20] Liu, H., Abraham, A., Sn´aˇsel, V., McLoone, S., 2012. Swarm scheduling approaches for work-flow applications with security constraints in distributed data-intensive computing environments. Information Sciences 192, 228–243. doi:10.1016/j.ins.2011.12.032. [21] Mahfoud, S.W., Goldberg, D.E., 1995. Parallel recombinative simulated annealing: A genetic algorithm. Parallel Computing 21, 1 – 28. doi:10.1016/0167-8191(94)00071-H. [22] Maniezzo, V., St¨ utzle, T., Voß, S., 2009. Matheuristics: Hybridizing metaheuristics and mathematical programming. volume 10. Springer. doi:10.1007/978-1-4419-1306-7. [23] Nissen, V., 1994. Solving the quadratic assignment problem with clues from nature. Neural Networks, IEEE Transactions on 5, 66–72. doi:10.1109/72.265961. ´ [24] Olafsson, S., Kim, J., 2002. Simulation optimization: Simulation optimization, in: Proceedings of the 34th Conference on Winter Simulation: Exploring New Frontiers, Winter Simulation Conference. pp. 79– 84. URL: http://dl.acm.org/citation.cfm?id=1030453.1030467. [25] Paul, R.J., Chanev, T.S., 1998. Simulation optimisation using a genetic algorithm. Simulation Practice and Theory 6, 601–611. doi:10.1016/S0928-4869(98)00007-X.

48

[26] Pedersen, M.B., Crainic, T.G., Madsen, O.B., 2009. Models and tabu search metaheuristics for service network design with asset-balance requirements. Transportation Science 43, 158–177. doi:10.1287/trsc.1080.0234. [27] Pinedo, M., 2012. Scheduling theory, algorithms, and systems. Springer, New York. doi:10.1007/978-1-4614-2361-4. [28] Renna, P., 2013. Capacity investment decision by monte carlo approach in a cooperation network. International Journal of Production Research 51, 6455–6469. doi:10.1080/00207543.2013.824628. [29] Sandeman, T., Fricke, C., Bodon, P., Stanford, C., 2010. Integrating optimization and simulation-a comparison of two case studies in mine planning, in: Simulation Conference (WSC), Proceedings of the 2010 Winter, IEEE. pp. 1898–1910. doi:10.1109/WSC.2010.5678882. [30] Shapiro, J., 2001. Modeling the supply chain. Brooks/Cole-Thomson Learning, Pacific Grove, CA. [31] Singh, G., Sier, D., Ernst, A.T., Gavriliouk, O., Oyston, R., Giles, T., Welgama, P., 2012. A mixed integer programming model for long term capacity expansion planning: A case study from the hunter valley coal chain. European Journal of Operational Research 220, 210–224. doi:10.1016/j.ejor.2012.01.012. [32] Solberg, J.J., 1981. Capacity planning with a stochastic workflow model. AIIE Transactions 13, 116–122. doi:10.1080/05695558108974543. [33] Swaminathan, J.M., 2000. Tool capacity planning for semiconductor fabrication facilities under demand uncertainty. European Journal of Operational Research 120, 545–558. doi:10.1016/S0377-2217(98)003890. [34] Talbi, E.G., 2013. A unified taxonomy of hybrid metaheuristics with mathematical programming, constraint programming and machine learning, in: Talbi, E.G. (Ed.), Hybrid Metaheuristics. Springer Berlin Heidelberg. volume 434 of Studies in Computational Intelligence, pp. 3–76. doi:10.1007/978-3-642-30671-6 1.

49

[35] Ting, C.K., Li, S.T., Lee, C., 2003. On the harmonious mating strategy through tabu search. Information Sciences 156, 189–214. doi:10.1016/S0020-0255(03)00176-2. [36] Welgama, P., Oyston, R., 2003. Study of options to increase the throughput of the hunter valley coal chain, in: Proceedings of MODSIM, pp. 1841–1846.

50