A resource constrained scheduling problem with multiple independent producers and a single linking constraint: A coal supply chain example

European Journal of Operational Research 236 (2014) 946–956

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

A resource constrained scheduling problem with multiple independent producers and a single linking constraint: A coal supply chain example Anu Thomas a,⇑, Jayendran Venkateswaran b, Gaurav Singh c, Mohan Krishnamoorthy a a

IITB-Monash Research Academy, IIT Bombay, Powai, Mumbai 400076, India Industrial Engineering and Operations Research, IIT Bombay, Mumbai 400076, India c CSIRO Computational Informatics, Clayton, VIC 3168, Australia b

a r t i c l e

i n f o

Article history: Available online 12 October 2013 Keywords: Distributed decision making Coordinated scheduling Resource constrained scheduling Column generation Job scheduling

a b s t r a c t This paper examines a resource constrained production planning and scheduling problem motivated by the coal supply chain. In this problem, multiple independent producers are connected with a resource availability (or, linking) constraint. A general description of such problems is provided, before decomposing the problem into two levels. In the first level, we deal with production planning and in the second level, we deal with tactical resource scheduling. A real-world coal supply chain example is presented to anchor the approach. The overall problem can be formulated as an integrated mixed integer programming model which, in several cases, struggles to find even a feasible solution in reasonable amount of time. This paper discusses a distributed decision making approach based on column generation (CG). Computational experiments show that, the CG scheme has significant advantages over the integrated model and a Lagrangian relaxation scheme proposed by Thomas et al. (2013). This paper concludes with detailed discussions on the results and future research directions. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction In this paper, we consider a resource constrained planning and scheduling problem in a supply chain in which many independent producers (or service providers) coordinate with a resource manager who provides resources such as trucks, trains or infrastructure, to transport products to their destinations. In such supply chains, producers independently make planning decisions and the resource manager (an independent decision maker) makes scheduling decisions that collectively affect the efficiency of the supply chain. The supply chain’s objective is to ensure a feasible allocation of resources that minimises total system cost. There are several real-life supply chains of this sort. In coal supply chains, multiple independent mines (producers) are often connected to a shipping terminal via a shared rail network. A single rail operator (resource manager) often coordinates rail schedules with the mines, to schedule trains (or, trips) between the mines and the shipping terminal in order to meet shipping orders (or, demand) of each mine (see, for example, Singh et al., 2012 and Thomas et al., 2013). A similar problem is also observed in wine supply chains (see Singh et al., 2005) in which, each grower has contracts with many wineries. Each winery needs certain varieties ⇑ Corresponding author. Tel.: +91 22 2576 7533. E-mail addresses: [email protected] (A. Thomas), [email protected] (J. Venkateswaran), [email protected] (G. Singh), [email protected] (M. Krishnamoorthy). 0377-2217/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ejor.2013.10.006

of grapes based on their production requirements. These need to be delivered to them in a preferred sequence. The winery needs to coordinate harvesting and transport of grapes to minimise wastage and waiting time (reducing perishability) and to maximise resource utilisation (of crushers, fermenters). See Singh et al. (2005) for a detailed analysis of coordinated scheduling problems wine supply chains in Australia. In a general problem, each independent producer receives a set of orders from their customers. These orders can be satisfied with different combinations of shared/common resources. The resources are often be grouped into classes based on their properties. The resource manager has a finite number of resources in each class. Therefore, there will be a conflict if the producers’ requests exceed the availability of these resources. At the production planning level, producers are not concerned about resource availability in each class. For producers, fulfilment of orders of all customers is of far greater importance. An optimal allocation and scheduling of resources is subsequently necessary to improve the overall performance, which can be achieved via either a cost minimisation, profit maximisation or a makespan or tardiness objective. It is possible to formulate (and maybe solve) the above supply chain problem in an integrated manner. In such a model, we could consider a single decision-making problem that includes all producers and the resource manager. Such an integrated model is likely to be large and complex. Most scheduling problems are NP-hard (see Brucker and Knust, 2006; Lenstra et al., 1977). The

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956

integrated scheduling problem adds to this complexity because the producers and the resource manager are independent decision making entities who are often reluctant to coordinate fully and transparently, to share their competitive information. In this paper, we develop an approach, using a customised column generation (CG) scheme, to solve resource allocation and resource scheduling problems. In the production planning stage, each producer defines a possible allocation of resources. This means that resource allocation and its scheduling can be carried out sequentially. Thus, the original problem is decomposed into smaller solvable sub-problems. Essentially, we develop a distributed decision making framework which is solved using a decomposition approach. Even though we explain the proposed solution approaches in the context of a coal supply chain, our method can be extended to other similar supply chains too. Section 2 provides a summary of related work. Section 3 describes the problem and its features. Section 4 explains the production planning problem for each producer. In Sections 5 and 6 we present a CG algorithm and its improvements, respectively. Section 7 gives an overview of the datasets and performance measures. Section 8 summarises detailed results and comparisons of the different modelling approaches. 2. Related work Integrated models are often solved using decomposition techniques such as Lagrangian relaxation, column generation and Benders decomposition. The vast literature available in this topic highlights its advantages and practical applications. For example, decomposition approaches are able to handle the complexity of planning and scheduling problems (see Wang et al., 2008). Given the inherent complexity of distributed decision making problems, decomposition methods can be employed for solving these problems. Bracken and McGill (1973) introduced hierarchical optimisation as a generalisation of mathematical programming, where a series of optimisation problems are solved in a predetermined sequence. In this, independent decision makers are considered separately, thereby reducing the problem size and complexity of the integrated model. Decomposed optimisation techniques have been extended to job-shop scheduling problems (see Brandimarte and Calderini, 1995; Zribi et al., 2007), resource constrained scheduling problems (RCSP) (see Hans et al., 2007; Kelly and Zyngier, 2008), manufacturing problems (see Wu and Ierapetritou, 2007; Ebadian et al., 2009; Timm and Blecken, 2011). Pochet and Wolsey (2006) discuss different production planning models and mixed integer programming formulations. They suggest many practical reformulation techniques and polyhedral results to strengthen the model and solve it effectively. Dantzig and Wolfe (1960) introduced an efficient decomposition algorithm for a special class of problems. Today, Dantzig– Wolfe decomposition (DWD) is a widely applied solution approach for solving large optimisation problems. Column generation (CG) methods have been incorporated along with the DWD to successfully solve many problems such as the cutting-stock problem, vehicle routing problem, crew scheduling problem, p-median problem, and graph colouring problem (see Borndorfer et al., 2005; Lubbecke and Desrosiers, 2005; Vanderbeck and Wolsey, 2010). Read Desrosiers and Lubbecke (2005) for an introduction to CG and its key features. Choi and Tcha (2007) proposed a CG approach to address a heterogeneous-fleet vehicle routing problem. Capone et al. (2010) compare two variants of CG schemes developed for a resource allocation problem in wireless mesh networks. All of these show that the CG scheme is effective for solving large problems. Wentges (1997) introduced a weighted DWD and presented some results on the convergence of this scheme under certain conditions. Pessoa et al. (2010) discuss a stronger convergence of

