# Optimizing Production Scheduling of Steel Plate Hot Rolling for Economic Load Dispatch under Time-of-Use Electricity Pricing

Authors: Mao Tan; Hua-li Yang; Bin Duan; Yong-xin Su; Feng He
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Time-of-Use (TOU) electricity pricing provides an opportunity for industrial users to cut electricity costs. Although many methods for Economic Load Dispatch (ELD) under TOU pricing in continuous industrial processing have been proposed, there are still difficulties in batch-type processing, since power load units are not directly adjustable and nonlinearly depend on production planning and scheduling. In this paper, for hot rolling, a typical batch-type and energy intensive process in steel industry, a production scheduling optimization model for ELD is proposed under TOU pricing, in which the objective is to minimize electricity costs while considering penalties caused by jumps between adjacent slabs. A NSGA-II based multiobjective production scheduling algorithm is developed to obtain Pareto optimal solutions, and then TOPSIS based multicriteria decision-making is performed to recommend an optimal solution to facilitate field operation. Experimental results and analyses show that the proposed method cuts electricity costs in production, especially in case of allowance for penalty score increase in a certain range. Further analyses show that the proposed method has effect on peak load regulation of power grid.


---

## Body

## 1. Introduction

Time-of-Use (TOU) electricity pricing, a practical demand response program implemented by many power suppliers to improve the peak load regulation ability of power grid, provides an opportunity for electricity users to implement Economic Load Dispatch (ELD), that is, cut electricity costs by reducing power loads during on-peak periods and shifting loads from on-peak to off-peak periods.Unlike conventional energy conservation to reduce absolute energy consumption, optimizing electricity costs under TOU pricing means that industrial users adjust their production schedule to avoid on-peak time periods, which will have significant effect on cutting electricity costs. In recent years, ELD under TOU pricing has become a hot area. Shrouf et al. [1] proposed a single machine scheduling problem, in which each time period has an associated price and the objective is to minimize electricity costs while considering traditional scheduling performance measures. Fang et al. [2] also considered job scheduling on a single machine to minimize total electricity costs under TOU pricing and proposed the algorithms for uniform-speed and speed-scalable machine environments, respectively. Mitra et al. [3] formulated mixed integer linear programming for continuous industrial processing, which allows optimal production planning, and provided a case study for time horizon of one week and hourly changing electricity prices. Furthermore, they improved the model with integration of operational and strategic decision-making [4]. Ashok [5] presented a theoretical model for batch-type load processing and proposed an integer programming method to reschedule their operations to reduce electricity costs under time-varying electricity price, but the model is an abstract theoretical model and is difficult to be applied to production directly. Wang et al. [6] proposed an optimization model to minimize electricity costs for steel plant, in which both power generation scheduling and batch production scheduling were considered; although the model has been believed to be effective under TOU pricing, the results cannot always be optimal because the production load units are determined by fixed production planning and scheduling.The above analyses motivate the potential for more benefits by ELD under TOU pricing in hot rolling production scheduling. Until now, most of the related literatures focused on specific part of the problem or the abstract simplified problem; thus there are still difficulties, since the power load units are not directly adjustable and nonlinearly depend on the results of production planning and scheduling.Hot rolling, a typical batch-type and energy intensive process in steel production with characteristics of strong schedulability, has become an important aspect of production organization and energy saving [7]. The general process flow of hot rolling production is illustrated in Figure 1. Hot rolling is mainly organized and carried out by batch scheduling program in steel mill, the primary task of which is arranging and sequencing slabs into rolling units to smooth jumps in width, gauge, and hardness between adjacent slabs; all of these will directly affect product quality. Hot rolling production scheduling has attracted attention from academia and industry for a long time. An early method proposed by Kosiba et al. treated steel production scheduling as a discrete event sequencing problem and thus formulated it as a traveling salesman problem [8]. Lopez et al. [9] formulated the problem as a generalized prize collecting traveling salesman problem with multiple conflicting objectives and constraints and proposed a heuristic tabu search method to determine good approximate solutions. Tang and Wang [10] proposed a modified genetic algorithm based on the multiple travelling salesman problem. Chen et al. [11] formulated the problem as a nonlinear integer programming model, and later it is corrected by Kim [12] and changed to a linear programming model. Furthermore, Alidaee and Wang [13] proposed a corrected integer programming formulation and reduced the quantity of variables. Nevertheless, most of proposed models are single objective or transformed models based on weighted-sum approach. Jia et al. [14] formulated the problem as a multiobjective vehicle routing problem with double time windows and proposed a decomposition-based hierarchical optimization algorithm to solve it. Soon after, they proposed a P-MMAS algorithm to solve the problem; multicriteria decision-making is performed to recommend the optimal solution from the Pareto frontier [15]. Moon et al. [16] proposed a production scheduling model with time-dependent and machine-dependent electricity cost, in which makespan was considered by using the weighted-sum objective but batch sizing was not considered, which is obviously simpler than batch scheduling problem. Because of complexity of batch sizing problem, Sarakhsi et al. [17] proposed a hybrid algorithm of scatter search and Nelder-Mead algorithms to improve the performance of solving algorithm.Figure 1
A process flow diagram of the hot rolling production procedure.Due to high energy consumption and rising energy costs in hot rolling production [18], energy saving has also been considered combined with the traditional objective mentioned above. As is shown in Figure 1, slabs are heated to high temperature before being rolled; the total energy consumed in heating is affected by batch schedule. Since Direct Hot Charge Rolling (DHCR) has significant benefits on energy cost, great efforts have been made to improve the ratio of DHCR while performing batch scheduling [19, 20]. Besides that, optimization of rolling schedule by adjusting thickness reduction ratio of slabs between the rolling passes, another way to reduce power consumption that used to drive rolling motor, has also been proposed [21–23].As mentioned previously, most methods of hot rolling production scheduling concentrate on internal production organization. Although some technical means have been proposed and applied to achieve energy conservation, their potential would be exhausted due to equipment and technology constraints. In this context, methods utilizing favorable external environments should be explored for energy saving. TOU pricing provides an opportunity to reduce electricity costs, but until now there are few published papers to implement ELD under TOU pricing for hot rolling production.This paper considers the Hot Rolling Production Scheduling Problem (HRPSP) as a mixture of batch scheduling problem and time-dependent job-shop scheduling problem. The rolling units, modeled as power load units, are planned and scheduled according to TOU prices. Primary objective of the proposed model is to minimize electricity costs while considering the traditional objective to minimize penalties caused by jumps between adjacent slabs. A multiobjective optimization model and corresponding solving algorithm are additionally proposed.The rest of this paper is classified as follows: in Section2, characteristics of the problem and opportunities under TOU pricing are presented, and a mathematical model with objective to minimize electricity costs in production is formulated. A multiobjective optimization algorithm is developed in Section 3 to solve the problem. Section 4 is dedicated to the experimental procedure and results to evaluate the proposed method; also the peak load regulation effect and robustness of the proposed method are further discussed. Finally, conclusion and future research planning are given in Section 5.


## 2. Problem Description and Formulation

HRPSP is an extremely complex problem which has significant influence on product quality, production efficiency, and energy consumption. In this paper, we study the Hot Rolling Batch Scheduling Problem (HRBSP) combined with the job-shop scheduling problem (JSP), where HRBSP focuses on how rolling units are organized and the JSP concentrates on when the rolling units are processed.


### 2.1. Problem Description

Hot rolling batch scheduling is a key process in hot rolling. The task of HRBSP, as is depicted in Figure2, is to select, group, and sequence slabs into rolling units with the constraints of production capacity and rolling rules. Each rolling unit has a coffin-shaped width profile consisting of a warming-up section and a coming-down section. In the previous section slabs are arranged from narrow to wide to warm up the rolls, and in the later section slabs are scheduled with decreasing width to avoid marking the coils surface. The major part of a rolling unit is the coming-down section, in which the quality of rolling mainly depends on the sequence of slabs. In most cases, the warming-up section is trivial and can be determined manually.Figure 2
Diagrammatic sketch of batch scheduling: (a) rolling unit and (b) rolling batch.
(a)
(b)Several constraints restrict the scheduling, the most important one of which is to smooth jumps in width, gauge, and hardness between adjacent slabs. Other constraints, such as cumulative rolling length of slabs in a rolling unit and continuous rolling length of slabs with same width, are also considered to ensure product quality and production capability.Because hot rolling is a key energy intensive process in steel industry, many approaches, such as optimization of batch scheduling with the objective of improving DHCR ratio and optimization of reduction schedule, have been proposed to achieve energy saving. In smart grid, TOU electricity pricing, which is one of the most commonly implemented demand response programs [24], provides a new opportunity for steel mill to achieve ELD in hot rolling production, which means cutting costs by shifting loads according to the electricity price.As is shown in Figure3, a whole day is partitioned into four types of periods based on the price of electricity: on-peak, mid-peak, flat-peak, and off-peak periods. We can see that the power cost for each rolling unit, which not only is determined by the quantity of power demand but also is dependent on the corresponding electricity pricing, should be accumulated piecewise during the processing time.Figure 3
Relationship between production scheduling and electricity costs under TOU pricing.Compared with flat electricity pricing, the objective of ELD under TOU pricing is to minimize total power cost, including charges for power consumed from shifting loads. In this paper, we assume that rolling units can be scheduled freely; therefore no operating costs from load shifting are included. Consequently, rolling production is encouraged during off-peak periods and discouraged during on-peak periods. In addition, we should know that the scheduling on fixed jobs is not always optimal, so the scheduled jobs (i.e., the rolling units obtained by hot rolling batch scheduling) should be created and associated with their operation time. Finally, the problem is turned into optimal scheduling for minimizing the electricity costs determined by batch scheduling solution and job-shop scheduling solution under specified electricity pricing, while the traditional objective that smoothing changes between adjacent slabs should not be ignored to ensure product quality.


