Energy Saving for Tissue Paper Mills by Energy-Efﬁciency Scheduling under Time-of-Use Electricity Tariffs

: Environmental concerns and soaring energy prices have brought huge pressure of energy saving and emission reduction to tissue paper mills. Electricity is one of the main energy sources of tissue paper mills. The production characteristics of tissue paper mills make it easy to decrease energy cost by using time-of-use (TOU) electricity tariffs. This study investigates the bi-objective energy-efﬁciency scheduling of tissue paper mills under time-of-use electricity tariffs, the objectives of which are makespan and energy cost. First, considering the processing energy cost, setup energy cost, and transportation energy cost, an energy cost model of a tissue paper mill under TOU electricity tariffs is established. Second, the energy-efﬁciency scheduling model under TOU electricity tariffs is built based on the energy cost model. Finally, on the basis of decomposition and teaching– learning optimization, this study proposes a novel multi-objective evolutionary algorithm and further combined with the variable neighborhood search to solve the problem. The case study results demonstrate that our study of tissue paper mill energy saving is feasible, and the proposed method has better performance than the existing methods.


Introduction
The paper industry is an energy-intensive industry. In 2019, its energy consumption reached 1483.2 trillion British thermal units (Btu), of which electricity consumption was 205.3 trillion Btu. [1]. At the same time, the carbon dioxide emissions of the paper industry reached 48 million tons in 2019, accounting for 4.45% of the total carbon dioxide emissions of the manufacturing industry [2]. Under the pressure from environmental degradation and rising energy costs, the papermaking industry is faced with huge energy-saving and emission reducing pressure. Studying the energy-efficiency scheduling under the time-ofuse (TOU) electricity tariffs is of great significance for papermaking enterprises to improve energy efficiency. In general, one day can be divided into three periods based on the electricity demand, namely, the on-peak period, the mid-peak period, and the non-peak period. Due to the imbalance in the demand for electricity, waste will be generated in the off-peak or mid-peak periods, while the demand will exceed the supply during the on-peak period. To deal with this issue, the TOU pricing scheme is used to encourage users to use electricity in the mid-peak period and off-peak period, thus balancing the demand of electricity. In the production process, the electricity consumption of each job is different, indicating that the operation with high electricity consumption is arranged in off-peak or mid-peak periods and the operation with low electricity consumption is arranged in an on-peak period, so that the electricity cost can be decreased. The energy-efficiency scheduling under TOU electricity tariffs is based on above idea, which has attracted more attention in recent years. composed of a pulper, a refiner, and a chest. The pulper is the first equipment used in the pulping and papermaking stage. The paperboard is first mixed with water and then disintegrated by the pulper to obtain the pulp with a certain concentration. The obtained pulp is stored in the chest (there are usually multiple chests). However, the quality of the pulp disintegrated by the pulper cannot meet the requirements of papermaking. Therefore, it is necessary to further improve the quality of the pulp using the refiner. The pulp obtained by the refiner will be stored in other chests. The pulper and the refiner operate intermittently, which mainly depends on the liquid level of the chest. Afterwards, the pulp is sent to the forming sector through the approach flow system. The forming sector removes part of the water, and finally, the paper is formed. Next, the paper is further dehydrated by the press section. The dryness of the paper from the press section falls far short of the requirements. Paper is dried in the dryer section to make the dryness of the paper meet the requirements. The last step of the pulping and papermaking stage is the winding section. The paper is reeled by the reel cylinder, thus forming the parent roll. The parent roll will be further processed in the conversion stage.
Processes 2021, 9, x FOR PEER REVIEW 3 of 24 which is composed of a pulper, a refiner, and a chest. The pulper is the first equipment used in the pulping and papermaking stage. The paperboard is first mixed with water and then disintegrated by the pulper to obtain the pulp with a certain concentration. The obtained pulp is stored in the chest (there are usually multiple chests). However, the quality of the pulp disintegrated by the pulper cannot meet the requirements of papermaking. Therefore, it is necessary to further improve the quality of the pulp using the refiner. The pulp obtained by the refiner will be stored in other chests. The pulper and the refiner operate intermittently, which mainly depends on the liquid level of the chest. Afterwards, the pulp is sent to the forming sector through the approach flow system. The forming sector removes part of the water, and finally, the paper is formed. Next, the paper is further dehydrated by the press section. The dryness of the paper from the press section falls far short of the requirements. Paper is dried in the dryer section to make the dryness of the paper meet the requirements. The last step of the pulping and papermaking stage is the winding section. The paper is reeled by the reel cylinder, thus forming the parent roll.
The parent roll will be further processed in the conversion stage. In the conversion stage, a certain number of parent rolls are simultaneously rolled into reels using the winder, and the number of the parent rolls depends on the numbers of paper layers in the finished product. The reel paper is cut into small reels with a cutter. Finally, the small reels are packaged into the finished product by the packaging machine. Discrete production mainly occurs in the pulping process and the processing between the pulping and papermaking stage and the conversion stage. The processing is continuous in the papermaking or conversion stage. The interval time between two stages is uncertain. Different from other flow shops, only when enough parent rolls are produced in the papermaking stage can the conversion stage be started, rather than after finishing a roll of In the conversion stage, a certain number of parent rolls are simultaneously rolled into reels using the winder, and the number of the parent rolls depends on the numbers of paper layers in the finished product. The reel paper is cut into small reels with a cutter. Finally, the small reels are packaged into the finished product by the packaging machine. Discrete production mainly occurs in the pulping process and the processing between the pulping and papermaking stage and the conversion stage. The processing is continuous in the papermaking or conversion stage. The interval time between two stages is uncertain. Different from other flow shops, only when enough parent rolls are produced in the papermaking stage can the conversion stage be started, rather than after finishing a roll of parent roll in the papermaking stage. The interval time between two stages is related to the product characteristics and the machine speed. Too long of an interval time will increase the pressure on work-in-process inventory. However, a short interval time is likely to cause frequent downtime of the machine in the conversion stage. Therefore, it is important to determine the optimal interval time in the production scheduling.
The scheduling of a tissue paper mill, which is two-stage flexible flow shop scheduling, is composed of the pulping and papermaking stage and the conversion stage. Each job is processed by a production line in the pulping and papermaking stage and a production line in the conversion stage. The energy consumption in the processing stage accounts for the vast majority of energy consumption in tissue paper mills. Product characteristics can determine the processing energy consumption. Different product characteristics mean different process parameters, and the energy consumed by machines under different process parameters varies significantly. For example, the products with high grammage will consume more energy in the papermaking stage, and vice versa. Under TOU electricity tariffs, the energy cost can be reduced in two ways: (i) arranging the job to the production line with low processing energy consumption; and (ii) arranging the job with high grammage production in the off-peak or mid-peak period and the job with low grammage production in the on-peak period. When there are two different consecutive jobs, production changeover will occur in tissue paper mills. The parameters of the machine should be adjusted to process the new later job, which will consume a lot of energy in the setup stage. In general, the setup energy cost can be reduced by arranging jobs with lower setup energy consumption processed consecutively. Under TOU electricity tariffs, the setup energy cost can be further reduced by arranging the jobs with higher setup energy consumption processed in the off-peak or mid-peak period and the jobs with lower setup energy consumption processed in the on-peak period. The parent rolls produced in the pulping and papermaking stage should be transported to the conversion stage for further processing. The transport energy consumption is related to the distance between the papermaking line and the conversion line. Similarly, the transportation distance can be shortened to decrease the transport energy cost. Moreover, the transport energy cost can be further reduced by arranging the jobs with high transportation energy consumption in the off-peak or mid-peak period and the jobs with low transportation energy consumption in the on-peak period. Tissue paper mills can decrease the energy cost by utilizing the TOU electricity pricing scheme in job scheduling. Table 1 shows the TOU electricity tariffs in Guangdong, China.