947

this algorithm. The major issues concerned with DWD are slow convergence and stabilisation of dual prices. Ben Amor et al. (2009) and Vanderbeck and Wolsey (2010) analyse some of the latest (dual price) stabilisation techniques. Thomas et al. (2013) propose a decomposition algorithm based on the Lagrangian relaxation (LR) to solve a similar problem. The main contributions of Thomas et al. (2013) are: (i) the mathematical models for production planning; (ii) distributed algorithm based on the LR and improvements with the Volume algorithm and the Wedelin algorithm; (iii) heuristics to compute the upper bound; and (iv) extensive computational experiments. The authors highlight the advantages of a distributed decision making approach over an integrated approach. Even though the LR algorithm performed better than the integrated model, it was time consuming for larger problems. Hence CG algorithms are explored to solve the integrated problem. A combination of multiple stabilisation methods are explored to improve the performance of the CG algorithm. Rather than directly applying the CG algorithm, we discuss and customise a few techniques that improve the computational efficiency. 3. Problem description Consider a model where a set of independent producers are connected with a single shared resource constraint. The constraints of each producer have some influence on the utilisation of the shared resource, which, in reality, is the independent decision of the resource manager. Since the resources are finite, these independent producers have to coordinate (implicitly) through the resource manager to achieve a better decision outcome for themselves. There are multiple ways to model and solve this problem. The first option is to model and solve this as a single large mixed integer programming problem. This model has to include all the constraints from the producers and the resource manager. Since this model integrates all the independent decision makers, we refer to it as the integrated model. Practical instances of such problems tend to be very large-sized. The integrated model has many solvable production planning sub-problems, and one resource scheduling problem linking the single resource manager and the (many) producers. Thus, we decompose the problem into two parts:  Production planning: Each producer plans its production based on its priorities and objective, and places a set of requests to the resource manager for a certain number of resources. In other words, each producer defines a set of requests (jobs) to meet its orders.  Resource scheduling: After receiving the requests (jobs) from the producers, the resource manager prepares a schedule based on resource availability. This problem is equivalent to a job scheduling problem. As a result of this decomposition, decisions are separated and distributed across the many independent partners in the system. Therefore, in such distributed models, it is not necessary to share too much information between the independent partners; while this is, of course, necessary in the integrated model. To understand the formulation and solution approaches better, we explain the problem with the help of an example from a coal supply chain. 3.1. A coal supply chain example In a coal supply chain, many mines produce/mine coal based on the market demand and use a single rail operator to transport the coal to a shipping terminal. The coal supply chain has multiple

948

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956

independent mines (or, a collection of mines owned by different owners). For the sake of simplicity, in this paper, we assume that each mine belongs to an owner. Coal is transported from the mines (producers) to the terminal, mainly using the rail network. Therefore the common rail operator (resource manager) has to coordinate the independent mines to avoid conflicts in train scheduling. In this paper, we ignore all terminal operations. Each mine receives a set of orders from their customers (or the terminal). Each order has a due-date and an order quantity associated with it. Each mine incurs an inventory cost at the mines, a stocking cost if the coal reaches the terminal before the due-date, and a demurrage cost for any late deliveries. A single rail operator provides the trains to move the coal from the mine to the terminal. Also, multiple train trips maybe required to fulfil any one order. Each of train trip is considered as a job. The orders at the mines need not have a one-to-one correspondence with jobs. Usually, more than one job is needed to complete an order. It is also possible that one job may satisfy two separate orders. In other words, mines do not know, a priori, which jobs are required for them; this is a decision that needs to be made. Each mine has to schedule production (subject to many constraints) and meet its orders at minimum cost. It is assumed that all trains are fully loaded at the mines. Partial loading or mixing loads from different mines are not permitted due to operational reasons. For example, in Australia, mine sites can be almost 400–500 kilometers1 away from the ports. Distances between mines and terminals are vast, coal trains are some 2–4 kilometers long, and the transport of material requires excessive labour and fuel costs. Therefore it does not makes sense to have trains only ‘partially-loaded’. We assume that at most only one order can be tardy at the mines. There are two due-dates in this problem: (i) the due-date of an order associated with a mine, and (ii) the due-date of a job, which is an internal due-date, between the mine and the rail operator. Similarly, tardiness is also defined with respect to the job as well as the order. The due-date/tardiness of a mine and the rail operator are computed with respect to the order and the job respectively. The rail operator has a pool of trains of different classes, classified based on the rake capacity and other features. Individual mines typically plan for and request trains of a particular class at appropriate times (so as to minimise their operational costs), irrespective of the demands from other mines. A mine is not necessarily concerned about the specific train that is allocated to it. The scheduling of specific trains that satisfy a set of jobs is a separate and non-trivial decision that needs to be made by the rail operator. The challenge here is that practical instances of such problems involve two sets of decisions: job-selection decisions and job-schedule decisions. In the job selection-scheduling problem, multiple mines could make simultaneous requests for the same class of trains thereby making the problem difficult. This can be viewed as a multi-resource constrained scheduling problem. The objective and constraints of this job scheduling problem are elaborated in Section 4.

3.2. Motivation Solution approaches that are developed for resource scheduling in a coal supply chain could be extended to other industries (airline, automobile manufacturing and the service industries). Hence it is important to study resource scheduling in a coal supply chain. The scheduling problem with a single resource class itself is a complicated problem. It does not have analytical solutions for most of its variants (see Lenstra et al., 1977). In this paper, we consider a multi-resource constrained scheduling problem with the addi1

http://www.hvccc.com.au.

