THE INTEGRATED CUTTING AND PACKING HETEROGENEOUS PRECAST BEAMS MULTIPERIOD PRODUCTION PLANNING PROBLEM

We introduce a novel variant of cutting production planning problems named Integrated Cutting and Packing Heterogeneous Precast Beams Multiperiod Production Planning (ICP-HPBMPP). We propose an integer linear programming model for the ICP-HPBMPP, as well as a lower bound for its optimal objective function value, which is empirically shown to be closer to the optimal solution value than the bound obtained from the linear relaxation of the model. We also propose a genetic algorithm approach for the ICP-HPBMPP as an alternative solution method. We discuss computational experiments and propose a parameterization for the genetic algorithm using D-optimal experimental design. We observe good performance of the exact approach when solving small-sized instances, although there are difficulties in finding optimal solutions for medium and large-sized problems, or even in finding feasible solutions for large instances. On the other hand, the genetic algorithm is shown to typically find good-quality solutions for large-sized instances within short computing times. Mathematics Subject Classification. 90C27, 90B30, 90C59, 62P30. Received March 8, 2021. Accepted July 20, 2021.


Literature review
To the best of our knowledge the ICP-HPBMPP is not defined in the literature, even though the problem has similarities with one-dimensional cutting stock problems (1DCSP) and one-dimensional packing problems (1DPP). On the order hand, 1DCSP, 1DPP, and their variants have been substantially studied in the literature.
As far as one-dimensional cutting and packing problems (C&P) are concerned, the studies of [10,11] proposed a column generation algorithm to solve the linear relaxation of large instances of 1DCSP. Such studies served as basis for a number of subsequent works. In [18], a heuristic for the 1DCSP was proposed based on the solution of the linear relaxation supplemented by a one-pass branching up procedure. The authors validated the proposed heuristic approach, testing on benchmark instances and on a case study. In [8], the authors introduced a typology of C&P problems, unifying notions in the literature to guide further research on particular types of those problems. Two branch-and-price approaches were proposed in [21] to find optimal solutions for the 1DCSP. In [24], a new typology was presented to categorize the types of C&P problems in the literature between the years 1995 and 2004, introducing new categorization criteria. A model was proposed in [20] for the multiperiod one-dimensional cutting stock problems (M1DCSP), considering the use of objects/leftovers in stock. In [16], an integer linear model for the M1DCSP was proposed, a column generation procedure was implemented to solve the linear relaxation, and two rounding heuristics were developed for finding integer solutions to the problem. In [13], a mathematical model for the general integrated lot-sizing and cutting stock problem was proposed; additionally, a vast classification of the literature on that problem was performed, providing directions for future research.
Regarding the C&P problems and optimization approaches in precast production, De Castilho et al. [7] described the problem of minimizing production costs for slabs of precast prestressed concrete joists and introduced a genetic algorithm to solve it. An integer linear programming model for multiperiod production planning of precast concrete beams was proposed in [3], which can be seen as a special case of the HPBMPP. In [2], the authors introduced a mathematical model for the cutting stock/leftover problem and suggested a column generation technique for finding the problem's linear relaxation solution. In [22], a mathematical model was proposed based on the multiperiod cutting stock problem for the production planning problem of joists in trusses slabs industries. The authors suggested a solution method based on column generation to solve the linear relaxation of the problem. In [1], several integer linear programming models for the Heterogeneous Prestressed Precast Beams Multiperiod Production Planning Problem were proposed. The authros established the NP-hardness of the problem and suggested a constraint programming model for generating cutting patterns for the problem. The authors also carried out computational experiments to validate the performance of the integer linear programming models. In [23], the authors introduced a two-hierarchy simulation-genetic algorithm hybrid model for precast production to ensure the on-time delivery of precast components minimizing the production cost, while simultaneously optimizing the resource waste under uncertainty in the processing time of each operation. The model was validated by means of a case study.
The problem which we study in this work is the integration of the cutting stock/leftover problem proposed in [2] and the HPBMPP introduced in [1]. We explore its solution via exact methods and heuristics methods in the case where instances cannot be solved by the state-of-art solvers.

Problem statement
In this section, we formally define the ICP-HPBMPP and propose an integer linear programming model for its solution based on the models proposed in [2] for the Cutting Stock/Leftover Problem (CSLP) and in [1] for the Heterogeneous Prestressed Precast Beams Multiperiod Production Planning Problem (HPPBMPP).
The ICP-HPBMPP consists in finding a feasible production planning to cast certain quantities of prestressed precast concrete beams, possibly of different types, while minimizing the total length of pieces of bars that cannot be used as leftover. A leftover is understood here as a piece of bar that can be cut or overlapped in the future to meet new demands and is not considered waste. The beam factory has a fixed amount of bars and bar leftovers with standard lengths in stock that can be used within a given time horizon.
Each mold can only be used to cast one type of beam at a time. It is possible, however, to simultaneously cast beams of different lengths in the same mold, as long as they are of the same type. The total length of the beams produced during a given period in a given mold cannot be greater than the mold's capacity, and the total number of days required to complete the entire production cannot be greater than a given time horizon. After the process of cutting bars is finished, they are packed in molds in order to produce the beams. Note that different beam types can demand different numbers of bars. For this reason, we refer to this problem as a Cutting and Packing problem. The ICP-HPBMPP process can be seen in Figure 1.
As input of the problem we have a deterministic, static demand of beams, with their respective types and lengths, stock of bars and stock of bars leftovers, with their respective lengths. The cutting planning of bars is made for the entire time horizon, resulting in more bars leftovers (which can be used in another production planning) and, possibly, incurring in bar loss. The bars that are cut are then packed in the molds for the beam production, along the given time horizon. After the production of all demanded beams is met, there will usually be concrete waste of the beams and additional loss of bars.