### 2.2. Mathematical Formulation

We interpret the basic model of the HRBSP as a vehicle routing problem (VRP), which is a classical combinatorial optimization problem. In the model, it can be considered that each rolling unit is a vehicle within limited capacity and each slab is a customer that should be visited at most once. Suppose that there aren slabs to be scheduled into m rolling units; the objective of the problem is to determine m routes (rolling units) to minimize the total distance traveled (penalties caused by jumps between adjacent slabs).The variables used in formulation are listed as follows:N: a set of slabs; N={1,2,…,n}M: a set of rolling units; M={1,2,…,m}T: a set of time periods; T={1,2,…,t}πj: electricity price during time period jWi: power demand of slab i during rolling procedurelj: rolling length of slab jpi: processing time for slab iPij: the penalty for rolling slab j immediately after slab i, where Pij=pijw+pijg+pijh; pijw, pijg, and pijh, respectively, represent the contribution due to width, gauge, and hardnesssij: binary variable with value of 1 if the widths of slab i and j are the same; otherwise it is 0tsi: processing start time of slab iL: lower bound of the cumulative length of slabs scheduled in a single rolling unitU: upper bound of the cumulative length of slabs scheduled in a single rolling unitR: upper bound of the cumulative length of slabs with the same width in a single rolling unitTS: total time that can be allocated for productionFive decision expressions are defined to identify the scheduling solution as follows:(1)xijk=1if slabjis immediately after slabiin unitk,0otherwise,yik=1if slabiis scheduled in rolling unitk,0otherwise,rijk=1if slabjis rolled after slabiin rolling unitk,0otherwise,dij=1if slabiis processed in time periodsj,0otherwise.v
i, a positive integer or 0, is a variable to indicate the idle time allocated to rolling unit i before production.Note that production efficiency may not always be the only one target in engineering, especially when production capacity is abundant; then the target of our model is to minimize electricity costs on the premise of processing all products in given time horizon. According to basic VRP model combined with consideration of relationship between slab processing sequence and processing time as shown in Figure3, we formulate the hot rolling production optimization problem as(2)minf1=∑k∈M∑i∈N∑j∈NPij·xijk(3)minf2=∑j∈Tπj·∑i∈NWi·dij(4)s.t.∑i∈Nxijk=yjk,j∈N,k∈M(5)∑j∈Nxijk=yik,i∈N,k∈M(6)∑k∈Myik=1,i∈N(7)∑i∈Nrijk·sij·lj≤R,j∈N,k∈M(8)L≤∑i∈Nyik·li≤U,k∈M(9)0≤∑i∈Mvi≤TS-∑i∈Npi(10)∑k∈Mrijk≤1,i∈N,j∈N(11)xijk≤rijk,i∈N,j∈N,k∈M(12)rijk≤yik,i∈N,j∈N,k∈M(13)rijk≤yjk,i∈N,j∈N,k∈M,where objective f1 is the traditional objective to ensure product quality, which means to minimize the total penalties caused by jumps between adjacent slabs, and objective f2 means to minimize the total electricity costs in hot rolling production, in which dij can be further formulated as(14)dij=1if∑α<jλα≤tsi<∑α≤jλα,0otherwise,where the condition correspond to dij=1 means that slab i is processed in time period j. Note that variable tsi not only is determined by which rolling units the slab is scheduled in but also depended on the processing time of previous slabs and the idle time allocated for rolling units; then it can be expressed as(15)tsi=∑δ∈Myiδ·∑β<δ∑α∈Nyαβ·pα+∑α≤δvα+∑β=δ∑α∈Nrαiβ·pα,where δ is a traversal variable to search the rolling unit that slab i is allocated in; expression in brackets indicates the cumulative time before processing slab i. If slab i is not allocated in rolling unit δ, the expression in brackets would be ignored because yiδ=0.Constraints (4) and (5) specify the sequence of slabs in a rolling unit. Constraint (6) ensures that each slab can be scheduled only once. Constraint (7) restricts the cumulative length of continuously rolled slabs with the same width in each rolling unit. Constraint (8) indicates rolling mill production capacity, which restricts the lower and upper bounds of cumulative length of slabs in each rolling unit. Constraint (9) means that the total idle time allocated for rolling units cannot be greater than margin of production capability. Constraints (10)–(13) restrict the value of rijk, xijk, and yik according to their logical relationship.


## 2.1. Problem Description

Hot rolling batch scheduling is a key process in hot rolling. The task of HRBSP, as is depicted in Figure2, is to select, group, and sequence slabs into rolling units with the constraints of production capacity and rolling rules. Each rolling unit has a coffin-shaped width profile consisting of a warming-up section and a coming-down section. In the previous section slabs are arranged from narrow to wide to warm up the rolls, and in the later section slabs are scheduled with decreasing width to avoid marking the coils surface. The major part of a rolling unit is the coming-down section, in which the quality of rolling mainly depends on the sequence of slabs. In most cases, the warming-up section is trivial and can be determined manually.Figure 2
Diagrammatic sketch of batch scheduling: (a) rolling unit and (b) rolling batch.
(a)
(b)Several constraints restrict the scheduling, the most important one of which is to smooth jumps in width, gauge, and hardness between adjacent slabs. Other constraints, such as cumulative rolling length of slabs in a rolling unit and continuous rolling length of slabs with same width, are also considered to ensure product quality and production capability.Because hot rolling is a key energy intensive process in steel industry, many approaches, such as optimization of batch scheduling with the objective of improving DHCR ratio and optimization of reduction schedule, have been proposed to achieve energy saving. In smart grid, TOU electricity pricing, which is one of the most commonly implemented demand response programs [24], provides a new opportunity for steel mill to achieve ELD in hot rolling production, which means cutting costs by shifting loads according to the electricity price.As is shown in Figure3, a whole day is partitioned into four types of periods based on the price of electricity: on-peak, mid-peak, flat-peak, and off-peak periods. We can see that the power cost for each rolling unit, which not only is determined by the quantity of power demand but also is dependent on the corresponding electricity pricing, should be accumulated piecewise during the processing time.Figure 3
Relationship between production scheduling and electricity costs under TOU pricing.Compared with flat electricity pricing, the objective of ELD under TOU pricing is to minimize total power cost, including charges for power consumed from shifting loads. In this paper, we assume that rolling units can be scheduled freely; therefore no operating costs from load shifting are included. Consequently, rolling production is encouraged during off-peak periods and discouraged during on-peak periods. In addition, we should know that the scheduling on fixed jobs is not always optimal, so the scheduled jobs (i.e., the rolling units obtained by hot rolling batch scheduling) should be created and associated with their operation time. Finally, the problem is turned into optimal scheduling for minimizing the electricity costs determined by batch scheduling solution and job-shop scheduling solution under specified electricity pricing, while the traditional objective that smoothing changes between adjacent slabs should not be ignored to ensure product quality.


## 2.2. Mathematical Formulation