tional complexity of defining jobs based on the production planning decisions. The dependency between the orders and jobs makes this problem challenging. It cannot be modelled as a batch scheduling or a makespan minimisation problem. All of these relations can be captured in a single integrated model. But that will be too large and complex and would make it impossible to solve a single integrated model for large problem instances within a reasonable amount of computational time. Hence the need for a computationally viable distributed decision making approach that is based on decomposition methods. There are two issues to consider while developing such an approach: (i) the right algorithm, and (ii) the right customisation. We develop a CG algorithm that significantly outperforms another popular decomposition algorithm based on Lagrangian relaxation. Multiple stabilisation techniques are implemented to overcome the weakness of the CG algorithm. A decentralised decision making approach is an alternative. In such models we have (i) multiple partners, (ii) information asymmetry, (iii) conflict in objectives, and (iv) a negotiation protocol. We model a situation with multiple mines that do not interact/ share-information with other mines but, constantly interact with the rail operator. The mines and the rail operator possess different sets of information. Hence our decision model has many players with an inherent information asymmetry and meets the first two requirements. The objectives of (each of) the mines and the rail operator are not same. Each mine expects the rail operator to send a train to it to load the coal and deliver it on-time at the terminal. However, the rail operator may have different priorities and constraints. Hence there is a conflict in objectives. We do not, however, incorporate the fourth requirement of negotiation protocols. In that sense, ours is not a truly decentralised decision making framework. A negotiation protocol is only partially implemented using the CG algorithm which is (in a sense) ‘controlled’ by the ‘honest broker’ rail operator. We made this assumption just so that we could converge to a solution to a very difficult real-world problem. The core objective of this paper is to analyse the interaction between the production planning and resource scheduling. The experimental data sets assume that the base time unit is 1 hour and the length of the planning horizon is either one week or one fortnight. This framework can be extended to include other complexities such as rolling horizon, different time units and longer time horizons. 4. Production planning model The decomposition scheme allows us to solve the production planning problems separately. The production planning sub-problem described in this section has to be solved in every iteration of the CG algorithm. We may use the nouns mine and producer, train and resource or rail operator and resource manger interchangeably based on the context. Let i 2 {1, 2, . . . , I}, t 2 {0, 1, . . . , T}, u 2 {1, 2, . . . , U} and w 2 {1, 2, . . . , W} be the indexes of mines, time periods, orders and train (resource) classes respectively. Let SPi denote the production planning sub-problem of ith mine. The solution space corresponding to SPi is denoted as S i . Parameters Atw cost incurred for requesting a train of class w at time t

Ct Ot Ht

tardiness/demurrage cost of an order at time t overstock cost per unit quantity at the terminal at time t inventory holding cost per unit quantity at the mine at time t (continued on next page)

949

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956

Fu Qu Dt P B Vw Kw Sw Lw Rw Ew

due date of order u uth order quantity cumulative demand at time t production capacity at the mine inventory holding capacity at the mine capacity of the trains in class w number of trains in class w forward travel time required for the trains in class w to reach a mine from the terminal loading time required for the trains in class w return travel time required for the trains in class w, from a mine to the terminal = Lw + Rw, time required for loading and return travel to the terminal

Decision variables ht Inventory level at the mine by the end of time t rt Amount produced in period t xt Over stock level of the product at the terminal by the end of time t wt Binary variable, equals to 1 if any order is tardy at time t. Zero indicates no order is tardy gtw Total number of trains from class w requested on or before time t

X

w gtE V w  D t 6 xt 8 t w

(u  1)th order must be delivered before the due-date of the uth order, Fu. This also implies that at any time t one order for a mine can be at the most tardy. In other words, cumulative supply at time Fu should be not less than the demand at time Fu1.

X

gFwu Ew V w P DFu1 8 u

min Z i ¼

" X

t

t

t

t

t

t

h H þx O þw C þ

t

X

t w

t1 w

g g



# Atw

ð1Þ

w

Constraints: Eqs. (2)–(12) completely define the constraints of the model. Total number of trains requested on or before time t is non-decreasing.

gtw P gt1 8 w; t w

ð2Þ

The amount produced at the mine at any time t is not more than its production capacity.

rt 6 P 8 t

ð3Þ

ð7Þ

w

If the cumulative supply at time t is less than its demand then the mine needs to pay demurrage cost.

wt P 1  1=Dt 

X

w gtE V w 8 t P F1 w

ð8Þ

w

All orders must be met by the end of the planning horizon.

X

gTw V w P DT

ð9Þ

w

Only one train can load at the mine at a given time period.

X



w gtw  gtL 61 8t w

ð10Þ

w

The total number of trains of a class w running at any time t cannot exceed the total number of trains in that class.

 A delivery order u from the terminal can be seen as an ordered pair of due date (Fu) and quantity (Qu). Without loss of generality we assume that Fu < Fu+1 and uth order must be met before satisfying the demand for (u + 1)th order. The cumulative demand for a P mine at time t is defined as Dt ¼ fujF u 6tg Q u . Objective: The overall objective, Zi, is to minimise total system cost. This includes the cost of holding inventory at the mines, the overstock cost at the terminal, demurrage cost and total cost of requesting trains.

ð6Þ

w



w w gtþS  gtE 6 K w 8 w; t w w

ð11Þ

Boundary conditions and scope.

g0w ¼ r0 ¼ h0 ¼ x0 ¼ wTu ¼ 0; ht ; xt ; rt P 0; wtu 2 f0; 1g; gtw 2 Zþ :

ð12Þ

Additional Constraints: The above constraints are necessary to represent the model. A few more additional integer cuts and bounds have been generated and used to tighten the Mixed Integer Programming (MIP) formulation. These additional constraints are redundant in the MIP formulation, but, active in the linear programming relaxation of the MIP. Hence, an MIP solver should take less time to solve a problem that contains these additional constraints. The cumulative supply at time t cannot be more than the total P production, w gtw V w 6 tP. Due to holding costs at the mine and terminal, the mine does not want to supply more than DT. Since (i) partial loading and order mixing are not permitted and (ii) DT need not be a linear combination of train sizes, we need to allow P t T a margin, Vmax, in the final supply. Hence w gw V w 6 D þ V max , where Vmax = maxw{Vw}. From these two, we define an upper bound Xt for the cumulative supply.

X

gtw V w 6 X t 8 t; where X t ¼ minftP; DT þ V max g

ð13Þ

w

Since the supply is bounded, overstock also has to be bounded. The inventory balancing constraint: Inventory at time t can be computed by subtracting the supply at time t from the sum of previous inventory and production at time t.

ht ¼ ht1 þ rt 

X

 gtw  gt1 Vw w

8 w; t

ð4Þ

w

The inventory at end of period t cannot be more than the holding capacity of the mine.

ht 6 B 8 t

ð5Þ

If a train of class w has to reach the terminal at time t, then it should reach the mine at t  Ew. Hence the cumulative supply at P w time t can be expressed as w gtE V w . If it is more than the mine’s w demand, then it has an overstock.

xt 6 X t  Dt 8 t

ð14Þ

If a mine overstocks at the terminal then it cannot incur demurrage.

wt 6 1  xt =ðX t  Dt Þ 8 t

ð15Þ

Constraint (13) is similar to a Knapsack constraint and therefore can be further tightened with the integer cuts proposed by Atamtürk (2005). Let Kw be the upper bound of gtw which satisfies (13) after fixing all other variables to zero. The set C  {1, 2, . . . , N} P t is a ‘cover’ if w2C ðKw V w  X Þ P 0. Then a stronger cut can be derived as

X

 Kw  gtw P dX t =Ve 8 t

w2C

where V ¼ maxfw2Cg V w and Kw = bXt/Vwc.

ð16Þ

950

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956

Similarly we can use an integer cut from constraint (8). Since both w and g are integer variables, this cut will help us to improve the relaxed solution.

ð1  wt ÞdDt =V min e 6

V MP ¼ min

X V ic xic

s:t:

X

X xic ¼ 1 8 i

where D ¼

P

fujF u 6tg Q u

vticw xic 6 K w 8 w; t

8 u;

and Vmin = minw{Vw}.

F u < t < F uþ1 :