Integer linear programming model
In order to define a model for the ICP-HPBMPP, we make use of the same parameters defined in [1], as follows: -M : number of molds in which the beams are produced; -T : number of available periods to complete the production; -C: number of beam types; -q c : number of distinct lengths of beams of type c, with c = 1, . . . , C; -l(c, 1), . . . , l(c, q c ): real numbers corresponding to the actual lengths of beams of type c, with c = 1, . . . , C; -d(c, k): demand for beams of type c and length l(c, k), with c = 1, . . . , C and k = 1, . . . , q c ; -t c : integer number corresponding to the curing time (in terms of periods) of beams of type c, for c = 1, . . . , C; -L m : real number corresponding to the capacity of mold m, with m = 1, . . . , M ; : packing pattern, where c i stands for the beam type associated with pattern P i and a i 1 , . . . , a i qc i represent the quantity of each beam of length l(c i , 1), . . . , l(c, q ci ) in patterns P i , with i = 1, . . . , r, c i = 1, . . . , C. Note that r represents the number of packing patterns; -P 0 : special pattern, which is used to denote that a mold is currently being used for the casting of a pattern that began in a previous period and whose production extends at least up to the current period.
Note that an idle mold (in other words, a mold that is not being used during a specific period) is not assigned the pattern P 0 . In fact, it has no pattern assigned to it.
In order to refer to specific information on a given pattern P i = c i , (a 1 , . . . , a qc i ) , we define the following notation: -N i (c, k): number of beams of type c and length l(c, k) that pattern P i includes. If c = c i , then N i (c, k) = a k , with k ∈ {1, . . . , q ci }; otherwise, N i (c, k) = 0, for any k.
-E i : number of periods required to produce the beams in P i , with i = 1, . . . , r. This number equals the quantity of consecutive periods in which P i remains occupying a mold and is precisely the curing time of beams of type c i , given by t ci .
Given a set of patterns P = {P 1 , . . . , P r }, not including P 0 , we define some important sets as follows: -Q(m): set containing the indices of the patterns in P whose capacity does not exceed the capacity of the In what follows, we present the parameters that concern bars and bars leftover: -W : number of different bar lengths; -V : number of different bar leftover lengths; -H: number of cutting patterns; -O: number of overlapping patterns; -Γ: number of different mold lengths; -b 1 , . . . , b W : bar lengths; -b W +1 , . . . , b W +V : bar leftover lengths allowed. Note that this data narrows the types of cutting, and overlapping patterns; -L 1 , . . . , L Γ : mold lengths. Note that this data narrows the types of cutting, and overlapping patterns; -G(L γ ) = set of molds which are of length L γ , γ = 1, . . . , Γ; -H w : set of cutting patterns for bar of length b w that do not include leftovers.
-H w (v): set of cutting patterns for bar type w that include leftovers of length b W +v ; -O: set of overlapping patterns; -O(γ): set of overlapping patterns that produce bars of length L γ .
. . , a µ V )): overlapping pattern that generates a bar of length L γµ , with γ µ = 1, . . . , Γ and µ = 1, . . . , O. Note that a µ 1 , . . . , a µ V are the number of bars of lengths b W +1 , . . . , b W +V ; -D ci = number of bars that a pattern P i with beam type c i demands; -e w = number of bars of length b w in stock, leftover or otherwise, with w = 1, . . . , W + V ; -a v,µ = number of leftovers of length b W +v in overlapping pattern O µ , with µ = 1, . . . , O.
-a γ,h,w = number of objects of length L γ cut from a bar of length b w following a cutting pattern I h that generates no leftover, with with w = 1, . . . , W + V ; -a γ,h,w,v = number of objects of length L γ cut from a bar of length b w following a cutting pattern I h that generates a leftover of length b W +v , with w = 1, . . . , W and v = 1, . . . , V . -f h,w = waste resulting from using a cutting pattern I h to cut a bar of length b w generating no leftover, with w = 1, . . . , W + V . -f h,w,v = waste resulting from using a cutting pattern I h to cut a bar of length b w generating a leftover of length b W +v , with w = 1, . . . , W and v = 1, . . . , V . -f µ = waste of bar produced by overlapping pattern O µ , with µ = 1, . . . , O.
We present the decision variables below: 1, if the packing pattern P i starts to be used in mold m at period t (and its usage, naturally, lasts for E i periods); 0, otherwise. z t = 1, if as least one mold is used at period t, for t = 1, . . . , T ; 0, otherwise.
y h,w : number of bars of length b w cut following a cutting pattern I h ∈ H w . y h,w,v : number of bars of length w cut following a cutting pattern I h ∈ H w (v) generating a leftover of length b W +v . o µ : number of times the overlapping pattern O µ was used, µ ∈ O.
Note that variables y h,w , y h,w,v , and o µ are nonnegative integer decision variables. We present the integer linear programming model proposed for the ICP-HPBMPP as follows: The objective function (3.1) is divided into 4 terms. The first term is the makespan value. The second term defines the waste related to the use of new bars to produce the demand of bars. The third term describes the waste associated to the use of new bars to produce the bars required by beam production while creating new leftovers. Finally, the fourth term specifies the waste corresponding to the bar leftovers in stock that are used to produce the amount of bars required. Note that each term of (3.1) could alternatively be regarded as an independent objective function to be minimized. We obtain (3.1) using the weighted sum method, in which the parameters λ i ∈ R + , with i = 1, . . . , 4, indicate the weight of each objective function term. A solution that minimizes (3.1) is, therefore, a Pareto optimum [12].
Constraints (3.2) ensure that at most one pattern must be assigned to mold m at period t, with the possibility of this pattern being P 0 . Constraint set (3.3) requires that all demands must be satisfied. Constraints (3.4) force that, if pattern P i is initiated at period t, then the next E i − 1 periods shall have the pattern P 0 assigned to them (the right-hand side of the constraint remains unconstrained, in case x m,t i = 0). Constraint sets (3.5) and (3.6) establish that P 0 shall only be used in mold m if there is some pattern associated with a previous period in the same mold, whose production has not yet been completed.
Each constraint in set (3.7) ensures that variable z t must be 1 if period t is used to produce beams. Constraints (3.8) force that there is no inactive period during beam production in the molds. This means that the production is continuous, i.e., if a mold is used it will be used with no interruption; in other words, if the production stops at a given mold and period, it will not resume in that mold at a subsequent period.
Constraints (3.9) establish that the number of bar leftovers cut plus the number of leftover bars used to produced bars via overlapping does not exceed the stock, note that the cutting of a leftover does not generate leftovers. Constraint set (3.10) ensures that the number of bars cut does not exceed the stock. Constraints (3.11) force that the amount of bars necessary to produce the beams is achieved, assuming that the required amount of bars is the number of bars used by the forms in the entire time horizon. Constraints (3.12)-(3.16) define the domains of the decision variables.
The model (ICP) has O(M T r + W V H + O) variables and O(q + M T r + V + W + Γ) constraints, with q = C i=1 q c . Thus, depending on the total number of possible packing, cutting, and overlapping patterns, there may be an excessive number of variables and constraints in the model. We choose to limit the number of packing patterns, which are the more numerous type of pattern, in practice, by using only maximal packing patterns, used successfully by [1,21]. We say that a pattern P i contains a pattern P j if c i = c j and a i k ≥ a j k , with k = 1, . . . , q ci . Proposition 3.1. Restricting the model (ICP) to using only maximal packing patterns does not modify its set of optimal solutions.
Proof. The proof can be found in [1].