We interpret the basic model of the HRBSP as a vehicle routing problem (VRP), which is a classical combinatorial optimization problem. In the model, it can be considered that each rolling unit is a vehicle within limited capacity and each slab is a customer that should be visited at most once. Suppose that there aren slabs to be scheduled into m rolling units; the objective of the problem is to determine m routes (rolling units) to minimize the total distance traveled (penalties caused by jumps between adjacent slabs).The variables used in formulation are listed as follows:N: a set of slabs; N={1,2,…,n}M: a set of rolling units; M={1,2,…,m}T: a set of time periods; T={1,2,…,t}πj: electricity price during time period jWi: power demand of slab i during rolling procedurelj: rolling length of slab jpi: processing time for slab iPij: the penalty for rolling slab j immediately after slab i, where Pij=pijw+pijg+pijh; pijw, pijg, and pijh, respectively, represent the contribution due to width, gauge, and hardnesssij: binary variable with value of 1 if the widths of slab i and j are the same; otherwise it is 0tsi: processing start time of slab iL: lower bound of the cumulative length of slabs scheduled in a single rolling unitU: upper bound of the cumulative length of slabs scheduled in a single rolling unitR: upper bound of the cumulative length of slabs with the same width in a single rolling unitTS: total time that can be allocated for productionFive decision expressions are defined to identify the scheduling solution as follows:(1)xijk=1if slabjis immediately after slabiin unitk,0otherwise,yik=1if slabiis scheduled in rolling unitk,0otherwise,rijk=1if slabjis rolled after slabiin rolling unitk,0otherwise,dij=1if slabiis processed in time periodsj,0otherwise.v
i, a positive integer or 0, is a variable to indicate the idle time allocated to rolling unit i before production.Note that production efficiency may not always be the only one target in engineering, especially when production capacity is abundant; then the target of our model is to minimize electricity costs on the premise of processing all products in given time horizon. According to basic VRP model combined with consideration of relationship between slab processing sequence and processing time as shown in Figure3, we formulate the hot rolling production optimization problem as(2)minf1=∑k∈M∑i∈N∑j∈NPij·xijk(3)minf2=∑j∈Tπj·∑i∈NWi·dij(4)s.t.∑i∈Nxijk=yjk,j∈N,k∈M(5)∑j∈Nxijk=yik,i∈N,k∈M(6)∑k∈Myik=1,i∈N(7)∑i∈Nrijk·sij·lj≤R,j∈N,k∈M(8)L≤∑i∈Nyik·li≤U,k∈M(9)0≤∑i∈Mvi≤TS-∑i∈Npi(10)∑k∈Mrijk≤1,i∈N,j∈N(11)xijk≤rijk,i∈N,j∈N,k∈M(12)rijk≤yik,i∈N,j∈N,k∈M(13)rijk≤yjk,i∈N,j∈N,k∈M,where objective f1 is the traditional objective to ensure product quality, which means to minimize the total penalties caused by jumps between adjacent slabs, and objective f2 means to minimize the total electricity costs in hot rolling production, in which dij can be further formulated as(14)dij=1if∑α<jλα≤tsi<∑α≤jλα,0otherwise,where the condition correspond to dij=1 means that slab i is processed in time period j. Note that variable tsi not only is determined by which rolling units the slab is scheduled in but also depended on the processing time of previous slabs and the idle time allocated for rolling units; then it can be expressed as(15)tsi=∑δ∈Myiδ·∑β<δ∑α∈Nyαβ·pα+∑α≤δvα+∑β=δ∑α∈Nrαiβ·pα,where δ is a traversal variable to search the rolling unit that slab i is allocated in; expression in brackets indicates the cumulative time before processing slab i. If slab i is not allocated in rolling unit δ, the expression in brackets would be ignored because yiδ=0.Constraints (4) and (5) specify the sequence of slabs in a rolling unit. Constraint (6) ensures that each slab can be scheduled only once. Constraint (7) restricts the cumulative length of continuously rolled slabs with the same width in each rolling unit. Constraint (8) indicates rolling mill production capacity, which restricts the lower and upper bounds of cumulative length of slabs in each rolling unit. Constraint (9) means that the total idle time allocated for rolling units cannot be greater than margin of production capability. Constraints (10)–(13) restrict the value of rijk, xijk, and yik according to their logical relationship.


## 3. Production Scheduling Optimization Method

As known that VRP is a classical NP-hard problem, it is hard to find the optimal solution for large scale problem. Since there are a large number of slabs in the day-ahead scheduling problem combined with complex objective functions, such asf2 with quadratic equation (15), it is difficult to find the exact optimal solution, even a feasible solution. In this paper, the production scheduling method consists of two stages. In the first stage, objectives shown in (2)-(3) are optimized simultaneously, and a set of Pareto optimal solutions is generated by the multiobjective optimization algorithm. In the second stage, TOPSIS based multicriteria decision-making is performed to recommend an optimal solution to facilitate field operation.


### 3.1. NSGA-II Based Multiobjective Optimization

Recently, many swarm intelligence algorithms are introduced to solve complex optimization problem, in which Nondominated Sorting Genetic Algorithm with Elitism (NSGA-II) proposed by Deb [25] is a typical method to solve multiobjective problem. NSGA-II has been widely used to solve combinatorial optimization problems in engineering, such as hydrothermal power scheduling problem [26], job sequencing problem [27], and flow-shop scheduling problem [28]. In this paper, a NSGA-II based Multiobjective Production Scheduling Algorithm (MOPSA) is developed to solve the HRPSP; some personalized changes are made to instantiate the algorithm, in which the most important things are designing customized chromosome code and genetic operators to adapt specific problem.


#### 3.1.1. Chromosome Encoding

In order to contain information of both batch scheduling and job-shop scheduling, a hybrid chromosome code consisting of two sections as shown in Figure4 is designed. The first section is a natural number sequence C that can be transformed to a two-dimensional matrix B through a code mapping procedure, where B represents a batch scheduling solution and element bij in B is the original sequence of slab j in rolling unit i. For each i, if the minimal j is found while bij=0, it can be resolved that the last slab in rolling unit i is bi,j-1. The second section is a floating number sequence V that represents the idle time allocated during job-shop scheduling, where job means production of rolling units.Figure 4
Relationship between production scheduling and electricity costs under TOU pricing.According to above description, the hybrid chromosome codeG can be expressed as(16)G=C,V,C=c1,c2,…,cm×n,V=v1,v2,…,vm,where element ci in C is a natural number that ranged from 1 to m×n, m is the quantity of rolling units, n is the quantity of slabs to be scheduled, any two numbers ci and cj are assigned to different values, and vi in V represents the idle time allocated to rolling unit i before rolling production.Detailed steps of the code mapping procedure as mentioned previously are listed as follows.Step 1.
Setfi(i=1,2,…,n) to 0, where fi is a flag and fi=1 represents the fact that slab i has been scheduled into rolling units; for rolling unit k(k=1,2,…,m), set numk=0, where numk means the slab quantity in rolling unit k; set dk=0, where dk is the accumulative rolling length in rolling unit k; set qk=0, where qk means the continuously rolled length of slabs with same width in rolling unit k; define a loop variable j and set j=1.Step 2.
Confirm the variabless and k in accordance with natural number cj, by which slab s scheduled in rolling unit k can be determined. s and k can be calculated by(17)s=cj-cj-1m×m,k=cj-1n+1.Step 3.
Check if conditionfs=0 is satisfied:(i)
If it is satisfied, it means that slabs is an unscheduled slab. Then if ws≠wk′, set qk=0, where ws is the width of the slab s and wk′ is the width of the latest appended slab in rolling unit k. Furthermore, if dk+ls≤U and qk+ls≤R, put slab s into rolling unit k and update matrix B(=[bij]) by bk,numk=s, set numk=numk+1, dk=dk+ls, qk=qk+ls, and fs=1.(ii)
Otherwise, go to Step  4.Step 4.
Updatej=j+1; go to Step  2 to repeat the above operations until j=m×n+1.Step 5.
Check iffi=1(i=1,2,…,n) and dk≥L(k=1,2,…,m) are all satisfied:(i)
If they are satisfied, it means that all slabs are scheduled into rolling units with subjection to given constraints. Perform idle time allocation procedure to generate sequenceV of chromosome code; then G=(C,V) represent a feasible solution of the model in this paper.(ii)
Otherwise, a large number should be assigned to functionsf1 and f2 to avoid chromosome being selected into new population in next selection operator.The detailed steps of idle time allocation as mentioned above based on sectionC of chromosome code are listed as follows.Step 1.
Confirm the electricity priceπis corresponding to each rolling unit i when the production starts and price πie when the production completes.Step 2.
Create a random floating number sequenceV that represents the idle time allocated for rolling units before production. For elements in V, the constraint as (9) in Section 2.2 must be satisfied.Step 3.
Sort time periods in descending order based on electricity price; after that, a new set of time periodsT′=t1′,t2′,…,tt′ is generated, in which the price associated with tk′ is πk′; define a loop variable j and set j=1.Step 4.
Adjust the idle time allocation for rolling units. For each rolling uniti that started from time period tj′, if πie<πj′ and vi+1>0, set vi+1=0; vi=vi+vi+1; for rolling unit i that completed in time period tj′, if πis<πj′ and vi>0, set vi=0; vi+1=vi+1+vi.Step 5.
Update variablej=j+1; go to Step  4 to repeat the above operation for the left time periods until j=t, which represent the fact that adjustment of idle time allocation is completed.The benefits of hybrid encoding and mapping procedure above are not only containing complete information of production scheduling but also handling constraints. From Step5, we can see that all constraints from (4)–(8) in Section 2.2 are satisfied in accepted feasible solution, which is helpful to reduce the difficulty of problem solving.