xic 2 f0; 1g

ð18Þ

The constraints (1)–(18) are applicable for individual mines. Hence, the decision variables and the parameters corresponding to mine i are used with a subscript of i. For example P i ; Hi ; gtiw , are used in the context of the integrated problem. Since the train is utilised across the mines. Hence the constraint (11) can be further tightened as

X



w w gtþS  gtE 6 K w 8 w; t iw iw

ð22Þ

i;c

The demand changes only at the due-date of orders. Once the demand for an order was met in an interval, it stays as ‘met’ till the due-date of the next order. This gives

wt 6 wt1

ð21Þ

c X

w gtE dV w =V min e 8 u; F u 6 t < F uþ1 ð17Þ w

w t

ð20Þ

i;c

ð19Þ

i

ð23Þ

The objective (20) minimises the total value of selected columns. Constraint (21) ensures that exactly one column will be selected for each mine. All the selected columns must satisfy constraint (22). Each of the columns satisfies the corresponding mine’s constraints. Therefore the optimal solution of the master problem will be feasible for all mines and optimal for the system. In the relaxed master problem (RMP), the integer constraint (23) is relaxed as

0 6 xic 6 1:

ð24Þ

At an iteration k, the number of columns will be finite and restricted to a subset of all the feasible columns. The restricted relaxed master problem at iteration k is represented as RMPk. 5.2. An iterative scheme based on column generation

5. A column generation approach Thomas et al. (2013) provide a Lagrangian relaxation (LR) approach for solving integrated decision problems. But, as discussed in Section 2, this has a few shortcomings. The LR approach does not solve large problems quickly. Moreover, the Lagrangian multipliers in any iteration relies only on the information that is available from the sub-problems at that iteration. This results in a slow convergence of the Lagrangian multipliers. Hence, there is a need for a more efficient approach to solve practical instances of such problems. A CG algorithm has two components: master problem (MP) and sub-problems. We consider the train scheduling problem without any production considerations, as the MP. The production planning problem, described in Section 4, is considered in the sub-problems. Solutions (schedules) of each mine form the columns of the MP, from which a globally feasible solution is obtained. However, this MP can be very large. Therefore, we consider a restricted MP which has a smaller subset of the columns. The sub-problems iteratively generate new columns for the MP using the dual prices from the restricted MP. The different solutions to the MP imply different allocations of trains to the mines. 5.1. Master problem The master problem selects the best column for each mine which does not violate the resource availability constraint and minimises the total system cost. Index c denotes the columns. Define, Vic

vticw xic

Value (cost) of column c of ith mine. For each column c, there exists a solution in S i . Hence Vic can be computed with the objective function (1) Number of servicing trains of class w for the ith mine with respect to column c at time t 1 if the column c of mine i is used, 0 otherwise

The parameters Vic an vticw are precomputed for each column. xic is the only decision variable of the MP. Then the master problem is

Column generation algorithm solves the relaxed master problem iteratively. Initially it starts with a finite set of feasible columns and iteratively adds more columns. The penalty for violating the resource constraint is communicated to the sub-problems through the dual prices computed from the master problem. Let ktw be the dual price of resource constraint for class w at time t, then the sub-problems are solved with the modified LR objective function:

Z i ðkÞ ¼ min

X

hti Hti þ xti Oti þ wti C ti þ

t

where Ati;w ¼ Ati;w þ

X



!

t gti;w  gt1 i;w Ai;w

ð25Þ

w tþE w 1 X

kuw :

u¼tSw

Hence LBðkÞ ¼

the lower bound at kth  P P  ðkÞ ðkÞ . i Z i ðk Þ  w K w  kw

iteration

is

defined

as

The objective of the CG algorithm is to compute a better lower bound for the relaxed master problem. The lower bound from the integrated model and the CG will be the same if the integer relaxation of the problem was being solved, but in this case it is an MILP problem and therefore the quality of bound also depends on the cuts and how deep the branch and bound search is able to reach. We show in this paper through our decomposition techniques that our method is much better in achieving this. In every iteration k, LBk is less than or equal to V RMP . Once both values are equal, sub-problems will not add any new columns to the solution pool. Therefore the algorithm will not guarantee good integer solutions. Section 5.3 illustrates heuristic procedures to achieve a better upper bound by generating more numbers of globally feasible columns from the solutions of RMPk. 5.3. Heuristics to generate feasible columns In this section, we propose two heuristics procedures: (i) Job scheduling model, where the rail operator attempts to derive a feasible solution to the problem with inputs from all the mines, (ii) Leader–follower mechanism, where each mine attempts to make a better schedule after fixing certain columns (partial schedules) for all the other mines. In these heuristics, importance is given to

951

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956

global feasibility rather than cost minimisation. Heuristic 1 and Heuristic 2 receive inputs from the solutions of RMP, which consists of many partial schedules for each mine. Since these heuristics are customised to generate global feasible columns, it improves the upper bound UB⁄. 5.3.1. Job scheduling A combination of the solutions from the sub-problems need not satisfy the resource constraint (22) of the master problem. Therefore the rail operator uses the available information from the subproblem solution and attempts to prepare a globally feasible schedule which minimises deviation from expected schedules. The rail operator, upon receiving the requests from the mines, is aware of the train class combination and due-dates. Hence it defines the train trips (jobs) required by each class at this stage. Each job has three stages: (i) forward travel from the terminal to the mine; (ii) loading; and (iii) return journey with full load to the terminal. We define additional variables required for the job scheduling model: Parameters j index of job Mj mine associated with job j Tj train class associated with job j Sj forward travel time required for the job j Lj loading time required for the job j Rj return travel time required for the job j Ej = Lj + Rj, time required for loading and return travel for the job j Fj due date of job j T W j tardiness weights of the job j, which reflect the importance of the job. Higher weight implies higher importance W Ej earliness weights of the job j C tj

weighted tardiness of job j at time

Otj

t; C tj ¼ W Tj maxft  F j ; 0g for all j, t weighted earliness of job j at time

Kw

t; ¼ maxfF j  t; 0g for all j, t number of trains with class w

Otj

j

tþEj

zj

V j 6 tP i

8 i; t

ð27Þ

jjM j ¼i

At any time, there should not be more than one job for loading at any production location.

 X  tþEj tþR zj  zj j 6 1 8 i; t

ð28Þ

jjM j ¼i

Once a job is completed it stays as completed.

ztj P zt1 j

8 j; t

ð29Þ

All jobs must be completed at the end of the planning horizon.

zTj ¼ 1 8 j 2 J

ð30Þ

The number of trains in any class is limited. This is equivalent to constraints (19) and (22).

 X  tþSj þEj zj  ztj 6 K w

8 w; t

ð31Þ

jjT j ¼w

Every mine must satisfy (u  1)th order before the due-date of the uth order.

X

F

zj i;u V j P DF i;u1

8 i; u

ð32Þ

jjM j ¼i

Scope and boundary conditions.

z0j ¼ 0;

