Very Large-Scale Neighborhood Search for Steel Hot Rolling Scheduling Problem With Slab Stack Shuffling Considerations

We study a steel hot rolling scheduling (HRS) problem considering slab stack shuffling (SSS), which is a kind of key issues in the steel strip production. The HRS problem is to select suitable slabs for a predefined sequence of hot rolling slots, from their respective candidate slab sets so that the number of shuffling operations of slabs is minimized. Different from the most researches, we consider a situation where there are intersections among candidate slab sets and shuffled slabs will not return original stacks. Basing on a special index method of slabs, we present an integer linear programming (ILP) to compute approximation solutions of the problem. In order to improve the solutions obtained by the ILP model, we propose a very large-scale neighborhood (VLSN) search algorithm. In the VLSN search process, the solutions of HRS problems are improved by an iteration procedure that partial solutions are interactively destroyed and repaired. In each iteration, partitioned hot rolling slots correspond to a sub-problem. In the repair operations, we propose an efficient branch-and-bound algorithm for solving sub-problem. The computational results on a number of different scale simulated instances show that the VLSN search algorithm is efficient for solving this kind of HRS problems.


I. INTRODUCTION
The steel hot rolling scheduling (HRS) is one of key tasks in operation management of steel production processes, which has an important impact on product quality and production efficiency of hot rolling processes. The HRS problem is a complex combinatorial optimization problem. It rises from an industrial problem where many practical constraints must be considered and thus hard to solve. Therefore, it has been studied extensively.
In early studies, the HRS problem was mainly converted into a kind of traveling salesman problems (TSP) or vehicle routing problems (VRP) to solve. The travel cost between nodes in the TSP/VRP corresponds to the penalty value on changes for the product width, thickness and hardness between pieces continuously rolled in the The associate editor coordinating the review of this manuscript and approving it for publication was Hisao Ishibuchi . hot rolling plan [1]- [3]. Pan et al. [4] considered multiobjective hot-rolling scheduling problems contain minimizing virtual sheet-strips and change times of the thickness at same time maximal changes in thickness between adjacent sheet-strips. Pan and Yang [5] took into account order delivery time, production capacity besides penalty values and proposed a variant of column generation algorithm to solve the problem. Witt and Voß [6] established a nonlinear optimization model for the HRS problem with minimizing makespan as optimization objective, and designed a parallel genetic algorithm to solve the problem. Lyu et al. [7] studied continuous casting and HRS problems, established a mixed integer programming model and proposed a heuristic algorithm to solve it. Zhao et al. [8]- [10] consider hot rolling scheduling problems for wire rod and bar products and propose two-stage decomposition method [8], iterated greedy algorithms [9] and memetic algorithms [10] to handle them. Zhang et al. [11] and Chen et al. [12] both studied hot rolling batch scheduling problems. Zhang et al. determine a sequence of the sheet strips with the objective of minimizing average thickness change and proposed a hybrid variable neighborhood search algorithm to solved the problem. Chen et al. decomposed the problem into a two-stage problem to solve. In the above researches, shuffling cost of slabs is not taken into account when slabs are charged into the furnace.
The operation efficiency of slab yard has a great impact on the production cost of hot rolling process. It is necessary to optimize the operation management that through hot rolling scheduling to reduce shuffling cost of slabs. In recent researches, some scholars begin to study HRS optimization considering slab stack shuffling (SSS). Tang et al. [13] established a nonlinear integer programming model for an HRS problem considering SSS and proposed a genetic algorithm. They assumed an ideal case that there is no intersection among candidate slab sets. Sing et al. [14] put forward an improved parallel genetic algorithm on the same problem. Wang et al. [15] proposed a two-stage heuristic for an HRS problem with SSS. Ren and Tang [16] established a nonlinear integer programming model for HRS problems with SSS considering crane capacity and converted the nonlinear model into a linear one. Different from existing studies, in this paper we study HRS problems considering SSS where there are intersections between candidate slab sets, and it is no need to put the shuffled slab back in place. To solve the problem, we establish an optimization model and propose a very large-scale neighborhood (VLSN) search algorithm.
The VLSN search algorithm was first proposed by Shaw et al. [17]. As a kind of metaheuristic, the VLSN search algorithm improves a solution by interactively destroying and repairing current feasible solutions. The VLSN search algorithm has been successfully applied to multiple combinational optimization problems. Breunig et al. [18] and Chen et al. [19] used the VLSN search algorithm to solve VRP. Sinclair et al. [20] implemented the VLSN search algorithm for an integrated aircraft and passenger recovery problem. Yu et al. [21] and Yagiura et al. [22] used this algorithm to solve assignment problems. Majid et al. [23] used this algorithm to solve supply chain network design problems. Brueggemann and Hurink [24] used this algorithm to solve a parallel machine scheduling problem. Guo et al. [25] and Fanjul-Peyro and Ruiz [26] used this algorithm to solve parallel machines scheduling problem. M. Zenker et al. [27] and Guo et al. [28] applied the VLSN search algorithm to railroad container terminals. In this paper, we apply the VLSN search algorithm to solve an HRS problem.
The main contributions of this paper include the following two parts: 1) We propose a special index method of slabs. By using it, the number of shuffles for selected slabs can be minimized by indirectly minimizing the sum of their serial numbers. In this way a hard procedure for calculating the number of shuffles can be avoided. 2) Based on the above idea, we establish a simplified integer linear programming (ILP) model that can approximately solve an HRS problem. This model can be optimally solved very fast by solving its linear relaxation since its constraint matrix is a totally unimodular matrix. 3) Taking the solution obtained by ILP model as an initial one, we propose a VLSN search algorithm to improve it in which a branch and bound (BAB) procedure is applied as a repairing process. We conduct extensive experiments to evaluate the proposed algorithm. The results show that the algorithm has great performance. Its high accuracy and fast speed imply its good potential to be used in practice.
The remaining part of this paper is organized as follows. Section 2 provides notations and an integer linear programming (ILP) model of the HRS problem. In Section 3, we propose a local search (LS) algorithm and a VLSN search algorithm. Section 4 presents the computational experiment results. In Section5, we summarize the research.

