Late Acceptance Hill-Climbing Matheuristic for the General Lot Sizing and Scheduling Problem with Rich Constraints

This paper considers the general lot sizing and scheduling problem with rich constraints exemplified by means of rework and lifetime constraints for defective items (GLSP-RP), which finds numerous applications in industrial settings, for example, the food processing industry and the pharmaceutical industry. To address this problem, we propose the Late Acceptance Hill-climbing Matheuristic (LAHCM) as a novel solution framework that exploits and integrates the late acceptance hill climbing algorithm and exact approaches for speeding up the solution process in comparison to solving the problem by means of a general solver. The computational results show the benefits of incorporating exact approaches within the LAHCM template leading to high-quality solutions within short computational times.


Introduction
In many production processes, a fraction of the manufactured products can be defective due to an unreliable production system. Often, defective items incorporate substantial value and it is possible to rework them into an "as-good-as-new" state. Incorporating this feature results in a more ecological friendly production process as disposal reduces. Furthermore, reworking of defective items can only take place within a specific time span. For instance, in the pharmaceutical industry rework often has to take place in a short amount of time due to likely increasing contamination. In this context, we formulate a related problem and address a model formulation of the general lot sizing and scheduling problem with rework termed as GLSP-RP, which was first presented by [1]. This problem also considers lifetime constraints for defective items. As the authors report, the GLSP-RP model is hardly solvable by a general-purpose solver in an adequate amount of time. Even for small problem instances of non-practical size, it requires more than 1800 s. This leads to the necessity of faster approaches providing solutions within reasonable execution times.
Matheuristics can be defined as those solution approaches obtained from the interoperation between metaheuristics and mathematical programming techniques [2]. Within matheuristics, decomposition approaches belong to one of the most relevant algorithms (some recent examples can be consulted in, for example, [3][4][5]). For instance, the principal idea behind proposing the matheuristic version of POPMUSIC is to decompose the problems into smaller (and likely more tractable) ones and solve them by means of an exact technique [6]. As remarked by [7] to successfully apply decomposition matheuristic approaches, one may develop an iterated procedure where the information is interchanged along the process. Adapting those ideas to well-known single point metaheuristics such as late acceptance hill-climbing (LAHC, [8]) permits, within a hill-climbing framework, to substitute the inherent method for generating candidate solutions by an exact method solving a given sub-problem. Furthermore, recently, [9] showed by means of the unrelated parallel machine scheduling problem the importance of including advanced local search methods to improve the solution quality in LAHC. This, thus, led to our research question about how LAHC would perform if we use an exact approach instead of a local search solver.
Based on the previous discussion, the contribution of this paper is twofold: • We investigate the relationship between the original formulation of the general lot sizing and scheduling problem (GLSP) and the GLSP-RP by providing a comprehensive discussion of both problems.