ztj 2 f0; 1g 8 j; t

ð33Þ

In general, this model will give us a feasible solution by advancing the jobs, unless the due-dates are very close and hence constraint (32) is violated. Once the job schedules are computed, the corresponding columns for each mine and their values can be computed easily.



Each job j can be represented as a pair of the mine, Mj, and the corresponding train class, Tj, that is job(j) = (Mj, Tj). The properties such as Gj, Sj, Lj, Rj are inherited from the corresponding train class. In other words, Gj ¼ GT j and so on. Suppose a mine makes a request 0 t 0 1 for a train from the class w at time t0 , then gtw  gw ¼ 1. Therefore the corresponding job is expected to arrive at time t0 . Hence its due date is Fj = t0 + Ej. The weights W Ej ; W Tj and costs C tj ; Otj are decided by the rail operator. In this section, the properties such as due-date, earliness, tardiness and weight are discussed only with respect to the jobs. These properties differ from the similar ones of the orders received from the terminal by the mines. Objective: The objective of this job scheduling model is to minimise the total weighted tardiness and earliness. This guides the model to schedule the jobs as close to their due-dates as possible. If the earliness cost is greater than the tardiness cost, jobs will be advanced only if it is required to maintain feasibility.

  i XXh t  þ Otj ztj  zt1 C j ztj  zt1 j j

X

5.3.2. Leader–follower mechanism As stated earlier, the resource constraint (22) is the only link between the mines. If we know a complete (or partial) schedule of all the other mines, then each mine can have a better understanding about the resource availability. Hence a particular mine can customise their schedule to obtain a globally feasible solution. In this case, constraint (11) can be modified as

W Ej

Decision Variables  ztj 1 if the job j is completed on or before time t; = 0 otherwise:

min Z R ¼

Cumulative supply from any mine should be less than its total possible production.

ð26Þ

t>0

Constraints: Eqs. (27)–(33) denote the constraints of the model.



w w gtþS  gtE 6 K tw 8 w; t w w

ð34Þ

K tw

where ¼ K w  ðTrains reserved by other mines at time tÞ. Solving the sub-problem with the above constraint will restrict the solution space and will increase the possibility of getting globally feasible columns. To achieve this, columns are randomly selected from the solution pool of RMPk for a few mines and solve the job scheduling for others. Hence it will be restricted in a small neighbourhood of feasible solutions of the RMP. Adding these extra columns back to the pool will improve the UB⁄. 5.4. Integer solutions We may get good integer solutions directly from the heuristics described in Section 5.3. Another option is to solve the integer version of the MP where the model will try to find optimal columns from the available pool. There is no guarantee that the optimal integer solution will be one among the optimal solutions of the RMP; some non-basic columns of RMP can also be in the optimal integer solution. To maintain a reasonable problem size, the CG algorithm removes a few non-contributing columns from the solution pool. Hence, there are chances of missing some good candidates for the optimal integer solution. To avoid this problem, we can have two versions of the solution pool S; one for the RMP and one for the

952

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956

integer MP. Unless there are memory restrictions, all columns from every iteration can be stored. But solving the integer MP with all columns is not advisable because of the size of the solution pool. Once the prices are stable, new columns will not be added to the solution pool. Then, we can filter the columns for the integer MP based on the current bounds such that they are within a neighbourhood of the relaxed solution. The basic column generation algorithm was implemented and tested on several randomly generated datasets (see Section 7). Even though CG showed significant improvement over the centralised model, there are multiple strengthening methods that can and must be explored. 6. Strengthening techniques for the CG algorithm The traditional CG algorithm can be improved through stabilisation and improvement techniques.

kÞ is added to the objective of the RMP as a conjugate of Dðk   stability term (see Ben Amor et al., 2009). Therefore the RMP can be updated as

V RMP ¼ min

X V ic xic  k  y þ D ðyÞ

ð36Þ

c

X

vticw xic  ytw 6 K w 8 w; t ðdual variables kÞ ð37Þ

i;c

0 6 xic 6 1; 0 6 ytw

ð38Þ

  where y ¼ ytw corresponds to the additional train provided for the class w at time t and the stability term is defined as

8 if 0 6 m 6 l; > < mM 1 D ðyÞ ¼ gðytw Þ; where gðmÞ ¼ lM 1 þ M2 ðm  lÞ if m > l; > : w;t 0 otherwise: X

ð39Þ

6.1. Stabilisation The major drawbacks of the basic CG scheme, as mentioned by Ben Amor et al. (2009) are (i) slow convergence and (ii) high oscillation of dual prices. Many stabilisation methods are proposed in the literature to overcome these problems. The reader can refer to Lubbecke and Desrosiers (2005), Ben Amor et al. (2009) and Vanderbeck and Wolsey (2010) for recent reviews. Stabilisation of dual price is centred around what is referred to as stability centre. In any iteration, the stability centre is the current best known dual price. Hence, we try to get dual prices close to the current stability centre. The concept of stability centre will minimise very high deviations in the dual prices and will eventually improve the quality of the solution. In our proposed scheme, we use the weighted Dantzig–Wolfe decomposition method discussed in Desrosiers and Lubbecke (2005) and Pessoa et al. (2010) along with a twopiece stability term motivated from Ben Amor et al. (2009). The weighted scheme stabilises the dual prices at the sub-problem level. The stability term added to the objective helps to stabilise the master problem. By introducing both methods we aim to minimise the instability of the scheme and achieve the optimal state in fewer number of iterations. k be the stability centre. To stabilise the dual price, we Let  penalise the difference between new dual prices k and  k. Fenchel’s

ð35Þ

i;c

X s:t: xic ¼ 1 8 i

and M1, M2 are constants such that M1 < M2. 6.2. Multiple columns In general, finding all possible solutions of each sub-problem is hard. The traditional CG algorithm considers only one optimal solution from the sub-problems. However, Huisman et al. (2005) suggest that faster convergence could be achieved by adding more columns. The complexity of MP and the additional time and computational resource required are the bottlenecks in using this option directly. A commercial solver, CPLEX 12.1, is used to perform the computational experiments. Instead of storing a single solution from each iteration, multiple solutions are recorded. CPLEX has the ability to generate more solutions, but at a higher cost and memory for each of the sub-problems and the MP. Hence the advantage of iterative improvement is lost. Therefore all computational experiments were performed with default settings. 6.3. Improved CG algorithm The basic CG algorithm is strengthened with the stabilisation techniques and the improvement features discussed in Section 6.2. An improved version of it is presented in Algorithm 1.

Algorithm 1. Improved CG algorithm 1: Initialise  k ¼ 0; LB ¼ 0; UB ¼ 1; k ¼ 0; a ¼ 0:5; Time Limit TL ¼ 3600;  ¼ 105 ; l ¼ 0:5; eTime ¼ elapsed time in seconds. 2: Initialise the solution pool S 0 with a feasible solution. 3: while (eTime < TL) do 4: 5: 6: 7: 8:

