Order Acceptance and Scheduling Problem with Carbon Emission Reduction and Electricity Tariffs on a Single Machine

Order acceptance and scheduling (OAS) problems are realistic for enterprises. They have to select the appropriate orders according to their capacity limitations and profit consideration, and then complete these orders by their due dates or no later than their deadlines. OAS problems have attracted significant attention in supply chain management. However, there is an issue that has not been studied well. To our best knowledge, no prior research examines the carbon emission cost and the time-of-use electricity cost in the OAS problems. The carbon emission during the on-peak hours is lower than the one in mid-peak and off-peak hours. However, the electricity cost during the on-peak hours is higher than the one during mid-peak and off-peak hours when time-of-use electricity (TOU) tariff is used. There is a trade-off between sustainable scheduling and the electricity cost. To calculate the objective value, a carbon tax and carbon dioxide emission factor are included when we evaluate the carbon emission cost. The objective function is to maximize the total revenue of the accepted orders and then subtract the carbon emission cost and the electricity cost under different time intervals on a single machine with sequence-dependent setup times and release date. This research proposes a mixed-integer linear programming model (MILP) and a relaxation method of MILP model to solve this problem. It is of importance because the OAS problems are practical in industry. This paper could attract the attention of academic researchers as well as the practitioners.


Introduction
Numerous traditional scheduling problems assume a set of n orders to be accepted and scheduled in a specified production environment. An order i typically quotes a due date d i and a deadline d i . The due date is less than or equal to deadline. After an order i is completed at time C i , the company receives a profit e i if C i is before d i and may get weighted penalties in the case of late delivery [1]. This profit calculation is depicted in Figure 1. A firm, nonetheless, must often choose the appropriate orders according to its available capacity, market focus, competitive advantage, or a combination of these [2]. This scheduling problem has transformed into an order acceptance and scheduling (OAS) problem involving the decision of accepting or rejecting orders and then scheduling the accepted orders on machines. Engels et al. [3] defined the OAS problem to be a combination of the knapsack and scheduling problems. The general objective is to maximize the total revenue Q i minus the tardiness penalties. Because OAS is a practical problem for enterprises in the supply chain, it has gained increasing research attention in recent decades [3][4][5][6][7][8]. OAS problems may involve a single machine [3,9,10], parallel machines [6,[11][12][13], or a permutation flow shop [5,8,14]. To address these three machine types, some studies have proposed exact algorithms, such as the branch-and-bound algorithm [9], mixed integer programming [15], and dynamic programming [3,16]. Oguz et al. [1] found some useful upper bound equations which reduce the computational effort in the branch-and-bound algorithm. However, if the number of job is larger than 15, branch-and-bound algorithm cannot solve them in 1 h. Nobibon and Leus [9] proposed the time-indexed formulations in a branch-and-bound algorithm with up to 50 jobs in a reasonable time.
For large-size and NP-hard problems [17], no algorithm can solve the problem in polynomial time. Many researchers have employed metaheuristic methods, such as genetic algorithms [2,[18][19][20], tabu search [4], estimation of distribution algorithm (EDAs) [21], simulated annealing [15], and the artificial bee colony algorithm [8]. Cesaret et al. [4] proposed a tabu search algorithm, which was effective compared with m-ATCS and ISFAN based on the deviation from the UBs provided by MILP. The work of Wattanapornprom et al. [21] might be the first research to solve the OAS problem by EDAs. Their node based EDAs gain better results compared to GA with local search. In addition, Wattanapornprom et al. [21] found the consideration of OAS which increased the utilization rate in the supply chain. Finally, Slotnick [7] wrote an extensive review of the OAS problem, which is a good reference to understand the background.
Even though OAS problems have been studied extensively [7], prior OAS studied have not adopted the concept of green manufacturing that has received increasing attention among enterprises [22][23][24][25]. In particular, some recent works studied the carbon emission in the scheduling problems [26][27][28]. This paper focuses on the electricity power and excludes fuel consumption. Electricity power suppliers use gas-fired generation plants to supplement coal-fired generation plants during on-peak hours [29]. Because gas-fired generation plants usually emit less carbon per kilowatt hour of electricity produced than coal-fired generation plants, carbon emission during on-peak hours is lower than that during mid-peak and off-peak hours [30,31].
By contrast, the electricity cost during the on-peak hours is higher than that during mid-peak and off-peak hours if time-of-use electricity (TOU) tariff is used. The electricity cost of coal-fired generation plants is lower than that of gas-fired generation plants. There is an apparent trade-off between green scheduling and the electricity cost. To the best of our knowledge, although Zhang et al. [31] studied the trade-off for flowshop scheduling problems, no study on OAS has investigated sustainable scheduling, TOU costs, and a combination of both. The present study aims to solve the OAS problem and achieve carbon footprint reduction and minimization of electricity fees, which are the major contributions of this study.
To address the trade-off between carbon footprint reduction and TOU electricity costs, the carbon emission tax was introduced. In some countries, the carbon emission tax is defined per tonne (https://en.wikipedia.org/wiki/Carbon_tax). This study examined the carbon emission tax of selected countries, and they are listed in Table 1. For example, in Japan, the carbon emission tax per tonne is 2400 Yen (US$21.852 as of May 2019). Sweden charges the highest carbon emission tax. The average carbon emission tax per ton is $26.73155 (Table 1). Given the average carbon emission tax, the revenue of each accepted order is subtracted from the TOU electricity cost and carbon emission tax according to the power consumption of each order. This paper proposes a new OAS problem considering the carbon emission cost and the TOU tariff. To solve this problem, a mathematical model of mixed-integer linear programming (MILP) was formulated and solved using IBM CPLEX 12.9. MILP normally requires a higher computational time to obtain an optimal solution. Prior studies investigated linear programming (LP) relaxation and applying less binary variables [32], preprocessing [33], objective scaling [34], Lagrangian Relaxation [12], MILP-based heuristics, an evolutionary algorithm within MILP [35,36], executing primal or the dual simplex method at the child nodes and different types of cuts [33,37], and so on. In this research, we use LP to obtain an initial solution which was input to the MILP model. This proposed method is MIPStart. We compared the performance of the original MILP model and that with an LP and MIPStart method.
The remainder of the paper is organized as follows. Section 2 shows the related works of order acceptance scheduling problems on a single machine. Section 3 presents the MILP formulations of the proposed problem. Section 3 presents LP relaxation and MIPStart method. Section 4 describes the LP relaxation and data generation procedures. Section 4 presents the experimental result obtained using the original MILP model and that with an initial solution. Section 5 concludes the paper and presents future research directions.

Literature Review
Regarding the OAS problem on a single machine, numerous works can be found. Stemming from Bartal et al. [16], Engels et al. [3] studied the objective function of minimizing the total weighted completion time and the penalty of the rejected jobs, ∑ i∈S w i C i + ∑ i∈S e i , instead of the makespan and sum of the penalty of the rejected orders, C max + ∑ i∈S e i . This paper proves that the studied problem is weak NP-Complete. Most importantly, this is the first research that proposes the idea of a rejection machine. The rejected orders are sent to the rejection machine. That is, a single machine scheduling problem becomes a parallel machine scheduling problem. Due to this important characteristic, they demonstrate a dynamic programming algorithm that also yields a factor approximation algorithm (FPTAS). Because this paper does not conduct any experiments on the dynamic programming algorithm and the FPTAS, we do not know whether the two algorithms could solve the maximum jobs. Rom and Slotnick [2] studied a genetic algorithm solving the OAS with tardiness penalties. The genetic algorithm employs the random key approach and determines the rejection orders when they evaluate the completion time of the orders. Additionally, some techniques are used to increase the population diversity, including clone removal, mutation, immigration, and changing the population size. Two local search operators are also used in the GAs. In their design-of-experiment (DOE), they suggested setting the population size to 200, randomly interchanging successive jobs when they perform the mutation, using the clone removal, and not applying immigration in the GAs. They compared a heuristic with the GAs. The results show the GAs outperform the myopic heuristic.
Ogus et al. [1] examined the release time, due day, deadlines, and the sequence-dependent setup times on a single machine. The objective function is to maximize the total revenue from a function of the total tardiness and deadline. They solved this problem by mixed-integer linear programming (MILP) in the small size problems up to 15 orders. Ogus et al. [1] also developed three heuristics to solve larger size problems up to 300 orders.
Unlike Engels et al. [3], when the single machine scheduling problem is transformed into a parallel machine scheduling problem, Tanaka [10] converted the sum of job completion cost into the general scheduling problem that modifies the job completion costs. The primary approach is to remove jobs i with the completion cost f (Ĉ S i ) being more than the revenue Q i . Furthermore, the rejected jobs are moved to the end of the schedule. Tanaka [10] integrated this approach in the branch-and-bound algorithm for solving six sets with different objective functions of the single machine scheduling problems. Their experimental results show the significant efficiency when compared with the branch-and-bound algorithm in Nobibon and Leus [9]. In [10], 50 orders may take less than 1 s on average; however, Nobibon and Leus [9] required at least 24 s. We could see the significant difference of their branch-and-bound algorithm even though the latter considers some firm planned orders and the computer performance is not the same either. The only two conditions this method currently does not solve are when jobs are give deadlines and the idle times are forbidden.
Cesaret et al. [4] maximized total revenue on a single machine with sequence-dependent setup times and release date. The objective function is to maximize the total revenue, ∑ i∈S (x i Q i − w i T i ). A tabu search algorithm and two heuristics are compared on the instances up to 100 jobs. Because they supply the testing instances of 10, 15, 25, 50, and 100, this research adopts their instances to test our unified approach.
Chen et al. [38] proposed improved GAs with a local search for solving the single machine order acceptance problems. Later, Chen et al. [39] presented a new diversity metric to control population diversity. Both papers apply a same-site-copy-first principle [19] in their GAs. However, the performance improvement due to this operator is not reported when compared with other crossover operators. That is why we to apply the same-site-copy-first principle in our proposed crossover operator.
For a detailed review, Slotnick [7] wrote an excellent survey of the OAS problems. This paper points out five taxonomies of the related works. According to this review and the above referred articles, there is no research study on the time-of-tariff energy cost and CO 2 emission cost. Hence, the two new considerations are the significant contribution of this research.

Mathematical Formulations
Let n be a set of orders in the OAS problem for the three machine types. An order i has its own due date d i , deadline d i , processing time p i , and arrival time r i . When order i is accepted and processed, if the completion time C i is before d i , the company receives a revenue e i . When tardiness occurs, T i must be equal to max(0, C i − d i ) and C i must be less than d i . If not, the order must be rejected without processing.
To model the formulations of the OAS problems, we employed the model provided by Oguz et al. [1], as shown in Section 3.1. Then, we added constraints for calculating the carbon emission cost and electricity tariff (Section 3.2). In the model of Oguz et al. [1], dummy order 0 and n+1 were added, denoted as O 0 and O n+1 . O 0 is the first position in the sequence, and O n+1 occupies the last position. The two dummy orders are available at time zero. In addition, the processing time, due date, deadline, and revenue are also 0. d n+1 is set as the maximum due date among all the orders. We adopted the aforementioned settings proposed by Oguz et al. [1]. We assumed that a factory works 24 h a day and the total manufacturing hours are no more than 24 h. Please access the variables used in MILP model in the abbreviation.
Finally, the time unit is the minutes for the arrival time, processing time, due day, and deadline. If an order arrives at 08:00, the value of r i is 240 min because the starting time is 0 at the midnight. The unit of profit, TOU cost and tax are in US dollars.

Basic Model of Order Acceptance Problem
Oguz et al. [1] defined two decision variables. a i represents whether the order is accepted or not. u ij is the decision variable indicating whether order i is processed by order j.
The following equations are the constraints for the OAS problems. A variable R i is used to evaluate the final profit of each order. The objective function in Equation (1) is used to maximize the total profit of each order.
∑ n j=0,j =i u ji = a i ∀i = 1, . . . , n + 1 (3) Equation (2) reveals that if order i is accepted, the order i precedes only one order, and Equation (3) shows that order i is preceded only by one order. If the order is rejected, it is not included in the sequence. Both constraints also forbid the preemption of orders. The constraint in Equation (4) indicates that, if order j is preceded by order i, order j has a greater completion time than order i; in addition, the sequence-dependent setup time s ij and processing time p j are higher. If order j is not preceded by order i, C j ≥ 0 is the only constraint. The constraint in Equation (4) ensures this result, which prevents subtour in the sequence.
Equation (5) specifies that, if the order j is accepted and begins with sequence number i in the timetable, then the completion time of order j should be at least the sum of the release date, processing time, and sequence-dependent setup cost between orders i and j. If order i does not precede order j, the completion time of order j has a relaxed limit. The constraint in Equation (6) ensures that any orders completed after the due date are not accepted. If order j is not accepted, the limit is reduced to C j ≥ 0. These constraints allow us to calculate the correct completion time of the order.
The constraints in Equations (7)-(9) define the tardiness for each order. Equations (11) and (12) set the completion time of dummy orders 0 and n+1, and the constraint in Equation (13) defines that both dummy orders are accepted.
Oguz et al. [1] explored the possibility of using an effective inequality to improve the upper bound of LP relaxation in Equations (14) and (15). These inequalities ensure that the model does not break the variable u ij into too many small pieces and cause some orders to be rejected.

Mathematical Formulations for Considering the TOU Cost and Carbon Emission
Previous studies on OAS problems have considered that the completion time C i is associated with due day d i . The present research studied the duration of each job in different TOU electricity hours and CO 2 emission intervals. To solve the problem more effectively, we introduce the starting time ST i of order i by considering TOU and CO 2 emission hours. Equation (16) can be used to calculate the starting time of order i based on its processing time, setup cost, and acceptance status a i . The constraint in Equation (17) indicates that the starting time of job i must be greater than or equal to its arrival time. Equations (16) and (17) provide the correct calculation of ST i .
Using Equations (18)- (20), we calculate the x ik according to the time-of-use electricity hours. Equation (18) indicates that, when some jobs are allocated within the period b k−1 to b k , the total processing time cannot exceed the duration of this period. Equation (19) indicates that the total processing time of job i across different electricity hours must be greater than or equal to the completion time of the job minus the starting time. This equation keeps the integrity of the processing time of each job. Equation (20) indicates that, when the job cross-interval processing occurs, it is necessary to separately calculate the processing time in which the current interval is allocated and the processing time to be allocated in the span. The parameter x i,k during this process is defined by the time spent between ST i and C i . If the job i starts and ends in the same interval, x ik is reduced to C i − ST i in Equation (20). By using Equation (20), we assume that the processing time of each order would not exceed more than two TOU emission periods.
Similarly, we use Equations (21)- (23) to calculate the production time of job i in the carbon emission period k. For most cases, we use Equation (22) to calculate z i,k because only a few orders would span different time intervals. Equation (22) also keeps the integrity of z i,k . If an order is processed across two intervals, Equation (23) can be used to calculate the first and second parts of the processing time. We also assume that the processing time of each order does not span more than two CO 2 emission periods.
In Equation (24), the first term provides the accepted order i minus the tardiness cost. The tardiness penalty is w i T i , where w i is the customer weight represented as a proportional lateness discount [2]. The second part of Equation (24) indicates the corresponding electricity bill generated by each accepted order according to the amount of time and power consumption during TOU periods. Finally, the last part is used to calculate the carbon emission cost of each order.

LP Relaxation and MIPStart Method
Although we can calculate an optimal solution in the proposed MILP optimization model, obtaining such solutions for large problems would be a time-consuming process. Relaxing the conditions for the integer decision variables u ij and a i might yield the upper-bound values in a short time. We can designate such a relaxed model as our LP model to represent the obtained results from this approach. Because u ij and a i are no longer integer values, the objective function values or upper-bound values might be higher than those derived in the MILP version. Specifically, our LP model may generate infeasible solutions. The reason is some u ij and a i are no longer 0, which increase the total revenue. However, the deadline of LP solution may not meet the real case.
In addition to our LP relaxation method, we propose another approach named MIPStart, which takes advantage of CPLEX. In this approach, the solution u ij obtained by LP relaxation is injected into the MILP as the starting solution. Because the value of u ij is a real number between 0 and 1, we convert the u ij into the integer value. If u ij is greater than or equal to 0.5, we set the u ij as 1. Otherwise, the revised u ij is 0. After this transformation, u ij is turned into the one-dimensional array required by CPLEX. CPLEX does not support an initial solution with two-dimensions.
When we inject the array, a i is obtained immediately. Because of the existence of infeasible solutions, we use an automatic model in CPLEX, which can convert an infeasible solution into a feasible solution. It is a built-in characteristic provided by CPLEX. We conducted evaluations to determine whether we could produce solutions of relatively high quality or reduce the computation time by using the mentioned models and approaches.

Testing Instances and Electricity Data
Cesaret et al. [4] supplied the testing OAS instances on a single machine. The order sizes are 10, 15, 25, 50, and 100. The problem comprises τ and R parameters of values 1, 3, 5, 7, and 9. Each problem combination has 10 instance replications. The total number of instances is 1250. The objective is to maximize the total revenue on a single machine with sequence-dependent setup times and release date. Because no study has investigated the TOU tariff for OAS problems, we generated the corresponding energy consumption requirements for the above instances. The power consumption (kW) of each order Ω i is in the range U(1, 20). The generated instances are available at Github (https: //github.com/worldstar/OpenGA/tree/master/instances/SingleMachineOASWithTOU). In our experiments, we selected τ and R with values 1, 5, and 9. In addition, the first instance replication was used because the MILP model takes a long time to yield results.
We calculated the TOU tariff and CO 2 emission cost by using data provided by Zhang et al. [31]. They provided the TOU tariff and CO 2 emissions in the summer season data shown in Tables 2 and 3, respectively. Figures 2 and 3 represent the results of the two tables. A summer day contains five TOU periods from Monday to Saturday and nine CO 2 emission periods. The carbon dioxide emission factors (q k ) for electricity are from 0.682 to 0.725. It is quite similar to the one of Ireland. Liu et al. [40] adopted the average carbon dioxide emission factors for 1991-2006, which was 0.785 kgCO 2 /kWh, while q i is 0.468 in the report for 2016 (https://www.seai.ie/resources/publications/Dublin-City-Baseline-Report. pdf). On Sunday, the off-peak is all day, hence, the electricity price does not influence the production sequence. Therefore, we only considered the TOU cost from Monday to Saturday.   The two tables clearly show that the carbon emission during the off-peak periods is higher than that during the on-peak period. We employed these electricity data to evaluate the objective function in Equation (24). In addition, based on the average carbon tax in Table 1, the average tax of carbon emission is $0.02673155 per kg.

Experiment Results
This section presents the results obtained by our proposed MILP mathematics model, LP relaxation model, and MIPStart method; these methods were run in IBM CPLEX 12.9 on an iMac 2017 with an Intel i5 3.4 GHz CPU and 16 GB of RAM. We applied the default setting of optimal gap value, which is 1e-06. Because problems with orders of magnitude larger than 25 are associated with high total computation times, each run of the program was limited to 1 h. Moreover, because we obtained the LP results before running any MIPStart experiments, we determined that the computation time of a problem with a 100-job order is 1 h instead of a few minutes. Hence, we set the computation time for both MILP and LP to 0.5 h when solving the 100-job instances.
In Table 4, the second column presents the minimum result of UB Min obtained from MILP UB , LP UB , and MIPStart UB for orders with 10, 15, and 20 jobs. UB Min ensures that we have a tight upper-bound collection in the three methods. The abbreviation Obj represents the objective functions of the three algorithms within 1 h of CPU time. The percentage difference between the objective of each algorithm and UB Min is presented in the Gap column. The gap of MILP can be calculated through Equation (25). Obj is the value of the objective function of the MILP model. Consider, for example, the instance 10-Tao1R1_1; the percentage difference between the objective value and the global upper bound UB Min is 0.01%. All the small instances can be solved in 1 h, except for 20-Tao5R1_1 and 20-Tao5R9_1. This thus explains why the Gap of MILP is close to 0. Both 20-Tao5R1_1 and 20-Tao5R9_1 may need additional CPU time to obtain optimal solutions.
Because LP relaxation usually yields infeasible solutions, the objective value or upper-bound value is larger than UB Min . Hence, we can use Equation (26) to calculate the error ratio between the objective value of LP relaxation and UB Min . Of 27 instances, 15 have error ratios that are higher than 3%, despite the maximum CPU time being less than 0.2 s. The last four columns in the table present the results obtained for MIPStart; among 27 instances, 15 (e.g., 10-Tao1R1_1, 10-Tao1R9_1, and 10-Tao5R1_1) are associated with lower upper bounds than those provided by MILP. (The MILP UB value observed for 10-Tao1R1_1 is larger than the MIPStart UB value by 0.0001168.) The error ratio of MIPStart can also be derived using Equation (25). Table 5 presents the results of medium to large problems. Only four instances can be solved in 1 h. Our LP model can solve 25-job and 50-job problems in a few seconds; nevertheless, six orders require at least 1 h of computation time. Therefore, when solving the 100-job problem, we limited the CPU time for the LP and MILP optimization processes to 0.5 h for MIPStart. We observed that, for 16 instances, MILP provides a lower upper bound than that provided by MIPStart. For 11 instances, MIPStart provides superior upper bounds to those provided by MILP. As indicated in Tables 4 and 5, the average error ratios derived for MILP, LP, and MIPStart are 5.97%, 9.59%, and 7.08%, respectively. The LP model has the highest error ratio among the three methods and generates some infeasible solutions; accordingly, LP optimization is not suitable for our studied problem. The MIPStart method has a higher error ratio (higher by 1.1%) than does the MILP optimization algorithm. The reason is that some infeasible solutions generated by LP are injected into MIPStart, and those injected solutions are not useful for pruning nodes in CPLEX 12.9. Hence, future research could employ some metaheuristic algorithms, such as genetic algorithms [2,18,20], Tabu search [4], and simulated annealing [15], to solve this new problem first and then inject the new solution into MILP.

Conclusions and Future Research
Even though OAS problem has been studied extensively [4,7,9,33], no previous research has studied the effects of the carbon footprint reduction and the TOU tariff. As a result, the major contribution of this paper is to solve the new OAS problem considering the two effects. Because a trade-off exists between carbon emission and electricity cost, a practical implication for electricity policy maker is that they cannot promote one direction alone. To calculate the CO 2 emission cost, a per-kilogram carbon tax was introduced so that the carbon emissions were transformed into the cost function. We extended the MILP formulations of Oguz et al. [1] to solve this problem with CPLEX. If the conditions for two decision variables were relaxed, this MILP optimization would be converted into an LP model.
LP model usually solves a problem in a short time compared to MILP model. We further use LP relaxation to obtain an initial solution and then inject that initial solution into the MILP model, thus yielding a third method, namely MIPStart. When we compared the performance of the three methods, we determined that MILP is often significantly superior to LP because LP cannot feasibly solve many solutions. MIPStart would obtain few tighter upper bounds, which might be not useful to improve computational burden. Original structure of MILP is slightly superior to MIPStart in terms of the error ratio. The result is similar to the one of Alemany et al. [32]. We recommend that practitioners first use an algorithm to derive a feasible constructed solution in a short time and then adopt this constructed solution in MILP optimization.
For future research, we suggest that scholars employ some mathematical methods that could enhance the solution quality of the MILP model [9,41], including time-indexed formulations [9], preprocessing, cover cuts [33], and Lagrangian Relaxation [12]. Nobibon and Leus [9] employed time-indexed formulations, and solved problems with up to 50 jobs in a reasonable time. Other promising techniques include preprocessing and cover cuts in commercial solvers [33] and Lagrangian Relaxation [12]. Moreover, we could apply real-world requirements to simplify the proposed problem, which considers an OAS with carbon footprint reduction or with a TOU tariff separately. Metaheuristic algorithms, such as Tabu search [4], simulated annealing [15], and genetic algorithms [2,20], have also been proven effective for OAS problems. These algorithms are suitable for generating suitable initial solutions in a short time; such solutions can be injected into an MILP model. Deriving lower UBs and calculating an optimal solution in a reasonable time would be valuable. appreciate Oguz who provided the testing instances and the detailed explanation of their math model. Finally, we thank IBM for providing the academic version of CPLEX 12.9 for free, which enabled us to solve the MILP model in this research.

Conflicts of Interest:
The authors declare no conflict of interest. Starting time of carbon emission period k, ∀ 1 ≤ k ≤ h q k CO 2 emission per kWh during period k, ∀ 1 ≤ k ≤ m z i,k Amount of time of order i in the CO 2 emission period k Tax CO 2 emission tax per kg m

Abbreviations
Number of periods of TOU tariff of a day b k Starting time of TOU tariff period k, ∀ 1 ≤ k ≤ m EC k Electricity cost during period k, ∀ 1 ≤ k ≤ m x i,k Amount of time of order i in the TOU tariff period k