NP-hardness
To argue the ICP-HPBMPP hardness note that for instances where D c = 0, for all c = 1, . . . , C, constraints (3.9)-(3.11) are naturally fulfilled and all variables y h,w , y h,w,v and o µ are set to zero, reducing an instance of ICP-HPBMPP to an HPPMBPP instance involving the minimization of the makespan, up to a constant multiplicative factor. Consequently, the ICP-HPBMPP is a generalization of HPPMBPP, which is already known to be NP-hard [1].

Objective function lower bound
Since the ICP-HPBMPP is a NP-hard problem, a lower bound for the optimal objective function value may help in evaluating the quality of feasible solutions in heuristic and exact methods. In order to simplify the presentation of our proposed lower bound for objective function (3.1) optimal value, we present the following notation. For a given γ ∈ {1, . . . , Γ} we define the following sets: An upper bound on the optimal value of model (ICP) is given by equation (3.17). l(c, k) · d(c, k) L γ · min{Ĉ γ } . (3.17) The first part of equation (3.17) corresponds to a lower bound for the makespan, while the second part stands for the minimum waste resulting from using molds of some fixed length L γ .

Patterns generation
Instead of carrying out exhaustive enumerations, we generated the desirable packing, cutting, and overlapping patterns for a given instance using constraint programming models, which are described in the remainder of this section.

Packing patterns generation
Consider the following notation, in addition to the notation presented in Section 3: -K: the largest number of different lengths among beam types, i.e., max q c with c = 1, . . . , C. For example, in an instance with 2 beam types, in which type 1 has 6 distinct beam lengths and type 2 has 4 distinct beam lengths, we have K = 6. -v i ∈ {1, . . . , C}: a decision variable that corresponds to the type of beam used by the pattern P i . -γ i ∈ {1, . . . , Γ}: auxiliary decision variable for generating patterns that will be maximal in at least one mold of the problem. It defines in which mold capacity the generated pattern P i is maximal.
-A i ∈ Z K : a vector of decision variables, with A j representing the number of beams of the length (v, j), for all j ∈ {1, . . . , K}. Given a pattern P i of type v, the nonzero components of vector A i correspond to : the generated pattern. For the generation of a packing pattern P i we propose the model, which is adapted from [1].
Constraint (4.1) implies that the pattern type has domain ∈ {1, . . . , C}. Constraint (4.2) defines the length of the molds in which the generated pattern should be maximal. Constraint set (4.3) implies that if the generated pattern is of type v then it includes no beam of size l(v, j), such that j > q v . Constraint set (4.4) imposes that the capacity used by the generated pattern is simultaneously larger than the mold length minus the shortest beam length from its type and no larger than the length of the actual mold. The empty pattern is, therefore, not generated and has to be manually included in the final set of patterns. We utilized the solver CPLEX CP Optimizer to enumerate all the solutions of model (4.1)-(4.5).