•
For solving the GLSP-RP, we propose the Late Acceptance Hill-climbing Matheuristic (LAHCM) as a general solution framework, which is inspired by and integrates the LAHC strategy with the exact approaches, i.e., the solution of optimization models for reduced problems configured in the spirit of the LAHC. To assess its performance, we evaluate it over a set of well-defined instances.
The computational results indicate that by means of our approach we are able to improve the objective function values for several instances where a general-purpose solver such as CPLEX (https://www.ibm.com/nl-en/analytics/cplex-optimizer) is not able to find optimal solutions in the given time limit of 1800 s.
The outline of this paper is as follows. In Section 2, those works related to the GLSP as well as its respective formulation are discussed. The general lot sizing and scheduling problem considering rework and lifetime constraints for defective items is explained in Section 3. Later, our matheuristic solution approach based on late acceptance hill-climbing is described in Section 4. The extensive computational experience is reported in Section 5. Finally, some conclusions and further research lines are exposed in Section 6.

Literature Review
To jointly decide on lot sizing and scheduling has attracted many researchers in recent years producing a rich body of literature on different models and problem formulations including, for example, the general lot sizing and scheduling problem (GLSP). In this paper, we present a novel solution approach for an extended version of the GLSP and focus on this specific problem while referring to [10], and [11] for extensive literature reviews on different simultaneous lot sizing and scheduling problems. In Section 2.1, we reappraise the original model formulation of the GLSP as it will later be used in the algorithm for creating an initial solution. In Section 2.2, we review solution approaches and model extensions to the GLSP.

The General Lot Sizing and Scheduling Problem
The GLSP was originally introduced by [12]. The authors presented two versions of the model, namely the GLSP-CS (Conservation of Setup state) where the setup state is preserved over idle periods and the GLSP-LS (Loss of Setup state). In this section, with the aim of making the paper self-contained, we present an extended mixed-integer programming (MIP) formulation for the GLSP-CS that includes sequence-dependent setup times as presented by [13] and a revised version of the minimum lot size constraint introduced by [14] and use the name GLSP for this formulation in the further course of this paper. The GLSP features integrated production planning and scheduling for capacitated single facility environments with multiple products. The setups for these products are sequence-dependent in times and costs, implying that setting up the production of the products in a certain sequence can be more favorable than another in terms of time and costs; see also [15,16]. To simultaneously determine the production sequence and the quantity of production lots, the GLSP includes a two-level time structure. Products j = 1, . . . , J are scheduled over a finite planning horizon consisting of macro-periods t = 1, . . . , T with given length. Micro-periods m = 1, . . . , L T (with L T denoting the last micro-period of the last macro-period T) are used to model the production sequence of products within the macro-periods. Each macro-period t consists of a predefined non-overlapping sequence of micro-periods with |M t | being an advance fixed total number of micro-periods within a macro-period t. The length of each macro-period is given by the input parameter for the resource capacity b t ; the lengths of micro-periods are decision variables expressed by the number of productions in the micro-period. A production lot consists of a sequence of consecutive micro-periods where the same product is produced and may continue over a series of micro-periods and macro-periods. As described, this model presents an extended version of the GLSP-CS where the setup state is preserved over idle periods, i.e., if after an idle micro-period the same product is produced as before, no additional setup is necessary.
The objective is to find an optimal production quantity and a sequence that minimizes overall production costs, consisting of setup changeover and inventory-holding costs. Table 1 presents the variables, parameters and indices of the GLSP model. Table 1. Variables, parameters and indices used in the model.

i, j
Products with i, j = 1, . . . , J L t Denoting the last micro-period of a macro-period t m Micro-periods m = 1, . . . , Demand of product j in period t f i,j Sequence-dependent setup costs for a change over from product i to j h j Inventory holding cost factor for product j κ j Minimum lot size amount of product j |M t | Fixed number of micro-periods within a macro-period t st i,j setup times for a change over from product i to j tp j Process time per unit of product j M big number, e.g., max(b) t∈T Variables s j,t Inventory of product j in period t x j,m Binary setup variable, 1 if product j is set up for production in micro-period m, 0 otherwise y j,m Production amount of product j in micro-period m z i,j,m Changeover variable, 1 if production is changed from product i to j in micro-period m, 0 otherwise The problem is formulated as follows: The objective Function (1) minimizes the overall costs consisting of the holding costs for produced items and the sequence-dependent setup costs. Equation (2) depicts the inventory balancing constraint. It states that the inventory in macro-period t is composed of the inventory of the previous macro-period t − 1, the production amount of all micro-periods belonging to t subtracted by the demand of period t. Inequality (3) states that the maximum available capacity (b t ) of macro-period t reduced by the sum of sequence-dependent setup times incurred by setup changes from product i to product j with m ∈ M t must be greater or equal to the sum of processing times in all micro-periods belonging to macro-period t. Constraints (4) and (5) ensure that product j can only be produced in micro-period m, if a setup for product j is performed in m, i.e., if the binary setup variable (x j,m ) is set to one and that only a single product can be produced at a time in a micro-period m. The authors of [12] pointed out that minimum lot sizes have to be included, otherwise wrong calculations of setup costs could occur as a result of setup state changes without product changes in cases where the setup cost matrix does not satisfy the triangle inequality (i.e., f i,k + f k,j ≥ f ij with i, j, k ∈ 1, . . . , J). As shown by [14], this originally introduced minimum lot size constraint was strictly limited as a production state between two consecutive periods; it is only conserved if the available capacity exceeds the minimum production quantity. Therefore, the authors adjusted Constraint (6) and introduced Constraint (7) to enforce the production continuity across macro-periods. Constraint (8) presents the changeover logic and Constraints (9) and (10) enforce non-negative and binary variables, respectively.
Many researchers have worked on this specific model and developed extensions, solution approaches, and applications for industrial settings, which will be reviewed in the next section.

Solution Approaches and Model Extensions
The authors of [12] present a solution procedure for the GLSP-based on the local search algorithm threshold accepting. The procedure starts from a configuration defined by a fixed setup pattern and calculates the corresponding lot sizes and costs through backward oriented heuristics. New candidates are created by neighborhood operations that change the old configuration and are accepted as a new configuration if they provide an improved solution. To escape local optima, new candidates are also accepted if their solution quality is not worse than a specific threshold, while the convergence of the algorithm is ensured by lowering this threshold. In [13,17], an improved version of this procedure is used that determines lot sizes optimally for a fixed pattern through "dual reoptimization" instead of heuristically for solving two extended versions of the GLSP, i.e., the GLSP with setup times and the GLSP for parallel lines, respectively. While this approach proves to be efficient for the single-line GLSP, the results of the multi-line formulation are insufficient. Therefore, Meyr and Mann (2013) adjust the approach by aggregating the original multi-line "master-problem" in order to decompose it into a set of isolated single-line "sub-problems", which are solved by the heuristic provided in [13]. A model extension of the single-stage GLSP to multiple production stages is presented in [18]. The authors of [19] adjust this model with respect to the synchronization of the machines and present several reformulations of the problem based on variable redefinitions. In addition, three heuristic approaches are presented to solve the problem. In [20], the authors present a GLSP-based model formulation called the synchronized and integrated two-level lot sizing and scheduling problem (SITLSP), which is motivated by industrial settings, mainly of soft-drink companies where the production process involves two interdependent levels with decisions concerning raw material storage and soft drink bottling. The authors of [21] solve the SITLSP by a genetic algorithm (GA) and study single-population and multi-population approaches. Additional solution approaches for the SITLSP are studied by [22]. They apply a tabu search (TS), threshold accepting (TA), genetic algorithm (GA) and hybrid genetic algorithms using a combination of GA with TS and TA, respectively, to the problem. The authors observe that for the hybrid approaches, an execution of TA and TS over the best individual of each population results in high computational times and that a strategy where these strategies are only executed when the best individual of all populations is updated performs best. The authors of [23] propose a multi-population GA with individuals structured in ternary trees to solve data provided by a soft-drink company. A simplified model formulation of the SITLSP called two-stage multi-machine lot scheduling model (P2SMM) is proposed by [24] where each filling line can only be connected to a single tank at the same point in time. As a relaxation approach, they first solve a single-stage version of the model and afterward determine the setup variables of the P2SMM. In addition, they propose fix-and-relax strategies and combinations of them with a relaxation approach to solve the P2SMM. Reference [25] study different relax-and-fix strategies to solve a model formulation similar to the single-stage P2SMM with only one production line and the production bottleneck in the bottling stage. Reference [26] propose four different alternative single-stage formulations of the P2SMM where the first two are based on the single-stage GLSP model with sequence-dependent setup times and costs, while the other two are asymmetric traveling salesman problem (ATSP)-based formulations with different subtour elimination constraints. Numerical experiments show that solution strategies based on the ATSP formulations perform best for the P2SMM. The authors of [27] apply a combined genetic algorithm and mathematical programming approach to solve the P2SMM that outperforms the approaches presented in [26]. In [28], the authors adjust the P2SMM for an application in the brewery industry and solve the problem heuristically by using a constructive relax-and-fix heuristic and several improvement procedures based on fix-and-optimize strategies. Note that these approaches borrow ideas from the much older POPMUSIC approach of [29].
The authors of [30] present a GLSP-based model formulation for a multi-product multi-level job shop production where they subdivide the macro-periods into three types of micro-periods: one type for production, one for setups, and one for idle time. Due to the high complexity of their formulation, only problems with a small number of periods, products, and machines can be solved to optimality. In [31], the authors propose a formulation for lot sizing and scheduling in a flow shop environment with sequence-dependent setups and develop four fix-and-relax heuristics that rely on successive resolutions of reduced size MIP-formulations of the problem. The authors of [32] propose a formulation based on the model of [30] for a flexible flow shop environment with sequence-dependent setups. Additionally, the authors of [33] develop two algorithms based on a rolling horizon approach where first, all binary variables of the problem are relaxed. Then, relaxed binary variables of the current period are divided into the two groups so that the members of the first and second groups get value one or zero, respectively.
The authors of [34] present two MIP-formulations for joint lot sizing and scheduling motivated by an animal feed plant. While the first formulation sequences each period independently, the second model takes linking setup states between periods into account. The authors present three relax-and-fix (RF) heuristic approaches where they either apply the RF heuristic for relaxing and fixing the lot sizes or the periods. While the RF heuristic on the integer lot sizes performs best, the authors observed that a forward-oriented RF over the periods acts like a greedy heuristic resulting occasionally in infeasible solutions for instances known to be feasible. In [35], the author presents an extended formulation for the GLSP with multiple production stages where an embedded State-Task-Network is used to model product substitution and flexible Bill-of-Materials by introducing several tasks that produce the same product while using different input products. In [36], the authors consider a GLSP-based model with different machine configurations. Therefore, the authors introduce tools as an additional, capacitated resource in a single-stage, multi-machine problem and model tool changeovers instead of product changeovers. As a specific tool cannot be used on more than one machine at a time, continuous variables for the starting and ending times of tool use are introduced. The authors of [37] present a hybrid mixed-binary model based on the GLSP for a practical problem of the process industry and propose two transportation-based reformulations, namely the quantity-based transportation problem (QTP) and the proportional transportation problem (PTP). In their statistical analysis, they observe that the PTP reformulation performs best on average, but for cases of a high level of minimum production quantities, it is outperformed by the QTP reformulation. In [38], the authors propose three novel mathematical formulations for a two-stage GLSP where multiple products are produced on parallel, non-identical machines inspired by the textile industry and solve them with a MIP-solver. The authors of [39] adapt this model for usage in the yarn production and solve the problem by means of a hybrid method called hamming-oriented partition search, which is a branch-and-bound-based procedure that incorporates a fix-and-optimize improvement method.
In [40], the authors propose a GLSP-based model for a pulp and paper mill planning problem and solve the problem with a two-step approach with a MIP-based construction heuristic followed by a MIP-based improvement heuristic. The authors of [41] extend this model to the multi-stage scenario and include additional constraints with the objective to maximize the steam output used to generate electrical energy. To solve the problem, the authors develop a hybrid approach that combines the general variable neighborhood search heuristic, a specific heuristic called speeds constraint heuristic and an exact solver. In [42], the authors also extend the model of [40] to the case of parallel non-identical paper machines and solve the problem by a hybrid solution method composed of an exact method embedded in a genetic algorithm. The authors of [43] present a model formulation for a production environment with wear and tear of resources and maintenance activities based on a combination of the PLSP and the GLSP where exogenous macro-periods with given length and demand are combined with endogenous micro-periods of flexible length. The authors propose a fix-and-optimize type decomposition heuristic to solve this problem.
The authors of [44] present a GLSP formulation with multiple machines in a job shop environment and solve the problem with two MIP-based algorithms based on the rolling horizon iterative procedure. A multi-stage GLSP formulation for a flexible job-shop problem is presented by [45] and solved with a metaheuristic method that combines a GA with particle swarm optimization (PSO) and a local search heuristic. Deterioration and perishability for the GLSP are considered by [46]. The authors propose a model formulation that includes products that deteriorate after a maximum lifetime while being stored in an inventory. In this paper, we regard a model formulation based on this GLSP version, which was first presented by [1] and combines lifetime constraints for defective items with the ideas of an imperfect production process and rework, called GLSP-RP.
Reference [47] study a robust optimization and a scenario-based two-stage stochastic programming model for the General Lot-Sizing and Scheduling Problem (GLSP) under demand uncertainty. The authors propose an extensive simulation experiment based on Monte Carlo to evaluate different characteristics of the solutions, such as average costs, worst-case costs, and standard deviation. The authors of [48] focus on a general lot-sizing and scheduling problem, including perishable food products, and apply two mixed-integer programming based heuristics to the problem. The authors use relax-and-fix approaches as a constructive heuristic to find an initial feasible solution and afterward apply a fix-and-optimize heuristic to improve the obtained initial solution.

The GLSP with Rework and Lifetime Constraint for Defective Items
The GLSP-RP was introduced by [1] to model process environments where conflicting items are required to avoid contact as they are subject to contamination, for example, the food processing industry where an-allergenic food items are produced on the same machine as regular food items. The problem considers an imperfect production system where a fraction of the planned lot size is not of serviceable quality and goes to a separate rework inventory where it is stored till it is reworked and shipped to customers or disposed of, which is illustrated in Figure 1. Figure 1. Illustration of the defective production process (see [49]).
Note that, in this paper, no distinction will be made between the quality of good quality items and reworked items. Furthermore, we assume that the rework process always succeeds. For scheduling rework lots, a joint setup variable is considered, i.e., if the machine is set up for the original production in micro-period m (i.e., x j,m = 1) and followed by a rework lot for product j in m + 1 and vice versa, no additional setup is needed and, therefore, no additional setup costs occur. To address this in the model, the binary setup variable x j,m of the GLSP formulation is redefined as the binary joint setup-variable. In addition, the defectives have a fixed lifetime while being stored in the rework inventory, i.e., if the defective items are not reworked in a fixed time frame expressed by a number of consecutive micro-periods (Ω j ) they become scrap and have to be disposed, resulting in disposal costs (λ j ). Table 2 shows the additional variables and parameters of the GLSP-RP. Table 2. Additional variables and parameters for the general lot sizing and scheduling problem with defective items and lifetime constraints (GLSP-RP).
Parameters h R j Inventory holding cost factor for defective items of product j λ j Disposal cost for one perished item of product j Ω j Lifetime of product j in micro-periods Θ j,t Proportion of defective items 0 ≤ Θ j,t < 1 of product j in t tp R j Process time for reworking a defective item per unit of product j The GLSP-RP is formulated as follows: The adjusted objective Function (11) of the GLSP-RP additionally considers the holding costs of defective items and the disposal costs for defective items that perished while holding. As a positive rework inventory at the end of the planning horizon is allowed for feasibility reasons any positive amount of defectives stored in the rework inventory or incurred in the last micro-period of the planning horizon (m j,m = L T ) has to be disposed, leading to disposal costs that are included in the last term of the objective function. Equations (12) and (13) present the production process of defectives as described in [49]. Concerning the planned lot size amount y j,m of a given micro-period m, only the reduced amount of y S j,m is of serviceable quality and can be used for demand satisfaction. The amount of items that are defective (R j,m ) has to be stored in the rework inventory and later reworked or disposed. Consequently, the GLSP-RP regards two separate inventories for serviceable items and defective items. Equation (14) shows the amended inventory balancing constraint for serviceable items. The serviceable inventory in macro-period t is composed of the inventory of the previous period t − 1, the amount of serviceable and rework production in all micro-periods belonging to macro-period t subtracted by the demand of period t. Equation (15) considers the rework inventory. The rework inventory is micro-period based as it symbolizes an interim inventory for the defective items. It is composed of the rework inventory in m − 1, the number of defective items of micro-period m minus the reworked items in m and the number of items that passed their lifetime while being stored (s D j,m ). The interrelation between storage time and lifetime of the defective items is shown in Equation (16). It states that the amount of spoiled defectives for a given micro-period m is calculated by the defectives produced by the original production process from the beginning of the planning horizon until micro-period m deduced by the lifetime length Ω j of the defectives, less those items that are already reworked until micro-period m − 1 and defectives that were already disposed in previous micro-periods. Constraint (17) ensures that defectives of micro-period m − 1 can be reworked in m the earliest. The capacity constraints given in Inequalities (18) show that in addition to the sum of the processing times for the original production, the rework production process also has to be taken into account. Consequently, it has the maximum available capacity (b t ) of macro-period t minus the sum of sequence dependent setup times by setup changes from product i to product j not to be exceeded by the total sum of processing times for the original and the rework production occurring in micro-periods belonging to t. Equation (19) depicts the joint setup constraint. This constraint ensures that product j can only be produced or reworked, respectively, in micro-period m, if a setup for product j is performed in m, with UB symbolizing a large number that has to be chosen in a way that it does not restrict the production volume. Constraint (20) ensures that only a single product is produced and/or reworked in micro-period m. Constraints (21) and (22) are the adjusted minimum lot size constraints and ensure that after a setup change the minimum lot size amount (κ j ) has to be fulfilled by production and/or rework. Constraint (22) enforces the production continuity across macro-periods as described by [14], i.e., if the time interval of the last micro-period of a macro-period t is not sufficient to produce or rework κ j , and production is continued in the next macro-period (t + 1), then the sum of production and rework of micro-period m and m + 1 has to satisfy the minimum production quantity. Constraint (23) presents the changeover logic and enforces variables (z i,j,m ) for micro-period m to equal one if product i is produced or reworked in micro-period m − 1 and followed by a production or rework lot of another product j in micro-period m. Constraints (24) and (25) enforce non-negative and binary variables, respectively.

Late Acceptance Hill-Climbing Matheuristic Template
The late acceptance hill-climbing (LAHC) algorithm is a single-point iterative search algorithm proposed by [8] that accepts non-improving movements when the objective function value of a candidate has a better value than the one number of iterations before. Namely, while in the standard hill-climbing algorithm a candidate solution is compared to the direct predecessor, in LAHC the candidate is compared with the current solution several iterations before.
LAHC has been successfully applied to several optimization problems. In the following, some of them are briefly reviewed. The authors of [8] apply the LAHC to exam timetabling problems reporting a competitive performance. This problem is later addressed by [50] where they present a hybridization of an adaptive artificial bee colony (ABC) and LAHC. The latter is applied as a local search procedure during the exploitation of the ABC algorithm. The results indicate that this hybridization performs better than other modifications of the ABC algorithm. In [51], the authors apply the LAHC to the traveling purchaser problem showing a suitable performance for solving that problem. The authors of [52] use this algorithm for balancing two-sided assembly lines with multiple constraints. The LAHC, in this case, is able to provide the same solutions as those provided by the optimization model implemented using LINGO for the small-sized instances. For large-sized instances, the model is not able to provide an optimal solution in reasonable computational time, thus leaving LAHC the best option for those instances. The author of [53] applies LAHC to the liner shipping fleet reposition problem. The author compares it to simulated annealing (SA) reporting that LAHC performs well on small instances, but on larger instances, SA exhibits better performance. The authors of [54] experimentally examine their proposed method and compare it to search techniques that employ a cooling schedule, i.e., simulated annealing, threshold accepting and the great deluge algorithm for the traveling salesman problem. In [55], the authors present a cutoff time strategy based on the coupon collector's problem and apply it to the LAHC algorithm. The authors tested that approach on the Travelling Salesman Problem, the Quadratic Assignment Problem, and the Permutation Flow-shop Scheduling Problem. Their results show that the proposed strategy is a valid stopping method for a one-point stochastic local search algorithm that accepts worsening moves. The authors of [56] apply the LAHC for finding train composition causing greatest fatigue damage in railway bridges. In [9], the authors propose an LAHC approach to solve the unrelated parallel machine scheduling problem with sequence and machine-dependent setup times, the reported results indicated the relevancy of including advanced local search methods to improve LAHC performance. This raises the research question related to using an exact solver instead of a local search within LAHC.
Based on the previous discussion, there are two reasons for developing and proposing the Late Acceptance Hill-Climbing Matheuristic (LAHCM) in this paper. First, the LAHC as a metaheuristic reportedly performs well for the traveling purchaser problem and the traveling salesman problem. The traveling purchaser problem incorporates a choice option regarding the inclusion of nodes. As including rework into dynamic lot sizing shows a high similarity to this property and lot sizing models are frequently formulated as traveling salesman based formulations including the GLSP (see [26]), it is worth examining the application of LAHC for solving our GLSP model formulation. Second, considering the relevant performance of this strategy in the related literature and the benefits of exact approaches in terms of solution quality, the interplay between the two types of approaches is considered in LAHCM. In this sense, exact approaches are exploited to solve subproblems optimally and generate candidate solutions within LAHC. To the best of our knowledge, this is the first LAHC matheuristic template proposed so far and, therefore, we are contributing in this paper by presenting a LAHCM as a general solution framework.
The pseudo-code of LAHCM considering the GLSP-RP is reported in Algorithm 1. The parameter of the algorithm is a list length parameter l, which establishes the number of solutions to be stored. The initialization of the algorithm consists of creating an initial solution for the problem instance p (line 1). The current best solution is set to the initial solution (line 2), the number of iterations and the array for storing the solutions collected in the last l iterations are initialized in lines 3-6. The LAHCM search procedure is depicted in lines 7-15 and executed until a given stopping criterion is met. The index of the solution from l iterations before is stored in v (line 8). A candidate solution s is generated by solving a subproblem by means of an exact method (fixandsolve()) defined by the user (line 9). At this point, it should be noted that the main difference between a standard LAHC and LAHCM resides in the way the candidate solution is generated. Furthermore, if the candidate solution enhances the objective value of the current solution or the one from l iterations before, then the candidate solution is the new current solution (lines 10-11). Moreover, the candidate solution is compared to the best solution found so far s best in order to check if the candidate solutions improve the best one and thus replace it (lines 12-13). Finally, the solution is stored in the solution array and the iteration counter is incremented. In LAHCM, the search process is stopped when no improvement is obtained by means of the performance comparison shown in line 10.
As can be seen, the LAHCM template does not require a neighborhood structure to generate candidate solutions but an exact approach permitting to optimally generate candidate solutions for its mathematical programming neighborhood that can be defined by fixing some parameters or reducing the instance at hand. The previous permits, because of the size reduction, solving the subproblems with exact approaches in reasonable times to optimality.

Initial Solution
In order to generate an initial solution for the LAHCM algorithm, the procedure createinitialsol() (see Algorithm 1) is proposed. Its functioning consists of relaxing the problem by not considering the constraints concerning the rework and lifetime for defective items, that is to say, we solve the GLSP (see Section 2.1). Once that problem is solved, the solution is appropriately entered by means of fixing the sequence variables x j,m of the GLSP-RP. That is done by extending its corresponding model (Section 3) with the following constraints: where x(GLSP) j,m denote the variable values taken in the GLSP solution.

return s best
By means of this procedure, an initial solution is obtained in a reasonable computational time. It is worth pointing out that the solution entered to the GLSP-RP optimization model (Section 3) has the binary setup variables fixed. Since the initial solution comes from the GLSP, the selection of the variables is performed without losing the sense of feasibility. With the goal of facilitating the understanding of the initial solution generation process, an illustrative example considering three products (j = 1, . . . , 3), three macro-periods (t = 1, . . . , 3), and a fixed number of five micro-periods (M t = 5) per macro-period is presented. The available capacity for all macro-periods is constant with b t = 400. For all products, the processing times per product are tp j = 1, the holding costs are h j = 5, and the minimum lot size amount is κ j = 10. The other parameters of this example are set as follows: The first step for generating an initial solution is to solve the corresponding GLSP. This way, for the above-mentioned example the optimal solution is reported in Table 3. It should be noted that since the values of variables s j,t correspond to macro-periods, to avoid repetition only one number is reported for each row/block.
In the solution, we observe six setup changes with a total of 15.75 in setup changing costs. As the available capacity in t = 2 with b t = 400 is not sufficient to produce the total demand for all three products in t = 2, proportions of the demand of products j = 1 and j = 2 are produced and stored in t = 1 leading to storage costs, of 410. Therefore, the overall cost for the GLSP solution is 425.75.
As previously described, additional parameters are needed for the GLSP-RP. Consequently, the lifetime of defective items while stored in the rework inventory is set to three micro-periods for all products (Ω j,m = 3) and the disposal cost is λ j = 1000 per item. Regarding rework, we include rework holding costs of h R j = 1, rework processing costs of tp R j = 1 for all j. The fraction of defectives per product and macro-period is given as follows: The solution for the GLSP-RP with the fixed setup sequence is shown in Table 4. As in Table 3, variables s j,t only report for each row/block one value as they correspond to macro-periods.  The above solution is used as the starting solution for solving the GLSP-RP with the LAHCM approach. Its objective function value is 4478.75. As can be checked, the setup costs of 15.75 for the GLSP-RP are the same as for the GLSP due to the same setup sequence. The imperfect production process results in much higher storage costs with 1460 costs for the inventory storage and an additional three costs for the rework inventory. The greatest cost driver by far are the disposal costs. As a result of the fixed production sequence, not all defective items can be reworked in a timely manner and therefore disposal costs of 3000 occur.

Mathematical Programming-Based Neighborhood
As discussed in the previous section, LAHCM differentiates from the standard LAHC in the way the candidate solutions are generated. While in the LAHC a neighborhood structure is usually used, in LAHCM generation of the candidate solution is done by means of an exact algorithm. Thus, in order to produce a candidate solution for the GLSP-RP, the method fixandsolve() is proposed. By means of this function, an input solution is used for extracting and fixing decision variables with the aim of generating smaller sub-problems and solving them in an exact way, producing thus a candidate solution. In a similar way as proposed in [6] through the definition of reduced sub-problems from larger and more complex ones, we are able to address them using an exact algorithm in a small and reasonable computational time. This, as discussed in the relevant section, allows improving the quality of the provided solution while reducing the overall solving time. The input of our approach considers different strategies in order to obtain a reduced sub-problem. In the current work, we use the following strategies: Strategy 1 (α 1 ): The x j,m variables are fixed except for those regarding a product j selected at random. Strategy 2 (α 2 ):The x j,m variables are fixed except for those regarding two products j and j selected at random and where j = j . Strategy 3 (α 3 ):The x j,m variables are fixed except for those regarding three products j, j , and j selected at random and where j = j = j .
Henceforth, the resulting subset of non-fixed variables of x j,m is denoted as x f j,m . Moreover, the application of these strategies requires previous information about those variables that are to be fixed. That is done by considering previous information from a given solution s.
The pseudo-code of fixandsolve() is depicted in Algorithm 2. For a given problem instance of the GLSP-RP, the algorithm fixes some of the decisions associated to the variables x taking into account an input solution s. To determine the subset of variables to be free from x (i.e., x f ), a strategy is selected at random (line 2). Afterward, the optimization model of the GLSP-RP is solved for the ensuing sub-problem obtained from the input instance p with the decisions concerning x fixed except for x f (line 3). Hence, the model to be solved includes the following: The solution s obtained from solving the reduced problem is compared to the input one to check if an improvement has been achieved (line 4). If no improvement is achieved, the input solution s is returned (line 7).

Algorithm 2: Fix and Solve algorithm for the GLSP-RP
Data: Problem instance of the GLSP-RP, p Input: Initial solution s, strategies α = {α 1 , ..., α 3 } Output: Updated solution for the GLSP-RP With the goal of facilitating the comprehension of LAHCM, this subsection presents an example case of its functioning using the problem instance provided in Section 4.1 and considering a list length (l) of two.
Initial Solution: As detailed in Section 4.1, for creating a starting solution, the production sequence of the GLSP is given to the GLSP-RP, i.e., the complete sequence of x j,m of the GLSP is given to the GLSP-RP with no rescheduling enabled. Table 5 shows that the resulting sequence of the initial GLSP-RP solution matches the sequence of the GLSP. Table 5. Production sequence of the initial solution for the GLSP-RP. The objective function value of the initial solution (i.e., 4478.75) is stored in the list l = [4478.75] and the iterative procedure with relaxation of the sequence is started.
Iteration 1: The strategy selected at random is α 2 and the selected products within the strategy are j = 1 and 2.
Therefore, the binary setup variables for j = 1 and 2 are reset to zero for all micro-periods m, which is once exemplified in Table 6. This allows a rescheduling of the sequence of products 1 and 2. Table 6. Relaxed production sequence of the initial solution for the GLSP-RP. The above sequence is given to the GLSP-RP, producing thus a reduced problem, which is solved by means of an exact algorithm. The optimal solution of this problem results in the following new production sequence for iteration 1: Having the sequence of products j = 1 and j = 2 free while the remaining products are fixed leads to four rescheduling decisions (shown in bold in Table 7) and results in an objective function value of 998.25. This new solution does not include disposal costs. Moreover, this objective value is added to the list l = [4478.75; 998.25] and the updated production sequence will be the new input sequence for iteration 2. As the list is filled (with the predefined list length l = 2), the next solutions of the next iterations will be compared to the first entry of the list until the stopping criterion is met. Table 7. Output production sequence of Iteration 1 for the GLSP-RP.

Iteration 2:
The strategy selected at random is α 1 and the selected product within the strategy is j = 1.
All values of the production sequence of iteration 1 are fixed except for product j = 1. This is given as an input to the GLSP-RP. This reduced problem does not lead to any rescheduling and, therefore, results in the same objective function value of 998.25, which is lower than the first entry of the list (i.e., 998.25 < 4478.75). Consequently, the first entry of the list is dropped and the objective function value of iteration 2 is added: l = [998. 25; 998.25].
Iteration 3: The strategy selected at random is α 1 and the selected product within the strategy is j = 3. In this iteration of the algorithm, the products' sequence is fixed except for product j = 3. The optimal solution of the resulting reduced problem leads to a solution for the GLSP-RP with two rescheduling decisions, which are shown in Table 8. Table 8. Production sequence of Iteration 3 for the GLSP-RP. The new objective function value is 994.25, which is lower than the corresponding entry of the list (i.e., 994.25 < 998.25). Therefore, that entry of the list is dropped and the objective function value of iteration 2 is added to the end of the list: l = [998. 25; 994.25].

Iteration 4:
The strategy selected at random is α 1 and the selected product within the strategy is All values of the production sequence of iteration 3 are fixed except for product j = 2. Having the decision concerning this product free does not lead to any improvement and thus the objective function value of 994.25 is maintained. Nevertheless, this objective function value is lower than the corresponding one of the list (i.e., 994.25 < 998.25). Therefore, the entry of the list is dropped and the objective function value of iteration 4 is added to the end of the list: l = [994. 25; 994.25].
Iteration 5: The strategy selected at random is α 2 and the selected products within the strategy are j = 1 and 3.
Having the setup sequence of product j = 2 fixed and those with respect to products j = 1 and j = 3 unfixed does not enable any improving rescheduling, this results thus in the same objective function value of 994.25. This value is compared to the corresponding entry list. As both values are the same, the stopping criterion is met and the method ends with the final objective value of 994.25.

Computational Results
In this section, we provide a numerical analysis of the algorithm. We compare the algorithm to benchmark results for data classes that are generated with the same methodology as in [1], which is explained in more detail in the next section. The algorithm is programmed in JAVA and CPLEX 12.6 is used as the standard MIP-solver for solving the MIP formulations of the GLSP and iteratively the GLSP-RP with updated production sequences. CPLEX is is a high-performance general-purpose solver commonly used for solving linear and mixed-integer programming problems. All experiments are conducted on a computer equipped with an Intel Core i5-24510M CPU with 2.3 GHz and 6 GB RAM. We will present the parameters for the data classes in Section 5.1, and afterward, we will present the solutions of the algorithm and compare them to the CPLEX results in Section 5.2.

Data Set
All data classes consist of 50 instances. For generating the data classes we followed the same methodology as described in [1]. While the authors focused on a sensitivity analysis by varying specific parameters and test them for different numbers of micro-periods within the macro-periods we are going to test the performance of our algorithm for a single number of M t per class. We build three different data classes (A, B and C) where class A has the highest number of periodic decisions (i.e., 4 macro-periods and 7 micro-periods within a macro-period) and is therefore supposed to be the hardest class to solve. Table 9 summarizes the parameter settings for all data classes.

Algorithm Results
For testing our algorithm, we set a 100 s time limit on every iteration and a global time limit of 1800 s and compared the results to the best found solution by CPLEX 12.6 in a time limit of 1800 s. Concerning the algorithm parameters, we analyzed its performance for three relaxation strategies, i.e., α 1 , α 2 and α 3 (see Section 4.2), respectively, and tested the algorithm with regards to its main parameter l, i.e., for l = 1, 10, 20, and 50. The results for all instances of the three data classes are shown in Tables 10-12. Finally, at the end of the section a study on the implication of having a very large parameter l (i.e., 100) is conducted.
The tables show that while the solution times of the algorithm are lowest for the small list length of l = 1, the method tends to get stuck in a local optimum very fast for this list length. As a consequence, we can observe that only for a small number of instances better solutions could be found in comparison to the CPLEX results. In addition, we can see that the average objective function values are much higher than for the CPLEX results for all data classes. The effect of getting stuck in local optima is compensated for higher values of l resulting in much better objective function values but with also a small increase in solution times.
Furthermore, for the class A instances, CPLEX is not able to find an optimal solution for any instance in the given time limit of 1800 s. With l = 1, the algorithm was able to find better solutions for six instances with a reduction in the average solution time of 92.2%. While the average objective values for l = 1 for class A are higher than those from CPLEX, this changes for l = 10. For l = 10, the algorithm is able to find 37 better solutions and a reduction in the average objective values of 7.2%. Also, the number of iterations rise resulting in larger computation times; these times are still much lower (78.8%) than for CPLEX. The best results for all the data classes are for the list length of l = 50.
Analyzing the performance of the algorithm for each class of instances, we can observe that for the class A instances, the algorithm is able to find 48 better solutions in half of the computation time of CPLEX. Regarding the class B instances, 49 optimal solutions could be found by CPLEX in an average time of 212.23 s. For l = 50, the algorithm is able to find 48 of these solutions resulting in an average objective value that is only 0.1% higher than the CPLEX results while the computation time reduces by 20.8% compared to CPLEX. For the class C instances, it can be observed that 20 instances could be solved to optimality by CPLEX and the overall average CPLEX time is very high with 1354.50 s. While the algorithm is not able to find a better solution for any of the 50 instances for l = 1 this number increases to 6 instances for l = 10 and 12 instances for l = 20. For l = 50 our algorithm is able to find 40 better or same-quality solutions for the instances that could be solved by CPLEX. We can observe an average increase of 0.2% in the objective values due to the 10 instances where the algorithm was not able to find a better solution. The small increases in the objective function values are superseded by the significant reduction of the computational time, which is a decrease of 86.6% in this case in comparison to solving the optimization problem at once with a general solver.
Concluding this section, we can see several effects for the single parameter of LAHCM l. Reducing our approach to a hill-climbing matheuristic without late-acceptance criterion, i.e., setting l = 1 leads to solutions that get stuck in a local optimum soon for all data classes. Setting l = 10 is most effective for the solution quality of class A in terms of better solutions found. That is, for 37 instances better solutions can be found in comparison to the CPLEX results. However, this list size is not sufficient to overcome local optima for a large number of the instances for classes B and C. While for l = 20 the solution quality of the method increases, the best results can be seen for l = 50. As we can observe, a strict solution quality improvement for longer list sizes we can also observe that choosing a good value for l is critical regarding the computational time of the method.
In Table 13, we test LAHCM with l = 100 for the instances where we could not find better or optimal results for l = 50 compared to the CPLEX results. This will permit analyzing the impact on LAHCM when larger l are considered and thus better exemplify the trade-off between that parameter and solution quality. For classes A and B for all instances the method with l = 100 was able to find the same or better solutions compared to CPLEX in shorter computational times. For class B, we can observe the discussed trade-off, that is, as the solution quality improves for l = 100, the computational times are higher than for CPLEX. For class C, out of the 10 instances where the method could not provide better values for l = 50 (see Table 12), LAHCM was able to find eight equivalent results in less time than CPLEX for l = 100. In this regard, the average number of iterations regarding the instances presented in Table 13 for finding the best solution is 73.21 iterations, which also means that a list length of 100 should have been sufficient to find better or optimal results for all instances. Therefore, we can state that the selection of the strategy has a much higher influence than increasing the list length. On the one hand, strategy α 3 needs the longest computation time as it allows the most rescheduling decisions. On the other hand, this strategy is very effective to escape local optima while the other strategies are useful to create fast neighborhood solutions.      With: l = list length, Objective = Objective function value, Time = CPU time in seconds, Iterations = Number of iterations, # of better sol. found* = number of better results or same results for instances with optimality.

Conclusions
In this paper, we have proposed the Late Acceptance Hill-climbing Matheuristic (LAHCM) as a general framework for solving the optimization problem. It has been successfully applied to the general lot sizing and scheduling problem with defective items and lifetime constraints (GLSP-RP). Moreover, we have discussed the differences between the original GLSP and the GLSP-RP and gave an extensive overview of solution approaches and conducted research on several extensions of the GLSP. We exemplified the proceeding of the LAHCM on an example case and later applied the algorithm on several data classes. The overall results show that LAHCM performs very well on the GLSP-RP and is able to find optimal solutions or at least better solutions in much less time than a general purpose exact solver like CPLEX. While we were still not able to solve all instances to optimality, we were able to find 48 better or equal solutions for data classes A and B and for 40 instances for class C out of 50 instances. The results of the method also showed that increasing the single-point parameter l led to better solution values while the trade-off between better solution values and computational times has to be closely monitored. Moreover, the performance of LAHCM advises its use in practical environments where, due to rework and lifetime constraints, fast and reliable solutions are required in reasonable computational times.
Further research will be conducted towards studying the performance of LAHCM on other optimization problems, especially those where exact approaches are computationally limited because of the instance size. In that regard, a further research direction is to investigate how different solvers could be exploited using LAHCM for the same given problem. Moreover, steering mechanisms inside the LAHCM such as defining upper-level rules for determining the application of the fixing strategies will be investigated. In addition, research will be conducted in the area of different model formulations for the GLSP-RP, for example, in the form of a simple plant location model formulation.