II. PROBLEM DESCRIPTION AND MATHEMATICAL MODEL A. PROBLEM DESCRIPTION
In a practical production process, there are four linkage modes between steelmaking-continuous casting process and hot rolling process, which are continuous casting -cold charge rolling (CC-CCR) mode, continuous casting -hot charge rolling (CC-HCR) mode, continuous casting -direct hot charge rolling (CC-DHCR) mode, and continuous casting -hot direct rolling (CC-HDR) mode [29]. Except for CC-CCR linkage mode, in the other three linkage modes the slabs do not need to shuffling operation. Therefore, in this paper, we consider CC-CCR linkage mode only, means that all slabs in an HRS problem come from a slab yard.
An HRS problem consists of N ordered rolling slots. A slot i (i = 1, . . . , N) corresponds to a set C i of candidate slabs from which we can select one slab for slot i. Slabs at the same rolling slot are similar in physical properties and order delivery date. Note that C l ∩ C k = ∅, 1 ≤ ∃l, k ≤ N , l = k. An HRS problem is to select a slab for each slot i from C i so as to minimize the number of shuffled slabs during the execution of hot rolling schedule. For an HRS problem including N hot rolling slots, the number of feasible solutions is about N i=1 |C i |, which grows exponentially with size of problems.
Slabs pile on stacks whose locations are fixed in a slab yard. Let s (s = 1, · · · , S) denote a slab, j (j = 1, · · · , P) denote a stack, and N s be the set of slots in which slab s can be applied. The number of slabs piled on stack j is called height H j of stack j. We number the slabs as following rules. The serial numbers of the top-level slabs are from 1 to P, then the serial numbers of the second top level slabs are from P + 1 to 2P, and so on. The level of slabs in stackj counts from the top to the bottom, so that the level of top slab is 1, and the level of bottom slab is H j . Lets = f (j, h) be slab number whose level is h in stack j.
During execution of hot rolling schedule, selected slabs are charged into reheating furnaces according to the rolled order to be heated to a specific temperature. Only the top slabs in stacks are allowed to be charged. If a target slab is VOLUME 9, 2021 not at the top of stacks, all the blocking slabs above it must be moved to a buffer area. Removing a slab into the buffer is called one shuffle. Shuffling operation is unproductive and occupy production resources. Shuffling operation not only increases crane workloads, but also increases the time of hot rolling processes such that is a bottleneck of the process. So, minimizing the number of the shuffled slab is taken as the optimization objective when we solve the HRS problems. The slabs in the buffer area will not be moved back to their original stacks and can be directly charged into heating furnaces without shuffles.
Example: We consider a small instance with N = 8 hot rolling slots, P = 6 stacks, and S = 37 slabs, as shown in Fig. 1. In Fig. 1, each box represents a slab. Grey boxes indicate candidate slabs. White boxes indicate slabs that are not involved in the schedule but may be shuffled. Slash boxes indicate slabs that are not involved in the schedule and cannot affect execution of the schedule. For example, C 1 = {1, 9, 19} is the set of candidate slabs for hot rolling slot 1. N 9 = {1, 7} is the set of slots which slab 9 can be applied in. The height of stack 1 is H 1 = 4.  Table 1 presents a feasible solution of the instance. In Table 1, the row label ''slot'' represents hot rolling slot number, ''selected slab'' represents the selected slab number, and ''shuffles'' represents the number of shuffles for selected slabs. In this solution, we select slab f (1, 1) = 1 for slot 1, and it can be charged directly without shuffle. We select slab f (3, 3) = 15 for slot 2. To charge slab 15, it needs to 2 shuffles that slabs 3 and 9 must be moved to the buffer area firstly. We select slab f (3, 1) = 3 for slot 3, and it do not need shuffle anymore since slab 3 has been moved to the buffer area. By analogy, we calculate the number of shuffles for the remaining selected slabs. The total number of shuffles for executing the schedule is 2 + 4 + 5 = 11.