Cutting patterns generation
In this section, we propose a constraint programming model for cutting patterns generation. The decision variables are given below: -w h : index of the bar that will be cut in the generated cutting pattern I h ; -A h i : number of items of length L i cut in the pattern, for i ∈ {1, . . . , Γ}; -A h i : number of items of length b W +i cut in the pattern, for i ∈ {Γ + 1, . . . , Γ + V }; : the generated pattern.
The proposed constraint model for generating a cutting pattern H h is given by equations (4.6)-(4.10).
Constraint (4.6) defines the domain of each decision variables w h . Each w h variable determines defines the bar that will be cut in the current pattern to generate items. If 1 ≤ w ≤ W , the bar that will be cut is a new bar. If W + 1 ≤ w ≤ W + V , the bar that will be cut is a bar leftover. Constraint (4.7) states that the total length of items cut in the pattern must be shorter than the length of the bar used to cut such pattern, with expression element(w h , b) standing for the w h -th element of array b [4]. Constraint set (4.8) implies that a cutting pattern only generates one type of leftover. Constraint (4.9) implies that a leftover does not generate more leftovers. We utilized the CPLEX CP Optimizer to enumerate all the solutions of model (4.6)-(4.10).

Overlapping patterns
In order to enrich the problem by allowing the possibility of using overlapping bars, we recall that an overlapping pattern O µ is a tuple O µ = (γ µ , (a µ 1 , . . . , a µ V )). Note that γ is associated to the length of the bar that is generated in such pattern. Such length must be equal to the capacity of some mold, since we are only required to produce bars via overlapping that are used for beam production. A bar produced by overlapping is only produced from leftovers in stock.
In order to simplify the model's notation, consider the following decision variables: -A µ i : decision variable that represents the number of items b W +i used in the overlapping pattern, for i ∈ {1, . . . , V }.
-f ≥ 0: decision variable that expresses the waste of bar associated to the overlapping pattern to produce a bar of length L γ .
: the generated pattern. The following constraint programming model can be used to produce an overlapping pattern: Constraint (4.11) ensures that the length of the bar produced is one of the possible mold lengths. Constraint (4.12) forces that the total length of the chosen leftovers is greater than the length of the bar produced via overlapping plus a constant which is the loss of the bar resulting from the overlapping process. Constraint (4.13) defines that only 2 leftovers are used in the production of the bar made via overlapping. Constraint (4.14) defines the bar waste resulting from the overlapping pattern. The constraint programming model for overlapping pattern generation is sufficiently flexible to accommodate the production planner's necessities. In a more general setting, we could require that a bar made via overlapping can only be produced by using more than 2 and no more than a predefined number of leftovers and specify the value to be proportional to the number of leftovers used in such pattern.

Genetic algorithm for the ICP-HPBMPP
In this section, we propose a genetic algorithm to solve the ICP-HPBMPP, formalize the solution representation chosen, the solution fixing procedure, the selection, mutation, and crossover operators, as well as the initial population generation, population restart, and local search.

Solution representation
The solution representation consists of a 2-row matrix, in which each column j consists of the genes a j and x j , where a j is a pattern index and x j is the number of times the pattern represented by a j is used.   Table 1. Description of instance cwp000.
Instance cwp000 In order to illustrate the solution representation we first present, in Table 1, the cwp000 instance, which was randomly generated. Its respective packing, cutting, and overlapping patterns are presented in Tables 2, 3, and 4, respectively.
Note that ID is associated with the pattern indices. An optimal solution for the cwp000 instance is shown as the chromosome in Figure 3.
In the solution shown in Figure 3, we obtain an objective function value of 2.1, with makespan of 2 periods and bar waste of 0.1 m. Figure 4 shows that packing patterns with indices 2 and 6 were used 4 and 2 times, respectively. Due to the fact that we are restricted to using only maximal packing patterns in their respective molds and a given packing pattern is maximal with respect to only one particular length of mold, we infer that packing pattern 2 is associated with molds of length 5.95 m, and packing pattern 6 is associated with molds of length 11.95 m. Therefore, we need to produce a total number of 2 bars of length 5.95 m and 6 bars of length 11.95, since the beam type produced by each solution packing pattern requires only one bar. The cutting patterns used are those with indices 11, 15 and 16, and their frequencies are 2, 1, and 2, respectively. None of the overlapping patterns was selected in the solution.    The production planning consists of the specification of the exact quantity of bars required for the beam production, as long as the available stock of bars is not violated. Thus, the solution encoded in the chromosome in Figure 3 is feasible.

Initial population generation
Since we typically need a large quantity of individuals to generate a population, deterministic methods are not the best choice, despite the high-quality solutions produced by them. We propose a pseudorandom approach to generate a large quantity of solutions, which is described in Algorithm 1.
We call this method pseudorandom because we choose the patterns to add to the solution randomly, although each pattern frequency in the solution is computed in such a way as to respect stock and satisfy the demand. pacp ← random packing pattern that has not yet been selected. 4 if There is some beam in pacp whose demand is unfulfilled then 5 Increment the number of times that pacp is used in solution until all beams in pacp have their demands fulfilled. cutp ← random cutting pattern that generates bars of length Lγ that has yet not been selected. 12 bars needed ← number of bars of length Lγ required. 13 n ← number of times cutp can be added to solution without violating bars stock. ovep ← random overlapping pattern that generates a bar of length Lγ that has not yet been selected. 18 bars needed ← number of bars of length Lγ required. 19 n ← number of times ovep can be added to solution without violating bars stock. 20 Increment ovep frequency in solution by max(bars needed, n) times.