Energy Cost Modeling
According to the characteristics of energy consumption in tissue paper mills, this study mainly considers energy cost from three aspects: processing energy cost, setup energy cost, and transportation energy cost. Processing energy cost is related to processing energy consumption and electricity price. Under the TOU electricity pricing scheme, different electricity prices are set for different periods in a day. This study establishes a processing energy cost model combined with the TOU electricity pricing scheme and the processing energy consumption. The processing energy cost model is as below.
Processes 2021, 9, 274 where EP is the processing energy cost, and Equation (1) is the calculation formula of the processing energy cost. The processing energy cost is divided into four parts, namely EP1, EP2, EP3, and EP4. Equation (2) represents the processing energy cost from the start of the job to the next period. Equation (3) represents the processing energy cost from the next period of the start of the job to the previous period of the completion of the job. Equation (4) represents the processing energy cost from the previous period of the completion of the job to the completion of the job. If the job spans for several cycles, Equation (5) represents the processing energy cost of the job in several cycles. m and n are the numbers of production lines and jobs, respectively. O i,j denotes whether production line j processes job i. O i,j is set to one when production line j processes job i; otherwise, it is set to zero. CEIL is a function that rounds up to an integer. T i,j is the processing time of job i in production line j. S i,j is the start time of job i processed in production line j. Function FG has two parameters. When the first parameter is equal to the second parameter, FG is set to one; otherwise, it is set to zero. U i,j represents the power of job i in production line j, P k represents the electricity price in k period. FLOOR is a function that rounds down to an integer, K is the number of periods, MOD represents modulo operation, and FIX is the quotient operation. The setup time is up to several hours, especially in the pulping and papermaking stage, which cannot be ignored in tissue paper mills. Under the TOU electricity pricing scheme, the machine setup energy cost may be different when the setup occurs at different periods due to different energy prices in different periods. The off-peak or mid-peak period has lower setup energy cost compared with the on-peak period. Moreover, the setup may span multiple periods. The setup energy cost model is formulated as follows: where ES represents the setup energy cost, and Equation (6) is the calculation formula of the setup energy cost. The setup energy cost is divided into four parts, namely ES1, ES2, ES3, and ES4. Equation (7) represents the setup energy cost from the start of the setup to the next period. Equation (8) represents the setup energy cost from the next period of the start of the setup to the previous period of the completion of the setup. Equation (9) represents the setup energy cost from the previous period of the completion of the setup to the completion of the setup. If the setup spans for several cycles, Equation (10) represents the setup energy cost of the job in several cycles. N j represents the number of jobs processed in production line j, F i−1,j is the completion time of job i − 1 in production line j, TS i,i−1,j is the setup time in the conversion from job i − 1 to job i in production line j. US i,i−1,j is the power when job i − 1 is converted to job i in production line j.
In tissue paper mills, the parent roll produced in the papermaking line should be transported to the conversion line for further processing. However, the distance between papermaking lines and conversion lines is different. The longer the distance, the higher the energy cost of transporting the parent roll from the papermaking line to the conversion line. The transportation time in tissue paper mills is up to tens of minutes. In the production scheduling process, it is necessary to consider the processing route with short distance to reduce the transportation energy cost. Under the TOU electricity pricing scheme, the transportation energy cost varies in different time periods. Similarly, the transportation time may span multiple time periods. The transportation energy cost model of tissue paper mills is formulated as below: where ET is the transportation energy cost, Equations (12)- (14) calculate the proportions of the processing time of job i in the off-peak, mid-peak, and on-peak periods to the total processing time, respectively. Equation (15) represents the transportation energy cost from the start of the transportation to the next period. Equation (16) represents the transportation energy cost from the next period of the start of the transportation to the previous period of the complete of the transportation. Equation (17) represents the transportation energy cost from the previous period of the completion of the transportation to the completion of the transportation. If the transportation spans for several cycles, Equation (18) represents the transportation energy cost of the job in several cycles. m1 represents the number of production lines in the pulping and papermaking stage, while m2 is the number of production lines in the conversion stage. TP1, TP2, and TP3 represent the proportions of the processing time of job i in the off-peak, mid-peak, and on-peak periods to the total processing time, respectively. W i is the scale of job i. C ∑ m1 h=1 O i,h , l represents energy consumed by transporting job i from papermaking line h to conversion line l. The function Q(k, s) can judge whether the time period k belongs to period type s. The parameter s is the period type, and s values equal to 1, 2, and 3, respectively represent the off-peak, mid-peak, and on-peak period, respectively.