B. MATHEMATICAL MODEL
We have established a mathematical programming model for the problem. The parameters and decision variables are defined as follows.

D s
The number of slabs upon slab s in original slab yard ϕ(s) stack that slab s is in where, the objective function (1) minimizes the number of shuffles. Constraints (2) guarantee that each slab is assigned to at most a hot rolling slot. Constraints (3) guarantee that each hot rolling slot selects exactly one slab from its candidate slab set. Constraints (5) denote that if slab s is selected for the ith hot rolling slot, the number of shuffles is dependent on the selected slabs for the previous (i − 1) hot rolling slot which is in the same stack with slab s [13]. Constraints (4)-(6) provide value domain of variables. The model is complex and hard to solve.

C. THE ILP MODEL FOR APPROXIMATION OPTIMIZATION
According to the rule which we number slabs, the number for involved slabs in HRS problems satisfy that max Fig. 1. Basing on the index method of slabs, it can be seen that the closer slabs are to the top of stacks, the smaller slab numbers are, and the less slabs shuffled will be. Therefore, the number of shuffles for selected slabs is indirectly reduced by minimizing the sum of serial numbers for selected slabs. Based on this observation, we propose the following ILP model for approximatively 47858 VOLUME 9, 2021 solving the HRS problem, which can effectively avoid the complex shuffling constraints.
x is ∈ {0, 1}, i = 1, · · · , N, s = 1, · · · , S where, the objective function (1) minimizes the sum of serial numbers for selected slabs. Constraints (8)-(9) are same with constraints (2)-(3). The ILP model is a classic assignment problem whose constraint matrix is a totally unimodular matrix. So, optimal solutions of ILP model are the same as optimal solutions for linear relaxation of ILP model. Therefore, ILP model is easy to solve and can be used to quickly generate an approximate optimal solution. Table 2 presents a solution X 2 for the instance in Fig. 1 through solving ILP model with CPLEX. The objective value of solution X 2 is 6.

III. HEURISTIC ALGORITHMS
In this section, we propose two heuristics for the HRS problem. One is a local search (LS) algorithm, the other is a VLSN search algorithm.