Fitness function and selection operator
We use the objective function (3.1) from the mathematical model (ICP) as the fitness function to evaluate the solution quality of a given chromosome. The selection operator consists of the process of selecting the best distinct solutions with respect to their respective fitness function value, i.e., the individuals with the lowest fitness values.

Crossover operators
In this subsection, we propose two alternatives to use as crossover operators: types 1 and 2. Given two parents, both crossover types generate one offspring, which consists of a new solution (chromosome).
In crossover type 1, we preserve all pattern indices from both parents, but the number of times each pattern is used in the offspring corresponds to the mean of the number of times they are used by the parents rounded to the largest integer. For each gene there is a probability of mutation. When the mutation occurs the number of times that the current pattern is used in such gene is set to zero. After this crossover process, if the generated offspring results in an infeasible solution, an iterative procedure, shown in Algorithm 6, is applied for its correction. If some pattern from the current offspring is used zero times, the gene associated to it is removed from the chromosome.
In crossover type 2, we first initiate the offspring using all patterns that used in both parents with their respective frequencies set to zero. For the genes that have patterns that are part of both parents simultaneously, their respective frequencies are set as the mean of their frequencies in the parents rounded to the largest integer. For each remaining gene we have a probability of 50% of setting its respective frequency to be equal to the originating parent frequency or keeping it equal to zero. If the resulting offspring is not feasible, the fixing procedure, shown in Algorithm 6, is applied to it and all patterns with final frequencies equal to zero have their respective genes removed from the chromosome.

Mutation operator
The mutation of an individual consists of (a) choosing one pattern p 1 that is in the solution, and (b) adding one pattern p 2 , chosen randomly, that is not part of the solution. The number of times that p 2 is used becomes the number of times that p 1 is used, and the number of times that p 1 is used is set to zero. If the solution is infeasible after this procedure we apply a fixing phase to the solution. This process is frequently required in practice and is described in the next subsection.

Infeasible solution fixing
Since the proposed genetic operators of crossover and mutation can compromise the feasibility of solutions, we must define a procedure to fix infeasible solutions, turning them into feasible ones.
A chromosome may be an infeasible solution due to different reasons, as follows: (1) Infeasibility type 1, due to beam demand: the frequencies of packing patterns in the solution are not enough to fulfill the beam demands; (2) Infeasibility type 2, due to bar stock: the number of bars which are used in cutting and overlapping patterns exceeds the bar stock; (3) Infeasibility type 3, due to inconsistent number of bars produced and required: the number of necessary bars generated by cutting and overlapping patterns is different from the number of bars that beam production requires.
If we detect any of those kinds of infeasibility, we must apply the infeasible solution fixing phase, which consists of Algorithm 6. Each infeasibility type is treated in a particular procedure: Algorithms 3, 4, and 5 are used to fix infeasibility types 1, 2, and 3, respectively.
The unnecessary packing patterns procedure, shown in Algorithm 2, works like a solution treatment phase, which is not a necessary part of the solution fixing process, although applying such procedure may improve solution quality and simplify the fixing process, i.e., it would be less likely that the modified solutions could not be fixed. The procedure consists of decreasing the frequency of packing patterns after the beam demands are already fulfilled if there are beam surplus.
In Figure 5, we show an example of the crossover operators, with offspring 1 as the solution generated by crossover operator type 1, and offspring 2 as the solution created by crossover operator type 2. Note that the fixing procedure was applied to offspring 2 and not to offspring 1. In Figure 6, we show an example of the proposed mutation operator. The resulting chromosome is infeasible; therefore, the solution fixing procedure must be applied. If the application of the solution fixing procedure to a given chromosome could not turn it into a feasible solution, the chromosome is discarded.

Population restart
The population restart consists of the creation of a new population to compose the next generation after a predefined number of epochs. We apply a population restart after a given number of generations with no improvement of the best-fitness value. We divide such procedure into three parts, as follows: (1) selecting a certain number of the best-fitness individuals from the current population; (2) generating a number new pseudo-random individuals; (3) creating a new population with individuals from steps 1 and 2 and applying the selection operator to form the next population.

Local search
In order to improve the quality of final solutions, we apply a local search to every individual of the final population. For the local search we use the insert movement, which consists of, given two genes indices i and k, with i < k, inserting the gene i one position in front of k-th gene, i.e., all the genes between positions i and Algorithm 5: Fix chromosome with respect to infeasibility type 3.

Data: Infeasible chromosome
Result: Potentially feasible chromosome 1 Calculate the #bars generated by cutting and overlapping patterns; 2 Calculate the #bars that beam production requires according to the frequency of packing patterns; 3 for each bar γ generated do 4 if #bars γ generated ¿ #bars γ that beam production requires then 5 for each cutting pattern I h that generates only bars γ do 6 rt ← #bars γ generated − #bars γ that beam production requires; if #bars γ generated > #bars γ that beam production requires then 13 for each overlapping pattern Oµ that generates a bar γ do 14 rt ← #bars γ generated − #bars γ that beam production requires; if #bars γ generated < #bars γ that beam production requires then 28 for each overlapping pattern Oµ that generates a bar γ do 29 Increment frequency(Oµ) until (#bars γ generated ≥ #bars γ that beam production requires) or the stock is violated with new increment;