Solve the RMPk ( (35)–(38)) with solution pool S k and the corresponding dual variables be k(k). Set M2 = (k/20)LB⁄ and M1 = M2/10.0. ^ ¼ akðkÞ þ ð1  aÞ Set k k. k. For each i, solve SPi with ^  P P  ðkÞ ðkÞ Compute LB ¼ i Z i ðk Þ  w K w  kðkÞ w .

9: Update solution pool with all the new columns. 10: if (LBk > LB⁄) then 11: Set LB⁄ = LBk. ¼k ^ k /⁄ Update the stability centre⁄/ 12: 13: Call the heuristics explained in Section 5.3 to create globally feasible columns. 14: Remove columns which are non-basic and have no contribution to the optimal solution of the RMP for 10 consecutive iterations. 15: if ððV RMP ðS k Þ  LB Þ=V RMP ðS k Þ < ) then 16: Evaluate the integer version of the MP to identify the best known upper bound UB⁄ and integer solution. 17: break /⁄Since the bounds are close enough, there is scope to improve the dual prices⁄/ 18: k = k + 1 19: Report LB⁄,UB⁄ and gap = (UB⁄  LB⁄)/UB⁄.

953

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956 Table 1 Different CG algorithms and their properties. Version

Stabilisation

Solution (s)

Tag

1 2 3 4

No Yes No Yes

Single Single Multiple Multiple

CG CG CG CG

In Algorithm 1, we apply stabilisation in two stages: (i) Stabilisation at sub-problem level is implemented by using a linear combination of previous and current dual prices (see Step 6) and (ii) Stabilisation at the master problem level is provided by adding a stabilisation term to the objective of the MP (see Section 6.1). Other than these stabilisation techniques, we implemented a feature to store and use multiple solutions (columns) at each iteration (see Step 9). Column management discussed in Section 5.4 is implemented in Step 14. Considering these two factors, Table 1 describes four different versions of Algorithm 1. The abbreviations ‘SS’ and ‘MS’ represent single solution and multiple solutions respectively. 7. Computational experiments In this section we compare the CG algorithm with the integrated model on a set of realistic problem instances. For convenience, the iterative CG Algorithm 1 is represented as ‘CG’. The integrated model is a collection of all sub-problems with an additional linking constraint. It can be expressed as

½IM

min

( X i

 ) X  t Z i ðsi Þ viw ðsi Þ 6 K w ; 8 w; t; si 2 S i ; 8 i  i

ð40Þ

where si corresponds to a schedule (column) of each mine from the solution pool S i , which satisfies constraints (2)–(18) for all i. vtiw ðsi Þ represents the number of active resources (trains) from the class w at time t for the mine i. The above integrated model is referred to as ‘IM’. The IM is implemented and solved using CPLEX 12.1. Thomas et al. (2013) use 240 randomly generated instances used for a comparison between the LR algorithm and the IM. The same datasets are used for the results presented here. The IM and CG schemes were terminated when the optimal solution is found or when a CPU time limit of one hour is reached. An additional 0.1% tolerance was permitted in the final optimality gap for the IM and 0.005% for the lower bounds in CG (see Step 15 of Algorithm 1). However, sub-problems and job scheduling models were terminated at optimality. The lower bound of CG corresponds to the objective value of RMP and the upper bound is computed from the integer MP. Therefore we allow a 2% optimality gap in CG. The randomly generated instances were bundled in eight series with 30 instances per series. Each series represents a scenario with 5, 6, 7, 8, 9, 10, 12 or 15 mines (producers). The reader may refer to Thomas et al. (2013) for a detailed description of the datasets. The optimal objective value of the problem depends on randomly generated demand and due-dates. Therefore it is not easy to capture the impact of the solution approaches on different datasets by comparing objective function face values. Therefore, a normalised performance measure is developed based on the t-test. The relative performance ratio for comparing the lower bounds and upper bounds obtained from CG and IM are defined as

LBCG  LBIM LBR ðCG; IMÞ ¼ and maxðLBCG ; LBIM Þ UBCG  UBIM UBR ðCG; IMÞ ¼ : minðUBCG ; UBIM Þ A positive LB ratio (LBR) means that CG is better, as it has a higher lower bound. Similarly for the upper bound, a positive UB ratio

RMP Reference – – – –

Non-stabilised SS Stabilised SS Non-stabilised MS Stabilised MS

(20)–(22) and (24) (35)–(38) (20)–(22) and (24) (35)–(38)

(UBR) means that IM is better. In some cases UBR cannot be computed since the IM was not able to compute a single feasible solution within the time limit. However, the lower bound computation is very quick and hence, LBR is available for all problem instances. The estimated mean and 95% confidence interval (CI) for the same ratios are computed using the t-test. The CI boundaries and the mean are marked with thick lines.

8. Results and discussions The results are summarised in two sections: Section 8.1 shows the impact of stabilisation of CG algorithms and we choose the best one for further analysis. In Section 8.3, we compare the results from the best CG algorithm with the results from the IM and the LR algorithm discussed in Thomas et al. (2013).

8.1. Impact of stabilisation and improvements The following results compare the four different versions of the CG algorithm given in Table 1. Fig. 1 and Table 2 compare the relative gap and execution time of different CG versions computed on 240 problem instances. We analyse the impact of two features: (i) multiple-solutions (MS) and (ii) stabilisation. MS algorithms have smaller relative gap than their SS equivalent versions. For example, in 91.3% of data instances, non-stabilised SS scheme has less than 10% relative gap. However the equivalent MS scheme has 95.8% of data instances with less than 10% relative gap. Similarly 95.4% is improved to 96.7% in the stabilised version. Having multiple solutions helps to improve the upper bound (integer solution). Sometimes in SS, the CG based algorithm may converge to a common lower bound with a worse integer solution. In the long run both approaches may yield the same solution. However, in a given time limit stabilised versions have smaller gaps. The relative gap computed by any CG algorithm is significantly better than the gap obtained by the IM. When we compare the mean and median execution time, stabilised MS is relatively better than any other models. In Table 2, the clustering around t = 3600 s is due to one of the termination criteria on the maximum CPU run time. The best approach from the four versions is the multiple solutions version with the stabilised column generation scheme. A reduction of 1–5% in relative gap is observed by the application of stabilisation techniques in the 240 datasets. These problem sets have only 150/200 times periods and 3–4 orders per production unit. Therefore the advantage will be significantly higher in large problems.

8.2. LR based solution approach Thomas et al. (2013) proposed a distributed solution approach based on the Lagrangian relaxation to solve this problem. This algorithm was strengthened with the Wedelin algorithm (Wedelin, 1995) and Volume algorithm (Barahona and Anbil, 2000). This iterative scheme is referred to as ‘LR’ in the rest of the paper. The same datasets and termination criteria were used to test this algorithm.

954

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956

Fig. 1. Comparison of relative gap obtained from different CG algorithms.

Table 2 CPU run time comparison of different versions. Model

Non-stabilised