#### 3.1.2. Design of Genetic Operators

In order to instantiate the MOPSA algorithm, customized genetic operators are defined to match hybrid chromosome code; the most important operators for genetic algorithm are selection, crossover, and mutation.Selection operator, which means selecting individuals from population, is done based on the frontier rank of individuals by nondominated sorting. If many individuals have the same rank, the individual with the maximum crowded distance will be selected preferentially.Partially Mapped Crossover (PMX) mentioned in [29] and Scramble Sublist Mutation (SSM) mentioned in [30] are adopted to perform operation on section C of chromosome code. The PMX operator is performed on two parent chromosomes: randomly select two crossover points k1 and k2 and separate the chromosome code into three sections; swap the gene codes in range k1,k2; after that, replace the other gene codes out of range [k1,k2] according to mapping relationship determined by the middle section.Unlike the crossover operator, SSM mutation operator is performed on single parent chromosome: randomly select two positionsp1 and p2 that separated less than a fixed length in the chromosome code; then rearrange the gene codes within [p1,p2].After crossover or mutation operation, update sectionV of the chromosome code to allocate idle time for rolling units immediately.


#### 3.1.3. Decision Expressions and Fitness Function Calculation

We choose the objective functionsf1 and f2 to be the fitness functions in our genetic algorithm. f1 represents penalties and f2 represents electricity costs in production, which are both cost-oriented and need to find minimum value.In fitness function calculation, most needed variables and expressions are static and can be precomputed except the variablesxijk in f1 and dij in f2, so the key of fitness function calculation is to determine the values of xijk and dij based on chromosome code.According to the characteristics of chromosome code in this paper, we use matrixB=bij generated in chromosome code mapping procedure instead of part C to perform the following calculation. In order to determine the value of xijk, each row in matrix B should be traversed to search the adjacent elements that satisfy the following equation:(18)bk,j1=i,bk,j1+1=j,where the first equation means slab i is assigned in rolling units k and processed with the sequence j1 and the next equation indicates that slab j is allocated after slab i immediately in rolling unit k. xijk can be determined to be 1 if (18) is satisfied; otherwise, it is 0. For each rolling unit k, penalties between adjacent slabs are accumulated by Pij·xijk.Meanwhile, it should be noted that calculation ofdij in f2 would depend not only on B but also on sequence V that represents the allocated idle time for rolling units. According to (14)-(15) defined in Section 2.2, the determination of dij mainly depends on variables yik, rijk, and vi, in which the first two variables can be easily calculated on matrix B by a traversal procedure as done in determining xijk, and the last variable vi can be directly identified by the sequence V in chromosome code. Once dij is known, fitness function f2 can be accumulated by πj·Wi·dij for each time period.


### 3.2. TOPSIS Based Multicriteria Decision-Making

As MOPSA generate more than one Pareto optimal solution, in order to facilitate field operation, only a few solutions should be accepted. In this paper, Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) [31], a widely used multicriteria decision-making method to identify solutions from finite alternatives, is adopted as the method to select a recommended optimal solution.Detailed steps of the TOPSIS based multicriteria decision-making for HRPSP are listed as follows.Step 1.
The decision matrixX can be expressed as(19)X=x11x12x21x22⋮⋮xm1xm2,where X is a two-dimensional matrix with the size of m×n, which means that there are m solutions generated by the multiobjective algorithm and n objectives for the HRPSP, where n=2. The element xij in X is the value of the jth objective with respect to the ith solution. Then the normalized decision matrix Z=zij can be calculated according to(20)zij=xij∑i=1mxij2.Step 2.
Multiply the normalized decision matrix by its associated weights to calculate the weighted normalized decision matrixV=vij, in which vij is calculated as(21)vij=wj·zij,where wj is a weight factor associated with the jth objective. In our context, w1 and w2 are set to different values according to preference of two objectives.Step 3.
Identify the ideal solutions+ and the nadir solution s- of each objective according to the following equations:(22)s+=s1+,s2+,sj+=max1≤i≤mvijiffjis benefit-oriented,min1≤i≤mvijiffjis cost-oriented,s-=s1-,s2-,sj-=min1≤i≤mvijiffjis benefit-oriented,max1≤i≤mvijiffjis cost-oriented.
It should be known that both of the objectives in HRPSP are cost-oriented, which are said to find the minimum of objective functions.Step 4.
Measure the distancesdi+ and di- of the ith solution from the ideal solution s+ and the nadir solution s- by(23)di+=∑j=1nvij-sj+2,i=1,2,…,m,di-=∑j=1nvij-sj-2,i=1,2,…,m.Step 5.
CalculateCi∗ that represents the relative closeness of ith solution with respect to the ideal solution according to(24)Ci∗=di-di-+di+,i=1,2,…,m.
After completing the above steps, the decision-making can be finally performed on the Pareto optimal solutions according to the sequence determined byCi∗i=1,2,…,m in descending order; the solution that owns the maximal relative closeness will be selected as the recommended optimal solution.


## 3.1. NSGA-II Based Multiobjective Optimization

Recently, many swarm intelligence algorithms are introduced to solve complex optimization problem, in which Nondominated Sorting Genetic Algorithm with Elitism (NSGA-II) proposed by Deb [25] is a typical method to solve multiobjective problem. NSGA-II has been widely used to solve combinatorial optimization problems in engineering, such as hydrothermal power scheduling problem [26], job sequencing problem [27], and flow-shop scheduling problem [28]. In this paper, a NSGA-II based Multiobjective Production Scheduling Algorithm (MOPSA) is developed to solve the HRPSP; some personalized changes are made to instantiate the algorithm, in which the most important things are designing customized chromosome code and genetic operators to adapt specific problem.


### 3.1.1. Chromosome Encoding

In order to contain information of both batch scheduling and job-shop scheduling, a hybrid chromosome code consisting of two sections as shown in Figure4 is designed. The first section is a natural number sequence C that can be transformed to a two-dimensional matrix B through a code mapping procedure, where B represents a batch scheduling solution and element bij in B is the original sequence of slab j in rolling unit i. For each i, if the minimal j is found while bij=0, it can be resolved that the last slab in rolling unit i is bi,j-1. The second section is a floating number sequence V that represents the idle time allocated during job-shop scheduling, where job means production of rolling units.Figure 4
Relationship between production scheduling and electricity costs under TOU pricing.According to above description, the hybrid chromosome codeG can be expressed as(16)G=C,V,C=c1,c2,…,cm×n,V=v1,v2,…,vm,where element ci in C is a natural number that ranged from 1 to m×n, m is the quantity of rolling units, n is the quantity of slabs to be scheduled, any two numbers ci and cj are assigned to different values, and vi in V represents the idle time allocated to rolling unit i before rolling production.Detailed steps of the code mapping procedure as mentioned previously are listed as follows.Step 1.
Setfi(i=1,2,…,n) to 0, where fi is a flag and fi=1 represents the fact that slab i has been scheduled into rolling units; for rolling unit k(k=1,2,…,m), set numk=0, where numk means the slab quantity in rolling unit k; set dk=0, where dk is the accumulative rolling length in rolling unit k; set qk=0, where qk means the continuously rolled length of slabs with same width in rolling unit k; define a loop variable j and set j=1.Step 2.
Confirm the variabless and k in accordance with natural number cj, by which slab s scheduled in rolling unit k can be determined. s and k can be calculated by(17)s=cj-cj-1m×m,k=cj-1n+1.Step 3.
Check if conditionfs=0 is satisfied:(i)
If it is satisfied, it means that slabs is an unscheduled slab. Then if ws≠wk′, set qk=0, where ws is the width of the slab s and wk′ is the width of the latest appended slab in rolling unit k. Furthermore, if dk+ls≤U and qk+ls≤R, put slab s into rolling unit k and update matrix B(=[bij]) by bk,numk=s, set numk=numk+1, dk=dk+ls, qk=qk+ls, and fs=1.(ii)
Otherwise, go to Step  4.Step 4.
Updatej=j+1; go to Step  2 to repeat the above operations until j=m×n+1.Step 5.
Check iffi=1(i=1,2,…,n) and dk≥L(k=1,2,…,m) are all satisfied:(i)
If they are satisfied, it means that all slabs are scheduled into rolling units with subjection to given constraints. Perform idle time allocation procedure to generate sequenceV of chromosome code; then G=(C,V) represent a feasible solution of the model in this paper.(ii)
Otherwise, a large number should be assigned to functionsf1 and f2 to avoid chromosome being selected into new population in next selection operator.The detailed steps of idle time allocation as mentioned above based on sectionC of chromosome code are listed as follows.Step 1.
Confirm the electricity priceπis corresponding to each rolling unit i when the production starts and price πie when the production completes.Step 2.
Create a random floating number sequenceV that represents the idle time allocated for rolling units before production. For elements in V, the constraint as (9) in Section 2.2 must be satisfied.Step 3.
Sort time periods in descending order based on electricity price; after that, a new set of time periodsT′=t1′,t2′,…,tt′ is generated, in which the price associated with tk′ is πk′; define a loop variable j and set j=1.Step 4.
Adjust the idle time allocation for rolling units. For each rolling uniti that started from time period tj′, if πie<πj′ and vi+1>0, set vi+1=0; vi=vi+vi+1; for rolling unit i that completed in time period tj′, if πis<πj′ and vi>0, set vi=0; vi+1=vi+1+vi.Step 5.
Update variablej=j+1; go to Step  4 to repeat the above operation for the left time periods until j=t, which represent the fact that adjustment of idle time allocation is completed.The benefits of hybrid encoding and mapping procedure above are not only containing complete information of production scheduling but also handling constraints. From Step5, we can see that all constraints from (4)–(8) in Section 2.2 are satisfied in accepted feasible solution, which is helpful to reduce the difficulty of problem solving.


