# A Column Generation Algorithm for the Resource-Constrained Order Acceptance and Scheduling on Unrelated Parallel Machines

Authors: Yujian Song; Ming Xue; Changhua Hua; Wanli Wang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, we investigate the resource-constrained order acceptance and scheduling on unrelated parallel machines that arise in make-to-order systems. The objective of this problem is to simultaneously select a subset of orders to be processed and schedule the accepted orders on unrelated machines in such a way that the resources are not overutilized at any time. We first propose two formulations for the problem: mixed integer linear programming formulation and set partitioning. In view of the complexity of the problem, we then develop a column generation approach based on the set partitioning formulation. In the proposed column generation approach, a differential evolution algorithm is designed to solve subproblems efficiently. Extensive numerical experiments on different-sized instances are conducted, and the results demonstrate that the proposed column generation algorithm reports optimal or near-optimal solutions that are evidently better than the solutions obtained by solving the mixed integer linear programming formulation.


---

## Body

## 1. Introduction

In many industries, the increasing demand for product variety and customization has intensified competitions among manufacturers to such an extent that attracting and retaining customers is more important than before. The fierce competition has stimulated manufacturers to switch from make-to-stock production to make-to-order production. In the make-to-order production system, manufacturing of the product begins after receiving a customer order. Undoubtedly, accepting all potential orders may cause overloads, order delays, and customer dissatisfaction, due to the limited production capacity and tight delivery requirements. Therefore, the first and foremost decision for a make-to-order system is to decide simultaneously which orders to accept and how to schedule the accepted orders.Due to its application in diverse industries, the order acceptance and scheduling problem has drawn considerate attention from academics as well as practitioners. It has been considered in single machine, parallel machine, and flow shop environments. Various formulations and solution methods for this problem have been proposed. Although there is abundant literature on the order acceptance and scheduling problem, most of these studies treat machine as the only required resource. However, in real-world make-to-order production systems, other resources such as machine operators, tools, fixtures, or industrial robots are also needed and restricted in quantity. An example can be found in the plastic injection molding industry. Plastic injection molding is a production method used to produce plastic products by injecting heated plastic into a preformed mold. The molten plastic is injected at a predetermined temperature and pressure and left to cool. When it is solidified and adapted the shape of the mold, it is ejected as a formed product and the molding cycle is repeated. Figure1 illustrates the process of plastic injection molding. With no doubt, the injection machine can only start its processing when the mold suitable for the products to be produced is mounted on the injection machine. In general, the molds are quite expensive, and thus the number of molds available is limited. It is obvious that when a manufacturing system like this has limited additional resources besides machines, and thus is unable to meet all customers’ demand, it has to decide which orders are to be accepted for processing as well as their associated schedule. Therefore, it is of great practical significance to investigate the resource-constrained order acceptance and scheduling problem.Figure 1
The process of plastic injection molding (http://www.aireplatstic.com).In this paper, we try to fill the gap and address the resource-constrained order acceptance and scheduling problem on unrelated parallel machines. The objective of this problem is to simultaneously select a subset of orders to be processed and schedule the accepted orders on unrelated machines so as to ensure that the resources are not overutilized at any time. The main contributions of this research are threefold.(1)
We consider a new variant of the order acceptance and scheduling problem on unrelated machines, which takes the additional limited resources into consideration. We present an integer programming model and a set partitioning formulation for this problem.(2)
We develop a column generation algorithm based on the set partitioning formulation. In order to solve the pricing subproblems efficiently, we devise a differential evolution algorithm to rapidly price out promising columns.(3)
We conduct comprehensive numerical experiments on randomly generated instances to illustrate the efficacy of the proposed column generation algorithm.The remainder of the paper is organized as follows. Section2 reviews relevant literature. Problem definition and the integer linear programming model are provided in Section 3, followed by a description of the proposed column generation algorithm and its implementation in Section 4. In Section 5, we discuss our numerical results and validate the effectiveness and the efficacy of the proposed solution approach. Section 6 concludes the contribution and recommends directions for future research.


## 2. Literature Review

The scheduling problem addressed in this work is related to two fields of research: (i) order acceptance and scheduling and (ii) resource-constrained unrelated parallel machine scheduling. Existing studies on both topics is reviewed in this section.


### 2.1. The Order Acceptance and Scheduling on Unrelated Parallel Machines

Due to its wide applications in the real world, the order acceptance and scheduling problem and its variants have been extensively studied. Interested readers can read Slotnick’s study [1] for a comprehensive literature review of the order acceptance and scheduling problem. The following survey only includes a representative subset of papers addressing the order acceptance and scheduling problem on unrelated parallel machines.Fanjul-Peyro and Ruiz [2] investigated two generalizations of the unrelated parallel machine scheduling problem, one of which is concerned with the integrated order selection and scheduling decision making problem. A mixed integer programming mathematical formulation is presented, along with a constructive method specifically tailored for such a generalization. Comprehensive numerical and statistical experiments show that the proposed method significantly improves the results yielded by the CPLEX solver. Hsu and Chang [3] dedicatedly study the unrelated parallel-machine scheduling problem with deteriorating jobs and rejection. The objective of this problem is to reject a subset of jobs and then schedule the accepted jobs in such a way to minimize the total cost. Computational experiments show that the problem can be solved in polynomial time when the quantity of machines is fixed. Jiang and Tan [4] study the unrelated parallel machine scheduling model with job acceptance option. Each job order is either accepted and processed without interruption, or is declined at a penalty cost. In an attempt to minimize the makespan of all accepted jobs, the penalty cost of rejected jobs, and the production cost of accepted jobs, the authors first develop a polynomial time algorithm to derive a lower bound and then propose a heuristic by introducing various relaxation techniques. The worst-case ratio bound of the proposed approach is proven to be 2.Emami et al. [5] consider the simultaneous order acceptance and scheduling problem in an unrelated parallel machines environment with the assumption of uncertain revenue and job-processing times. The objective of this problem is to find a subset of orders and schedule them in such a way that the total profit is maximized. A Lagrangian relaxation algorithm enhanced by a cutting plane technique is implemented and capable of solving instances with up to 40 orders and 6 machines. The same problem is studied in Emami et al. [6] and solved by Benders decomposition approach. Wang and Wang [7] mainly focus on the biobjective optimization of order acceptance and scheduling on unrelated parallel machines with machine eligibility constraints, where a subset of jobs has to be selected to simultaneously maximize the profit level and to minimize the total completion time. The problem is formulated as a biobjective mixed integer linear programming. In view of its NP-hardness, a multiobjective parthenogenetic algorithm with parthenogenetic operators and Pareto-ranking and selection methods are implemented to obtain approximate Pareto front. Mor and Mosheiov [8] study the order acceptance and scheduling problem on unrelated parallel machines, which takes the position-dependent processing times into consideration. This problem aims at minimizing total absolute deviation of job completion times. It is shown that the problem is solved in polynomial time in the number of jobs, and a polynomial time solution approach based on solving a sequence of linear assignment problems is designed to solve this problem. Wang and Ye [9] investigate the order acceptance and scheduling problem on unrelated parallel machines. The purpose is to determine the orders to be accepted for processing and the processing sequence for the accepted orders to get the optimal profit. Two mixed integer programming formulations are presented, which are further enhanced by various valid inequalities. For the purpose of solving complicated instances, the authors propose a tailored branch-and-bound algorithm. The proposed branch-and-bound algorithm which follows the paradigm of “divide and conquer” branches on the number of accepted orders and recursively solves the reduced unrelated parallel machine scheduling problem at each child node. Computational experiments confirm its effectiveness.


### 2.2. The Resource-Constrained Parallel Machine Scheduling

Despite the plentiful achievements mentioned above, research on order acceptance and scheduling with explicit consideration of additional constrained resources is somewhat limited. Since the considered problem is closely related to the resource-constrained parallel machine scheduling, an insight into its recent advances has substantial implications in mathematical modeling and solution development for the studied problem. Therefore, a brief literature review on it is presented as follows.Ventura and Kim [10] investigate a resource-constrained parallel machine scheduling problem which arises in the burn-in process of semiconductor manufacturing. The problem is formulated as an integer programming model with the objective of minimizing the total absolute deviations of job completion times from due dates. The problem is proved to be solved in polynomial time when there is one single type of additional resource, and each job consumes one unit of additional resource at most. Ventura and Kim [11] extend the work of [10] by considering that the number of kinds of additional resources and the resource requirements per job are unrestricted. This more general problem is modelled as a 0-1 integer program, and the Lagrangian relaxation algorithm is applied to obtain high-quality solutions and tight lower bounds for large-scale problems. Chen [12] studies the resource-constrained unrelated parallel machine scheduling problem which originates from the job scheduling in the plastics-forming industry. Beside plastic moulding machines, scarce and costly dyes are also indispensable equipment and thus treated as secondary resource constraints for plastic forming operations. A heuristic is developed to minimize the makespan. The computational results show that it outperforms the simulation annealing algorithm in terms of efficiency and effectiveness. This work is further improved by Chen and Wu [13], who propose an effective tabu search algorithm. Ruiz-Torres et al. [14] investigate the problem of scheduling jobs on uniform parallel machines, where the speed of the machine is determined by the quantity of the assigned additional resource. The objective aim is to minimize the total number of delayed jobs. Two different versions of the problem are discussed. The first one makes the assumption that all jobs are preallocated to the machines, while the second one considers simultaneous optimization of job scheduling and resource allocation. An integer programming formulation is presented for the first case and five heuristics for the second. Edis et al. [15] comprehensively survey the related literature on resource-constrained parallel machine scheduling. They classify the studies with respect to the characteristics of machine environments, objective functions, additional resource types, complexity of the problem and solution approaches, and other important issues. Strengths and limitations of the existing literature as well as open areas for future studies are also highlighted in this paper. Afzalirad and Rezaeian [16] deal with a resource-constrained unrelated parallel machine scheduling problem which considers various realistic features, such as release dates, sequence-dependent setup times, machine eligibility, and precedence restrictions. An integer programming model is proposed to represent the complicated problem mathematically. Two different heuristic algorithms including genetic algorithm and artificial immune algorithm are designed to identify high-quality solutions. Comparative experiments are performed on a set of randomly generated test instances and demonstrate that both algorithms perform rather well for small-scale instances and the artificial immune algorithm is suggested to solve large-scale instances. Li et al. [17] delicately address a job-scheduling problem arising in the manufacturing systems exposed to both machine and worker availability constraints. A branch population genetic algorithm is proposed to minimize the makespan and cost. The effectiveness of the proposed algorithm is validated by comparative experiments and an actual case study, with the conclusion that the proposed algorithm is very useful under the realistic conditions. Fanjul-Peyro et al. [18] investigate a resource-constrained parallel machine scheduling problem where the processing of jobs on the machines demands a number of units of a scarce resource. Two mixed integer linear programming (MILP) models are presented to solve small-sized instances. Besides, three meta-heuristic algorithms based on the MILP models are developed to solve large-sized instances. Vallada et al. [19] mainly focus on the solution development for the resource-constrained unrelated parallel machine scheduling problem. They propose a basic scatter search algorithm, an enriched scatter search algorithm and an enriched iterated greedy algorithm, to solve this problem. In all the proposed methods, nonfeasible solutions with violation of resource limitation are allowed and further fixed to a feasible solution using a repairing mechanism. Different local search procedures based on insertion, swap, and restricted neighborhoods are implemented to enhance the performance of the proposed approaches. Nguyen et al. [20] deal with a resource-constrained job scheduling problem abstracted from mineral transportation. A parallel differential evolution algorithm is implemented in conjunction with iterated greedy search and column generation improvement strategies to solve this problem. They conduct an exhaustive numerical evaluation of the proposed and existing methods using different-sized instances. The computational results confirm the efficiency and effectiveness of the proposed algorithm and advantage of the adopted improvement strategies. More recently, Yepes-Borrero et al. [21] address the unrelated parallel machine scheduling problem with setup times and additional resources in the setups. A MILP formulation and three meta-heuristics are proposed for this problem.In summary, there are a number of studies on order acceptance and scheduling as well as resource-constrained parallel machine scheduling, but this field of research still requires further development due to its widespread applications. The key extension of our paper which differentiates ours from existing studies is the consideration of order acceptance and resource-constrained parallel machine scheduling in an integrated manner.


## 2.1. The Order Acceptance and Scheduling on Unrelated Parallel Machines

Due to its wide applications in the real world, the order acceptance and scheduling problem and its variants have been extensively studied. Interested readers can read Slotnick’s study [1] for a comprehensive literature review of the order acceptance and scheduling problem. The following survey only includes a representative subset of papers addressing the order acceptance and scheduling problem on unrelated parallel machines.Fanjul-Peyro and Ruiz [2] investigated two generalizations of the unrelated parallel machine scheduling problem, one of which is concerned with the integrated order selection and scheduling decision making problem. A mixed integer programming mathematical formulation is presented, along with a constructive method specifically tailored for such a generalization. Comprehensive numerical and statistical experiments show that the proposed method significantly improves the results yielded by the CPLEX solver. Hsu and Chang [3] dedicatedly study the unrelated parallel-machine scheduling problem with deteriorating jobs and rejection. The objective of this problem is to reject a subset of jobs and then schedule the accepted jobs in such a way to minimize the total cost. Computational experiments show that the problem can be solved in polynomial time when the quantity of machines is fixed. Jiang and Tan [4] study the unrelated parallel machine scheduling model with job acceptance option. Each job order is either accepted and processed without interruption, or is declined at a penalty cost. In an attempt to minimize the makespan of all accepted jobs, the penalty cost of rejected jobs, and the production cost of accepted jobs, the authors first develop a polynomial time algorithm to derive a lower bound and then propose a heuristic by introducing various relaxation techniques. The worst-case ratio bound of the proposed approach is proven to be 2.Emami et al. [5] consider the simultaneous order acceptance and scheduling problem in an unrelated parallel machines environment with the assumption of uncertain revenue and job-processing times. The objective of this problem is to find a subset of orders and schedule them in such a way that the total profit is maximized. A Lagrangian relaxation algorithm enhanced by a cutting plane technique is implemented and capable of solving instances with up to 40 orders and 6 machines. The same problem is studied in Emami et al. [6] and solved by Benders decomposition approach. Wang and Wang [7] mainly focus on the biobjective optimization of order acceptance and scheduling on unrelated parallel machines with machine eligibility constraints, where a subset of jobs has to be selected to simultaneously maximize the profit level and to minimize the total completion time. The problem is formulated as a biobjective mixed integer linear programming. In view of its NP-hardness, a multiobjective parthenogenetic algorithm with parthenogenetic operators and Pareto-ranking and selection methods are implemented to obtain approximate Pareto front. Mor and Mosheiov [8] study the order acceptance and scheduling problem on unrelated parallel machines, which takes the position-dependent processing times into consideration. This problem aims at minimizing total absolute deviation of job completion times. It is shown that the problem is solved in polynomial time in the number of jobs, and a polynomial time solution approach based on solving a sequence of linear assignment problems is designed to solve this problem. Wang and Ye [9] investigate the order acceptance and scheduling problem on unrelated parallel machines. The purpose is to determine the orders to be accepted for processing and the processing sequence for the accepted orders to get the optimal profit. Two mixed integer programming formulations are presented, which are further enhanced by various valid inequalities. For the purpose of solving complicated instances, the authors propose a tailored branch-and-bound algorithm. The proposed branch-and-bound algorithm which follows the paradigm of “divide and conquer” branches on the number of accepted orders and recursively solves the reduced unrelated parallel machine scheduling problem at each child node. Computational experiments confirm its effectiveness.


## 2.2. The Resource-Constrained Parallel Machine Scheduling

Despite the plentiful achievements mentioned above, research on order acceptance and scheduling with explicit consideration of additional constrained resources is somewhat limited. Since the considered problem is closely related to the resource-constrained parallel machine scheduling, an insight into its recent advances has substantial implications in mathematical modeling and solution development for the studied problem. Therefore, a brief literature review on it is presented as follows.Ventura and Kim [10] investigate a resource-constrained parallel machine scheduling problem which arises in the burn-in process of semiconductor manufacturing. The problem is formulated as an integer programming model with the objective of minimizing the total absolute deviations of job completion times from due dates. The problem is proved to be solved in polynomial time when there is one single type of additional resource, and each job consumes one unit of additional resource at most. Ventura and Kim [11] extend the work of [10] by considering that the number of kinds of additional resources and the resource requirements per job are unrestricted. This more general problem is modelled as a 0-1 integer program, and the Lagrangian relaxation algorithm is applied to obtain high-quality solutions and tight lower bounds for large-scale problems. Chen [12] studies the resource-constrained unrelated parallel machine scheduling problem which originates from the job scheduling in the plastics-forming industry. Beside plastic moulding machines, scarce and costly dyes are also indispensable equipment and thus treated as secondary resource constraints for plastic forming operations. A heuristic is developed to minimize the makespan. The computational results show that it outperforms the simulation annealing algorithm in terms of efficiency and effectiveness. This work is further improved by Chen and Wu [13], who propose an effective tabu search algorithm. Ruiz-Torres et al. [14] investigate the problem of scheduling jobs on uniform parallel machines, where the speed of the machine is determined by the quantity of the assigned additional resource. The objective aim is to minimize the total number of delayed jobs. Two different versions of the problem are discussed. The first one makes the assumption that all jobs are preallocated to the machines, while the second one considers simultaneous optimization of job scheduling and resource allocation. An integer programming formulation is presented for the first case and five heuristics for the second. Edis et al. [15] comprehensively survey the related literature on resource-constrained parallel machine scheduling. They classify the studies with respect to the characteristics of machine environments, objective functions, additional resource types, complexity of the problem and solution approaches, and other important issues. Strengths and limitations of the existing literature as well as open areas for future studies are also highlighted in this paper. Afzalirad and Rezaeian [16] deal with a resource-constrained unrelated parallel machine scheduling problem which considers various realistic features, such as release dates, sequence-dependent setup times, machine eligibility, and precedence restrictions. An integer programming model is proposed to represent the complicated problem mathematically. Two different heuristic algorithms including genetic algorithm and artificial immune algorithm are designed to identify high-quality solutions. Comparative experiments are performed on a set of randomly generated test instances and demonstrate that both algorithms perform rather well for small-scale instances and the artificial immune algorithm is suggested to solve large-scale instances. Li et al. [17] delicately address a job-scheduling problem arising in the manufacturing systems exposed to both machine and worker availability constraints. A branch population genetic algorithm is proposed to minimize the makespan and cost. The effectiveness of the proposed algorithm is validated by comparative experiments and an actual case study, with the conclusion that the proposed algorithm is very useful under the realistic conditions. Fanjul-Peyro et al. [18] investigate a resource-constrained parallel machine scheduling problem where the processing of jobs on the machines demands a number of units of a scarce resource. Two mixed integer linear programming (MILP) models are presented to solve small-sized instances. Besides, three meta-heuristic algorithms based on the MILP models are developed to solve large-sized instances. Vallada et al. [19] mainly focus on the solution development for the resource-constrained unrelated parallel machine scheduling problem. They propose a basic scatter search algorithm, an enriched scatter search algorithm and an enriched iterated greedy algorithm, to solve this problem. In all the proposed methods, nonfeasible solutions with violation of resource limitation are allowed and further fixed to a feasible solution using a repairing mechanism. Different local search procedures based on insertion, swap, and restricted neighborhoods are implemented to enhance the performance of the proposed approaches. Nguyen et al. [20] deal with a resource-constrained job scheduling problem abstracted from mineral transportation. A parallel differential evolution algorithm is implemented in conjunction with iterated greedy search and column generation improvement strategies to solve this problem. They conduct an exhaustive numerical evaluation of the proposed and existing methods using different-sized instances. The computational results confirm the efficiency and effectiveness of the proposed algorithm and advantage of the adopted improvement strategies. More recently, Yepes-Borrero et al. [21] address the unrelated parallel machine scheduling problem with setup times and additional resources in the setups. A MILP formulation and three meta-heuristics are proposed for this problem.In summary, there are a number of studies on order acceptance and scheduling as well as resource-constrained parallel machine scheduling, but this field of research still requires further development due to its widespread applications. The key extension of our paper which differentiates ours from existing studies is the consideration of order acceptance and resource-constrained parallel machine scheduling in an integrated manner.


## 3. Problem Definition and Mathematical Model

In this section, we mainly elaborate on the problem statement and provide an integer programming formulation.


### 3.1. Problem Statement and Integer Programming Formulation

This paper considers a setting where a manufacture receives orders from its customers but is confronted by additional limited resources. Specifically, letM = {1, ..., |M|} be the set of unrelated machines and N = {1, ..., |N|} be the set of orders. The maximum amount of resource k {k∈K} is constrained to be Qmaxk. Each machine can process only one order at any time and is always available within the planning horizon T = {0, 1, ..., |T|}. Each order can be processed on atmost one machine without interruption and requires a constant number of the additional resources during its entire process. Thus, each order can only be processed when one machine and sufficient additional resources are available during its processing time.Associated with orderi are as follows: (1) release date ri; (2) processing time pim, corresponding to the time required to process order i on machine m; (3) due date di; (4) resource requirement of order i to resource k, qik; (5) revenue ei gained from its customer if it is accepted; and (6) penalty wi for each unit time of tardiness. Since tardiness penalties will incur a loss of revenue, the objective of this problem is to select a subset of orders and schedule them in such a way that all constraints are satisfied and the total net revenue is maximized. To formulate the problem, the following notations are defined.Parameters:T: set of time periods in the scheduling horizon, T = {0, 1, ..., |T|}M: set of unrelated machines, M = {1, ..., |M|}N: set of orders, N = {1, ..., |N|}K: set of resource types, N = {1, ..., |K|}ri: release time of order ipim: processing time of order i on machine mdi: due time of order iqik: the amount of resource k required for the processing of order iei: revenue of order iwi: unit penalty cost of order iQmaxk: maximum amount of the additional limited resource kDecision variables:Xi: a binary variable, equal to 1 if order i is accepted and 0 otherwiseYitm: a binary variable, equal to 1 if order i completes its processing at time t on machine m and 0 otherwiseCi: completion time of order iTi: tardiness of order iThe problem described above can be mathematically formulated as the following MILP model.(1)max∑i∈Nei∗Xi−wi∗Ti,(2)s.t.∑i∈N∑s∈t,t+pim−1Yism≤1,∀m∈M,t∈T,(3)∑m∈M∑t∈ri+pim,TYitm≤1,∀i∈N,(4)∑m∈M∑i∈N∑s∈t,t+pim−1qik∗Yism≤Qmaxk,∀t∈T,k∈K,(5)∑m∈M∑t∈0,ri+pim−1Yitm=0,∀i∈N,(6)∑m∈M∑t∈Tt∗Yitm=Ci,∀i∈N,(7)∑m∈M∑t∈TYitm=Xi,∀i∈N,(8)Ti≥Ci−di,∀i∈N,(9)Ci,Ti≥0,Xi,Sitm∈0,1,∀i∈N,t∈T,m∈M.The objective function (1) seeks to maximize the total net revenue of accepted orders, which is computed by the total revenue minus the total tardiness penalty. Constraint (2) ensures that not more than one order can be assigned to any machine at any time period. Constraint (3) specifies that an order is either accepted or rejected. Constraint (4) restricts that the total amount of the additional resource k consumed at any time should not exceed the available amount Qmaxk. Constraint (5) restricts that an order cannot be processed before its release time. Constraint (6) computes the completion time for each order. Constraint (7) links variable Yitm and variable Xi and is introduced to determine if an order is accepted or not. Constraint (8) calculates the tardiness for each order. Constraint (9) specifies the domains of decision variables.


### 3.2. An Illustrative Example

To verify the proposed MILP formulation and intuitively understand this complicated scheduling problem, we consider a simple example involving 6 jobs, 2 machines, and one additional resource type. The maximum amount of the additional resource is 6 units. The parameters related to the example are provided in Table1. The mathematical formulation of this example is solved to optimality using CPLEX 12.10 software. The net revenue of the optimal solution is 48. In Figure 2, the optimal schedule is visualized in the Gantt chart. The horizontal axis of the chart shows the time when each order is under processing, while the vertical axis indicates the name of each machine. Each block on the chart corresponds to an order, which is marked by its order name and the amount of resource required for processing.Table 1
Example problem data.
OrderriqidieiwipimMachine 1Machine 2Order 112510232Order 223612323Order 334610433Order 43379131Order 544811212Order 643910266Figure 2
The optimal solution to the illustrative example. (a) The Gantt chart of the optimal schedule. (b) Resource consumption.
(a)(b)It can be seen that all orders except order 3 are accepted. We can also observe that machine 2 is idle from time 4 to time 5 since order 5 and order 6 cannot be processed at the same time due to limited resources. From this example, we can know that the considered problem is difficult to solve because its solution requires not only selecting a subset of orders but also scheduling them on machines in such a way that the resources are not overused at any time.


### 3.3. Complexity of the Problem

Theorem 1.
The resource-constrained order acceptance and scheduling on unrelated parallel machines is strongly NP-hard.Proof.
When we set the maximum amount of the additional resource to sufficiently large and introduce a dummy machine to which all rejected orders are assigned with zero processing time, zero revenue, and zero penalty, the studied problem can be reduced to the classical unrelated parallel machine scheduling problem which is known to be strongly NP-hard (Liaw et al. [22]). Because the classical unrelated parallel machine scheduling problem is a special case of the considered problem, it follows that the considered problem is also strongly NP-hard.


## 3.1. Problem Statement and Integer Programming Formulation

This paper considers a setting where a manufacture receives orders from its customers but is confronted by additional limited resources. Specifically, letM = {1, ..., |M|} be the set of unrelated machines and N = {1, ..., |N|} be the set of orders. The maximum amount of resource k {k∈K} is constrained to be Qmaxk. Each machine can process only one order at any time and is always available within the planning horizon T = {0, 1, ..., |T|}. Each order can be processed on atmost one machine without interruption and requires a constant number of the additional resources during its entire process. Thus, each order can only be processed when one machine and sufficient additional resources are available during its processing time.Associated with orderi are as follows: (1) release date ri; (2) processing time pim, corresponding to the time required to process order i on machine m; (3) due date di; (4) resource requirement of order i to resource k, qik; (5) revenue ei gained from its customer if it is accepted; and (6) penalty wi for each unit time of tardiness. Since tardiness penalties will incur a loss of revenue, the objective of this problem is to select a subset of orders and schedule them in such a way that all constraints are satisfied and the total net revenue is maximized. To formulate the problem, the following notations are defined.Parameters:T: set of time periods in the scheduling horizon, T = {0, 1, ..., |T|}M: set of unrelated machines, M = {1, ..., |M|}N: set of orders, N = {1, ..., |N|}K: set of resource types, N = {1, ..., |K|}ri: release time of order ipim: processing time of order i on machine mdi: due time of order iqik: the amount of resource k required for the processing of order iei: revenue of order iwi: unit penalty cost of order iQmaxk: maximum amount of the additional limited resource kDecision variables:Xi: a binary variable, equal to 1 if order i is accepted and 0 otherwiseYitm: a binary variable, equal to 1 if order i completes its processing at time t on machine m and 0 otherwiseCi: completion time of order iTi: tardiness of order iThe problem described above can be mathematically formulated as the following MILP model.(1)max∑i∈Nei∗Xi−wi∗Ti,(2)s.t.∑i∈N∑s∈t,t+pim−1Yism≤1,∀m∈M,t∈T,(3)∑m∈M∑t∈ri+pim,TYitm≤1,∀i∈N,(4)∑m∈M∑i∈N∑s∈t,t+pim−1qik∗Yism≤Qmaxk,∀t∈T,k∈K,(5)∑m∈M∑t∈0,ri+pim−1Yitm=0,∀i∈N,(6)∑m∈M∑t∈Tt∗Yitm=Ci,∀i∈N,(7)∑m∈M∑t∈TYitm=Xi,∀i∈N,(8)Ti≥Ci−di,∀i∈N,(9)Ci,Ti≥0,Xi,Sitm∈0,1,∀i∈N,t∈T,m∈M.The objective function (1) seeks to maximize the total net revenue of accepted orders, which is computed by the total revenue minus the total tardiness penalty. Constraint (2) ensures that not more than one order can be assigned to any machine at any time period. Constraint (3) specifies that an order is either accepted or rejected. Constraint (4) restricts that the total amount of the additional resource k consumed at any time should not exceed the available amount Qmaxk. Constraint (5) restricts that an order cannot be processed before its release time. Constraint (6) computes the completion time for each order. Constraint (7) links variable Yitm and variable Xi and is introduced to determine if an order is accepted or not. Constraint (8) calculates the tardiness for each order. Constraint (9) specifies the domains of decision variables.


## 3.2. An Illustrative Example

To verify the proposed MILP formulation and intuitively understand this complicated scheduling problem, we consider a simple example involving 6 jobs, 2 machines, and one additional resource type. The maximum amount of the additional resource is 6 units. The parameters related to the example are provided in Table1. The mathematical formulation of this example is solved to optimality using CPLEX 12.10 software. The net revenue of the optimal solution is 48. In Figure 2, the optimal schedule is visualized in the Gantt chart. The horizontal axis of the chart shows the time when each order is under processing, while the vertical axis indicates the name of each machine. Each block on the chart corresponds to an order, which is marked by its order name and the amount of resource required for processing.Table 1
Example problem data.
OrderriqidieiwipimMachine 1Machine 2Order 112510232Order 223612323Order 334610433Order 43379131Order 544811212Order 643910266Figure 2
The optimal solution to the illustrative example. (a) The Gantt chart of the optimal schedule. (b) Resource consumption.
(a)(b)It can be seen that all orders except order 3 are accepted. We can also observe that machine 2 is idle from time 4 to time 5 since order 5 and order 6 cannot be processed at the same time due to limited resources. From this example, we can know that the considered problem is difficult to solve because its solution requires not only selecting a subset of orders but also scheduling them on machines in such a way that the resources are not overused at any time.


## 3.3. Complexity of the Problem

Theorem 1.
The resource-constrained order acceptance and scheduling on unrelated parallel machines is strongly NP-hard.Proof.
When we set the maximum amount of the additional resource to sufficiently large and introduce a dummy machine to which all rejected orders are assigned with zero processing time, zero revenue, and zero penalty, the studied problem can be reduced to the classical unrelated parallel machine scheduling problem which is known to be strongly NP-hard (Liaw et al. [22]). Because the classical unrelated parallel machine scheduling problem is a special case of the considered problem, it follows that the considered problem is also strongly NP-hard.


## 4. Column Generation Approach

The MILP formulation presented in Section3 is capable of solving small instances to optimality in reasonable time, but it becomes computationally expensive with the increase of the problem size, as shown in Section 5. This fact necessitates the development of efficient solution methodologies. Column generation is an effective and efficient methodology which has been well studied and successfully applied to a wide range of problems. In this section, we provide an alternative set partitioning formulation and then develop a tailored column generation approach to solve the considered problem.


### 4.1. An Overview of the Proposed Column Generation Algorithm

Column generation is an efficient approach for solving linear programming problems with a large number of variables (columns). Its main advantage is that not all columns need to be enumerated. Instead, the problem is first formulated as a restricted linear master problem (RLMP). This RLMP considers only a subset of columns. And then, one or more subproblems with respect to current dual solution to the RLMP are solved to identify favorable new columns which can be incorporated into the RLMP to improve the objective value. The RLMP and subproblems are solved iteratively until no promising columns can be further added to the RLMP. This termination condition implies that the RLMP is optimal. The RLMP probably returns a fractional solution. In this paper, in order to obtain an integer solution, the RLMP with all columns which have ever been generated as decision variables is converted to an IP model and then solved by an IP solver.In an attempt to enhance the performance of the column generation, we also propose a differential evolution (DE) algorithm to solve the subproblems efficiently. Besides, we develop a constructive heuristic to obtain the initial solution for the RLMP. The procedures of the proposed column generation approach are shown in Algorithm1. Next, we present a detailed description of each step.Algorithm 1: The procedure of the proposed column generation.
Step 1: implement the proposed constructive heuristic to obtain a feasible solutionStep 2: initialize the RLMP with the solution obtained by the constructive heuristicStep3: set the number of consecutive iterations of failing to find a column with positive reduced costIternon to be 0, i.e., Iternon=0Step 4: whileIternon is less than a specified number IternonCG, doStep 4.1: solve the RLMP and get the dual solutionStep 4.2: define the subproblems by the current dual solution and solve them by the proposed DE algorithmStep 4.3: If there exist columns with positive reduced cost, doAdd new columns to the RLMP;Iternon=0ElseIternon=Iternon+1EndEndStep 5: convert the RLMP model to an IP modelStep 6: solve the IP model by CPLEXStep 7: output the solution


### 4.2. Master Problem: A Set Partitioning Formulation

In this section, we reformulate the considered problem as a set partitioning formulation which serves as the master problem in the column generation framework. The set partitioning formulation which is a special class of linear integer programs has been utilized to formulate a wide range of practical problems in fields such as crew scheduling, vehicle routing, and timetable scheduling. Based on the special structure of the set partitioning formulation, efficient approaches techniques, for example, column generation, have been proposed to solve large-sized problems, making it possible to deal with many real-life applications. Therefore, we attempt to propose an alternative mathematical formulation based on the set partitioning formulation.In the set partitioning formulation, there is a set of columns (assignments) denoted by Ω that composes the set partitioning formulation matrix. Each column corresponds to a feasible assignment of a single order to a machine in accordance with the time-related constraints. This matrix contains three submatrices, denoted asA, B, and C, all containing |Ω| columns. Matrix A = (aiω) contains a row for each order, and aiω = 1, if and only if column ω represents a feasible assignment of order i ∈ N. In matrix B = (bmtω), m ∈ M; t ∈ T, a single row corresponds to a particular (machine, time) position. If machine m is in use at time t in assignment ω, then bmtω = 1 and 0 otherwise. In matrix C = (ctωk), a single row corresponds to a time interval, each element of which represents the number of the additional resource assigned in time interval t to process the order which column ω represents. Each column ω is associated with a profit Rω which is defined as the net revenue of the respective assignment of a specific order to a machine. The profit Rω of column ω can be easily computed based on its corresponding order index and the (machine, time) position. The decision variable is denoted by xω which equals 1 if column ω is used in the solution, and 0, otherwise.Therefore, the set partitioning formulation is mathematically represented as follows:(10)max∑ω∈ΩRωxω,(11)s.t.∑ω∈Ωaiωxω≤1,∀i∈N,(12)∑ω∈Ωbmtωxω≤1,∀t∈T,m∈M,(13)∑ω∈Ωctωkxω≤Qmaxk,∀t∈T,k∈K,(14)xω∈0,1,∀ω∈Ω.The objective function (10) is to maximize the sum of revenue of selected assignments. Constraint (11) indicates that any order can be accepted atmost once. Constraint (12) limits that a given (machine, time) position can only process atmost one order. Constraint (13) ensures that the resources are not overutilized at any time. Constraint (14) defines the domains of variables.In practice, the number of feasible assignments (columns) grows exponentially with the problem size. Therefore, it is computationally intractable to completely enumerate all feasible columns, making the explicit resolution of the set partitioning formulation impractical for most real-world applications. As an alternative, we solve the set partitioning formulation problem in a column generation framework.The underlying idea of column generation is not to explicitly list all of the columns of the set partitioning formulation (10)–(14), but rather to dynamically generate promising columns. To be more specific, column generation starts with a “simple” formulation including a few columns ω¯∈Ω and relaxes the integrality restrictions on the variables, i.e., xω¯≤0, which is referred to as a restricted linear master problem (RLMP). Subproblems with respect to current dual variables of the RLMP are solved to identify favorable new columns which can be incorporated into the RLMP to improve the objective. The RLMP and subproblems are solved iteratively until no promising columns can be further identified.


### 4.3. The Initial Solution to the RLMP

As mentioned before, column generation starts with only a small subset of columns and then generates necessary columns by an iterative procedure between the restricted master problem and subproblems. Therefore, we should initialize the RLMP. In this section, we design a constructive heuristic to obtain the initial solution for the RLMP. The heuristic works as follows. In the first step, the most urgent order is selected (i.e., the one with the earliest release time). Then, in the second step, it is tentatively inserted into the end of each machine’s current partial schedule. The order is rejected if its every possible insertion incurs a decrease in the net revenue. Otherwise, the machine which leads to the largest positive growth in the net revenue is selected to process the order, and the order is placed at the end of this machine’s partial schedule. These two steps iterate until all orders have been checked. The constructive heuristic can always yield a feasible scheduling plan. Afterwards, the scheduling plan is transformed into a series of assignments of orders to machines, which are mathematically represented by columns and construct our initial RLMP.


### 4.4. Pricing Subproblems and the DE Algorithm

#### 4.4.1. Mathematical Formulation of the Pricing Subproblem

We solve the pricing subproblems to identify promising columns with positive reduced cost, whose inclusion in the restricted linear master problem will improve the objective value. Because machines are assumed to be unrelated, there are |M| pricing subproblems in total, one for each machine. That is, for each machine m, we try to find an assignment ω which has a positive reduced cost calculated by C¯mω=Rω−∑i∈Naiωβi−∑t∈Tbmtωσmt−∑t∈Tctωγt where βi, σmt, and γkt denote the optimal dual solutions associated with constraints (11)–(13), respectively, in the current iteration of the restricted master problem. The net revenue Rω can also be represented in terms of the variables defined in the original integer programming model (1)–(9), i.e., Rω=∑i∈Nei∗Xi−wi∗Ti. Similarly, we have aiω=Xi, bmtω=∑i∈N∑t¯∈maxt−pim+1,0,tSit¯m, and ctωk=∑i∈N∑t¯∈maxt−pim+1,0,tqik∗Sit¯m, where Sitm is equal to 1 if order i starts its processing at time t on machine m and 0 otherwise. Therefore, the pricing subproblem for machine m can be mathematically formulated as follows:(15)max∑i∈Nei∗Xi−wi∗Ti−∑i∈NXiβi−∑t∈T∑i∈N∑t¯∈maxt−pim+1,0,tSit¯mσmt−∑t∈T∑k∈K∑i∈N∑t¯∈maxt−pim+1,0,tqik∗Sit¯mγkt,(16)s.t.∑i∈NXi≤1,∀i∈N,(17)∑t∈TSitm=Xi,∀i∈N,(18)Sitm=0,∀i∈N,t<r,(19)∑t∈Tt+pim∗Sitm=Ci,∀i∈N,(20)Ti≥Ci−di,∀i∈N,(21)Ci,Ti≥0,∀i∈N,(22)Xi,Sitm∈0,1,∀i∈N,t∈T.In this model, the objective function (15) maximizes the reduced costs. Constraint (16) ensures that a column is associated with atmost one order. Constraint (17) ensures that an order can start its processing at time t on the given machine m only if the order is assigned to machine m. Constraint (18) guarantees that no order is processed before its release time. Constraint (19) computes the completion time of each order. Constraint (20) calculates the tardiness of each order. Constraints (21) and (22) define the regions of decision variables.Ideally, the solution approach for solving the subproblems should be efficient because it has to be executed many times during the column generation procedure. Undoubtedly, the subproblems of large-sized instances will become computationally challenging for generic state-of-the-art solvers, such as CPLEX. Alternatively, we next introduce a DE algorithm for the pricing subproblems, which can solve the subproblems more efficiently.


#### 4.4.2. DE Algorithm for Pricing Subproblems

Algorithm 2: The procedure of the DE algorithm.
Step 1: initialize a population ofNP individualsStep 2: evaluate the fitness values of the NP individualsStep 3:While the preset stopping condition is not satisfied DOStep 3.1: Mutation—create a mutant population using a mutation operatorStep 3.2: Crossover—recombine the parent and the mutant to generate a trial populationStep 3.3: Selection—compare the parent and trial individuals and select the better ones for the next generationEndStep 4: output the best solutionThe DE algorithm is a simple yet effective metaheuristic algorithm which was first introduced by Storn and Price [23]. Similar to conventional evolutionary algorithms, DE maintains and evolves a population. Specifically, DE starts with a randomly initialized population and then iteratively proceeds from the current population to a new better-performing one via three evolutionary operators, including mutation, crossover, and selection operators. Mutation operator is responsible for generating new mutant individuals, while crossover operator is responsible for exchanging some subcomponents between parent and mutant individuals to produce new trial individuals. Based on a fitness metric, the selection operator compares the parent individuals with the trial ones and the better ones are admitted to the next generation. The three operators iterate until a termination condition is satisfied. The basic evolution process of the DE algorithm is shown in Algorithm 2.Solution Representation. In order to apply the DE algorithm to solve the subproblems, we need to develop an appropriate encoding approach to represent a feasible solution of a given subproblem. As stated in Section 4.3, the purpose of solving the subproblems is to identify new columns with positive reduced cost, and each column corresponds to a feasible assignment of a single order to a machine in accordance with the time-related constraints. In view of this, we represent each individual of the population by a vector, whose length is 1 + |N|. The first element represents the candidate order, which takes value within [1, |N| ]. The remaining |N| elements corresponds to the starting times of |N| orders, each of which is sampled from its release time to the latest starting time. To better illustrate the solution representation scheme, a simple example for the case mentioned in Section 3.2 is presented in Figure 3.Figure 3
An example for solution representation.In this example, the candidate order is order 3, and its starting time is 4. Suppose that the dual variable is given and we are solving the subproblem for machine 1, we can easily obtain a feasible assignment of order 3 to machine 1 and then compute its corresponding reduced cost (the objective function of this individual).Initialization. The DE population maintains a population of NP vectors. The initial population is randomly generated in the feasible solution domain. For example, the i-th individual can be initialized as follows: (23)Xi0=LB+rand∗UB−LB,where rand is a random number that is sampled from a uniform distribution within (0,1). LB and UB is the lower bound vector and upper bound vector of all elements, respectively.Mutant Operator. At each generation, the mutation operator creates a new mutant individual Vit with respect to each individual (target individual Xr1t) of the current population by the mutation operator. There are various mutation strategies. We adopt the DE/rand/1 scheme, which is expressed as follows: (24)Vit=Xr1t+F∗Xr2t−Xr3t,where Xr1t, Xr2t, and Xr3t are mutually different individuals randomly selected from the current population. The scaling factor F controls the mutation scale and is generally distributed within [0, 1].If the components of the obtained mutant individual exceed the corresponding lower or upper bounds, they can be reinitialized randomly and uniformly within the given range. SinceF is a fractional number, Vit may also contain fractional elements. Therefore, if it happens, each fractional element in Vit is rounded to its nearest integer.Crossover Operator. After mutation, in order to increase the diversity of the population, each pair of target individual Xit and its corresponding mutant individual Vit is crossed to produce a new trial individual Uit. In the DE algorithm, the most commonly used crossover strategy is the binomial crossover scheme, which is defined as follows: (25)Uij,t=Vij,t,ifrandij≤CRorj=jrand,Xij,t,otherwise,j=1,2,…,1+N,where rand is a function returning a randomly generated number within [0, 1]. CR is the crossover rate which is distributed within [0, 1] and mainly controls how many elements of the trial individual originates from the mutant individual. jrand is an integer number randomly selected from {1, 2, 3, ..., 1 + |N|} to enforce that the trial individual differs from its corresponding parent individual in at least one element.Selection Operator. Finally, to keep a constant population size in the following generations, the selection operator is executed to decide whether the trial or the parent individual survives to the next generation. In this paper, we employ an one-to-one selection operator. To be specific, each trial individual Uit is evaluated, and then its objective function value (reduced cost) fUit is compared with the objective function value of the corresponding parent individual fXit in the current population. If the objective function value of the trial individual is greater than or equal to that of the corresponding parent individual, then the parent individual is replaced by the trial one. Otherwise, the parent individual is preserved for the next generation. The selection operator can be described as follows. (26)Xi,t+1=Ui,t,iffUi,t≥fXi,t,Xi,t,otherwise.Stopping Criteria of the DE Algorithm. In the proposed DE algorithm, two stopping conditions are adopted. If any of these stopping conditions is met, the DE algorithm terminates. These two stopping conditions are as follows:(1)
Maximum number of iterations:Itermax iterations(2)
Maximum number of nonimprovement iterations:Itermax−non iterationsAfter termination of the DE algorithm, all columns with positive reduced cost that have been found in the final DE population are added to the RLMP.


### 4.5. Stopping Condition of the Column Generation Process

The column generation process which is an iterative procedure between the restricted master problem and subproblems aims to identify columns with positive reduced cost. Therefore, we terminate the column generation process when it fails to find columns with positive reduced cost for a specified number of consecutive iterations.


### 4.6. Finding an Integer Solution

Column generation is used to solve the LP relaxation of the restricted master problem and thus may return a fractional solution. Therefore, when column generation reports a fractional solution, we adopt a primal heuristic to obtain a feasible integer solution based on the columns already enumerated in the RLMP. The heuristic simply imposes integrality to the master variables, converts the RLMP model to an IP model, and then calls an IP solver to solve the IP model. Since we are solving the IP model on a subset of columns, this process does not necessarily obtain an optimal integer solution to the considered problem. However, it has been shown to be useful in practice, especially for the problem we address in this paper, as we will validate in Section5.3. This heuristic has also been widely used in literature related to column generation-based heuristics (Hauge et al. [24]; Song et al. [25]; Liang et al. [26]).


## 4.1. An Overview of the Proposed Column Generation Algorithm

Column generation is an efficient approach for solving linear programming problems with a large number of variables (columns). Its main advantage is that not all columns need to be enumerated. Instead, the problem is first formulated as a restricted linear master problem (RLMP). This RLMP considers only a subset of columns. And then, one or more subproblems with respect to current dual solution to the RLMP are solved to identify favorable new columns which can be incorporated into the RLMP to improve the objective value. The RLMP and subproblems are solved iteratively until no promising columns can be further added to the RLMP. This termination condition implies that the RLMP is optimal. The RLMP probably returns a fractional solution. In this paper, in order to obtain an integer solution, the RLMP with all columns which have ever been generated as decision variables is converted to an IP model and then solved by an IP solver.In an attempt to enhance the performance of the column generation, we also propose a differential evolution (DE) algorithm to solve the subproblems efficiently. Besides, we develop a constructive heuristic to obtain the initial solution for the RLMP. The procedures of the proposed column generation approach are shown in Algorithm1. Next, we present a detailed description of each step.Algorithm 1: The procedure of the proposed column generation.
Step 1: implement the proposed constructive heuristic to obtain a feasible solutionStep 2: initialize the RLMP with the solution obtained by the constructive heuristicStep3: set the number of consecutive iterations of failing to find a column with positive reduced costIternon to be 0, i.e., Iternon=0Step 4: whileIternon is less than a specified number IternonCG, doStep 4.1: solve the RLMP and get the dual solutionStep 4.2: define the subproblems by the current dual solution and solve them by the proposed DE algorithmStep 4.3: If there exist columns with positive reduced cost, doAdd new columns to the RLMP;Iternon=0ElseIternon=Iternon+1EndEndStep 5: convert the RLMP model to an IP modelStep 6: solve the IP model by CPLEXStep 7: output the solution


## 4.2. Master Problem: A Set Partitioning Formulation

In this section, we reformulate the considered problem as a set partitioning formulation which serves as the master problem in the column generation framework. The set partitioning formulation which is a special class of linear integer programs has been utilized to formulate a wide range of practical problems in fields such as crew scheduling, vehicle routing, and timetable scheduling. Based on the special structure of the set partitioning formulation, efficient approaches techniques, for example, column generation, have been proposed to solve large-sized problems, making it possible to deal with many real-life applications. Therefore, we attempt to propose an alternative mathematical formulation based on the set partitioning formulation.In the set partitioning formulation, there is a set of columns (assignments) denoted by Ω that composes the set partitioning formulation matrix. Each column corresponds to a feasible assignment of a single order to a machine in accordance with the time-related constraints. This matrix contains three submatrices, denoted asA, B, and C, all containing |Ω| columns. Matrix A = (aiω) contains a row for each order, and aiω = 1, if and only if column ω represents a feasible assignment of order i ∈ N. In matrix B = (bmtω), m ∈ M; t ∈ T, a single row corresponds to a particular (machine, time) position. If machine m is in use at time t in assignment ω, then bmtω = 1 and 0 otherwise. In matrix C = (ctωk), a single row corresponds to a time interval, each element of which represents the number of the additional resource assigned in time interval t to process the order which column ω represents. Each column ω is associated with a profit Rω which is defined as the net revenue of the respective assignment of a specific order to a machine. The profit Rω of column ω can be easily computed based on its corresponding order index and the (machine, time) position. The decision variable is denoted by xω which equals 1 if column ω is used in the solution, and 0, otherwise.Therefore, the set partitioning formulation is mathematically represented as follows:(10)max∑ω∈ΩRωxω,(11)s.t.∑ω∈Ωaiωxω≤1,∀i∈N,(12)∑ω∈Ωbmtωxω≤1,∀t∈T,m∈M,(13)∑ω∈Ωctωkxω≤Qmaxk,∀t∈T,k∈K,(14)xω∈0,1,∀ω∈Ω.The objective function (10) is to maximize the sum of revenue of selected assignments. Constraint (11) indicates that any order can be accepted atmost once. Constraint (12) limits that a given (machine, time) position can only process atmost one order. Constraint (13) ensures that the resources are not overutilized at any time. Constraint (14) defines the domains of variables.In practice, the number of feasible assignments (columns) grows exponentially with the problem size. Therefore, it is computationally intractable to completely enumerate all feasible columns, making the explicit resolution of the set partitioning formulation impractical for most real-world applications. As an alternative, we solve the set partitioning formulation problem in a column generation framework.The underlying idea of column generation is not to explicitly list all of the columns of the set partitioning formulation (10)–(14), but rather to dynamically generate promising columns. To be more specific, column generation starts with a “simple” formulation including a few columns ω¯∈Ω and relaxes the integrality restrictions on the variables, i.e., xω¯≤0, which is referred to as a restricted linear master problem (RLMP). Subproblems with respect to current dual variables of the RLMP are solved to identify favorable new columns which can be incorporated into the RLMP to improve the objective. The RLMP and subproblems are solved iteratively until no promising columns can be further identified.


## 4.3. The Initial Solution to the RLMP

As mentioned before, column generation starts with only a small subset of columns and then generates necessary columns by an iterative procedure between the restricted master problem and subproblems. Therefore, we should initialize the RLMP. In this section, we design a constructive heuristic to obtain the initial solution for the RLMP. The heuristic works as follows. In the first step, the most urgent order is selected (i.e., the one with the earliest release time). Then, in the second step, it is tentatively inserted into the end of each machine’s current partial schedule. The order is rejected if its every possible insertion incurs a decrease in the net revenue. Otherwise, the machine which leads to the largest positive growth in the net revenue is selected to process the order, and the order is placed at the end of this machine’s partial schedule. These two steps iterate until all orders have been checked. The constructive heuristic can always yield a feasible scheduling plan. Afterwards, the scheduling plan is transformed into a series of assignments of orders to machines, which are mathematically represented by columns and construct our initial RLMP.


## 4.4. Pricing Subproblems and the DE Algorithm

### 4.4.1. Mathematical Formulation of the Pricing Subproblem

We solve the pricing subproblems to identify promising columns with positive reduced cost, whose inclusion in the restricted linear master problem will improve the objective value. Because machines are assumed to be unrelated, there are |M| pricing subproblems in total, one for each machine. That is, for each machine m, we try to find an assignment ω which has a positive reduced cost calculated by C¯mω=Rω−∑i∈Naiωβi−∑t∈Tbmtωσmt−∑t∈Tctωγt where βi, σmt, and γkt denote the optimal dual solutions associated with constraints (11)–(13), respectively, in the current iteration of the restricted master problem. The net revenue Rω can also be represented in terms of the variables defined in the original integer programming model (1)–(9), i.e., Rω=∑i∈Nei∗Xi−wi∗Ti. Similarly, we have aiω=Xi, bmtω=∑i∈N∑t¯∈maxt−pim+1,0,tSit¯m, and ctωk=∑i∈N∑t¯∈maxt−pim+1,0,tqik∗Sit¯m, where Sitm is equal to 1 if order i starts its processing at time t on machine m and 0 otherwise. Therefore, the pricing subproblem for machine m can be mathematically formulated as follows:(15)max∑i∈Nei∗Xi−wi∗Ti−∑i∈NXiβi−∑t∈T∑i∈N∑t¯∈maxt−pim+1,0,tSit¯mσmt−∑t∈T∑k∈K∑i∈N∑t¯∈maxt−pim+1,0,tqik∗Sit¯mγkt,(16)s.t.∑i∈NXi≤1,∀i∈N,(17)∑t∈TSitm=Xi,∀i∈N,(18)Sitm=0,∀i∈N,t<r,(19)∑t∈Tt+pim∗Sitm=Ci,∀i∈N,(20)Ti≥Ci−di,∀i∈N,(21)Ci,Ti≥0,∀i∈N,(22)Xi,Sitm∈0,1,∀i∈N,t∈T.In this model, the objective function (15) maximizes the reduced costs. Constraint (16) ensures that a column is associated with atmost one order. Constraint (17) ensures that an order can start its processing at time t on the given machine m only if the order is assigned to machine m. Constraint (18) guarantees that no order is processed before its release time. Constraint (19) computes the completion time of each order. Constraint (20) calculates the tardiness of each order. Constraints (21) and (22) define the regions of decision variables.Ideally, the solution approach for solving the subproblems should be efficient because it has to be executed many times during the column generation procedure. Undoubtedly, the subproblems of large-sized instances will become computationally challenging for generic state-of-the-art solvers, such as CPLEX. Alternatively, we next introduce a DE algorithm for the pricing subproblems, which can solve the subproblems more efficiently.


### 4.4.2. DE Algorithm for Pricing Subproblems

Algorithm 2: The procedure of the DE algorithm.
Step 1: initialize a population ofNP individualsStep 2: evaluate the fitness values of the NP individualsStep 3:While the preset stopping condition is not satisfied DOStep 3.1: Mutation—create a mutant population using a mutation operatorStep 3.2: Crossover—recombine the parent and the mutant to generate a trial populationStep 3.3: Selection—compare the parent and trial individuals and select the better ones for the next generationEndStep 4: output the best solutionThe DE algorithm is a simple yet effective metaheuristic algorithm which was first introduced by Storn and Price [23]. Similar to conventional evolutionary algorithms, DE maintains and evolves a population. Specifically, DE starts with a randomly initialized population and then iteratively proceeds from the current population to a new better-performing one via three evolutionary operators, including mutation, crossover, and selection operators. Mutation operator is responsible for generating new mutant individuals, while crossover operator is responsible for exchanging some subcomponents between parent and mutant individuals to produce new trial individuals. Based on a fitness metric, the selection operator compares the parent individuals with the trial ones and the better ones are admitted to the next generation. The three operators iterate until a termination condition is satisfied. The basic evolution process of the DE algorithm is shown in Algorithm 2.Solution Representation. In order to apply the DE algorithm to solve the subproblems, we need to develop an appropriate encoding approach to represent a feasible solution of a given subproblem. As stated in Section 4.3, the purpose of solving the subproblems is to identify new columns with positive reduced cost, and each column corresponds to a feasible assignment of a single order to a machine in accordance with the time-related constraints. In view of this, we represent each individual of the population by a vector, whose length is 1 + |N|. The first element represents the candidate order, which takes value within [1, |N| ]. The remaining |N| elements corresponds to the starting times of |N| orders, each of which is sampled from its release time to the latest starting time. To better illustrate the solution representation scheme, a simple example for the case mentioned in Section 3.2 is presented in Figure 3.Figure 3
An example for solution representation.In this example, the candidate order is order 3, and its starting time is 4. Suppose that the dual variable is given and we are solving the subproblem for machine 1, we can easily obtain a feasible assignment of order 3 to machine 1 and then compute its corresponding reduced cost (the objective function of this individual).Initialization. The DE population maintains a population of NP vectors. The initial population is randomly generated in the feasible solution domain. For example, the i-th individual can be initialized as follows: (23)Xi0=LB+rand∗UB−LB,where rand is a random number that is sampled from a uniform distribution within (0,1). LB and UB is the lower bound vector and upper bound vector of all elements, respectively.Mutant Operator. At each generation, the mutation operator creates a new mutant individual Vit with respect to each individual (target individual Xr1t) of the current population by the mutation operator. There are various mutation strategies. We adopt the DE/rand/1 scheme, which is expressed as follows: (24)Vit=Xr1t+F∗Xr2t−Xr3t,where Xr1t, Xr2t, and Xr3t are mutually different individuals randomly selected from the current population. The scaling factor F controls the mutation scale and is generally distributed within [0, 1].If the components of the obtained mutant individual exceed the corresponding lower or upper bounds, they can be reinitialized randomly and uniformly within the given range. SinceF is a fractional number, Vit may also contain fractional elements. Therefore, if it happens, each fractional element in Vit is rounded to its nearest integer.Crossover Operator. After mutation, in order to increase the diversity of the population, each pair of target individual Xit and its corresponding mutant individual Vit is crossed to produce a new trial individual Uit. In the DE algorithm, the most commonly used crossover strategy is the binomial crossover scheme, which is defined as follows: (25)Uij,t=Vij,t,ifrandij≤CRorj=jrand,Xij,t,otherwise,j=1,2,…,1+N,where rand is a function returning a randomly generated number within [0, 1]. CR is the crossover rate which is distributed within [0, 1] and mainly controls how many elements of the trial individual originates from the mutant individual. jrand is an integer number randomly selected from {1, 2, 3, ..., 1 + |N|} to enforce that the trial individual differs from its corresponding parent individual in at least one element.Selection Operator. Finally, to keep a constant population size in the following generations, the selection operator is executed to decide whether the trial or the parent individual survives to the next generation. In this paper, we employ an one-to-one selection operator. To be specific, each trial individual Uit is evaluated, and then its objective function value (reduced cost) fUit is compared with the objective function value of the corresponding parent individual fXit in the current population. If the objective function value of the trial individual is greater than or equal to that of the corresponding parent individual, then the parent individual is replaced by the trial one. Otherwise, the parent individual is preserved for the next generation. The selection operator can be described as follows. (26)Xi,t+1=Ui,t,iffUi,t≥fXi,t,Xi,t,otherwise.Stopping Criteria of the DE Algorithm. In the proposed DE algorithm, two stopping conditions are adopted. If any of these stopping conditions is met, the DE algorithm terminates. These two stopping conditions are as follows:(1)
Maximum number of iterations:Itermax iterations(2)
Maximum number of nonimprovement iterations:Itermax−non iterationsAfter termination of the DE algorithm, all columns with positive reduced cost that have been found in the final DE population are added to the RLMP.


## 4.4.1. Mathematical Formulation of the Pricing Subproblem

We solve the pricing subproblems to identify promising columns with positive reduced cost, whose inclusion in the restricted linear master problem will improve the objective value. Because machines are assumed to be unrelated, there are |M| pricing subproblems in total, one for each machine. That is, for each machine m, we try to find an assignment ω which has a positive reduced cost calculated by C¯mω=Rω−∑i∈Naiωβi−∑t∈Tbmtωσmt−∑t∈Tctωγt where βi, σmt, and γkt denote the optimal dual solutions associated with constraints (11)–(13), respectively, in the current iteration of the restricted master problem. The net revenue Rω can also be represented in terms of the variables defined in the original integer programming model (1)–(9), i.e., Rω=∑i∈Nei∗Xi−wi∗Ti. Similarly, we have aiω=Xi, bmtω=∑i∈N∑t¯∈maxt−pim+1,0,tSit¯m, and ctωk=∑i∈N∑t¯∈maxt−pim+1,0,tqik∗Sit¯m, where Sitm is equal to 1 if order i starts its processing at time t on machine m and 0 otherwise. Therefore, the pricing subproblem for machine m can be mathematically formulated as follows:(15)max∑i∈Nei∗Xi−wi∗Ti−∑i∈NXiβi−∑t∈T∑i∈N∑t¯∈maxt−pim+1,0,tSit¯mσmt−∑t∈T∑k∈K∑i∈N∑t¯∈maxt−pim+1,0,tqik∗Sit¯mγkt,(16)s.t.∑i∈NXi≤1,∀i∈N,(17)∑t∈TSitm=Xi,∀i∈N,(18)Sitm=0,∀i∈N,t<r,(19)∑t∈Tt+pim∗Sitm=Ci,∀i∈N,(20)Ti≥Ci−di,∀i∈N,(21)Ci,Ti≥0,∀i∈N,(22)Xi,Sitm∈0,1,∀i∈N,t∈T.In this model, the objective function (15) maximizes the reduced costs. Constraint (16) ensures that a column is associated with atmost one order. Constraint (17) ensures that an order can start its processing at time t on the given machine m only if the order is assigned to machine m. Constraint (18) guarantees that no order is processed before its release time. Constraint (19) computes the completion time of each order. Constraint (20) calculates the tardiness of each order. Constraints (21) and (22) define the regions of decision variables.Ideally, the solution approach for solving the subproblems should be efficient because it has to be executed many times during the column generation procedure. Undoubtedly, the subproblems of large-sized instances will become computationally challenging for generic state-of-the-art solvers, such as CPLEX. Alternatively, we next introduce a DE algorithm for the pricing subproblems, which can solve the subproblems more efficiently.


## 4.4.2. DE Algorithm for Pricing Subproblems

Algorithm 2: The procedure of the DE algorithm.
Step 1: initialize a population ofNP individualsStep 2: evaluate the fitness values of the NP individualsStep 3:While the preset stopping condition is not satisfied DOStep 3.1: Mutation—create a mutant population using a mutation operatorStep 3.2: Crossover—recombine the parent and the mutant to generate a trial populationStep 3.3: Selection—compare the parent and trial individuals and select the better ones for the next generationEndStep 4: output the best solutionThe DE algorithm is a simple yet effective metaheuristic algorithm which was first introduced by Storn and Price [23]. Similar to conventional evolutionary algorithms, DE maintains and evolves a population. Specifically, DE starts with a randomly initialized population and then iteratively proceeds from the current population to a new better-performing one via three evolutionary operators, including mutation, crossover, and selection operators. Mutation operator is responsible for generating new mutant individuals, while crossover operator is responsible for exchanging some subcomponents between parent and mutant individuals to produce new trial individuals. Based on a fitness metric, the selection operator compares the parent individuals with the trial ones and the better ones are admitted to the next generation. The three operators iterate until a termination condition is satisfied. The basic evolution process of the DE algorithm is shown in Algorithm 2.Solution Representation. In order to apply the DE algorithm to solve the subproblems, we need to develop an appropriate encoding approach to represent a feasible solution of a given subproblem. As stated in Section 4.3, the purpose of solving the subproblems is to identify new columns with positive reduced cost, and each column corresponds to a feasible assignment of a single order to a machine in accordance with the time-related constraints. In view of this, we represent each individual of the population by a vector, whose length is 1 + |N|. The first element represents the candidate order, which takes value within [1, |N| ]. The remaining |N| elements corresponds to the starting times of |N| orders, each of which is sampled from its release time to the latest starting time. To better illustrate the solution representation scheme, a simple example for the case mentioned in Section 3.2 is presented in Figure 3.Figure 3
An example for solution representation.In this example, the candidate order is order 3, and its starting time is 4. Suppose that the dual variable is given and we are solving the subproblem for machine 1, we can easily obtain a feasible assignment of order 3 to machine 1 and then compute its corresponding reduced cost (the objective function of this individual).Initialization. The DE population maintains a population of NP vectors. The initial population is randomly generated in the feasible solution domain. For example, the i-th individual can be initialized as follows: (23)Xi0=LB+rand∗UB−LB,where rand is a random number that is sampled from a uniform distribution within (0,1). LB and UB is the lower bound vector and upper bound vector of all elements, respectively.Mutant Operator. At each generation, the mutation operator creates a new mutant individual Vit with respect to each individual (target individual Xr1t) of the current population by the mutation operator. There are various mutation strategies. We adopt the DE/rand/1 scheme, which is expressed as follows: (24)Vit=Xr1t+F∗Xr2t−Xr3t,where Xr1t, Xr2t, and Xr3t are mutually different individuals randomly selected from the current population. The scaling factor F controls the mutation scale and is generally distributed within [0, 1].If the components of the obtained mutant individual exceed the corresponding lower or upper bounds, they can be reinitialized randomly and uniformly within the given range. SinceF is a fractional number, Vit may also contain fractional elements. Therefore, if it happens, each fractional element in Vit is rounded to its nearest integer.Crossover Operator. After mutation, in order to increase the diversity of the population, each pair of target individual Xit and its corresponding mutant individual Vit is crossed to produce a new trial individual Uit. In the DE algorithm, the most commonly used crossover strategy is the binomial crossover scheme, which is defined as follows: (25)Uij,t=Vij,t,ifrandij≤CRorj=jrand,Xij,t,otherwise,j=1,2,…,1+N,where rand is a function returning a randomly generated number within [0, 1]. CR is the crossover rate which is distributed within [0, 1] and mainly controls how many elements of the trial individual originates from the mutant individual. jrand is an integer number randomly selected from {1, 2, 3, ..., 1 + |N|} to enforce that the trial individual differs from its corresponding parent individual in at least one element.Selection Operator. Finally, to keep a constant population size in the following generations, the selection operator is executed to decide whether the trial or the parent individual survives to the next generation. In this paper, we employ an one-to-one selection operator. To be specific, each trial individual Uit is evaluated, and then its objective function value (reduced cost) fUit is compared with the objective function value of the corresponding parent individual fXit in the current population. If the objective function value of the trial individual is greater than or equal to that of the corresponding parent individual, then the parent individual is replaced by the trial one. Otherwise, the parent individual is preserved for the next generation. The selection operator can be described as follows. (26)Xi,t+1=Ui,t,iffUi,t≥fXi,t,Xi,t,otherwise.Stopping Criteria of the DE Algorithm. In the proposed DE algorithm, two stopping conditions are adopted. If any of these stopping conditions is met, the DE algorithm terminates. These two stopping conditions are as follows:(1)
Maximum number of iterations:Itermax iterations(2)
Maximum number of nonimprovement iterations:Itermax−non iterationsAfter termination of the DE algorithm, all columns with positive reduced cost that have been found in the final DE population are added to the RLMP.


## 4.5. Stopping Condition of the Column Generation Process

The column generation process which is an iterative procedure between the restricted master problem and subproblems aims to identify columns with positive reduced cost. Therefore, we terminate the column generation process when it fails to find columns with positive reduced cost for a specified number of consecutive iterations.


## 4.6. Finding an Integer Solution

Column generation is used to solve the LP relaxation of the restricted master problem and thus may return a fractional solution. Therefore, when column generation reports a fractional solution, we adopt a primal heuristic to obtain a feasible integer solution based on the columns already enumerated in the RLMP. The heuristic simply imposes integrality to the master variables, converts the RLMP model to an IP model, and then calls an IP solver to solve the IP model. Since we are solving the IP model on a subset of columns, this process does not necessarily obtain an optimal integer solution to the considered problem. However, it has been shown to be useful in practice, especially for the problem we address in this paper, as we will validate in Section5.3. This heuristic has also been widely used in literature related to column generation-based heuristics (Hauge et al. [24]; Song et al. [25]; Liang et al. [26]).


## 5. Numerical Experiments

In this section, we conduct extensive computational experiments to verify the applicability and effectiveness of column generation based on randomly generated instances. Column generation is implemented in MATLAB 2016a and compared against the original MILP formulation solved by CPLEX 12.10. All runs are performed on a Windows 10 professional 64 bit operating system with an Intel Core i5-8250U CPU with 1.60 GHz and 8 GB RAM laptop. For each run, we enforce a time limit of 3600 s of CPU time.


### 5.1. Data Generation

There are no publicly available benchmark instances reflecting the setting of the studied problem. Therefore, we randomly generate instances in varying parameter and problem size via the generation schemes of Cesaret et al. [27] and Vallada et al. [19] with some adaptations to the characteristics of our problem. We generate three different sets: small-sized with |N| = {15, 20} and |M| = {2, 3}, medium-sized with |N| = {40, 60} and |M| = {2, 3}, and large-sized with |N| = {90, 100} and |M| = {3, 4}.Each |N|∗|M| combination contains 9 subclasses, and each subclass comprises 5 randomly generated instances. Each subclass is characterized by parameters TF and R. More specifically, the parameter TF is the tardiness factor, which controls average due date values and parameter R determines the relative range of due dates. Parameters TF and R all take values from the set {0.2, 0.6, 1.0}, yielding the aforementioned 9 subclasses. Other parameters are generated as follows. Processing times are randomly generated from the interval [1, 100], i.e., pim ∼ U[1, 100], i∈N, m∈M. The due dates of jobs di are generated from the formula di=ri+maxslack,pim. slack is drawn from a discrete uniform distribution in the interval p¯1−TF−R/2,p¯1−TF+R/2, where p¯=∑i∈N∑m∈Mpim/M2, similar to the study by Cesaret et al. [27] and Wang and Ye [9]. The release time ri is generated from the interval [0, TF∗p¯]. The number of resource types is set to be 2. For each type of resources, the resource consumption of each job is randomly generated in the interval [1, 9]. And, the total number of each type of resources is set to be 6∗|M|. Penalties and revenues are randomly generated from a discrete uniform distribution in the interval [1, 10].


### 5.2. Parameter Settings and Performance Metrics

#### 5.2.1. Parameter Settings

The parameter settings affect the quality of the solutions obtained by the proposed algorithm. In this paper, the Taguchi method is applied to tune the parameters. The Taguchi method which is essentially a fractional factorial experiment approach uses orthogonal arrays to drastically reduce the number of experiments to be performed. The Taguchi method supposes two types of factors affect a process, i.e., controllable factors and uncontrollable noise factors. Since removal of the noise factors is impractical, the Taguchi method seeks to determine the optimal setting of controllable factors and to minimize the effect of noise, in order to achieve a robust performance. The responses at each setting of controllable parameters are considered as a measure, indicative of not only the mean of the performance but also its variance. The mean and the variance are transformed to a single performance measure known as the signal-to-noise (S/N) ratio. The term “signal” signifies the desirable value (mean response variable) and “noise” signifies the undesirable value (standard deviation), so the S/N ratio represents the amount of variation presents in the response variable. Thus, the goal is to maximize the S/N ratio.This paper uses the Taguchi method as a reliable method of setting the parameters for the proposed algorithm. The DE algorithm has five parameters: (1) population sizeNP; (2) scaling factor F; (3) crossover rate CR; (4) maximum number of iterations Itermax; and (5) maximum number of nonimprovement iterations Itermax−non. And, the stopping condition of the column generation process involves one parameter IternonCG, i.e., the number of consecutive iterations of failing to find a column with positive reduced cost. Five levels are considered for each parameter, as shown in Table 2. The orthogonal array L25(65) shown in Table 3 is designed as the experimental layout, which displays 25 different combinations of factor levels for the experiments.Table 2
Levels of parameters.
Parameter (factor)NPFCRItermaxItermax−nonIternonCGLevel1100.20.251001052150.40.5020025103200.60.7530050204250.80.9540075305301150010040Table 3
The orthogonal array.
Experiment no.ParameterNPFCRItermaxItermax−nonIternonCG11111112122222313333341444445155555621234572234518234512924512310251234113135241232413513335241143413521535241316414253174253141843142519442531204531422151543222521543235321542454321525554321For each of three different-sized problems, five instances are randomly selected and then tested on the 25 parameter combinations. The obtained computational data are processed by Minitab 19 statistical software. The S/N ratio trends of parameter levels for different-sized instances are displayed in Figures4–6. As indicated in the three figures, we set NP = 20, F = 0.6, CR = 0.75, Itermax = 100, Itermax−non = 50, and IternonCG = 10 for small-sized and medium-sized instances. And, for large-sized instances, we set NP = 30, F = 0.6, CR = 0.75, Itermax = 500, Itermax−non = 100, and IternonCG = 20, respectively.Figure 4
The S/N ratio trend for small-sized problems.Figure 5
The S/N ratio trend for medium-sized problems.Figure 6
The S/N ratio trend for large-sized problems.


#### 5.2.2. Performance Metrics

We use the gap between the revenue and the upper bound as the performance metric. The revenue is obtained from a given solution method, including either the proposed column generation approach or the integer programming formulation solved by CPLEX solver. The upper bound is obtained by solving the MILP formulation in Section3 with a time limit of 3600 s. Furthermore, the number of optimal solutions out of 5 instances of each subclass and the average computation time are recorded as auxiliary metrics.


### 5.3. Computational Results and Analysis

#### 5.3.1. Results of Small-Sized Instances

We first verify the performance of the proposed column generation approach (CG for short) for solving small-sized instances. Tables4 and 5 display the comparison results between the MILP formulation and the proposed column generation approach. The first column represents the |N| × |M| combinations. Column 2 and column 3 show parameter TF and parameter R, respectively. We report in columns 4–9 the minimum (Min), average (Avg) and maximum (Max) gaps between the upper bound and the best found solution obtained by the respective method for every subclass. In columns 10 and 11, we display the number of optimal solutions found by the given solution method. The following two columns with the same heading “Total” report the average computation time consumed by solution approaches. The remaining columns show the time spent on each step of column generation. Specifically, columns “CH,” “LP,” “MILP,” and “DE” report the time spent by the constructive heuristic, LP solver, MILP solver, and the proposed DE algorithm, respectively.Table 4
Computation results of small-sized instances with |N| = 15.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU (s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE15 × 20.20.20000005550.803.740.030.821.501.380.60000005578.683.460.031.330.511.5910000005531.003.810.021.061.201.520.60.20000005574.302.200.030.670.251.250.60000005582.756.170.12.571.212.2910000005521.006.420.023.500.452.4510.2000000559.001.790.020.470.081.220.60000005567.522.940.021.090.051.7810000005511.803.570.021.290.142.1315 × 30.20.20000005515.213.250.031.090.591.540.6000000559.532.040.030.550.231.2310000005510.861.790.030.490.161.110.60.20000005512.321.390.020.310.080.970.60000005511.065.740.032.640.192.8810000005518.621.840.020.480.191.1410.20000005513.681.340.020.240.081.000.60000005510.221.750.020.380.101.2510000005511.851.230.020.220.050.94Average——000000——30.013.030.031.070.391.54Table 5
Computation results of small-sized instances with |N| = 20.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE20 × 20.20.200000055119.143.240.031.080.681.450.60000005568.634.490.031.910.641.9110000005532.273.580.031.300.281.970.60.20000005545.836.110.041.902.251.920.60000005522.234.880.022.460.771.6310000005542.655.050.022.040.462.5310.20000005520.606.550.022.460.054.020.60000005519.867.190.022.740.424.0110000005529.805.830.022.230.323.2620 × 30.20.20000005545.215.240.071.530.543.100.60000005539.035.280.121.740.323.1010000005571.076.270.031.950.613.680.60.20000005548.036.030.031.950.783.270.60000005557.816.750.022.050.514.1710000005545.025.710.021.580.573.5410.20000005536.206.820.022.100.154.550.60000005537.693.190.020.680.112.3810000005542.033.320.130.610.182.40Average——000000——45.735.310.041.800.542.94As indicated by the gap values in Tables4 and 5, both the MILP formulation and column generation solve all instances to optimality. In comparison with the MILP formulation, column generation consumes significantly less running time to find the optimal solutions. That is, the MILP formulation needs computing time increasing from few seconds to hundreds of seconds, whereas column generation is able to find the optimal solution in seconds. For example, the average CPU time reduces from 119.14 s with the MILP formulation to 3.24 s with the proposed algorithm for instances with |N| = 20, |M| = 2, TF = 0.2, and R = 0.2.To sum up, the numerical experiments on small-sized instances show that the proposed algorithm significantly outperforms the MILP formulation in terms of computational efficiency.


#### 5.3.2. Results of Medium-Sized Instances

We further examine the performances of the proposed method and the MILP formulation for solving medium-sized instances. Tables6 and 7 present the experiment results.Table 6
Computation results of medium-sized instances with |N| = 40.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE40 × 20.20.20.002.1310.66000451610.2030.180.113.2125.071.790.600000055134.4710.180.073.554.422.14100000055213.2010.250.123.484.022.630.60.2014.6833.1901.333.91222204.82358.080.142.96347.327.660.603.7917.4500.350.8833447.6787.960.072.5918.7966.521010.3722.0200.661.96132907.4052.630.074.1816.7031.6810.200000055758.6123.040.0311.890.4410.680.600000055182.6516.580.035.291.1310.13100000055110.0014.020.202.940.7810.1040 × 30.20.20000005596.619.870.084.092.912.790.60000005592.255.580.092.670.682.14100000055130.822.380.070.730.181.400.60.200000055107.2010.610.054.303.322.940.600000055239.688.020.043.181.962.8410000005575.2514.000.057.951.784.2210.20000005577.217.160.042.990.323.810.60000005596.6352.560.046.001.2545.2710000005592.203.460.031.280.301.85Average——01.724.6300.130.38——532.0539.810.074.0723.9711.70Table 7
Computation results of medium-sized instances with |N| = 60.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE60 × 20.20.24.6712.9518.48000053600103.260.146.2893.932.910.66.3213.3920.34000053600121.100.256.70110.343.81110.3312.0613.1500005360080.230.275.9570.303.710.60.215.6818.9524.5400.882.11023600180.370.124.9583.4491.860.612.5715.9720.8500.412.07043600144.270.084.41125.4214.3619.8614.3117.8200005360020.020.129.545.095.2710.200000055417.89107.970.0569.936.6731.320.6021.2831.6101.032.5011360059.80.0512.193.9143.6517.0316.5525.4600.521.3203360023.330.0810.405.517.3460 × 30.20.200000055971.8111.610.175.532.313.600.605.648.69000152961.2212.630.176.701.873.89106.6217.39000252287.8610.530.114.772.033.620.60.200000055221.2128.930.055.6819.883.320.6013.7126.03000152926.0710.190.085.141.733.241014.2833.12000252328.0312.190.056.461.793.8910.200000055228.0918.750.029.630.698.410.600000055235.2110.960.035.461.244.23102.7513.7700.261.3344935.258.390.043.830.803.72Average——3.699.3615.0700.170.52——2350.7053.590.1010.2029.8313.45As shown in Tables6 and 7, we can observe that the proposed column generation algorithm is superior to the MILP formulation in terms of efficiency and effectiveness. For instances with |N| = 40, column generation finds 83 optimal solutions out of 90 instances with 0.13% deviation on average. In contrast, the MILP formulation identifies the optimal solution in 80 out of 90 instances and reports 1.72% deviation on average. As expected, the performance of the MILP formulation degenerates since the solution space exponentially grows with the increase in problem size. When |N| equals 60, the MILP formulation only reports 36 instances out of 90 to optimality within 3600 s, whereas column generation finds 79 optimal solutions. The average CPU time of the MILP formulation increases to 2350.70 s, while column generation runs in 53.59 s on an average.Based on above comparison analysis, we can come to the conclusion that column generation still performs significantly better than the MILP formulation with regards to efficiency and effectiveness and is an effective approach for solving medium-sized instances.


#### 5.3.3. Results of Large-Sized Instances

We further report computational results obtained for large-sized instances in Tables8 and 9. These instances become more challenging for the MILP formulation, which can be indicated by computation time. We can observe that column generation performs significantly better than the MILP formulation in terms of number of instances solved to optimality and computation times. Specifically, the average percentage deviation of the MILP formulation is consistently larger than that of column generation. We can also observe that the MILP formulation fails to find a feasible solution to instances with TF = 0.2, R = 1 and |N| = 100. In contrast, the proposed column generation algorithm is able to prove the optimality of 177 instances and to provide very low gaps for the remaining, always below 1%. Column generation can generally solve each instance within 800 s, which further demonstrate its efficiency.Table 8
Computation results of large-sized instances with |N| = 90.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE90 × 30.20.2∗1.608.5611.68000052082.23357.670.21301.6812.6643.120.6∗2.945.186.30000051148.00296.320.18239.8410.7145.591 ∗0.395.8311.97000052100.61218.420.19160.607.1950.440.60.2011.7418.28000152880.21248.820.07192.1219.9036.730.68.3312.8716.63000053600223.360.06177.009.5136.791 ∗7.7314.4820.38000052964.35228.280.07173.027.5747.6210.200.110.5500.060.30441080.89121.080.0463.7710.4846.790.606.6919.25000351848.10107.610.0562.123.7141.731 ∗11.1913.1315.63000052412.80133.220.0477.589.3046.3090 × 40.20.205.8214.11000053600316.040.13254.717.0654.140.61.303.6811.07000053600274.330.14214.287.7652.151∗---00005414.81225.310.12164.484.5456.170.60.206.7615.99000252376.00308.200.05243.179.7355.250.605.0013.70000151762.24271.050.04200.519.4061.101017.5178.60000151139.31189.580.04125.935.6857.9310.200000055302.1267.800.0429.611.5336.620.600000055408.0888.780.0343.623.0842.05104.4911.29000252423.28112.080.0458.964.4348.65Average——1.977.1715.6100.000.02——2007.95210.440.09154.618.0147.73“-”indicates that the CPLEX fails to find feasible solutions for this subclass. “∗” indicates that the CPLEX runs of memory when solving some instances of this subclass.Table 9
Computation results of large-sized instances with |N| = 100.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE100 × 30.20.2∗2.248.0013.9200005571.84374.270.25208.5810.50154.940.6∗5.387.318.6700005627.31247.860.24147.8010.5689.261∗---00005436.57200.050.26129.568.9061.330.60.210.0415.2419.45000053600359.940.1192.67122.7144.460.610.3715.5720.41000053600254.400.07185.3712.2756.691∗14.2415.0516.00000051711.34248.250.07165.4811.8770.8310.2∗00.411.3500.320.86331198.63147.540.0673.9512.4261.110.6∗08.5418.24000152454.21160.960.0988.2112.6060.061∗10.0210.8112.25000051055.10176.530.0599.087.4769.93100 × 40.20.2∗3.975.577.04000051751.80488.650.32389.7012.0086.630.6∗0.824.798.7500005519.20489.110.25376.2414.4198.211∗---00005457.40402.910.35264.6612.91124.990.60.207.4913.13000153000793.430.06619.8624.58148.930.6∗9.5112.1217.90000052977.87562.540.09395.9020.20146.351∗7.547.547.5400005438.96337.990.08142.0021.61174.3010.200000055410.21121.410.0639.766.2775.320.601.286.42000451164.28141.240.0673.006.2261.961∗0.555.629.39000051073.99155.210.0778.846.3469.96Average——4.677.8311.2800.020.05——1502.71314.570.14203.9318.5591.96“-” indicates that the CPLEX fails to find feasible solutions for this subclass. “∗” indicates that the CPLEX runs of memory when solving some instances of this subclass.The encouraging computational results of the proposed algorithm may be attributed to two aspects. Firstly, the column generation approach keeps most columns of the set partitioning formulation implicit and only generates necessary columns that can improve the objective value. Therefore, a considerate number of times can be saved. Secondly, the DE algorithm is an efficient heuristic. It is used to solve the subproblems and helps to identify promising columns in a very short time.In summary, the capability of the proposed algorithm to consistently find optimal or near-optimal solutions confirms its effectiveness and superiority against the MILP formulation. Our experimental results clearly show that column generation is a computationally viable solution approach for practical situations.


### 5.4. Effect of ParametersTF and R on the Performance of the Methods

In this section, we further investigate the impact of the instance generation parameters on the performance of the methods. As stated in Section5.1, parameter TF determines due date tightness while parameter R controls the due date range. Additionally, parameter TF also controls the dispersion degree of the release dates of the jobs, in a way that the higher the TF is, the more dispersed the release dates are.As indicated by the gap and computation time values in Tables4–9, we can notice that given that all other specifications are the same, instances with large values of TF are generally easier to solve than those with small values. A reasonable explanation for this behaviour could be as follows. A small TF value corresponds to a situation where the release dates are distributed in the very beginning of the time horizon and then are not very restrictive. And, a very large TF value leads to a situation where the release dates are very dispersed. Hence, as TF increases, we receive incoming orders released at dispersed time points with tight due dates. Therefore, more orders are expected to be rejected, which in turn implies the reduction of the size of the solution space. Consequently, the instances with large TF values are less challenging. We can observe from the tables that the number of instances solved to optimality by the MILP formulation increases and the average CPU time spent by the MILP formulation tends to decrease as the TF gets large when it comes to medium- and large-sized instances. A similar observation can be found when instances are solved by the proposed column generation algorithm. That is, two subclasses are all exactly solved by the proposed algorithm, but the subclass with TF = 1 evidently consumes less CPU time than the subclass with TF = 0.2.Given theTF value, the due date range increases with the increase of the parameter R. Results in Tables 4–9 show that given the TF value, some subclasses with large R value are more difficult to solve than those with small R value, but we can also find some exceptions. This coincides with observations reported by Wang and Ye [9]. Therefore, there is no clear trend in correlation between the R value and the hardness of the problem.


## 5.1. Data Generation

There are no publicly available benchmark instances reflecting the setting of the studied problem. Therefore, we randomly generate instances in varying parameter and problem size via the generation schemes of Cesaret et al. [27] and Vallada et al. [19] with some adaptations to the characteristics of our problem. We generate three different sets: small-sized with |N| = {15, 20} and |M| = {2, 3}, medium-sized with |N| = {40, 60} and |M| = {2, 3}, and large-sized with |N| = {90, 100} and |M| = {3, 4}.Each |N|∗|M| combination contains 9 subclasses, and each subclass comprises 5 randomly generated instances. Each subclass is characterized by parameters TF and R. More specifically, the parameter TF is the tardiness factor, which controls average due date values and parameter R determines the relative range of due dates. Parameters TF and R all take values from the set {0.2, 0.6, 1.0}, yielding the aforementioned 9 subclasses. Other parameters are generated as follows. Processing times are randomly generated from the interval [1, 100], i.e., pim ∼ U[1, 100], i∈N, m∈M. The due dates of jobs di are generated from the formula di=ri+maxslack,pim. slack is drawn from a discrete uniform distribution in the interval p¯1−TF−R/2,p¯1−TF+R/2, where p¯=∑i∈N∑m∈Mpim/M2, similar to the study by Cesaret et al. [27] and Wang and Ye [9]. The release time ri is generated from the interval [0, TF∗p¯]. The number of resource types is set to be 2. For each type of resources, the resource consumption of each job is randomly generated in the interval [1, 9]. And, the total number of each type of resources is set to be 6∗|M|. Penalties and revenues are randomly generated from a discrete uniform distribution in the interval [1, 10].


## 5.2. Parameter Settings and Performance Metrics

### 5.2.1. Parameter Settings

The parameter settings affect the quality of the solutions obtained by the proposed algorithm. In this paper, the Taguchi method is applied to tune the parameters. The Taguchi method which is essentially a fractional factorial experiment approach uses orthogonal arrays to drastically reduce the number of experiments to be performed. The Taguchi method supposes two types of factors affect a process, i.e., controllable factors and uncontrollable noise factors. Since removal of the noise factors is impractical, the Taguchi method seeks to determine the optimal setting of controllable factors and to minimize the effect of noise, in order to achieve a robust performance. The responses at each setting of controllable parameters are considered as a measure, indicative of not only the mean of the performance but also its variance. The mean and the variance are transformed to a single performance measure known as the signal-to-noise (S/N) ratio. The term “signal” signifies the desirable value (mean response variable) and “noise” signifies the undesirable value (standard deviation), so the S/N ratio represents the amount of variation presents in the response variable. Thus, the goal is to maximize the S/N ratio.This paper uses the Taguchi method as a reliable method of setting the parameters for the proposed algorithm. The DE algorithm has five parameters: (1) population sizeNP; (2) scaling factor F; (3) crossover rate CR; (4) maximum number of iterations Itermax; and (5) maximum number of nonimprovement iterations Itermax−non. And, the stopping condition of the column generation process involves one parameter IternonCG, i.e., the number of consecutive iterations of failing to find a column with positive reduced cost. Five levels are considered for each parameter, as shown in Table 2. The orthogonal array L25(65) shown in Table 3 is designed as the experimental layout, which displays 25 different combinations of factor levels for the experiments.Table 2
Levels of parameters.
Parameter (factor)NPFCRItermaxItermax−nonIternonCGLevel1100.20.251001052150.40.5020025103200.60.7530050204250.80.9540075305301150010040Table 3
The orthogonal array.
Experiment no.ParameterNPFCRItermaxItermax−nonIternonCG11111112122222313333341444445155555621234572234518234512924512310251234113135241232413513335241143413521535241316414253174253141843142519442531204531422151543222521543235321542454321525554321For each of three different-sized problems, five instances are randomly selected and then tested on the 25 parameter combinations. The obtained computational data are processed by Minitab 19 statistical software. The S/N ratio trends of parameter levels for different-sized instances are displayed in Figures4–6. As indicated in the three figures, we set NP = 20, F = 0.6, CR = 0.75, Itermax = 100, Itermax−non = 50, and IternonCG = 10 for small-sized and medium-sized instances. And, for large-sized instances, we set NP = 30, F = 0.6, CR = 0.75, Itermax = 500, Itermax−non = 100, and IternonCG = 20, respectively.Figure 4
The S/N ratio trend for small-sized problems.Figure 5
The S/N ratio trend for medium-sized problems.Figure 6
The S/N ratio trend for large-sized problems.


### 5.2.2. Performance Metrics

We use the gap between the revenue and the upper bound as the performance metric. The revenue is obtained from a given solution method, including either the proposed column generation approach or the integer programming formulation solved by CPLEX solver. The upper bound is obtained by solving the MILP formulation in Section3 with a time limit of 3600 s. Furthermore, the number of optimal solutions out of 5 instances of each subclass and the average computation time are recorded as auxiliary metrics.


## 5.2.1. Parameter Settings

The parameter settings affect the quality of the solutions obtained by the proposed algorithm. In this paper, the Taguchi method is applied to tune the parameters. The Taguchi method which is essentially a fractional factorial experiment approach uses orthogonal arrays to drastically reduce the number of experiments to be performed. The Taguchi method supposes two types of factors affect a process, i.e., controllable factors and uncontrollable noise factors. Since removal of the noise factors is impractical, the Taguchi method seeks to determine the optimal setting of controllable factors and to minimize the effect of noise, in order to achieve a robust performance. The responses at each setting of controllable parameters are considered as a measure, indicative of not only the mean of the performance but also its variance. The mean and the variance are transformed to a single performance measure known as the signal-to-noise (S/N) ratio. The term “signal” signifies the desirable value (mean response variable) and “noise” signifies the undesirable value (standard deviation), so the S/N ratio represents the amount of variation presents in the response variable. Thus, the goal is to maximize the S/N ratio.This paper uses the Taguchi method as a reliable method of setting the parameters for the proposed algorithm. The DE algorithm has five parameters: (1) population sizeNP; (2) scaling factor F; (3) crossover rate CR; (4) maximum number of iterations Itermax; and (5) maximum number of nonimprovement iterations Itermax−non. And, the stopping condition of the column generation process involves one parameter IternonCG, i.e., the number of consecutive iterations of failing to find a column with positive reduced cost. Five levels are considered for each parameter, as shown in Table 2. The orthogonal array L25(65) shown in Table 3 is designed as the experimental layout, which displays 25 different combinations of factor levels for the experiments.Table 2
Levels of parameters.
Parameter (factor)NPFCRItermaxItermax−nonIternonCGLevel1100.20.251001052150.40.5020025103200.60.7530050204250.80.9540075305301150010040Table 3
The orthogonal array.
Experiment no.ParameterNPFCRItermaxItermax−nonIternonCG11111112122222313333341444445155555621234572234518234512924512310251234113135241232413513335241143413521535241316414253174253141843142519442531204531422151543222521543235321542454321525554321For each of three different-sized problems, five instances are randomly selected and then tested on the 25 parameter combinations. The obtained computational data are processed by Minitab 19 statistical software. The S/N ratio trends of parameter levels for different-sized instances are displayed in Figures4–6. As indicated in the three figures, we set NP = 20, F = 0.6, CR = 0.75, Itermax = 100, Itermax−non = 50, and IternonCG = 10 for small-sized and medium-sized instances. And, for large-sized instances, we set NP = 30, F = 0.6, CR = 0.75, Itermax = 500, Itermax−non = 100, and IternonCG = 20, respectively.Figure 4
The S/N ratio trend for small-sized problems.Figure 5
The S/N ratio trend for medium-sized problems.Figure 6
The S/N ratio trend for large-sized problems.


## 5.2.2. Performance Metrics

We use the gap between the revenue and the upper bound as the performance metric. The revenue is obtained from a given solution method, including either the proposed column generation approach or the integer programming formulation solved by CPLEX solver. The upper bound is obtained by solving the MILP formulation in Section3 with a time limit of 3600 s. Furthermore, the number of optimal solutions out of 5 instances of each subclass and the average computation time are recorded as auxiliary metrics.


## 5.3. Computational Results and Analysis

### 5.3.1. Results of Small-Sized Instances

We first verify the performance of the proposed column generation approach (CG for short) for solving small-sized instances. Tables4 and 5 display the comparison results between the MILP formulation and the proposed column generation approach. The first column represents the |N| × |M| combinations. Column 2 and column 3 show parameter TF and parameter R, respectively. We report in columns 4–9 the minimum (Min), average (Avg) and maximum (Max) gaps between the upper bound and the best found solution obtained by the respective method for every subclass. In columns 10 and 11, we display the number of optimal solutions found by the given solution method. The following two columns with the same heading “Total” report the average computation time consumed by solution approaches. The remaining columns show the time spent on each step of column generation. Specifically, columns “CH,” “LP,” “MILP,” and “DE” report the time spent by the constructive heuristic, LP solver, MILP solver, and the proposed DE algorithm, respectively.Table 4
Computation results of small-sized instances with |N| = 15.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU (s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE15 × 20.20.20000005550.803.740.030.821.501.380.60000005578.683.460.031.330.511.5910000005531.003.810.021.061.201.520.60.20000005574.302.200.030.670.251.250.60000005582.756.170.12.571.212.2910000005521.006.420.023.500.452.4510.2000000559.001.790.020.470.081.220.60000005567.522.940.021.090.051.7810000005511.803.570.021.290.142.1315 × 30.20.20000005515.213.250.031.090.591.540.6000000559.532.040.030.550.231.2310000005510.861.790.030.490.161.110.60.20000005512.321.390.020.310.080.970.60000005511.065.740.032.640.192.8810000005518.621.840.020.480.191.1410.20000005513.681.340.020.240.081.000.60000005510.221.750.020.380.101.2510000005511.851.230.020.220.050.94Average——000000——30.013.030.031.070.391.54Table 5
Computation results of small-sized instances with |N| = 20.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE20 × 20.20.200000055119.143.240.031.080.681.450.60000005568.634.490.031.910.641.9110000005532.273.580.031.300.281.970.60.20000005545.836.110.041.902.251.920.60000005522.234.880.022.460.771.6310000005542.655.050.022.040.462.5310.20000005520.606.550.022.460.054.020.60000005519.867.190.022.740.424.0110000005529.805.830.022.230.323.2620 × 30.20.20000005545.215.240.071.530.543.100.60000005539.035.280.121.740.323.1010000005571.076.270.031.950.613.680.60.20000005548.036.030.031.950.783.270.60000005557.816.750.022.050.514.1710000005545.025.710.021.580.573.5410.20000005536.206.820.022.100.154.550.60000005537.693.190.020.680.112.3810000005542.033.320.130.610.182.40Average——000000——45.735.310.041.800.542.94As indicated by the gap values in Tables4 and 5, both the MILP formulation and column generation solve all instances to optimality. In comparison with the MILP formulation, column generation consumes significantly less running time to find the optimal solutions. That is, the MILP formulation needs computing time increasing from few seconds to hundreds of seconds, whereas column generation is able to find the optimal solution in seconds. For example, the average CPU time reduces from 119.14 s with the MILP formulation to 3.24 s with the proposed algorithm for instances with |N| = 20, |M| = 2, TF = 0.2, and R = 0.2.To sum up, the numerical experiments on small-sized instances show that the proposed algorithm significantly outperforms the MILP formulation in terms of computational efficiency.


### 5.3.2. Results of Medium-Sized Instances

We further examine the performances of the proposed method and the MILP formulation for solving medium-sized instances. Tables6 and 7 present the experiment results.Table 6
Computation results of medium-sized instances with |N| = 40.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE40 × 20.20.20.002.1310.66000451610.2030.180.113.2125.071.790.600000055134.4710.180.073.554.422.14100000055213.2010.250.123.484.022.630.60.2014.6833.1901.333.91222204.82358.080.142.96347.327.660.603.7917.4500.350.8833447.6787.960.072.5918.7966.521010.3722.0200.661.96132907.4052.630.074.1816.7031.6810.200000055758.6123.040.0311.890.4410.680.600000055182.6516.580.035.291.1310.13100000055110.0014.020.202.940.7810.1040 × 30.20.20000005596.619.870.084.092.912.790.60000005592.255.580.092.670.682.14100000055130.822.380.070.730.181.400.60.200000055107.2010.610.054.303.322.940.600000055239.688.020.043.181.962.8410000005575.2514.000.057.951.784.2210.20000005577.217.160.042.990.323.810.60000005596.6352.560.046.001.2545.2710000005592.203.460.031.280.301.85Average——01.724.6300.130.38——532.0539.810.074.0723.9711.70Table 7
Computation results of medium-sized instances with |N| = 60.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE60 × 20.20.24.6712.9518.48000053600103.260.146.2893.932.910.66.3213.3920.34000053600121.100.256.70110.343.81110.3312.0613.1500005360080.230.275.9570.303.710.60.215.6818.9524.5400.882.11023600180.370.124.9583.4491.860.612.5715.9720.8500.412.07043600144.270.084.41125.4214.3619.8614.3117.8200005360020.020.129.545.095.2710.200000055417.89107.970.0569.936.6731.320.6021.2831.6101.032.5011360059.80.0512.193.9143.6517.0316.5525.4600.521.3203360023.330.0810.405.517.3460 × 30.20.200000055971.8111.610.175.532.313.600.605.648.69000152961.2212.630.176.701.873.89106.6217.39000252287.8610.530.114.772.033.620.60.200000055221.2128.930.055.6819.883.320.6013.7126.03000152926.0710.190.085.141.733.241014.2833.12000252328.0312.190.056.461.793.8910.200000055228.0918.750.029.630.698.410.600000055235.2110.960.035.461.244.23102.7513.7700.261.3344935.258.390.043.830.803.72Average——3.699.3615.0700.170.52——2350.7053.590.1010.2029.8313.45As shown in Tables6 and 7, we can observe that the proposed column generation algorithm is superior to the MILP formulation in terms of efficiency and effectiveness. For instances with |N| = 40, column generation finds 83 optimal solutions out of 90 instances with 0.13% deviation on average. In contrast, the MILP formulation identifies the optimal solution in 80 out of 90 instances and reports 1.72% deviation on average. As expected, the performance of the MILP formulation degenerates since the solution space exponentially grows with the increase in problem size. When |N| equals 60, the MILP formulation only reports 36 instances out of 90 to optimality within 3600 s, whereas column generation finds 79 optimal solutions. The average CPU time of the MILP formulation increases to 2350.70 s, while column generation runs in 53.59 s on an average.Based on above comparison analysis, we can come to the conclusion that column generation still performs significantly better than the MILP formulation with regards to efficiency and effectiveness and is an effective approach for solving medium-sized instances.


### 5.3.3. Results of Large-Sized Instances

We further report computational results obtained for large-sized instances in Tables8 and 9. These instances become more challenging for the MILP formulation, which can be indicated by computation time. We can observe that column generation performs significantly better than the MILP formulation in terms of number of instances solved to optimality and computation times. Specifically, the average percentage deviation of the MILP formulation is consistently larger than that of column generation. We can also observe that the MILP formulation fails to find a feasible solution to instances with TF = 0.2, R = 1 and |N| = 100. In contrast, the proposed column generation algorithm is able to prove the optimality of 177 instances and to provide very low gaps for the remaining, always below 1%. Column generation can generally solve each instance within 800 s, which further demonstrate its efficiency.Table 8
Computation results of large-sized instances with |N| = 90.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE90 × 30.20.2∗1.608.5611.68000052082.23357.670.21301.6812.6643.120.6∗2.945.186.30000051148.00296.320.18239.8410.7145.591 ∗0.395.8311.97000052100.61218.420.19160.607.1950.440.60.2011.7418.28000152880.21248.820.07192.1219.9036.730.68.3312.8716.63000053600223.360.06177.009.5136.791 ∗7.7314.4820.38000052964.35228.280.07173.027.5747.6210.200.110.5500.060.30441080.89121.080.0463.7710.4846.790.606.6919.25000351848.10107.610.0562.123.7141.731 ∗11.1913.1315.63000052412.80133.220.0477.589.3046.3090 × 40.20.205.8214.11000053600316.040.13254.717.0654.140.61.303.6811.07000053600274.330.14214.287.7652.151∗---00005414.81225.310.12164.484.5456.170.60.206.7615.99000252376.00308.200.05243.179.7355.250.605.0013.70000151762.24271.050.04200.519.4061.101017.5178.60000151139.31189.580.04125.935.6857.9310.200000055302.1267.800.0429.611.5336.620.600000055408.0888.780.0343.623.0842.05104.4911.29000252423.28112.080.0458.964.4348.65Average——1.977.1715.6100.000.02——2007.95210.440.09154.618.0147.73“-”indicates that the CPLEX fails to find feasible solutions for this subclass. “∗” indicates that the CPLEX runs of memory when solving some instances of this subclass.Table 9
Computation results of large-sized instances with |N| = 100.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE100 × 30.20.2∗2.248.0013.9200005571.84374.270.25208.5810.50154.940.6∗5.387.318.6700005627.31247.860.24147.8010.5689.261∗---00005436.57200.050.26129.568.9061.330.60.210.0415.2419.45000053600359.940.1192.67122.7144.460.610.3715.5720.41000053600254.400.07185.3712.2756.691∗14.2415.0516.00000051711.34248.250.07165.4811.8770.8310.2∗00.411.3500.320.86331198.63147.540.0673.9512.4261.110.6∗08.5418.24000152454.21160.960.0988.2112.6060.061∗10.0210.8112.25000051055.10176.530.0599.087.4769.93100 × 40.20.2∗3.975.577.04000051751.80488.650.32389.7012.0086.630.6∗0.824.798.7500005519.20489.110.25376.2414.4198.211∗---00005457.40402.910.35264.6612.91124.990.60.207.4913.13000153000793.430.06619.8624.58148.930.6∗9.5112.1217.90000052977.87562.540.09395.9020.20146.351∗7.547.547.5400005438.96337.990.08142.0021.61174.3010.200000055410.21121.410.0639.766.2775.320.601.286.42000451164.28141.240.0673.006.2261.961∗0.555.629.39000051073.99155.210.0778.846.3469.96Average——4.677.8311.2800.020.05——1502.71314.570.14203.9318.5591.96“-” indicates that the CPLEX fails to find feasible solutions for this subclass. “∗” indicates that the CPLEX runs of memory when solving some instances of this subclass.The encouraging computational results of the proposed algorithm may be attributed to two aspects. Firstly, the column generation approach keeps most columns of the set partitioning formulation implicit and only generates necessary columns that can improve the objective value. Therefore, a considerate number of times can be saved. Secondly, the DE algorithm is an efficient heuristic. It is used to solve the subproblems and helps to identify promising columns in a very short time.In summary, the capability of the proposed algorithm to consistently find optimal or near-optimal solutions confirms its effectiveness and superiority against the MILP formulation. Our experimental results clearly show that column generation is a computationally viable solution approach for practical situations.


## 5.3.1. Results of Small-Sized Instances

We first verify the performance of the proposed column generation approach (CG for short) for solving small-sized instances. Tables4 and 5 display the comparison results between the MILP formulation and the proposed column generation approach. The first column represents the |N| × |M| combinations. Column 2 and column 3 show parameter TF and parameter R, respectively. We report in columns 4–9 the minimum (Min), average (Avg) and maximum (Max) gaps between the upper bound and the best found solution obtained by the respective method for every subclass. In columns 10 and 11, we display the number of optimal solutions found by the given solution method. The following two columns with the same heading “Total” report the average computation time consumed by solution approaches. The remaining columns show the time spent on each step of column generation. Specifically, columns “CH,” “LP,” “MILP,” and “DE” report the time spent by the constructive heuristic, LP solver, MILP solver, and the proposed DE algorithm, respectively.Table 4
Computation results of small-sized instances with |N| = 15.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU (s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE15 × 20.20.20000005550.803.740.030.821.501.380.60000005578.683.460.031.330.511.5910000005531.003.810.021.061.201.520.60.20000005574.302.200.030.670.251.250.60000005582.756.170.12.571.212.2910000005521.006.420.023.500.452.4510.2000000559.001.790.020.470.081.220.60000005567.522.940.021.090.051.7810000005511.803.570.021.290.142.1315 × 30.20.20000005515.213.250.031.090.591.540.6000000559.532.040.030.550.231.2310000005510.861.790.030.490.161.110.60.20000005512.321.390.020.310.080.970.60000005511.065.740.032.640.192.8810000005518.621.840.020.480.191.1410.20000005513.681.340.020.240.081.000.60000005510.221.750.020.380.101.2510000005511.851.230.020.220.050.94Average——000000——30.013.030.031.070.391.54Table 5
Computation results of small-sized instances with |N| = 20.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE20 × 20.20.200000055119.143.240.031.080.681.450.60000005568.634.490.031.910.641.9110000005532.273.580.031.300.281.970.60.20000005545.836.110.041.902.251.920.60000005522.234.880.022.460.771.6310000005542.655.050.022.040.462.5310.20000005520.606.550.022.460.054.020.60000005519.867.190.022.740.424.0110000005529.805.830.022.230.323.2620 × 30.20.20000005545.215.240.071.530.543.100.60000005539.035.280.121.740.323.1010000005571.076.270.031.950.613.680.60.20000005548.036.030.031.950.783.270.60000005557.816.750.022.050.514.1710000005545.025.710.021.580.573.5410.20000005536.206.820.022.100.154.550.60000005537.693.190.020.680.112.3810000005542.033.320.130.610.182.40Average——000000——45.735.310.041.800.542.94As indicated by the gap values in Tables4 and 5, both the MILP formulation and column generation solve all instances to optimality. In comparison with the MILP formulation, column generation consumes significantly less running time to find the optimal solutions. That is, the MILP formulation needs computing time increasing from few seconds to hundreds of seconds, whereas column generation is able to find the optimal solution in seconds. For example, the average CPU time reduces from 119.14 s with the MILP formulation to 3.24 s with the proposed algorithm for instances with |N| = 20, |M| = 2, TF = 0.2, and R = 0.2.To sum up, the numerical experiments on small-sized instances show that the proposed algorithm significantly outperforms the MILP formulation in terms of computational efficiency.


## 5.3.2. Results of Medium-Sized Instances

We further examine the performances of the proposed method and the MILP formulation for solving medium-sized instances. Tables6 and 7 present the experiment results.Table 6
Computation results of medium-sized instances with |N| = 40.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE40 × 20.20.20.002.1310.66000451610.2030.180.113.2125.071.790.600000055134.4710.180.073.554.422.14100000055213.2010.250.123.484.022.630.60.2014.6833.1901.333.91222204.82358.080.142.96347.327.660.603.7917.4500.350.8833447.6787.960.072.5918.7966.521010.3722.0200.661.96132907.4052.630.074.1816.7031.6810.200000055758.6123.040.0311.890.4410.680.600000055182.6516.580.035.291.1310.13100000055110.0014.020.202.940.7810.1040 × 30.20.20000005596.619.870.084.092.912.790.60000005592.255.580.092.670.682.14100000055130.822.380.070.730.181.400.60.200000055107.2010.610.054.303.322.940.600000055239.688.020.043.181.962.8410000005575.2514.000.057.951.784.2210.20000005577.217.160.042.990.323.810.60000005596.6352.560.046.001.2545.2710000005592.203.460.031.280.301.85Average——01.724.6300.130.38——532.0539.810.074.0723.9711.70Table 7
Computation results of medium-sized instances with |N| = 60.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE60 × 20.20.24.6712.9518.48000053600103.260.146.2893.932.910.66.3213.3920.34000053600121.100.256.70110.343.81110.3312.0613.1500005360080.230.275.9570.303.710.60.215.6818.9524.5400.882.11023600180.370.124.9583.4491.860.612.5715.9720.8500.412.07043600144.270.084.41125.4214.3619.8614.3117.8200005360020.020.129.545.095.2710.200000055417.89107.970.0569.936.6731.320.6021.2831.6101.032.5011360059.80.0512.193.9143.6517.0316.5525.4600.521.3203360023.330.0810.405.517.3460 × 30.20.200000055971.8111.610.175.532.313.600.605.648.69000152961.2212.630.176.701.873.89106.6217.39000252287.8610.530.114.772.033.620.60.200000055221.2128.930.055.6819.883.320.6013.7126.03000152926.0710.190.085.141.733.241014.2833.12000252328.0312.190.056.461.793.8910.200000055228.0918.750.029.630.698.410.600000055235.2110.960.035.461.244.23102.7513.7700.261.3344935.258.390.043.830.803.72Average——3.699.3615.0700.170.52——2350.7053.590.1010.2029.8313.45As shown in Tables6 and 7, we can observe that the proposed column generation algorithm is superior to the MILP formulation in terms of efficiency and effectiveness. For instances with |N| = 40, column generation finds 83 optimal solutions out of 90 instances with 0.13% deviation on average. In contrast, the MILP formulation identifies the optimal solution in 80 out of 90 instances and reports 1.72% deviation on average. As expected, the performance of the MILP formulation degenerates since the solution space exponentially grows with the increase in problem size. When |N| equals 60, the MILP formulation only reports 36 instances out of 90 to optimality within 3600 s, whereas column generation finds 79 optimal solutions. The average CPU time of the MILP formulation increases to 2350.70 s, while column generation runs in 53.59 s on an average.Based on above comparison analysis, we can come to the conclusion that column generation still performs significantly better than the MILP formulation with regards to efficiency and effectiveness and is an effective approach for solving medium-sized instances.


## 5.3.3. Results of Large-Sized Instances

We further report computational results obtained for large-sized instances in Tables8 and 9. These instances become more challenging for the MILP formulation, which can be indicated by computation time. We can observe that column generation performs significantly better than the MILP formulation in terms of number of instances solved to optimality and computation times. Specifically, the average percentage deviation of the MILP formulation is consistently larger than that of column generation. We can also observe that the MILP formulation fails to find a feasible solution to instances with TF = 0.2, R = 1 and |N| = 100. In contrast, the proposed column generation algorithm is able to prove the optimality of 177 instances and to provide very low gaps for the remaining, always below 1%. Column generation can generally solve each instance within 800 s, which further demonstrate its efficiency.Table 8
Computation results of large-sized instances with |N| = 90.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE90 × 30.20.2∗1.608.5611.68000052082.23357.670.21301.6812.6643.120.6∗2.945.186.30000051148.00296.320.18239.8410.7145.591 ∗0.395.8311.97000052100.61218.420.19160.607.1950.440.60.2011.7418.28000152880.21248.820.07192.1219.9036.730.68.3312.8716.63000053600223.360.06177.009.5136.791 ∗7.7314.4820.38000052964.35228.280.07173.027.5747.6210.200.110.5500.060.30441080.89121.080.0463.7710.4846.790.606.6919.25000351848.10107.610.0562.123.7141.731 ∗11.1913.1315.63000052412.80133.220.0477.589.3046.3090 × 40.20.205.8214.11000053600316.040.13254.717.0654.140.61.303.6811.07000053600274.330.14214.287.7652.151∗---00005414.81225.310.12164.484.5456.170.60.206.7615.99000252376.00308.200.05243.179.7355.250.605.0013.70000151762.24271.050.04200.519.4061.101017.5178.60000151139.31189.580.04125.935.6857.9310.200000055302.1267.800.0429.611.5336.620.600000055408.0888.780.0343.623.0842.05104.4911.29000252423.28112.080.0458.964.4348.65Average——1.977.1715.6100.000.02——2007.95210.440.09154.618.0147.73“-”indicates that the CPLEX fails to find feasible solutions for this subclass. “∗” indicates that the CPLEX runs of memory when solving some instances of this subclass.Table 9
Computation results of large-sized instances with |N| = 100.
|N| × |M|TFRGap (%)No. of optimal solutionsCPU(s)MILPCGMILPCGMinAveMaxMinAveMaxMILPCGTotalTotalCHLPMILPDE100 × 30.20.2∗2.248.0013.9200005571.84374.270.25208.5810.50154.940.6∗5.387.318.6700005627.31247.860.24147.8010.5689.261∗---00005436.57200.050.26129.568.9061.330.60.210.0415.2419.45000053600359.940.1192.67122.7144.460.610.3715.5720.41000053600254.400.07185.3712.2756.691∗14.2415.0516.00000051711.34248.250.07165.4811.8770.8310.2∗00.411.3500.320.86331198.63147.540.0673.9512.4261.110.6∗08.5418.24000152454.21160.960.0988.2112.6060.061∗10.0210.8112.25000051055.10176.530.0599.087.4769.93100 × 40.20.2∗3.975.577.04000051751.80488.650.32389.7012.0086.630.6∗0.824.798.7500005519.20489.110.25376.2414.4198.211∗---00005457.40402.910.35264.6612.91124.990.60.207.4913.13000153000793.430.06619.8624.58148.930.6∗9.5112.1217.90000052977.87562.540.09395.9020.20146.351∗7.547.547.5400005438.96337.990.08142.0021.61174.3010.200000055410.21121.410.0639.766.2775.320.601.286.42000451164.28141.240.0673.006.2261.961∗0.555.629.39000051073.99155.210.0778.846.3469.96Average——4.677.8311.2800.020.05——1502.71314.570.14203.9318.5591.96“-” indicates that the CPLEX fails to find feasible solutions for this subclass. “∗” indicates that the CPLEX runs of memory when solving some instances of this subclass.The encouraging computational results of the proposed algorithm may be attributed to two aspects. Firstly, the column generation approach keeps most columns of the set partitioning formulation implicit and only generates necessary columns that can improve the objective value. Therefore, a considerate number of times can be saved. Secondly, the DE algorithm is an efficient heuristic. It is used to solve the subproblems and helps to identify promising columns in a very short time.In summary, the capability of the proposed algorithm to consistently find optimal or near-optimal solutions confirms its effectiveness and superiority against the MILP formulation. Our experimental results clearly show that column generation is a computationally viable solution approach for practical situations.


## 5.4. Effect of ParametersTF and R on the Performance of the Methods

In this section, we further investigate the impact of the instance generation parameters on the performance of the methods. As stated in Section5.1, parameter TF determines due date tightness while parameter R controls the due date range. Additionally, parameter TF also controls the dispersion degree of the release dates of the jobs, in a way that the higher the TF is, the more dispersed the release dates are.As indicated by the gap and computation time values in Tables4–9, we can notice that given that all other specifications are the same, instances with large values of TF are generally easier to solve than those with small values. A reasonable explanation for this behaviour could be as follows. A small TF value corresponds to a situation where the release dates are distributed in the very beginning of the time horizon and then are not very restrictive. And, a very large TF value leads to a situation where the release dates are very dispersed. Hence, as TF increases, we receive incoming orders released at dispersed time points with tight due dates. Therefore, more orders are expected to be rejected, which in turn implies the reduction of the size of the solution space. Consequently, the instances with large TF values are less challenging. We can observe from the tables that the number of instances solved to optimality by the MILP formulation increases and the average CPU time spent by the MILP formulation tends to decrease as the TF gets large when it comes to medium- and large-sized instances. A similar observation can be found when instances are solved by the proposed column generation algorithm. That is, two subclasses are all exactly solved by the proposed algorithm, but the subclass with TF = 1 evidently consumes less CPU time than the subclass with TF = 0.2.Given theTF value, the due date range increases with the increase of the parameter R. Results in Tables 4–9 show that given the TF value, some subclasses with large R value are more difficult to solve than those with small R value, but we can also find some exceptions. This coincides with observations reported by Wang and Ye [9]. Therefore, there is no clear trend in correlation between the R value and the hardness of the problem.


## 6. Conclusion

In this paper, we study a resource-constrained order acceptance and scheduling on unrelated parallel machines. We present two different models for this problem, i.e., the MILP formulation and the set partitioning formulation. With the MILP formulation including all decision variables, we can directly solve problem instances via general-purpose commercial solvers without complicated programming. In order to solve large-scale instances, we develop a column generation algorithm based on the set partitioning formulation. In the proposed column generation algorithm, we propose a constructive heuristic to initialize the restricted linear master problem and design a differential evolution algorithm to solve the subproblems efficiently. Numerical experiments based on a set of randomly generated test instances are carried out to verify the performance of the proposed algorithm. The results demonstrate that the MILP formulation can identify optimal solutions for small-scale instances and report solutions with large optimality gaps at termination for medium- and large-sized instances. Column generation is advantageous over the MILP formulation in terms of effectiveness and efficiency and quite suitable for problems with different sizes. We also analyse the effect of parametersTF and R on the performance of the proposed algorithm and the MILP formulation.Further research will concentrate on problem extensions and algorithmic enhancements. One possible extension of the current work is to consider stochasticities in processing times and develop a stochastic and robust counterpart for the column generation algorithm. Considering some practical assumptions including machine breakdown, parallel batch processing, and order interruption are complex but extremely valuable variations to the studied problem. Furthermore, we also attempt to improve the performance of the proposed column generation algorithm by solving the subproblems via a dynamic algorithm in conjunction with heuristics.



---

*Source: 5566002-2021-10-20.xml*