Stabilised

Single-solution

Multiple-solution

Single-solution

Multiple-solution

# Producers

Mean

Median

Mean

Median

Mean

Median

Mean

Median

5 6 7 8 9 10 12 15

1819 2339 1844 2227 2286 2335 3154 2796

1783 2481 1698 2603 2699 3145 3582 3362

2022 2407 1045 1612 1513 1877 2604 2314

1704 2535 1044 1291 1274 1530 3229 2360

2231 2569 1840 2347 2738 2460 3273 2935

2859 3318 1750 2955 3470 2551 3602 3539

1795 2033 1228 1808 1806 2122 2843 2838

1433 1875 1251 1518 1490 1832 3303 3137

8.3. Comparison of the distributed models and the integrated model From the above runs, we identify that the stabilised multiple solutions method is the best among all CG versions. For convenience this version is henceforth referred to as ‘stabilised CG’ or in short CG, unless it is stated otherwise. The LBR (LR, IM) and UBR (LR, IM) were also defined similar to the ratios defined in Section 7. Some of the results and part of the discussions comparing the LR with IM are presented in Thomas et al. (2013). This section emphasise the new results and comparisons with the CG scheme. Figs. 2 and 3 compare LBR and UBR calculated with the stabilised CG, LR and IM. The observations are summarised below: 1. The UBR and LBR increase as the number of mines increases. The trend is identical in both, LR and CG, schemes. The highest LBR is observed for the series with 8 mines in both cases. The highest UBR is obtained for the series with 12 mines in both cases. 2. On an average, LBR for the CG is higher than that of the LR by 1–2%. However, the UBR varies up to 8% since the CG could get better upper bounds than that from the LR. 3. The lower bounds obtained by the CG are always at least 6% better than the lower bounds found by the IM. The difference was much more significant in instances with more than 7 mines; the quality of lower bounds found by the CG was at least 20% better. 4. For series with 5 or 6 mines, IM, LR and CG, were able to converge to a solution that is close to the optimal in the given time. Hence the UBR is close to zero for these series. 5. The CI for the UBR with 12 mines is relatively large. In 26 out of 30 instances; the IM could not find the solution. Therefore the CI has been computed from the remaining 4 values. That might be a reason for the larger interval.

Fig. 2. 95% confidence interval for the relative difference in lower bound.

Fig. 3. 95% confidence interval for the relative difference in upper bound.

955

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956 Table 3 Comparison of the mean LBR and UBR. # Mines ?

5

6

7

8

9

10

12

15

LBR

(LR, IM) (CG, IM)

0.08002 0.09361

0.14533 0.16364

0.14448 0.15273

0.22862 0.24237

0.20223 0.21289

0.21553 0.22873

0.21671 0.23734

0.21723 0.22846

UBR

(LR, IM) (CG, IM)

0.01502 0.01505

0.01124 0.02333

0.04828 0.06203

0.11903 0.13051

0.10914 0.12948

0.29863 0.31871

0.54168 0.62244

0.33725 0.33859

Fig. 4. Relative gap comparison of IM and CG.

Table 4 CPU run time comparison of IM, LR and CG (TL = 3600 seconds). # Mines

5 6 7 8 9 10 12 15

IM

LR

CG

Min

Max

Mean

Median

Min

Max

Mean

Median

Min

Max

Mean

Median

41 1666 369 TL TL TL TL TL

TL TL TL TL TL TL TL TL

2911 3498 3383 TL TL TL TL TL

TL TL TL TL TL TL TL TL

950 2978 TL 357 2058 2867 TL TL

TL TL TL TL TL TL TL TL

3444 TL TL 3424 3572 TL TL TL

TL TL TL TL TL TL TL TL

21 202 114 307 249 284 691 697

TL TL 3541 TL TL TL TL TL

1795 2033 1228 1808 1806 2122 2843 2838

1433 1875 1251 1518 1490 1832 3303 3137

The UBR and LBR increases as the number of mines increases. Figs. 2 and 3 show that the distributed algorithm, specially CG algorithm, is superior to the IM. Table 3 compares the mean LBR, UBR from the two decomposed approaches. As we can see from Table 3, the LBR from the CG is always higher than that from the LR. This implies that the CG scheme is able to get tighter lower bound than that of the LR. The same is true for the upper bound too. For the instances with greater than 7 mines, the mean LBR is greater than 20% and the mean UBR is greater than 10% for both distributed algorithms. It shows the advantage of the distributed model over the integrated model. Fig. 4 shows the histogram of relative gap observed for IM, LR and CG on 240 data sets. Some of the key observations are. 1. The CG outperforms the LB and IM in all cases with significant difference. The number of cases versus relative gap is summarised for all the approaches as Relative Gap (%)

CG

LR

IM

<2 <5 <10

155/240 216/240 232/240

80/240 134/240 198/240

15/240 17/240 24/240

2. When the CG achieves less than 5% gap in 89% of cases, the IM could not find a single solution in 111 (46% of) cases. Hence IM will be ineffective for larger problems. 3. The CG reaches 99% cumulative frequency with a relative gap less than 20% and the LR achieved the same with a relative gap less than 30%. 4. Fig. 4 shows that the CG scheme converges faster than the LR scheme. Table 4 compares the CPU run time of IM, LR and CG. It is evident that the CG has a clear advantage over the LR and IM in terms of run time. The clustering around t = 3600 seconds is due to one of the termination criteria of fixing the maximum CPU run time to one hour. 9. Conclusions We study a resource constrained scheduling problem (RCSP) with multiple producers and a single linking constraint. This is motivated from a real-world mining example. We present mathematical models for the integrated problem, which includes production planning and resource scheduling. The integrated problem is decomposed and solved with the column generation (CG) techniques and compared with a similar Lagrangian relaxation

956

A. Thomas et al. / European Journal of Operational Research 236 (2014) 946–956