30
Update the #bars γ generated; 31 end 32 end 33 end 34 return Chromosome k + 1 are moved one position to the right after the insertion of the k-th gene. In Figure 7 an insert movement neighbor is shown for a given solution after inserting 2nd gene in front of 5th gene.
Considering the function INSERT(solution, i, k) as the movement of insertion given indices i and k, we describe the local search procedure in the Algorithm 7.

Algorithm description
In order to describe the proposed genetic algorithm we define the following parameters: population size (TP), number of generations (NG), crossover type (CRS), number of pseudo-random solutions generated for the initial population and restart selections (AS), mutation probability (MUT), number of generations with  no fitness improvement to apply population restart (RST), and the number of individuals from the current population selected to be used in restart operator procedure (TER).

Computational experiments
In this section, we present computational experiments on a set of benchmark instances that were produced with the intent to mimic real-world scenarios, in order to evaluate the solution methods proposed in this study. The instances used in this work were randomly generated based on an existing order arising from a real-world production plant [1,3]. For privacy reasons, we are not allowed to use or provide here the actual data coming from the aforementioned instance.
The patterns corresponding to each test instance were generated using the constraint programming solver IBM ILOG CPLEX 12.8 CP Optimizer. For the integer programming model implementation we adopted the solver IBM ILOG CPLEX 12.8. Both solvers were used with Concert technology using the C++ programming language. The genetic algorithms were also developed with the C++ programming language.
We carried out every test in this paper on a Linux Ubuntu 18.04 64 bits machine with 8 GB of memory and Intel Core i5-3470 CPU 3.20 GHz ×4 processor. We compiled the created codes with the GNU GCC 7.3.0 compiler using Code::Blocks 17.12 IDE. Note that, for different values of λ i we can form the Pareto front and may have different behaviors of the proposed model and algorithms. However, for the purpose of the study, we did not approach the multi-objective nature of the problem and considered, for each test described in this section, λ i = 1, with i = 1, . . . , 4.

Test instances generation
In this subsection, we describe how we generate the set of benchmark instances used in this section. We introduce a set of instances that are based on data arising from a possible real-world scenario. The different instances represent a sample of the variability of the problem's parameters, such as number of beam types, number of molds, and mold lengths.
In Table 5, we present details about each test instance parameter. We can see that the number of packing patterns increases as the number of beam types increases. However, the number of cutting and overlapping patterns remains constant because of the fact that we expect that the possible distinct bar lengths are standardized in real-world scenarios and therefore do not lead to variability.
We consider mold capacities of 5.95 m and 11.95 m, while we take 1.12 m, 1.45 m, 2.35 m, 2.5 m, 2.65 m, 2.95 m, and 3.3 m as possible beam lengths. For each instance, the possible curing times may be 1, 2, or 3 periods, chosen randomly when instances have more than 3 beam types. In addition, if the instance has up to 3 beam types, we associate the curing time to the beam type index, for example the beam type 2 needs a curing time of 2 periods. With respect to the number of bars that some beam type demands, we choose randomly a value between 1 and 3 for each beam type. We choose the beam demands uniformly between 17 and 50. For total time horizon T , we calculate it as the ceiling of 150% of the optimal makespan lower bound, defined in equation (6.1).
For all instances, we consider an unique length of new bars as 12 m and the possible lengths of bar leftovers as 2 m, 5 m, 6 m, and 8 m. We do not vary such lengths along the test instances since, in practice, it is expected that they are standardized. For generating realistic bar stocks we introduce an upper bound for the number of bars needed to fulfill the beam demand as UB (see Eq. (6.2)).
We set the stock of new bars of length 12 m equal to UB, whilst we choose the stock of each leftover randomly between UB/5 and UB following an uniform distribution. We implemented the instance generator using the MATLAB programming language.

Computational experiments with the mathematical model
In this subsection, we discuss the results of the computational tests with the benchmark instance set that we generated following the scheme described in Section 6.1. In Table 6, we show the results of the computational experiments for the model (ICP) and its linear relaxation. The solution time was limited to 3600 s. As regards to notation in Table 6, we consider LB, IP, and LP standing for the optimal objective function value lower bound, best solution value by CPLEX for model (ICP), and its linear relaxation value, respectively. When we say gap we mean the relative percentage deviation between the best integer objective and the objective of the best node remaining in the CPLEX B&C tree, calculated like this: gap = 100 · |bestbound − bestinteger|/(1e − 10 + |bestinteger|) (0% means a proven optimal solution). We denote by "B&C nodes" the number of nodes generated in the branch-and-cut tree in the solution process, and t(s) as the solution time in seconds. Table 6 shows that the linear relaxation of all instances could be solved, with the average time of 53.21 s, and with 624.61 s being the longest time to get to the optimal solution. On the other hand, only 11 instances could be solved to optimality by the integer programming model (4 of them solved in the root node of the B&C tree). For 23 instances CPLEX could not even find a single feasible solution, a situation that we denote by "-". Moreover, for 36 instances feasible solutions were found by CPLEX, but the search was not completed within the time limit. The computational test results show that the larger the instance parameter values are, the harder it is to find solutions for the resulting problem. With increased values of the instance parameters, when solutions are found, the optimality gap tends to be worse, i.e., the solutions achieved within the time limit are further from the optimal solution, in terms of objective function value.
We compare the results of the integer linear model (ICP), its linear relaxation, and our lower bound, in equation (3.17), for the optimal value of objective function in the chart in Figure 9.
In Figure 9, the lower bound proposed in this work for the optimal objective function value was greater than the linear relaxation for all test instances and highly close to the objective function values obtained by CPLEX.