### 3.1.2. Design of Genetic Operators

In order to instantiate the MOPSA algorithm, customized genetic operators are defined to match hybrid chromosome code; the most important operators for genetic algorithm are selection, crossover, and mutation.Selection operator, which means selecting individuals from population, is done based on the frontier rank of individuals by nondominated sorting. If many individuals have the same rank, the individual with the maximum crowded distance will be selected preferentially.Partially Mapped Crossover (PMX) mentioned in [29] and Scramble Sublist Mutation (SSM) mentioned in [30] are adopted to perform operation on section C of chromosome code. The PMX operator is performed on two parent chromosomes: randomly select two crossover points k1 and k2 and separate the chromosome code into three sections; swap the gene codes in range k1,k2; after that, replace the other gene codes out of range [k1,k2] according to mapping relationship determined by the middle section.Unlike the crossover operator, SSM mutation operator is performed on single parent chromosome: randomly select two positionsp1 and p2 that separated less than a fixed length in the chromosome code; then rearrange the gene codes within [p1,p2].After crossover or mutation operation, update sectionV of the chromosome code to allocate idle time for rolling units immediately.


### 3.1.3. Decision Expressions and Fitness Function Calculation

We choose the objective functionsf1 and f2 to be the fitness functions in our genetic algorithm. f1 represents penalties and f2 represents electricity costs in production, which are both cost-oriented and need to find minimum value.In fitness function calculation, most needed variables and expressions are static and can be precomputed except the variablesxijk in f1 and dij in f2, so the key of fitness function calculation is to determine the values of xijk and dij based on chromosome code.According to the characteristics of chromosome code in this paper, we use matrixB=bij generated in chromosome code mapping procedure instead of part C to perform the following calculation. In order to determine the value of xijk, each row in matrix B should be traversed to search the adjacent elements that satisfy the following equation:(18)bk,j1=i,bk,j1+1=j,where the first equation means slab i is assigned in rolling units k and processed with the sequence j1 and the next equation indicates that slab j is allocated after slab i immediately in rolling unit k. xijk can be determined to be 1 if (18) is satisfied; otherwise, it is 0. For each rolling unit k, penalties between adjacent slabs are accumulated by Pij·xijk.Meanwhile, it should be noted that calculation ofdij in f2 would depend not only on B but also on sequence V that represents the allocated idle time for rolling units. According to (14)-(15) defined in Section 2.2, the determination of dij mainly depends on variables yik, rijk, and vi, in which the first two variables can be easily calculated on matrix B by a traversal procedure as done in determining xijk, and the last variable vi can be directly identified by the sequence V in chromosome code. Once dij is known, fitness function f2 can be accumulated by πj·Wi·dij for each time period.


## 3.1.1. Chromosome Encoding

In order to contain information of both batch scheduling and job-shop scheduling, a hybrid chromosome code consisting of two sections as shown in Figure4 is designed. The first section is a natural number sequence C that can be transformed to a two-dimensional matrix B through a code mapping procedure, where B represents a batch scheduling solution and element bij in B is the original sequence of slab j in rolling unit i. For each i, if the minimal j is found while bij=0, it can be resolved that the last slab in rolling unit i is bi,j-1. The second section is a floating number sequence V that represents the idle time allocated during job-shop scheduling, where job means production of rolling units.Figure 4
Relationship between production scheduling and electricity costs under TOU pricing.According to above description, the hybrid chromosome codeG can be expressed as(16)G=C,V,C=c1,c2,…,cm×n,V=v1,v2,…,vm,where element ci in C is a natural number that ranged from 1 to m×n, m is the quantity of rolling units, n is the quantity of slabs to be scheduled, any two numbers ci and cj are assigned to different values, and vi in V represents the idle time allocated to rolling unit i before rolling production.Detailed steps of the code mapping procedure as mentioned previously are listed as follows.Step 1.
Setfi(i=1,2,…,n) to 0, where fi is a flag and fi=1 represents the fact that slab i has been scheduled into rolling units; for rolling unit k(k=1,2,…,m), set numk=0, where numk means the slab quantity in rolling unit k; set dk=0, where dk is the accumulative rolling length in rolling unit k; set qk=0, where qk means the continuously rolled length of slabs with same width in rolling unit k; define a loop variable j and set j=1.Step 2.
Confirm the variabless and k in accordance with natural number cj, by which slab s scheduled in rolling unit k can be determined. s and k can be calculated by(17)s=cj-cj-1m×m,k=cj-1n+1.Step 3.
Check if conditionfs=0 is satisfied:(i)
If it is satisfied, it means that slabs is an unscheduled slab. Then if ws≠wk′, set qk=0, where ws is the width of the slab s and wk′ is the width of the latest appended slab in rolling unit k. Furthermore, if dk+ls≤U and qk+ls≤R, put slab s into rolling unit k and update matrix B(=[bij]) by bk,numk=s, set numk=numk+1, dk=dk+ls, qk=qk+ls, and fs=1.(ii)
Otherwise, go to Step  4.Step 4.
Updatej=j+1; go to Step  2 to repeat the above operations until j=m×n+1.Step 5.
Check iffi=1(i=1,2,…,n) and dk≥L(k=1,2,…,m) are all satisfied:(i)
If they are satisfied, it means that all slabs are scheduled into rolling units with subjection to given constraints. Perform idle time allocation procedure to generate sequenceV of chromosome code; then G=(C,V) represent a feasible solution of the model in this paper.(ii)
Otherwise, a large number should be assigned to functionsf1 and f2 to avoid chromosome being selected into new population in next selection operator.The detailed steps of idle time allocation as mentioned above based on sectionC of chromosome code are listed as follows.Step 1.
Confirm the electricity priceπis corresponding to each rolling unit i when the production starts and price πie when the production completes.Step 2.
Create a random floating number sequenceV that represents the idle time allocated for rolling units before production. For elements in V, the constraint as (9) in Section 2.2 must be satisfied.Step 3.
Sort time periods in descending order based on electricity price; after that, a new set of time periodsT′=t1′,t2′,…,tt′ is generated, in which the price associated with tk′ is πk′; define a loop variable j and set j=1.Step 4.
Adjust the idle time allocation for rolling units. For each rolling uniti that started from time period tj′, if πie<πj′ and vi+1>0, set vi+1=0; vi=vi+vi+1; for rolling unit i that completed in time period tj′, if πis<πj′ and vi>0, set vi=0; vi+1=vi+1+vi.Step 5.
Update variablej=j+1; go to Step  4 to repeat the above operation for the left time periods until j=t, which represent the fact that adjustment of idle time allocation is completed.The benefits of hybrid encoding and mapping procedure above are not only containing complete information of production scheduling but also handling constraints. From Step5, we can see that all constraints from (4)–(8) in Section 2.2 are satisfied in accepted feasible solution, which is helpful to reduce the difficulty of problem solving.


## 3.1.2. Design of Genetic Operators

In order to instantiate the MOPSA algorithm, customized genetic operators are defined to match hybrid chromosome code; the most important operators for genetic algorithm are selection, crossover, and mutation.Selection operator, which means selecting individuals from population, is done based on the frontier rank of individuals by nondominated sorting. If many individuals have the same rank, the individual with the maximum crowded distance will be selected preferentially.Partially Mapped Crossover (PMX) mentioned in [29] and Scramble Sublist Mutation (SSM) mentioned in [30] are adopted to perform operation on section C of chromosome code. The PMX operator is performed on two parent chromosomes: randomly select two crossover points k1 and k2 and separate the chromosome code into three sections; swap the gene codes in range k1,k2; after that, replace the other gene codes out of range [k1,k2] according to mapping relationship determined by the middle section.Unlike the crossover operator, SSM mutation operator is performed on single parent chromosome: randomly select two positionsp1 and p2 that separated less than a fixed length in the chromosome code; then rearrange the gene codes within [p1,p2].After crossover or mutation operation, update sectionV of the chromosome code to allocate idle time for rolling units immediately.