A. THE LS ALGORITHM
We use LS algorithm to solve the HRS problem, which is a common search method in engineering. In the LS process, let X denote a solution and X cur denote a current solution of the HRS problem, X [i] be the selected slab number for hot rolling slot i in X , and (X ) = {s | s = X [i] , i = 1, · · · , N }. Let cost (X ) be the objective value of solution X . Starting from the current solution X cur , each operation creates a neighbor solution X by reselecting a candidate slab for hot rolling slot i from the set C i . Compared cost (X cur ) and cost (X ), if an improving solution X is found, the current solution X cur is replaced by X . The search traverses candidate slab sets of all the hot rolling slots, stopping when none of the neighbors can improve the current solution. Fig. 2 shows a search operation for solution X 2 of the instance in Table 2. The procedure of the LS algorithm is as follows. Step 1: Get an initial solution X solving ILP model, X cur = X , let flag improved = true Step 3: Return X cur

B. THE VLSN SEARCH ALGORITHM
The VLSN search algorithm is based on the observation that searching a large neighborhood results in finding local optima of high quality. Starting from an initial feasible solution, the VLSN search algorithm iteratively improves the current solution by using destruction and reconstruction operations until the termination conditions are met. In the destruction operation, the algorithm releases some decision variable values in the current solution and retains the remaining variable assignments. In the reconstruction operation, the algorithm uses a specially designed procedure to reassign the remained variable values.
The proposed VLSN search algorithm starts from a feasible initial solution X . At every step, we release k hot rolling slots to reassign candidate slabs and the neighborhood is composed of all feasible assignments to released k hot rolling slots, while the selected slabs have been fixed for the remain N −k hot rolling slots. Then we make use of a BAB algorithm to optimally reassign slabs to these released hot rolling slots forming a neighboring solution X new . If the chosen neighboring solution has a lower cost of objective function than the current one, it becomes the new current solution. We set the number of iterations, stopping until exceeds the predefined limit. Fig. 3 shows a search operation for solution X 2 of instance in Table 2 and the neighborhood structure of the VLNS. Let sub k denote the sub problem involving k released hot rolling slots forming a set B sub k , and X sub k denote the optimal solution of sub k . All other symbols have the same meaning as described above. The procedure of the VLSN search algorithm is as follows.
Step 2: while (count ≤ maxUnsuccess) do a. Based on X cur , release k hot rolling; slot assignments to get a sub problem sub k . ; b. Solve sub k by the BAB algorithm to; get its optimal solution X sub k . ; c. Combine X cur and X sub k to get a; new complete solution X new .; If cost (X new ) < cost(X cur ) Then X cur := X new , count = 0 Else count++ Step 3: Return X cur In the following, we describe the destruction operators and reconstruction operators in details.

1) DESTRUCTION OPERATOR
For the HRS problem, destruction operator is to randomly release k hot rolling slots and reselect slabs for them. The operation is a key factor in the VLNS algorithm that the number of released hot rolling slots determines both the size of the neighborhood and the scale of subproblems.