Experimental design and computational experiments with the proposed genetic algorithm
In order to achieve a better parameterization for the robustness of the proposed genetic algorithm, we apply fractional factorial parameter design. The authors of [9] used Taguchi experimental design (see [15]) to achieve improved robustness of the genetic algorithm which they proposed. In this method the optimal parameter choice is found with the analysis of different level combinations of the control factors in an orthogonal array, with no necessity of testing all of the possible level combinations. Table 7 displays the proposed levels for the genetic algorithm parameters (control factors) introduced in Section 5.   We must have one degree of freedom for total mean, one degree of freedom for each factor with two levels, and two degrees of freedom for the factor with 3 levels amounting to a total of nine degrees of freedom (1 + 1 × 6 + 2 × 1 = 9). However, with the control factors and respective levels that we defined, there is no orthogonal array aside from the full factorial array. Thus, we are not able to use a classical Taguchi orthogonal array design. In such circumstances one alternative is to use a D-optimal design approach [6], which aims at minimizing the generalized variance of the estimated regression coefficients. Note that D-optimality is only one possible criterion to choose a particular design. We obtain the D-optimal design, by Fedorov algorithm [19] using R programming language for 9 trials for the chosen factors and their respective levels, illustrated in Table 8.
Furthermore, the effectiveness characteristic of the genetic algorithms proposed is the expected fitness value, which we seek to minimize, i.e., "the lower is better principle". Thus, for increased robustness of the algorithm we use the S/N (signal-to-noise) ratio, defined as follows. Note that the larger the value of S/N ratio the better.
Trial TP NG MUT RST AS CRS TER with i and j denoting index of trial and index of replication, while FIT stands for the best objective function value obtained by running the GA. We denote by "trial" a certain combination of the control factor levels.
We define a replication as one GA run of some trial for a given instance, and N as the number of test instances multiplied by the number of replications. Since we have an instance set of size 70 and we run each instance 10 times, we perform 700 replications for each trial.
Since CPLEX could not find optimal or even a feasible solution for most instances, we are unable to use the relative percentage deviations from the optimal solution as a performance indicator for the GA. Thus, we use the lower bound relative percentage deviations (LBD) of the fitness values for such purpose. Given a trial i and a replication j the LBD value is defined as follows: where LB j stands for the lower bound of the optimal objective function value for the test instance used in replication j. The LBD for a given trial i, denoted by LBD i , is the average LBD for all replications of instance set, as we can see in the following equation: The remainder of the experimental design procedure consists of three phases: (1) Evaluate the impacts of the control factors on the S/N ratios and on the LBD values.
(2) For each factor, which has significant impact on the S/N ratios values, we choose the level which increases the S/N ratios. (3) For the factors, which do not have significant impact on the S/N ratios and have significant impact the LBD values, we choose the level which better approximate the lower bound values. (4) For the factors, which have significant impact neither on the S/N ratios nor on the LBD values, we select the factor levels regarding the more economic manner, that is, we choose the levels which have less impact on the algorithm running time. Table 9 shows the results after carrying out the computational tests for each trial with the test instance set. The value of N is 700, which corresponds to 10 replications to each one of the 70 test instances.
In Figure 10, we can see the main effects plot for the control factors using S/N ratios as the response variable. In Figure 11, we show the boxplots for each factor also using S/N ratios as the response variable. The mean response is clearly influenced by the type of crossover, while it is not so clear to affirm whether or not the other factors influence the response variable (Tab. 10).     0 "***" "**" 0.01 "*" 0.05 "." 0.1 Adjusting the linear regression model for all seven factors and performing an ANOVA test, we observe that only CRS is statistically significant with P -value 0.0259. Then we remove, one by one, the factors whose Pvalue is the greatest and readjust the regression model until all factors are statistically significant obtaining the ANOVA results in Table 11.
The number of generations, mutation rate, restart, and type of crossover showed to be statistically significant, meaning that we chose the levels whose average S/N ratios are the greater. The parameter levels chosen as a result of the ANOVA test are 1000r generations, 0.05 of mutation rate, 0.2r generations with no improvement to apply restart, and crossover type 1.
As regards to the LBD as response variable to the linear regression model. We observe in the main effects plot in Figure 12 and in boxplots in Figure 13 that LBD have a similar behavior on the control factors. However, we note that, in this case, the lower the LBD value the better.
Adjusting the linear regression model for all the seven factors and performing an ANOVA test using the LDB as response variable, we conclude that only CRS is statistically significant with P -value 0.03219. Therefore, we remove from the regression model the variables, one by one, whose P -value is the greatest and readjust the model until all factors are statistically significant achieving the ANOVA results illustrated in Table 11.
Taking into consideration the LBD as response variable to the regression model, only the mutation rate, and type of crossover are statistically significant, meaning that we would choose the mutation rate 0.05, and crossover type 1. However, these variables were already fixed at the S/N ratios analysis, and no factors that were not statistically significant for the S/N ratios showed to be statistically significant with LBD values. This leads us to choose the levels that would spend less computational time, for the factors whose level was not selected yet. Therefore, the most robust parameterization of the levels for the proposed control factors is: population size 25, 1000r generations, 200r generations with no improvement to apply restart, 100r pseudo-random solutions  generated in the initial population and restarts, crossover type 1, 5 preserved individuals upon restart, and mutation rate of 0.05 (Tabs. 12 and 13).