## 3.1.3. Decision Expressions and Fitness Function Calculation

We choose the objective functionsf1 and f2 to be the fitness functions in our genetic algorithm. f1 represents penalties and f2 represents electricity costs in production, which are both cost-oriented and need to find minimum value.In fitness function calculation, most needed variables and expressions are static and can be precomputed except the variablesxijk in f1 and dij in f2, so the key of fitness function calculation is to determine the values of xijk and dij based on chromosome code.According to the characteristics of chromosome code in this paper, we use matrixB=bij generated in chromosome code mapping procedure instead of part C to perform the following calculation. In order to determine the value of xijk, each row in matrix B should be traversed to search the adjacent elements that satisfy the following equation:(18)bk,j1=i,bk,j1+1=j,where the first equation means slab i is assigned in rolling units k and processed with the sequence j1 and the next equation indicates that slab j is allocated after slab i immediately in rolling unit k. xijk can be determined to be 1 if (18) is satisfied; otherwise, it is 0. For each rolling unit k, penalties between adjacent slabs are accumulated by Pij·xijk.Meanwhile, it should be noted that calculation ofdij in f2 would depend not only on B but also on sequence V that represents the allocated idle time for rolling units. According to (14)-(15) defined in Section 2.2, the determination of dij mainly depends on variables yik, rijk, and vi, in which the first two variables can be easily calculated on matrix B by a traversal procedure as done in determining xijk, and the last variable vi can be directly identified by the sequence V in chromosome code. Once dij is known, fitness function f2 can be accumulated by πj·Wi·dij for each time period.


## 3.2. TOPSIS Based Multicriteria Decision-Making

As MOPSA generate more than one Pareto optimal solution, in order to facilitate field operation, only a few solutions should be accepted. In this paper, Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) [31], a widely used multicriteria decision-making method to identify solutions from finite alternatives, is adopted as the method to select a recommended optimal solution.Detailed steps of the TOPSIS based multicriteria decision-making for HRPSP are listed as follows.Step 1.
The decision matrixX can be expressed as(19)X=x11x12x21x22⋮⋮xm1xm2,where X is a two-dimensional matrix with the size of m×n, which means that there are m solutions generated by the multiobjective algorithm and n objectives for the HRPSP, where n=2. The element xij in X is the value of the jth objective with respect to the ith solution. Then the normalized decision matrix Z=zij can be calculated according to(20)zij=xij∑i=1mxij2.Step 2.
Multiply the normalized decision matrix by its associated weights to calculate the weighted normalized decision matrixV=vij, in which vij is calculated as(21)vij=wj·zij,where wj is a weight factor associated with the jth objective. In our context, w1 and w2 are set to different values according to preference of two objectives.Step 3.
Identify the ideal solutions+ and the nadir solution s- of each objective according to the following equations:(22)s+=s1+,s2+,sj+=max1≤i≤mvijiffjis benefit-oriented,min1≤i≤mvijiffjis cost-oriented,s-=s1-,s2-,sj-=min1≤i≤mvijiffjis benefit-oriented,max1≤i≤mvijiffjis cost-oriented.
It should be known that both of the objectives in HRPSP are cost-oriented, which are said to find the minimum of objective functions.Step 4.
Measure the distancesdi+ and di- of the ith solution from the ideal solution s+ and the nadir solution s- by(23)di+=∑j=1nvij-sj+2,i=1,2,…,m,di-=∑j=1nvij-sj-2,i=1,2,…,m.Step 5.
CalculateCi∗ that represents the relative closeness of ith solution with respect to the ideal solution according to(24)Ci∗=di-di-+di+,i=1,2,…,m.
After completing the above steps, the decision-making can be finally performed on the Pareto optimal solutions according to the sequence determined byCi∗i=1,2,…,m in descending order; the solution that owns the maximal relative closeness will be selected as the recommended optimal solution.


## 4. Experimental Results and Analyses

In this section, we perform a series of experiments to evaluate the effectiveness and performance of the proposed method in different scenario.


### 4.1. Experimental Procedure

In experimental procedure, four groups of production data as is shown in Table1 are collected from a steel mill for experimental procedure. For each group of production data, if there are many slab varieties in width, gauge, and hardness, the penalty score between adjacent slabs will be larger. At the same time, full production load means that the idle time for processing slabs will be short.Table 1
Production data description.
Group Id
Slab quantity
Rolling units quantity
Processing time/(min)
Characteristics
1
450
8
1421.75
Many varieties of slabs and full production load
2
415
8
1323.05
Many varieties of slabs and not full production load
3
450
8
1427.88
Few varieties of slabs and full production load
4
415
8
1318.33
Few varieties of slabs and not full production loadAccording to constraints of production equipment and capability, the lower and upper bounds of the cumulative length of slabs scheduled in a single rolling unit are, respectively, set to 5 and 10 kilometers, and the upper bound of the continuously rolled length of slabs with the same width is set to 1 kilometer. For specific slab, rolling length, processing time, and power consumption can be obtained by the hot rolling process control system in steel mill. The penalties caused by jumps between adjacent slabs in width, gauge, and hardness are adopted by referring to [8]. The data in Table 2 are actually performed TOU electricity tariffs in steel mill. According to daily power load distribution, a whole day is split into eight periods that contain four types of time periods; each type of time period is associated with corresponding price.Table 2
TOU electricity tariffs.
Time period
Time frame
Electricity price/(CNY⋅kWh−1)
On-peak
18:00–21:00
0.878
Mid-peak
08:00–11:00, 15:00–18:00
0.778
Flat-peak
07:00-08:00, 11:00–15:00, 21:00-22:00
0.628
Off-peak
00:00–07:00, 22:00–24:00
0.428In order to obtain excellent algorithm performance, the NSGA-II parameters are determined by parameter sensitivity analysis based on empirical value and a lot of tests. The probability of crossover and mutation are set to 0.4 and 0.6, respectively, the population size is set to 50, and the maximum iterations of algorithm are set to 5000. The production scheduling optimization algorithm and TOPSIS decision-making procedure are both implemented and performed in MATLAB.In experimental procedure, the proposed method (named as PM) is compared with two conventional methods to evaluate effectiveness and performance. Since exact algorithm for large scale HRBSP problem is too difficult to implement, genetic algorithm is often used for solving this problem. In this paper, a relatively new method in [32] with the traditional objective to minimize jump penalties is adopted as a comparison method (named as CM1), in which a hybrid evolutionary algorithm with integration of genetic algorithm and extremal optimization is designed to solve the hot rolling scheduling problem.Because electricity price during hot rolling changes over time, it is natural to allocate the processing sequence and the idle time of rolling units to avoid on-peak time periods; then the MILP method proposed by [6] is adopted as another comparison method (named as CM2) to find the low bound of electricity costs on the basis of solution obtained in CM1.Unlike single objective optimization, the result of multiobjective optimization is not a single solution but a set of Pareto optimal solutions; in order to facilitate field operation, we choose different values of objective weight factorswj in TOPSIS decision-making procedure to recommend solution with different preference of penalty score and electricity cost. In our experimental procedure, the objective factors w=w1,w2 of the proposed method are set to 0.9,0.1, 0.4,0.6, and [0.1,0.9], respectively.Optimization results obtained by different methods are provided in Table3, in which we can see that penalties obtained by PM with w=[0.1,0.9] are roughly the same as those obtained by CM1 and CM2 but electricity costs are cut down obviously. It is obvious that load shifting to reduce electricity cost inevitably results in an increasing of penalty score, and we would just like to point out that minimizing jump penalties is a guiding target but not a strictly rigid constraint in engineering. If there is allowance for penalty increase on electricity cost, more significant effect on electricity cost reduction is shown, which tells us that penalty relaxation can play an import role, while electricity costs are the key consideration in production; as a consequence, we can utilize objective weight factors in TOPSIS procedure to adjust preferences of the two objectives. In our cases, electricity cost obtained by PM with TOPSIS decision-making on each group of data is less than CM1; even compared to CM2, which includes load shifting on fixed rolling batches, the result is still better; this advantage is attributed to TOU pricing based batching to construct rolling units. Besides that, we can see that the optimization effect is more significant, while the production load is not full, that is, groups 2 and 4, which is caused by more idle time margin which existed to avoid on-peak time periods in such situation.Table 3
Scheduling results obtained by different methods.
Method
Group 1
Group 2
Group 3
Group 4
f
1
f
2
f
1
f
2
f
1
f
2
f
1
f
2
CM1
5035
313254
4528
296357
2957
315894
2659
299623
CM2
5035
309078
4528
276813
2957
312753
2659
278717
PM,w=[0.9,0.1]
5129
308281
4573
275898
3090
311114
2710
277214
PM,w=0.4,0.6
7493
305691
6905
274022
3445
309242
3022
275397
PM,w=0.1,0.9
7701
305680
7665
273729
3478
309234
3308
274994