Energy-Efficiency Scheduling Modeling
In the pulping and papermaking stage or the conversion stage, there are multiple parallel production lines with different speeds. Each job must be first processed in a production line in the first stage and then processed in a production line in the second stage. The pulping and papermaking line or the conversion line only processes one job at a time. In this study, the job cannot be stopped once being processed in the pulping and papermaking line or the conversion line. In addition, the time interval between two stages is uncertain. The scheduling problem is the variant of a two-stage flexible flow shop scheduling problem. According to the production scheduling characteristics of tissue paper mills and the above established energy cost model, the energy-efficiency scheduling model is as below. The model in this study contains two optimization objectives, which are multi-objective optimization problems. Multi-objective optimization problems are different from single-objective optimization problems. The quality of the solution needs to be evaluated through the dominance relationship. When all the objectives of solution S1 are better than all the objectives of solution S2, it can be said that S1 dominates S2. If S1 only has some objectives better than S2, then S1 and S2 are non-dominated. If a solution is not dominated by any other solution, the solution is called a non-dominated solution. The solution to the multi-objective optimization problem is a set of non-dominated solutions.
I j,l ≥ F j,l−1 + TS l,l−1,j ∀j ∈ 1, . . . , m; ∀l ∈ 2, . . . , N j (22) Formula (19) is the optimization objectives, and there are two optimization objectives, which are makespan and energy cost. C i is the completing time of the job i in conversion line. Constraint (20) guarantees every job can and must be processed once in the pulping and papermaking line. Constraint (21) guarantees that every job can and must be processed once in the conversion line. Constraint (22) guarantees that every papermaking line or converting line can only process one job at the same time, indicating that only when the previous job and the setup are complete can the next job be started, where I j,l is the beginning time of the l-th job in production line j. Constraint (23) guarantees there is a time interval between two stages. B i,1 is the beginning time of job i in the papermaking line, and B i,2 is the beginning time of job i in the conversion line. There is a time interval between B i,1 of the length of the parent roll multiplying the number of paper layers in finished products. Formulas (25) and (26) are equality constraints.