Relax-and-fix heuristics
Aiming to compare our proposed GA with other methods developed to closely related cutting and packing problems, we adapt the relax-and-fix (RF) matheuristic proposed by Signorini et al. [17] to the problem of one-dimensional multi-period cutting stock problem for the production of precast beams. In the RF approach, a group of integer decision variables in a mixed-integer linear programming model is partitioned into mutually exclusive sets in which the number of partitions determines the number of iterations of the heuristic [25].
In each iteration, only the variables of a given partition are defined as integer decision variables, and the other decision variables are relaxed. The resulting sub-model is then solved, and two situations can occur: infeasibility and feasibility. In the case of infeasibility, the execution is stopped. In the case of feasibility, the decision variables in the previous partition are fixed in their current values, and the process is repeated for the other partitions.
We extend three RF strategies, adapted from Signorini et al. [17], to solve the ICP-HPBMPP. We decompose variables x imt and z t together according to dimension of periods, since they are the most numerous set of

Experimental results and discussion
In order to illustrate the behavior of the genetic algorithm, we run the GA with instance cwp021. Figure 14 illustrates the best fitness and mean fitness of the populations along all generations. The x-axis of the Figure 14 is on logarithmic scale. The largest improvement takes place during the first generations of the GA, while in the last ones the best fitness is stagnant with some improvement upon the first restart.
In Figure 15, the best fitness values obtained by running the GA were better than CPLEX in five instances, while solutions were obtained for all instances which CPLEX could not solve.
In Figures 16 and 17, the time spent by the GA on solving each instance was significantly better than the CPLEX solution time on the large and medium-sized instances. Thus, CPLEX was faster than the GA in the small-sized instances. The y-axis Figure 17 is in logarithmic scale.
In Table 14, "nfsol", "nbsol", and "ave lbd" stand for number of feasible solutions, best solutions, and average lower bound deviation, respectively. Based on Table 15, we can observe that GA was the only method able to find feasible solutions for all the evaluated test instances. Aiming to provide a fair statistical comparison between the evaluated methods, we consider a lower bound deviation of 100% for a method that has not returned a feasible integer solution for a given test instance (Figs. 18 and 19).      In order to achieve a pairwise comparison between the evaluated methods (RF1, RF2, RF3, MILP, and GA), we performed an ANOVA procedure, followed by a Tukey's test [14]. Lower bound deviations of the five methods are statistically different with p-value of 0.001578 in the ANOVA test. Figure 20 illustrates the Tukey multiple comparisons of means with 95% family-wise confidence level. Based on the results obtained, we can observe that GA outperforms all the other methods under comparison. Furthermore, we different methods for each instance group using Borda counting [5], as illustrated on Table 16. We note that we use Borda's original counting for handling ties. Taking this indicator into account, we can observe that the proposed GA obtained the largest count of votes of all methods under comparison, pointing to the superiority of this method.

Final remarks
In this work, we proposed a novel variant of cutting sequencing problems called the integrated cutting and packing heterogeneous precast beams multiperiod production planning (ICP-HPBMPP), which, to the best of the authors' knowledge, has not yet been studied and may have a large impact on both real-world and theoretical studies. The ICP-HPBMPP consists in integrating the problem of production planning of precast beams with the problem of cutting the traction elements used in such production, while taking into consideration the generation of leftovers and bars generated via overlapping.
We argued that the problem is NP-hard and proposed an integer linear programming model for its solution, in addition to a lower bound on its optimal objective function value. We also showed that restricting the  formulation to using exclusively maximal packing patterns does not change the optimal solution set of the problem.
We also proposed three constraint programming models for generating distinct types of beam production patterns. Additionally, we introduced a set of benchmark instances and carried out computational experiments in order to evaluate the relative performance of the different solution methods studied.
The experiments showed that the integer programming model can be used to solve small size instances, while it typically does not reach optimality while solving medium-sized ones. In addition, the model usually does not find feasible solutions for large size instances. We introduced a genetic algorithm for solving the problem and fine tuned its parameters by means of a D-optimal experimental design to achieve improved robustness of the algorithm. The final genetic algorithm is an attractive alternative to the integer programming model, resulting in high-quality solutions in shorter solution times as compared with the exact model.
There are numerous opportunities for future work regarding the ICP-HPBMPP. In the domain of modeling, the problem can be modified to take into consideration distinct types of bars varying in matter of diameter or material, instead of only in matter of length. Also, dynamic demand could be considered, i.e., in each period a new demand of beams could be included, while not exceeding a prescribed stock of bars. Regarding solution approaches, multi-objective optimization algorithms can naturally be applied to the problem, since it involves preferences between makespan and bar waste. Decomposition approaches, such as column generation, or MIP heuristics, e.g., size-reduction heuristics, can also be interesting methods to be explored in conjunction with the proposed integer programming model.