### 4.2. Scheduling Results Analysis

In this section, group 1 of data is chosen to have a detailed analysis on job scheduling results firstly. Because the main idea of this paper is ELD, the proposed method PM with[w1,w2]=[0.1,0.9], which has the most significant effect on electricity cost reduction, is selected to compare with the conventional method CM1. Rolling parameters obtained by both methods are given in Table 4, from which we can see that the parameters are subjected to instantiated constraints, which represent the fact that the schedule is feasible solution. Then, we analyze the scheduling results from two aspects.Table 4
Detailed parameters of scheduling results for group 1 of production data.
Method
RUS
SQ
RL
PT
PD
APL
PST
PCT
AIT
PM,w=0.1,0.9
1
51
8.53
2.44
56.74
23.25
00:00
02:26
0
2
58
9.95
2.79
64.17
23.00
02:26
05:14
0
3
55
9.92
2.80
60.75
21.70
05:14
08:02
0
4
61
9.94
3.44
68.20
19.83
08:02
11:28
0
5
59
9.96
3.16
66.46
21.03
11:28
14:38
0
6
57
9.81
3.17
62.24
19.63
14:38
17:48
0
7
52
9.10
2.89
54.39
18.82
17:48
20:41
0
8
57
9.94
3.00
66.69
22.23
20:59
23:59
0.3
CM1
1
57
9.80
3.07
63.38
20.64
00:00
03:04
0
2
62
10.00
3.23
67.76
20.98
03:04
06:18
0
3
56
9.92
2.96
63.13
21.33
06:18
09:16
0
4
55
9.47
2.96
62.15
21.00
09:16
12:14
0
5
58
9.91
2.99
63.72
21.31
12:14
15:13
0
6
53
9.53
2.83
59.68
21.09
15:13
18:03
0
7
51
8.54
2.70
56.86
21.06
18:03
20:45
0
8
58
9.98
2.96
62.95
21.27
20:45
23:42
0
RUS: rolling unit sequence; SQ: slab quantity; RL: rolling length (km); PT: processing time (h); PD: power demand (MW⋅h); APL: average power load (MW); PST: processing start time (HH:mm); PCT: processing complete time (HH:mm); AIT: allocated idle time (hour).On one hand, rolling units in Table4 are considered as production jobs and are illustrated in Figure 5. As it can be seen, in any subfigure, heavy loads are allocated in off-peak and flat-peak periods by PM, while light loads are allocated in on-peak or mid-peak periods. In addition, idle time is allocated at 18:00 to 21:00 for our scenarios. Another phenomenon is that the power load difference between heavy load and light load in PM is greater than that in CM1 and CM2, which is due to the fact that rolling units in PM are organized by TOU electricity price and their processing time.Figure 5
Illustration of job scheduling results obtained by PM and CM1.
(a)
Group 1
(b)
Group 2
(c)
Group 3
(d)
Group 4On the other hand, average power load distribution among time periods is illustrated in Figure6. Compared to CM1, power load obtained by PM reduces greatly in the last on-peak periods and increases substantially in last off-peak period, especially for group 2 and group 4, which are characterized by not full production load. At the same time, power load in the first off-peak period increases in a certain extent. In addition, power load distribution obtained by PM is also better than that obtained by CM2 based on the principle of load shifting corresponding to TOU pricing, which confirms the effectiveness and advancement of the proposed method furthermore.Figure 6
Power load distribution among time periods.
(a)
Group 1
(b)
Group 2
(c)
Group 3
(d)
Group 4From above results and analyses, we know that the advantages of our proposed method on electricity cost reduction can be attributed to two aspects: one is load shifting to avoid on-peak time periods and the other one is TOU pricing based load planning.


### 4.3. Robustness of the Algorithm

It is well known that NSGA-II is a randomized algorithm; each run of the algorithm may get different results. For evaluating robustness of the algorithm, we use box plot to portray the convergence metric in repeated operation, which is represented by average value of minimum normalized Euclidean distance and indicates the disparity between approximate Pareto frontier and ideal Pareto frontier.Assume thatP∗=p1,p2,…,pP∗ are the optimal solutions evenly distributed on ideal Pareto frontier and A∗=a1,a2,…,aA∗ are the approximate solutions obtained in a single run of the proposed algorithm. For any ai, minimum normalized Euclidean distance di between ai and P∗ can be calculated by(25)di=minj=1P∗∑m=12fmai-fmajfmmax-fmmin2,where fmmax and fmmin, respectively, represent the maximum and minimum values of the mth objective function in P∗, and then the convergence metric C can be expressed as(26)CA∗≜∑i=1A∗diA∗.Note that the ideal Pareto frontier is always unknown in real problem; the algorithm proposed in this application is run 30 times, respectively, on each group of production data, and then a pseudo-Pareto frontier, which consists of all the solutions in 30 runs with removing dominated solutions, is constructed to compare with the approximate Pareto frontiers. For every run, box plots based on convergence metrics are illustrated in Figure7. In general, metric C in less than 10-2 means good statistical convergence performance in Pareto optimality based multiobjective optimization. The symbol “+” in Figure 7 refers to an outlier in box statistics; nevertheless, it can be seen that the outlier is very close to 10-2. Overall, we can see that the upper edges on different groups of data are all less than 10-2, except a slightly larger value on group 4 and an outlier on group 2. However, the 3rd quartile on group 4 is totally in the range of less than 10-2. The statistical results show that the proposed algorithm is stable in repeated run. On the whole, we can conclude that the proposed algorithm is robust and suitable for application in engineering.Figure 7
Box plots based on convergence metrics.


## 4.1. Experimental Procedure

In experimental procedure, four groups of production data as is shown in Table1 are collected from a steel mill for experimental procedure. For each group of production data, if there are many slab varieties in width, gauge, and hardness, the penalty score between adjacent slabs will be larger. At the same time, full production load means that the idle time for processing slabs will be short.Table 1
Production data description.
Group Id
Slab quantity
Rolling units quantity
Processing time/(min)
Characteristics
1
450
8
1421.75
Many varieties of slabs and full production load
2
415
8
1323.05
Many varieties of slabs and not full production load
3
450
8
1427.88
Few varieties of slabs and full production load
4
415
8
1318.33
Few varieties of slabs and not full production loadAccording to constraints of production equipment and capability, the lower and upper bounds of the cumulative length of slabs scheduled in a single rolling unit are, respectively, set to 5 and 10 kilometers, and the upper bound of the continuously rolled length of slabs with the same width is set to 1 kilometer. For specific slab, rolling length, processing time, and power consumption can be obtained by the hot rolling process control system in steel mill. The penalties caused by jumps between adjacent slabs in width, gauge, and hardness are adopted by referring to [8]. The data in Table 2 are actually performed TOU electricity tariffs in steel mill. According to daily power load distribution, a whole day is split into eight periods that contain four types of time periods; each type of time period is associated with corresponding price.Table 2
TOU electricity tariffs.
Time period
Time frame
Electricity price/(CNY⋅kWh−1)
On-peak
18:00–21:00
0.878
Mid-peak
08:00–11:00, 15:00–18:00
0.778
Flat-peak
07:00-08:00, 11:00–15:00, 21:00-22:00
0.628
Off-peak
00:00–07:00, 22:00–24:00
0.428In order to obtain excellent algorithm performance, the NSGA-II parameters are determined by parameter sensitivity analysis based on empirical value and a lot of tests. The probability of crossover and mutation are set to 0.4 and 0.6, respectively, the population size is set to 50, and the maximum iterations of algorithm are set to 5000. The production scheduling optimization algorithm and TOPSIS decision-making procedure are both implemented and performed in MATLAB.In experimental procedure, the proposed method (named as PM) is compared with two conventional methods to evaluate effectiveness and performance. Since exact algorithm for large scale HRBSP problem is too difficult to implement, genetic algorithm is often used for solving this problem. In this paper, a relatively new method in [32] with the traditional objective to minimize jump penalties is adopted as a comparison method (named as CM1), in which a hybrid evolutionary algorithm with integration of genetic algorithm and extremal optimization is designed to solve the hot rolling scheduling problem.Because electricity price during hot rolling changes over time, it is natural to allocate the processing sequence and the idle time of rolling units to avoid on-peak time periods; then the MILP method proposed by [6] is adopted as another comparison method (named as CM2) to find the low bound of electricity costs on the basis of solution obtained in CM1.Unlike single objective optimization, the result of multiobjective optimization is not a single solution but a set of Pareto optimal solutions; in order to facilitate field operation, we choose different values of objective weight factorswj in TOPSIS decision-making procedure to recommend solution with different preference of penalty score and electricity cost. In our experimental procedure, the objective factors w=w1,w2 of the proposed method are set to 0.9,0.1, 0.4,0.6, and [0.1,0.9], respectively.Optimization results obtained by different methods are provided in Table3, in which we can see that penalties obtained by PM with w=[0.1,0.9] are roughly the same as those obtained by CM1 and CM2 but electricity costs are cut down obviously. It is obvious that load shifting to reduce electricity cost inevitably results in an increasing of penalty score, and we would just like to point out that minimizing jump penalties is a guiding target but not a strictly rigid constraint in engineering. If there is allowance for penalty increase on electricity cost, more significant effect on electricity cost reduction is shown, which tells us that penalty relaxation can play an import role, while electricity costs are the key consideration in production; as a consequence, we can utilize objective weight factors in TOPSIS procedure to adjust preferences of the two objectives. In our cases, electricity cost obtained by PM with TOPSIS decision-making on each group of data is less than CM1; even compared to CM2, which includes load shifting on fixed rolling batches, the result is still better; this advantage is attributed to TOU pricing based batching to construct rolling units. Besides that, we can see that the optimization effect is more significant, while the production load is not full, that is, groups 2 and 4, which is caused by more idle time margin which existed to avoid on-peak time periods in such situation.Table 3
Scheduling results obtained by different methods.
Method
Group 1
Group 2
Group 3
Group 4
f
1
f
2
f
1
f
2
f
1
f
2
f
1
f
2
CM1
5035
313254
4528
296357
2957
315894
2659
299623
CM2
5035
309078
4528
276813
2957
312753
2659
278717
PM,w=[0.9,0.1]
5129
308281
4573
275898
3090
311114
2710
277214
PM,w=0.4,0.6
7493
305691
6905
274022
3445
309242
3022
275397
PM,w=0.1,0.9
7701
305680
7665
273729
3478
309234
3308
274994