2) RECONSTRUCTION OPERATOR
BAB algorithm is an exact method which is frequently used to solve combinatorial optimization problems [30], and we use it to solve a sub-problem sub k . The BAB algorithm adopts depth-first search strategy and the maximum depth of the tree is k. During the solving procedure for sub k , the objective   function value of nodes equals to the number of shuffles for reassigned slabs in the sub-problem and selected slabs in the main problem. Because if reassigned slabs for released hot rolling slots in sub-problem locate in same stack as selected slabs for retained hot rolling slots in the main problem, it will change the number of shuffles for the selected slabs. Since the current solution is feasible before the destruction operator is applied, we take the objective function value of the current solution as the upper bound (UB) of the BAB algorithm.
The BAB solution tree contains k + 1 levels. Level 0 represents the virtual root. Each level l(l = 1, 2 · · · , k) corresponds to a slot in sub k and each node denotes a candidate slab. Based on the index method of slabs in section 2.1, we state that the smaller the selected slab serial number is, the less number of shuffles for slabs. Therefore, the branched nodes are sorted for each node in non-decreasing order according to serial number of slab that the tree will be pruned earlier in the BAB process.
There are two rules to prune branches in the search tree. The one, comparing with UB, we prune the branches in which  the total number of shuffled slabs of assigned slots is bigger than UB. The other one is that branched slabs have been occupied by other slots, which means that we cannot find a feasible solution. For sub k , a k + 1-level tree can be formed containing i∈B sub k |C i | nodes at most by a BAB algorithm.
In each operator, both the scale of selected subproblem and the capacity of the candidate slab set for each hot rolling slot are small, therefore, the BAB algorithm is efficient for solving subproblems.
Example: We randomly select hot rolling slot 1, 3 and 8 to form a sub-problem of the instance in Section 2.1. Fig. 4 shows a BAB tree of the sub-problem. The current solution is X 1 in Table 1, namely UB = 11. In Fig. 2 each node h(h = 1, 2, · · · ) shows a branching status, the currentUB and the cost of the node B h , where ''i'' represents slot number and ''s'' represents slab number. There are 4 levels in the tree, where level 0 is virtual root as node 0 and the other three levels correspond to hot rolling slot 1, 3 and 8 in turn. The BAB algorithm starts from node 0, which is branched into three nodes, namely node 1, node 2 and node 3, each of which denotes a candidate slab for hot rolling slot 1 in level 1. According to selection strategy, we selected node 1 denoting the smallest number slab in set C 1 . By depth-first search strategy, since B 1 < UB, node 1 can further be branched into two nodes, namely node 4 and node 5. Since node 4 contains a smaller number slab in set C 3 and B 4 < UB, node 4 is branched into node 6, node 7 and node 8 in turn, all of which are leaf, and the corresponding objective function values are B 6 = 7, B 7 = 6 and B 8 = 10. The node 7 has the minimal objective function value B 7 = 6, so we update UB = B 7 = 6 and prune node 8. We traverse all the nodes forming the path 0−1−5, 0−2 and 0−3 in turn by backtrack method, all which are pruned since B 5 , B 2 , B 3 ≥ UB. As shown in Fig. 2, dotted line indicates the whole depth-first search procedure. Since the minimum objective value is 6(cost (X 1 ) > 6) in the sub -problem, we can conclude that slab 1, 3, 20 are respectively selected for slot 1, 3, 8 in sub-problem to repair X 1 forming a new solution X 3 as shown in Table 3 and the objective value is 6.

IV. COMPUTATIONAL RESULTS
In this section, we report our computational results of ILP model, the LS algorithm and the VLSN search algorithm. All the computational experiments were conducted on a desktop computer with 4 GB of RAM, the Windows 10 Professional 64-bit operating system, and an Intel Core i5 processor with four 3.3 GHz cores. Algorithms are implemented in C++ using CPLEX 12.61 and compiled with the Visual Studio 2013 C++ compiler. For all instances, the run time limit is set to 3,600 seconds for CPLEX solver.

A. INSTANCE GENERATION
In order to test the performance of the proposed model and algorithms, we generated 7 sets of instances simulating practical situations, each of which includes 10 instances. The HRS problems vary by the following parameters: the number of hot rolling slots(N ), the number of stacks (P), the height of stacks (H j , j = 1, . . . , P), and the capacity of candidate slab sets(|C i | , i = 1, . . . , N ). In our experiments, the parameters of each instance set are shown in Table 6.
In each instance set, there are three sizes of candidate slab sets, whose capacities are 2, 3, and 4 slabs, respectively. The quantitative proportion for three sizes of sets is 2:3:2. Each stack randomly stores 3-11 slabs and candidate slabs are randomly distributed in the stacks. The number of candidate slabs accounts for about 2/3 of the total in the yard.