The Proposed Optimization Algorithm
In this section, the MOEA/DTL is proposed based on two outstanding methods, namely teaching-learning-based optimization (TLBO) and MOEA/D. In order to further improve the quality of the solutions, a variable neighborhood search (VNS) algorithm is designed to be combined with the MOEA/DTL to improve the performance of the algorithm (IMOEA/DTL).

MOEA/DTL
MOEA/D uses multiple weight vectors to decompose the multi-objective optimization problem into multiple single-objective sub-optimization problems. When the weight vectors are distributed uniformly, the MOEA/D assumes Pareto optimal solutions with uniform distribution generated by combining the solutions of sub-optimization problems. Every individual is corresponding to a sub-optimization problem in MOEA/D. In each iteration of MOEA/D, each individual exchanges information with neighbors and coevolves [28]. The flow chart of the standard MOEA/D is shown in Figure 2. The input parameters are consistent with the input parameters of the proposed algorithm; refer to the input parameters of the proposed algorithm below. The initialization is The input parameters are consistent with the input parameters of the proposed algorithm; refer to the input parameters of the proposed algorithm below. The initialization is the same as the proposed algorithm; refer to the proposed algorithm below. g te x i λ i , z is the decomposition approach. There are many decomposition approaches; among them, the boundary intersection approach, Tchebycheff approach, and weighted sum approach are three popular decomposition approaches.
The basic idea of TLBO is to simulate the learning process of the class, which is completed in two phases: the teaching phase and the learning phase. Each individual learns from the teacher (the best individual is selected as the teacher) in the teaching phase. Individuals learn from each other (an individual is randomly selected) in the learning phase. The TLBO is originally designed to solve the continuous optimization problem because the formulas used to update individuals in the teaching and learning phase are continuous variable oriented. Some changes should be made when the scheduling problem is solved by the TLBO. This study modifies the TLBO by replacing the update formula by the discrete crossover operation. According to the idea of MOEA/D and TLBO, this study proposes a hybrid method (MOEA/DTL). The Figure 3 shows  The nearest T weight vectors of the i-th weight vector is calculated according to Euclidean distances. 8. Set is the neighbors of λ i with T weight vectors. 9. end for 10. Generate x 1 , x 2 , . . . , x N by a specific method. 11. FV i = F x i is set by the decomposition approach.
12. Set z = (z 1 , z 2 , . . . , z n ) T , where z j = min 1≤i≤N f j x i , f j x i is the j-th objective of solution x i . 13. while true do 14.
for i = 1 to N do 15.
if x i is not the best individual among its T neighbors 17.
An individual from its T neighbors is randomly set as teacher Th with a probability Ps, or an individual is randomly selected from EP as teacher Th with a probability 1-Ps. 18.
An individual is randomly selected from EP as Th.
A new individual y is generated by a crossover operation from x i and Th.

22.
Mutation operation is applied to y with probability Pm to generate a new individual y . 23.
Set z j = min z j , f j (y ) . If g te y λ i , z ≤ g te x i λ i , z , set x i = y and FV i = F(y ). 25.
end for 26.
for each j ∈ B(i) do 28.
if x j is not the best individual in the T neighbors of x i 29.
A new individual y is first generated by crossover operation from x i and x j .

30.
A new individual y is generated by applying the mutation operation to y with probability Pm. 31. else 32.
Continue to the next iteration. 33.
set z j = min z j , f j (y ) . Update EP: a new non-dominated set is selected from EP and the new population as the new EP, and get rid of the duplicate solution. If the size of EP is smaller than the setting value, some randomly generated individuals (EP ) are added to EP until the size of EP is equal to the setting value. 43.
for each i = 1 to size(EP) do 44.
An individual x i is randomly selected from EP − EP , and an individual x j is ranbdomly selected from EP − EP with a probability Pc, or an individual x j is selected from EP with a probability 1 − Pc.