## 4.2. Scheduling Results Analysis

In this section, group 1 of data is chosen to have a detailed analysis on job scheduling results firstly. Because the main idea of this paper is ELD, the proposed method PM with[w1,w2]=[0.1,0.9], which has the most significant effect on electricity cost reduction, is selected to compare with the conventional method CM1. Rolling parameters obtained by both methods are given in Table 4, from which we can see that the parameters are subjected to instantiated constraints, which represent the fact that the schedule is feasible solution. Then, we analyze the scheduling results from two aspects.Table 4
Detailed parameters of scheduling results for group 1 of production data.
Method
RUS
SQ
RL
PT
PD
APL
PST
PCT
AIT
PM,w=0.1,0.9
1
51
8.53
2.44
56.74
23.25
00:00
02:26
0
2
58
9.95
2.79
64.17
23.00
02:26
05:14
0
3
55
9.92
2.80
60.75
21.70
05:14
08:02
0
4
61
9.94
3.44
68.20
19.83
08:02
11:28
0
5
59
9.96
3.16
66.46
21.03
11:28
14:38
0
6
57
9.81
3.17
62.24
19.63
14:38
17:48
0
7
52
9.10
2.89
54.39
18.82
17:48
20:41
0
8
57
9.94
3.00
66.69
22.23
20:59
23:59
0.3
CM1
1
57
9.80
3.07
63.38
20.64
00:00
03:04
0
2
62
10.00
3.23
67.76
20.98
03:04
06:18
0
3
56
9.92
2.96
63.13
21.33
06:18
09:16
0
4
55
9.47
2.96
62.15
21.00
09:16
12:14
0
5
58
9.91
2.99
63.72
21.31
12:14
15:13
0
6
53
9.53
2.83
59.68
21.09
15:13
18:03
0
7
51
8.54
2.70
56.86
21.06
18:03
20:45
0
8
58
9.98
2.96
62.95
21.27
20:45
23:42
0
RUS: rolling unit sequence; SQ: slab quantity; RL: rolling length (km); PT: processing time (h); PD: power demand (MW⋅h); APL: average power load (MW); PST: processing start time (HH:mm); PCT: processing complete time (HH:mm); AIT: allocated idle time (hour).On one hand, rolling units in Table4 are considered as production jobs and are illustrated in Figure 5. As it can be seen, in any subfigure, heavy loads are allocated in off-peak and flat-peak periods by PM, while light loads are allocated in on-peak or mid-peak periods. In addition, idle time is allocated at 18:00 to 21:00 for our scenarios. Another phenomenon is that the power load difference between heavy load and light load in PM is greater than that in CM1 and CM2, which is due to the fact that rolling units in PM are organized by TOU electricity price and their processing time.Figure 5
Illustration of job scheduling results obtained by PM and CM1.
(a)
Group 1
(b)
Group 2
(c)
Group 3
(d)
Group 4On the other hand, average power load distribution among time periods is illustrated in Figure6. Compared to CM1, power load obtained by PM reduces greatly in the last on-peak periods and increases substantially in last off-peak period, especially for group 2 and group 4, which are characterized by not full production load. At the same time, power load in the first off-peak period increases in a certain extent. In addition, power load distribution obtained by PM is also better than that obtained by CM2 based on the principle of load shifting corresponding to TOU pricing, which confirms the effectiveness and advancement of the proposed method furthermore.Figure 6
Power load distribution among time periods.
(a)
Group 1
(b)
Group 2
(c)
Group 3
(d)
Group 4From above results and analyses, we know that the advantages of our proposed method on electricity cost reduction can be attributed to two aspects: one is load shifting to avoid on-peak time periods and the other one is TOU pricing based load planning.


## 4.3. Robustness of the Algorithm

It is well known that NSGA-II is a randomized algorithm; each run of the algorithm may get different results. For evaluating robustness of the algorithm, we use box plot to portray the convergence metric in repeated operation, which is represented by average value of minimum normalized Euclidean distance and indicates the disparity between approximate Pareto frontier and ideal Pareto frontier.Assume thatP∗=p1,p2,…,pP∗ are the optimal solutions evenly distributed on ideal Pareto frontier and A∗=a1,a2,…,aA∗ are the approximate solutions obtained in a single run of the proposed algorithm. For any ai, minimum normalized Euclidean distance di between ai and P∗ can be calculated by(25)di=minj=1P∗∑m=12fmai-fmajfmmax-fmmin2,where fmmax and fmmin, respectively, represent the maximum and minimum values of the mth objective function in P∗, and then the convergence metric C can be expressed as(26)CA∗≜∑i=1A∗diA∗.Note that the ideal Pareto frontier is always unknown in real problem; the algorithm proposed in this application is run 30 times, respectively, on each group of production data, and then a pseudo-Pareto frontier, which consists of all the solutions in 30 runs with removing dominated solutions, is constructed to compare with the approximate Pareto frontiers. For every run, box plots based on convergence metrics are illustrated in Figure7. In general, metric C in less than 10-2 means good statistical convergence performance in Pareto optimality based multiobjective optimization. The symbol “+” in Figure 7 refers to an outlier in box statistics; nevertheless, it can be seen that the outlier is very close to 10-2. Overall, we can see that the upper edges on different groups of data are all less than 10-2, except a slightly larger value on group 4 and an outlier on group 2. However, the 3rd quartile on group 4 is totally in the range of less than 10-2. The statistical results show that the proposed algorithm is stable in repeated run. On the whole, we can conclude that the proposed algorithm is robust and suitable for application in engineering.Figure 7
Box plots based on convergence metrics.


## 5. Conclusions

This paper presented the challenge of energy saving in hot rolling production and formulated a multiobjective optimization model of HRPSP under TOU electricity pricing. Objective of the model is to minimize electricity costs in production while considering penalties caused by jumps between adjacent slabs. Since exact algorithm is difficult to implement for solving the large scale problem, a NSGA-II based production scheduling algorithm was developed to obtain Pareto optimal solutions, and then TOPSIS decision-making method was adopted to recommend solution with different objective preferences. Experimental results and analyses showed that the proposed method cuts electricity costs in production, and the performance is better than load shifting on fixed production load. Consider that multiple production lines existed in most steel mills; HRPSP integrated multiple parallel machine job-shop scheduling will be the subject of further study, which is expected to have greater benefits. Besides that, multistage scheduling problem will also be our next work.



---

*Source: 1048081-2017-03-07.xml*