B. COMPUTATIONAL RESULTS
We solve each instance using VLSN search algorithms with the different size of sub-problems (VLSN-k)k ∈ {3, 4, 5, 6, 7}, taking solutions obtained by the ILP model as initial solutions. The computational results are shown in Table 5. In Table 5, column ''obj'' represents the average objective function values on 10 instances of each set. Column ''CPU'' represents the average solution time in second on 10 instances of each set.
Based on the results shown in Table 5 we note that as sizes of sub-problems increase, the objective function values of each set are improved, but more time is taken. Balancing time and objective function values, it is more reasonable to use VLNS-5 to solve HRS problems.
To the best of our knowledge, no the same or similar problems have been considered by existing papers. To evaluate the effectiveness and accuracy of LS and VLNS algorithms, we have compared our algorithms with a constructive heuristic method that is currently used in real-world steel plants. In the constructive heuristic (CH), we randomly select a slab for each hot rolling slot from its candidate slab set forming a feasible schedule and calculate the number of shuffling operations. To avoid the randomness of the constructive heuristic, we run it 10000 times and select the best solution among the found ones as a competitor. The results of objective function value and running time are shown in column ''CH'' of Table 6.
Besides, Table 6 compares the performance of the LS and VLNS search algorithms. Both the LS and VLNS search algorithms take solutions obtained by the ILP model as initial solutions. We take the objective function value of the initial solutions as the benchmark, and the improved ratio of the LS algorithm and the VLSN-5 algorithm are shown in column ''impr'' of Table 6. The improved ratio is computed as follows.
The results in Table 6 show that compared to the constructive heuristic, both the model ILP, LS and VLNS algorithms have significant advantages in terms of time and objective function value. Besides, the results show that starting from the same initial solutions, the improvement of the VLSN search algorithm is better than the LS algorithm. But the runtimes of the VLSN search algorithms are longer than that of the LS algorithm. Balancing runtime and quality of solutions, the VLSN search algorithm significantly improves objective function values in acceptable runtime ranges. This shows the VLSN algorithm's accuracy and efficiency. In order to test the influence of different initial solutions on performances of the LS and VLSN search algorithms, we take randomly generated 5 feasible solutions as initial solutions for LS and VLSN-5 algorithm to solve a stance in each set. The computational results are shown in Table 7. Table 7 shows that as the objective function values of initial solutions increase, the objective function values of the LS algorithm increase significantly, but that of VLNS search algorithm have little change. It can be concluded that, comparing with the LS algorithm, quality of initial solutions has little influence on performance of the VLNS search algorithm, namely, the VLNS search algorithm is more robust than the LS.

V. CONCLUSION
In summary, we study HRS problems considering SSS. We build an integer linear programming to compute approximation solutions. Besides we propose the VLSN search algorithm to solve the HRS problems, which can effectively improve the quality of solutions with different scale of problems in acceptable time. At the same time, the VLSN search algorithm is robust since that does not depend on the quality of initial solutions. So, the proposed VLSN search algorithm is effective in both theoretical and practical production for HRS problems. In the future, we tend to explore other intelligence algorithms, e.g. monarch butterfly optimization, earthworm optimization algorithm, elephant herding optimization and moth search algorithm to solve the problem to further improve the solution quality and computing time. Meanwhile, we will attempt to find a suitable combinational optimization problem to which the HRS can convert and use it to obtain optimal solution of the problem. Besides, we will tend to design effective methods to obtain good lower bounds of the HRS.
YARONG SHI received the B.S. degree in automation from the Hunan University of Technology, Zhuzhou, China, in 2014, and the M.S. degree in systems engineering from Northeastern University, Shenyang, China, in 2016. She is currently pursuing the Ph.D. degree in systems engineering with Northeastern University, Shenyang, China. Her research focuses on planning and scheduling of steel production, mathematical programming, and heuristics algorithm.
SHIXIN LIU (Member, IEEE) received the B.S. degree in mechanical engineering from Southwest Jiaotong University, Sichuan, China, in 1990, and the M.S. and Ph.D. degrees in systems engineering from Northeastern University, Shenyang, China, in 1993 and 2000, respectively. He is currently a Professor with the College of Information Science and Engineering, Northeastern University, Shenyang, China. His research interests are in intelligent manufacturing, industrial big data, intelligent decision-making, production planning and scheduling. He has over 100 publications including one book. VOLUME 9, 2021