45.
A new y is generated by crossover operation from x i and x j . The weight vectors can determine the quality of solution obtained by MOEA/DTL. With the uniformly distributed weight vectors, the solution obtained by MOEA/DTL is more evenly distributed and closer to the Pareto optimal solution. Thus, the method of generating some weight vectors with even distribution is of great importance to MOEA/DTL. In this study, a number of weight vectors are generated as follows: Step Step (2) Generate weight vectors λ = λ 1 , λ 2 , . . . , λ N by randomly selecting N vectors from γ.
The parameter Ps is one of the important factors affecting the convergence speed of the algorithm. If Ps is too small, it means that most of the teachers are from the global optimal individual, and it is easy to fall into the local optimal value. If Ps is too large, it means that most of the teachers come from the best individuals in the neighborhood, and the convergence speed will be slower. The value of Ps between 0.7 and 0.9 is reasonable. The decomposition approach is also important for MOEA/DTL. The makespan and the energy cost are converted into an objective using the weighted sum approach. The makespan and the energy cost are the two different scalar objectives. Therefore, the decomposition method must have the function of normalization. The decomposition approach is as below: where g i is the scalar optimization objective of solution x i , f max j is the maximum of the j-th objective, and f min j is the minimum of the j-th objective.

Variable Neighborhood Search
With the advantages of high efficiency and simple implementation, the Variable Neighborhood Search (VNS) algorithm has been successfully applied in many engineering fields [29,30]. The VNS algorithm uses the neighborhood structure of different actions to search alternately, thus achieving a good balance between concentration and evacuation. Firstly, the VNS is carried out in a small neighborhood. If the quality of the solution cannot be improved, VNS will be carried out in a larger neighborhood. When a better solution is found, the VNS is conducted in a smaller neighborhood again. The cycle is repeated until the end of the algorithm. The VNS algorithm is simple in principle, easy to implement, and has good optimization performance. The procedure of the VNS algorithm used in this paper is shown in Algorithm 2: for each j ∈ [1, 2, . . . , M1] do 6.
Based on the current solution, the first neighborhood structure S 1 is used to search for better solutions than the current solution, and it means that the solutions found can dominate the current solution. 7.
end for 8.
Based on the current solution, the second neighborhood structure S 2 is used to search for better solutions than the current solution, and it means that the solutions found can dominate the current solution. 10.
end for 11.
Based on the current solution, the third neighborhood structures S 3 is used to search for better solutions than the current solution, and it means that the solutions found can dominate the current solution. 13. end for 14.
If the stopping criteria are not satisfied, continue; otherwise, stop and output the best solution.

end for
The core of the VNS algorithm lies in the design of neighborhood structure set. In the production scheduling, the neighborhood structure of VNS can be obtained by exchanging the order of two jobs or inserting one job before another. This study designs three neighborhood structures, which are S 1 , S 2 , and S 3 . The neighborhood solutions of S 1 are generated as follows: select two jobs randomly and exchange their positions to generate a new neighborhood solution. The neighborhood solutions of S 2 are generated as follows: select three jobs randomly and then disrupt their order to generate five new neighborhood solutions. The neighborhood solutions of S 3 are generated as follows: select four jobs randomly and then disrupt their order to generate 23 new neighborhood solutions. The parameters M, M1, M2, and M3 affect the search range. A large value means a large search range, but at the same time, it reduces the operating efficiency of the algorithm. On the contrary, the operating efficiency of the algorithm will improve, but it means that the search range becomes smaller. In order to balance the two, it is reasonable to set their value within 10.

IMOEA/DTL
Despite the outstanding global search capability of MOEA/DTL, it is likely to miss the locally better solution because MOEA/DTL does not search the neighborhood of the solution in detail in each iteration process. The MOEA/DTL is combined with VNS (IMOEA/DTL) to enhance its local search ability and performance. In IMOEA/DTL, the solution obtained in each iteration of MOEA/DTL is searched locally using the VNS algorithm to find a better solution. This study carries out VNS after the first update of EP, and VNS is carried out as Algorithm 3.

Case Study
In order to demonstrate the feasibility of the energy-efficiency scheduling model under TOU electricity tariffs for tissue paper mills and the superiority of the proposed IMOEA/DTL, three popular multi-objective optimization algorithms, namely, MOEA/D-MR [31], Non-dominated Sorting Genetic Algorithm II (NSGA-II), and Strength Pareto Evolutionary Algorithm 2 (SPEA2), are compared with the IMOEA/DTL in eight scheduling problems. The energy-saving potential of the IMOEA/DTL is evaluated by comparing IMOEA/DTL with workers in practice. The whole experiments are implemented in MAT-LAB 2015b.