(LR) approach (see Thomas et al., 2013). The distributed modelling, solved using a decomposition approach has many advantages over the integrated modelling approach, in term of the quality of solutions and computational time. Computational experiments show that the CG outperforms LR and IM. Although the methods that have been developed in this paper are tested within the context of a coal supply chain, it can be extended to other general RCSPs in contexts such as airline, wine, automobile manufacturing and the service industry. Two major decisions are considered in our approach: production planning and resource scheduling. The upper bound of the distributed model can be further strengthened by the proposed heuristics. Different stability functions can also be explored to further strengthen the performance of the CG. In future research, we plan to develop new algorithms for cases with partial information sharing and also for cases where there are multiple objectives. Ongoing research extends this study by including a third level of decision related to procurement and warehousing. For example, the coal supply chain problem can be extended to truck operations and terminal management operations. It is necessary to test the efficiency of our method on extensions that include more than two sets of decision makers. An alternative approach is decentralised decision making which requires (i) multiple stakeholders, (ii) information asymmetry, (iii) conflict in objectives, and (iv) a negotiation protocol. Our model considers multiple mines that do not interact/share-information with other mines while they constantly interact with the rail operator. The mines and the rail operator possess different sets of information. Hence the decision model has an inherent information asymmetry. This means that the decision model meets the first two requirements. The objectives of (each of) the mines and the rail operator are not same. Every mine expects the rail operator to send a train to the mine to load the coal and deliver it on-time at the terminal. However, the rail operator may have different priorities and constraints. Hence there is a conflict in the objectives of the players. The fourth requirement of ‘negotiation mechanisms’ is not available in our model. In that sense, ours is not a decentralised decision making framework. In our approach, the negotiation protocol is only partially implemented using the CG algorithm which is (in a sense) ‘controlled’ by the ‘honest broker’ rail operator. We made this assumption just so that we could converge to a solution to that is in fact a very difficult real-world problem. In this paper, we have proposed a solution approach for a distributed supply chain which has information asymmetries, and conflict in objectives amongst the partners. The distributed framework can be upgraded to a decentralised framework by introducing a honest broker. A decentralised framework would also have a negotiation scheme. The honest broker who manages the negotiation decides whether a decision is a wrong or right and manages the compensatory schemes for bad solutions that any one producer receives. In our paper, the negotiation protocol is only partially implemented using the CG algorithm which is (in a sense) ‘controlled’ by the ‘honest broker’ rail operator. Further research will need to include proper negotiation mechanisms and ensuring that the rail operator is no longer a ‘de facto honest broker’. Acknowledgements We thank the anonymous reviewers and the editors for their valuable contributions and inputs, which have greatly enhanced the quality of the manuscript.

References Atamtürk, A. (2005). Cover and pack inequalities for (mixed) integer programming. Annals of Operations Research, 139, 21–38. Barahona, F., & Anbil, R. (2000). The volume algorithm: Producing primal solutions with a subgradient method. Mathematical Programming, 87, 385–399. Ben Amor, H., Desrosiers, J., & Frangioni, A. (2009). On the choice of explicit stabilizing terms in column generation. Discrete Applied Mathematics, 157, 1167–1184. Borndorfer, R., Schelten, U., Schlechte, T., & Weider, S. (2005). A column generation approach to airline crew scheduling. Technical Report 05-37 ZIP. Bracken, J., & McGill, J. T. (1973). Mathematical programs with optimization problems in the constraints. Operations Research, 21, 37–44. Brandimarte, P., & Calderini, M. (1995). A hierarchical bicriterion approach to integrated process plan selection and job shop scheduling. International Journal of Production Research, 33, 161. Brucker, P., & Knust, S. (2006). Complex scheduling. Berlin, Heidelberg: SpringerVerlag. Capone, A., Carello, G., Filippini, I., Gualandi, S., & Malucelli, F. (2010). Solving a resource allocation problem in wireless mesh networks: A comparison between a CP-based and a classical column generation. Networks, 55, 221–233. Choi, E., & Tcha, D.-W. (2007). A column generation approach to the heterogeneous fleet vehicle routing problem. Computers & Operations Research, 34, 2080–2095. Dantzig, G. B., & Wolfe, P. (1960). Decomposition principle for linear programs. Operations Research, 8, 101–111. Desrosiers, J., & Lubbecke, M. E. (2005). A primer in column generation. In G. Desaulniers, J. Desrosiers, & M. Solomon (Eds.), Column generation (pp. 1–32). US: Springer. Ebadian, M., Rabbani, M., Torabi, S. A., & Jolai, F. (2009). Hierarchical production planning and scheduling in make-to-order environments: reaching short and reliable delivery dates. International Journal of Production Research, 47, 5761–5789. Hans, E. W., Herroelen, W., Leus, R., & Wullink, G. (2007). A hierarchical approach to multi-project planning under uncertainty. Omega, 35, 563–577. Huisman, D., Jans, R., Peeters, M., & Wagelmans, A. P. M. (2005). Combining column generation and Lagrangian relaxation. In G. Desaulniers, J. Desrosiers, & M. Solomon (Eds.), Column generation (pp. 247–270). Kelly, J. D., & Zyngier, D. (2008). Hierarchical decomposition heuristic for scheduling: Coordinated reasoning for decentralized and distributed decisionmaking problems. Computers and Chemical Engineering, 32, 2684–2705. Lenstra, J. K., Kan, A. H. G. R., & Brucker, P. (1977). Complexity of machine scheduling problems. Annals of Discrete Mathematics, 1, 343–362. Lubbecke, M. E., & Desrosiers, J. (2005). Selected topics in column generation. Operations Research, 53, 1007–1023. Pessoa, A., Uchoa, E., de Aragao, M. P., & Rodrigues, R. (2010). Exact algorithm over an arc-time-indexed formulation for parallel machine scheduling problems. Mathematical Programming Computation, 2, 1–32. Pochet, Y., & Wolsey, L. A. (2006). Production planning by mixed integer programming. NY: Springer. Singh, G., Dunstall, S., Eitzen, G., Ernst, A. T., Horn, M., & Krishnamoorthy, M. (2005). On coordination and scheduling in wine supply chains. In 5th International conference on operational research for development (ICORD), Jamshedpur, India (pp. 219–229). Singh, G., Sier, D., Ernst, A. T., Gavriliouk, O., Oyston, R., Giles, T., et al. (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. Thomas, A., Singh, G., Krishnamoorthy, M., & Venkateswaran, J. (2013). Distributed optimisation method for multi-resource constrained scheduling in coal supply chains. International Journal of Production Research, 51, 2740–2759. Timm, T., & Blecken, A. (2011). A method for the hierarchical planning of the structure, dimension and material requirements of manufacturing systems. International Journal of Production Research, 49, 3431–3453. Vanderbeck, F., & Wolsey, L. A. (2010). Reformulation and decomposition of integer programs. In M. Junger, T. Liebling, D. Naddef, G. L. Nemhauser, W. R. Pulleyblank, & G. Reinelt, et al. (Eds.), 50 Years of integer programming 1958– 2008 (pp. 431–504). Berlin: Springer. Wang, G. L., Lin, L., & Zhong, S. S. (2008). Task-oriented hierarchical planning and scheduling for discrete manufacturing enterprise. Applied Mechanics and Materials, 114–120. Wedelin, D. (1995). An algorithm for large scale 0–1 integer programming with application to airline crew scheduling. Annals of Operations Research, 57, 283–301. Wentges, P. (1997). Weighted Dantzig–Wolfe decomposition for linear mixed-integer programming. International Transactions in Operational Research, 4, 151–162. Wu, D., & Ierapetritou, M. (2007). Hierarchical approach for production planning and scheduling under uncertainty. Chemical Engineering and Processing: Process Intensification, 46, 1129–1140. Zribi, N., Kacem, I., El Kamel, A., & Borne, P. (2007). Assignment and scheduling in flexible job-shops by hierarchical optimization. IEEE Transactions on Systems, Man, and Cybernetics, Part C: Applications and Reviews, 37, 652–661.