Experiment Data
The data are collected from a tissue paper mill in Guangdong, China. The data are collected from the Jiangmen branch of Vinda Paper (China) Co., Ltd., Jiangmen City, Guangdong Province, China. The data collection time is from 2017 to 2018. The scheduling data and product process data are mainly collected from the production scheduling department. These data come from the enterprise resource planning system. Machine data and energy consumption data are automatically collected by the manufacturing execution system. There are six production lines in the pulping and papermaking stage and seven production lines in the conversion stage. A total of 16 scheduling instances are generated randomly, the numbers and names of jobs in each scheduling problem are listed in Table 2. In each scheduling problem, the size of each job is generated randomly and obeys uniform distribution in the range of [10,000, 600,000]. In the tissue paper mill, the speed of the production line varies with the type of products, but each production line has a maximum speed. The maximum speeds in the pulping and papermaking stage and the conversion stage are shown as Table 3. Similarly, the power of each production line also varies with the product type. The maximum power of the pulping and papermaking stage and the conversion stage is shown Table 4. The setup energy consumption is dependent on the setup time and the setup power. In this study, the setup time of the pulping and papermaking stage and the conversion stage is randomly generated based on the real situation. The setup time of the pulping and papermaking stage shows uniform distribution in the range of [10,100], while the setup time of the conversion stage shows uniform distribution in the range of [10,60]. Table 5 lists the setup power of the pulping and papermaking stage and the conversion stage. The transportation energy consumption is related to the distance between production lines. The transportation energy consumption per unit product is shown in Table 6.

Performance Metrics
The performance of the IMOEA/DTL, MOEA/D-MR, NSGA-II, and SPEA2 is evaluated using two metrics, which are Hypervolume indicator (HV-metric) and Set Coverage (C-metric). C-metric reflects the dominance between two non-dominated solution sets.  C(B, A), then the quality of the solution A is better than that of B.
The HV-metric of a non-dominated set can be represented by the size of portion of objective space that is dominated by those solutions collectively and bounded above by a reference point. The HV-metric can not only assess the proximity of the non-dominated set but also can assess the diversity of the non-dominated set. The HV-metric needs to be given a reference point, where we use the maximum value of the solutions obtained by all algorithms. When the HV-metric has a larger value, the approximate non-dominated solution set is better.
Moreover, considering the randomness of the evolutionary algorithm, this study carries out the Wilcoxon signed-rank test to detect whether the results obtained by different algorithms are different significantly. The significance level is 0.05. A p-value of the Wilcoxon signed-rank test less than 0.05 indicates significant differences between the two algorithms. Table 7 [31][32][33][34][35]. Moreover, considering the randomness of the evolutionary algorithm, three algorithms are run ten times for each scheduling problem, and the result is the average of ten times.  Table 8 describes the mean C-metric of the solutions obtained by IMOEA/DTL, MOEA/D-MR, SPEA2, and NSGA-II. W, X, Y, and Z represent the solutions obtained by MOEA/D-MR, NSGA-II, SPEA2, and IMOEA/DTL, respectively. From Table 8, it can be seen that C(Z, W) is larger than C(W, Z), C(Z, X) is larger than C(X, Z), and C(Z, Y) is larger than C(Y, Z) in 16 scheduling problems. The values of C(Z, W), C(Z, X) and C(Z, Y) are close to one in the eight scheduling problems. Meanwhile, the values of C(W, Z), C(X, Z), and C(Y, Z) are almost equal to zero. The C-metric indicates that other methods dominate the few solutions obtained by IMOEA/DTL. However, some solutions of IMOEA/DTL dominate all solutions obtained by MOEA/D-MR, NSGA-II, and SPEA2. As shown in Table 8, the solutions obtained by IMOEA/DTL have better performance than those obtained by MOEA/D-MR, NSGA-II, and SPEA2 in C-metric. Table 9 shows the p-value in the Wilcoxon signed-rank test of the C-metric between two algorithms. P1 represents the p-value of the Wilcoxon signed-rank test of the C-metric between IMOEA/DTL and MOEA/D-MR, P2 represents the p-value of the Wilcoxon signed-rank test of the C-metric between IMOEA/DTL and NSGA-II, and P3 represents the p-value of the Wilcoxon signed-rank test of the C-metric between IMOEA/DTL and SPEA2. P1, P2, and P3 are less than 0.05 in 16 scheduling problems, indicating significant differences between IMOEA/DTL and MOEA/D-MR, IMOEA/DTL and NSGA-II, IMOEA/DTL and SPEA2. Tables 8 and 9 indicate that IMOEA/DTL is superior to MOEA/D-MR, NSGA-II, and SPEA2 in terms of C-metric.  Table 9. The p-value of the Wilcoxon signed-rank test of the C-metric between two algorithms.

P1 P2 P3
Job1 0.0020 0.0020 0.0020 Job2 0.0020 0.0020 0.0020 Job3 0.0020 0.0020 0.0020 Job4 0.0020 0.0020 0.0020 Job5 0.0020 0.0020 0.0020 Job6 0.0020 0.0020 0.0020 Job7 0.0020 0.0020 0.0020 Job8 0.0020 0.0020 0.0020 Job9 0.0020 0.0020 0.0020 Job10 0.0020 0.0020 0.0020 Job11 0.0020 0.0020 0.0020 Job12 0.0020 0.0020 0.0020 Job13 0.0020 0.0020 0.0020 Job14 0.0020 0.0020 0.0020 Job15 0.0020 0.0020 0.0020 Job16 0.0020 0.0020 0.0020 Table 10 lists the average value of HV-metric of the solutions obtained by IMOEA/DTL, MOEA/D-MR, SPEA2, and NSGA-II. The solutions obtained by IMOEA/DTL have a larger average HV-metric than those obtained by MOEA/D-MR, SPEA2, and NSGA-II in all scheduling problems, indicating that the proximity and diversity of the solutions of IMOEA/DTL are best in all solutions. As shown in Table 10, the solutions obtained by MOEA/D-MR, NSGA-II, and SPEA2 are inferior to that obtained by IMOEA/DTL in terms of HV-metric. The p-value of the Wilcoxon signed-rank test of the HV-metric between two algorithms is shown in Table 11. P1 represents the p-value of the Wilcoxon signed-rank test of the HV-metric between IMOEA/DTL and MOEA/D-MR, P2 represents the p-value of the Wilcoxon signed-rank test of the HV-metric between IMOEA/DTL and NSGA-II, and P3 represents the p-value of the Wilcoxon signed-rank test of the HV-metric between IMOEA/DTL and SPEA2. P1, P2, and P3 are less than 0.05, demonstrating the significant differences between IMOEA/DTL and MOEA/D-MR, IMOEA/DTL and NSGA-II, and IMOEA/DTL and SPEA2. Tables 10 and 11 indicate that IMOEA/DTL is superior to MOEA/D-MR, NSGA-II, and SPEA2 in terms of HV-metric.  In order to intuitively show the scheduling schemes obtained by different algorithms, some Gantt charts of scheduling schemes obtained by different algorithms are provided, as shown in Figures 4-6. Each scheduling scheme corresponds to a solution in the Pareto front. We only provide Gantt charts of Job1, Job8, and Job16, and the Gantt charts for other scheduling problems are similar to these three. In the Gantt charts below, we can see how the jobs are organized and processed for the scheduling scheme obtained by each algorithm, and it can be seen that the scheduling scheme obtained by IMOEA/DTL has a smaller makespan. Not only that, in all the following scheduling schemes, IMOEA/DTL has the smallest energy cost.

Job12
0.0020 0.0020 0.0020 Job13 0.0020 0.0020 0.0020 Job14 0.0020 0.0020 0.0020 Job15 0.0020 0.0020 0.0020 Job16 0.0020 0.0020 0.0020 In order to intuitively show the scheduling schemes obtained by different algorithms, some Gantt charts of scheduling schemes obtained by different algorithms are provided, as shown in Figures 4-6. Each scheduling scheme corresponds to a solution in the Pareto front. We only provide Gantt charts of Job1, Job8, and Job16, and the Gantt charts for other scheduling problems are similar to these three. In the Gantt charts below, we can see how the jobs are organized and processed for the scheduling scheme obtained by each algorithm, and it can be seen that the scheduling scheme obtained by IMOEA/DTL has a    Table 12 shows the average makespan and average energy cost obtained by MOEA/D-MR, NSGA-II, SPEA2, and IMOEA/DTL. Table 12 indicates that both the average makespan and average energy cost obtained by IMOEA/DTL are lower than those obtained by MOEA/D-MR, NSGA-II, and SPEA2. Therefore, IMOEA/DTL can provide the solution with lower makespan and energy cost than MOEA/D-MR, NSGA-II, and SPEA2. This study adopts a multi-objective algorithm, and a non-dominated solution set is obtained. The solutions in the non-dominated solution set have lower energy cost or makespan, but they are mutually exclusive. Decision makers can select an optimal solution on the basis of actual situations. For instance, the solution with a low makespan can be selected in the case of a heavy production task, which can improve the utilization efficiency of the workshop and ensure that the task can be completed on schedule. On the contrary, the decision-makers can choose the optimal solution with low energy cost to reduce the production cost.   Table 12 shows the average makespan and average energy cost obtained by MOEA/D-MR, NSGA-II, SPEA2, and IMOEA/DTL. Table 12 indicates that both the average makespan and average energy cost obtained by IMOEA/DTL are lower than those obtained by MOEA/D-MR, NSGA-II, and SPEA2. Therefore, IMOEA/DTL can provide the solution with lower makespan and energy cost than MOEA/D-MR, NSGA-II, and SPEA2. This study adopts a multi-objective algorithm, and a non-dominated solution set is obtained. The solutions in the non-dominated solution set have lower energy cost or makespan, but they are mutually exclusive. Decision makers can select an optimal solution on the basis of actual situations. For instance, the solution with a low makespan can be selected in the case of a heavy production task, which can improve the utilization efficiency of the workshop and ensure that the task can be completed on schedule. On the contrary, the decision-makers can choose the optimal solution with low energy cost to reduce the production cost. In China, most of the scheduling jobs of the tissue paper mills are completed by workers based on their experience. The quality of the scheduling depends on the experience of the workers. In this study, a real scheduling problem is solved by IMOEA/DTL, and the solution is compared with that obtained by workers. The solution obtained by workers has the makespan of 17,035 and the energy cost of 1,368,497. Table 13 indicates that the average makespan obtained by IMOEA/DTL is about 7.5% less than that obtained by the worker. The average energy cost obtained by IMOEA/DTL is about 1.5% less than that obtained by the worker. From Table 13, it can be seen that the minimum energy cost is 1,333,494, and the corresponding makespan is 16,754, which reduces by 2.6% and 1.6%, respectively. The minimum makespan is 15,367, and the corresponding energy cost is 1,361,331, which decrease by 9.8% and 0.5%, respectively. The largest energy cost is 1,375,361, which is 0.5% higher than the cost obtained by the worker, but the corresponding makespan is reduced by 9.5%. The largest makespan is 18,675, which is 8.8% higher than the cost obtained by the worker, but the corresponding energy cost is reduced by 2.5%. It can be concluded that IMOEA/DTL can reduce both the energy cost and the makespan. Although the maximum makespan and the maximum energy cost are larger than those obtained by the worker, the corresponding energy cost and makespan are smaller. Figure 7 shows the scatter plot for the solutions obtained by IMOEA/DTL and the worker. The energy cost obtained by IMOEA/DTL is lower than that obtained by the worker. The solutions of IMOEA/DTL with lower energy cost than that obtained by the worker account for about 98.56% of the total solutions of IMOEA/DTL. While the solutions of IMOEA/DTL with lower makespan than that obtained by the worker account for about 97.12% of the total solutions of IMOEA/DTL. About 95.68% of the total solutions obtained by IMOEA/DTL have lower energy cost and makespan than those obtained by the worker.   In practical applications, the construction of basic data is a particularly critical issue, especially setup-related data, which requires a lot of statistical data. We suggest that a part of the key data can be constructed through the manufacturing execution system, and through the manufacturing execution system, we can know the status of machine, personnel, and materials in real time, and thus make timely responses to emergencies in the production process. In addition, we recommend using commercial production scheduling software and then integrating the model into commercial software. Compared with developing a scheduling system from scratch, this has a higher success rate.

Conclusions
This study explored a novel energy-saving approach by scheduling under TOU electricity tariffs for tissue paper mills. First of all, an energy cost model was established based on TOU electricity tariffs. The energy cost consists of transportation energy cost, set-up energy cost, and processing energy cost. Second, an energy-efficiency scheduling model was built on the basis of the scheduling model and the energy cost model. Finally, the model was solved by proposing a novel IMOEA/DTL method. The experiment made a comparison between the proposed method and two popular multi-objective optimization methods, MOEA/D-MR, NSGA-II, and SPEA2. The experiment results demonstrate that the proposed method performs better than NSGA-II and SPEA2 in terms of C-metric and HV-metric. Moreover, the proposed method has been proven to have an energy-saving ratio of 2.6% through a real scheduling problem. Considering the similarity between the production processes in tissue paper mills and other papermaking mills, the results are of great reference significance for other papermaking mills.
However, this study implemented the algorithms by MATLAB, and the speed of the algorithms should be further improved. Therefore, the speed of the algorithms will be improved in future research.