REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

Keywords: FOD, electric arc, aircrafts, busbars

Current development of more electrical aircraft and future all electrical aircraft leads increased power densities and power on board. The increase in power densities take the form of increased voltage (230/400 Vac, ± 270 Vdc) and/or decreased distances between conductors. The addition of these evolutions leads to higher probability of electric arc, which might be potentially more destructive. In this paper electric arcs ignited by a Foreign Object Debris (FOD) that shorts out two busbars, are characterized. The electrical parameters meet the aeronautical standards, i.e. 230/400 Vac at 400 Hz. A statistical study of the impact of the FOD geometry on the arc size is carried out with the aim to give guidelines for the development of passive protections. The results show that there is no direct correlation between the arc dimensions and the FOD diameter or mass, although the object behavior was different. In a second part the electric arc was characterized in terms of dimensions over time and electrical signature to state on the time in which active protection should trip and to try to find a signature of the arc to detect. The results tend to prove that the dV/dt might be a relevant marker of the arc.

I. INTRODUCTION

The evolution toward more electrical aircrafts has major impact on today and tomorrow aircrafts architecture [START_REF] Roboam | New trends and challenges of electrical networks embedded in 'more electrical aircraft[END_REF]. Besides the replacement of parts of pneumatic and hydraulic energy vectors by the electrical one and the increase in electrical power available, this evolution is characterized by an increase of the voltages levels (230/400 Vac and ± 270 Vdc) and a densification of the systems with the aim of increasing power densities. While the goals and benefits of this strategy are well identified (lower the mass of the aircraft and the fuel burned), it comes with many challenges to tackle. One of these, directly related to the changes in the voltage amplitude and the reduction of the distances between conductors at different potentials, is the increased risk of electric arc ignition, potentially destructive (flight SR 111, R-100/61-ZR, A319 in 2009). In parallel to these changes, another impacting factor is the move toward variable frequency generators ranging from 360 to 800 Hz, whereas high frequency increase the risk of long arc duration. One major concern in the study of electric arc in aeronautical environment is the presence of FOD. Indeed, a metallic object is able to initiate an electric arc if it shorts out two conductors. This kind of event is very hard to quantify and rationalize as many factors can impact the phenomenon occurrence and behavior, such as the object and conductors geometry, nature, the object dynamics (dropped or static), the potential presence of other pieces in the vicinity of the arc and the electrical parameters (voltage, current and power magnitude and limits, signal frequency). With the addition of the electric arc chaotic behavior, its global study is very complex. The most important factors in aeronautical systems for the arc events and their impact are summarized in Fig. 1 Obviously, it is impossible to conclude in one study on how an arc will behave under all the possible combination of parameters presented above. In the present article, two parameters are studied, the object nature and geometry, all the others being fixed. First a statistical analysis is realized to determine whether a relationship exists between the object type and the object and arc behavior. The industrial purpose of this section is to give, if possible, guidelines for the design of mechanical protections. In a second time, the electric measurement are correlated to high speed camera images of the arc to correlate the events to electrical signatures. The aim of the second section is to provide reference signature for the electrical signals in the presence of an arc.

.[2]

II. EXPERIMENTAL SETUP

The experimental setup used to simulate a part of aircraft power distribution system is presented in Fig. 3. The left and right electrodes are connected between two phases of a threephase Pacific power source 65 kVA. The phase to phase voltage is set at 400 VRMS and the frequency at 400 Hz, in agreement with recent aircraft standards. It is important to note that this is not the most critical configuration from the electric arc point of view, the frequency being allowed to vary between 360 and 800 Hz. This variable frequency brings many challenges such as a diminution of the arc extinction probability at high frequency. The arc voltage and current are measured thanks to HZ-115 and Hioki 3275 probes respectively, connected to an oscilloscope Hioki MR8847-02. The video acquisition is achieved thanks to a Photron SA4 high speed camera. The framerate is set at 20 000 fps with a 1/30 000 s shutter, approximately 2.5 s can be recorded, enough to observe the entire arc duration (the arc is extinguished manually if it last too long, in aircraft condition it is very unlikely to observe an arc this long as circuit breaker would have tripped before). To avoid overexposing the camera sensor in the presence of the arc, 64 or 400 filters are used. All the acquisition devices are trigged on the contactor 2. Temperature and humidity are measured but not controlled. The choice of the objects used as FOD is rather complex. The question beeing: how to simulate a FOD when it can be almost anything (screw, wire, metal shavings, etc.)? As it is impossible to test every objects, conductive wire were selected as FOD. These objects present simple geometry that can be easily controlled and varied (in length, diameter and nature of the metal), allowing to study the impact of the object shape and nature. In the present study two kinds of wires were tested, cores of aeronautical cables (stranded, gauge from 20 to 08 AWG) and solid cables (diameter from 0.3 to 4 mm). Four different metal natures were selected, copper, Ni-platted copper, aluminum and stainless steel. The details of the FOD tested are summarized in the Table I

A. Statistical study of arc ignition

The first step of the study was a statistical study of the arc ignition for the different FOD used as defined in Table 1. Prior to the detailed description, a major observation was that in every case an electric arc occurred. This is very important because it means that, in this specific setup with the electrical parameters used, when a foreign object is present and static on the electrodes, it will always lead to an electric arc. However, all the arc observed were not equivalent in terms of duration, extinction and volume occupied by the plasma.

In order to attempt to understand the origin of these different behavior, an analysis of the FOD behavior is presented in the Fig. 4. Two cases are discriminated, when the FOD is ejected, when the FOD is welded on the electrode (either on one or two sides). The first observation from the figure above is that the FOD behavior is correlated to the object mass. This is likely due to the fact that the force opposed to the ejection is directly proportional to the object mass, making heavier objects harder to eject. The second observation is that stranded and solid wires present the same comportment. This means that the interface between the object and the electrodes have little impact on the arc behavior. As it shows no difference, to keep simpler explanations, for the rest of the study the analysis will be focused on solid wires. The third observation is that there is a distinct transition between FOD being ejected and welded when the mass increased. Following this observation, the initial hypothesis was that, when the FOD is ejected, it allowed the arc to expand more easily as the arcs would stay "attached" to the object (as in Fig. 2). If this assumption was verified, it would give a useful tool to sort the FOD as a function of the damaging potential of the arc they can initiate, the lighter objects leading to more damages. In order to assess this, the arc dimensions were computed thanks to the high speed camera acquisition. The raw and processed images for three different cases with increasing size can be observed in Fig. 6. 12 mm diameter, (center) aluminum 1.92 mm diameter and (right) aluminum 0.7 mm diameter. The processed images correspond to successive integrations of the area occupied by the plasma. The color scale correspond to the time information. The image on the center presents two very distinct color areas, this traduces the occurrence of two arcs, a first one small, leading to the partial ejection of the FOD, followed by a second larger when the FOD fell on the electrodes. On the picture on the right, straight lines can be observed on the top, this is not plasma but the path of ejected molten metal.

The processed images correspond to successive integration of the area occupied by the arc during its whole duration. Several information can be extracted from these images, the maximum width and height taken by the plasma, the total surface and also the area occupied by the plasma at given time thanks to the color scale. The time information is certainly a useful data for the development of active protection, to rule on the available time to extinguish the arc before it reaches a certain dimension and also to quantify the damages induced.The pictures in Fig. 6 confirm that all the arc are not equivalent in size, in order to assess the hypothesis of the correlation between the FOD mass and the arc behavior, the arc maximum width and height are plotted according to the objects mass in the Fig. 5.

As it is possible to observe in the figure above, the hypothesis expressed previously is not verified. Indeed, even if there is a slight decay in the arc width when increasing the object mass, the arc dimensions appear to have a chaotic behavior making it impossible to state on a link with the FOD mass. It is important to note that the same study was performed to analyze if there was a pattern linked to the FOD diameter or metal nature reaching the same conclusion although the metal nature might have an impact on the composition of the plasma [START_REF] André | Arc tracking energy balance for copper and aluminum aeronautic cables[END_REF].

Another important observation is that the dimensions of the plasma can reach very large values, width and height about 20 cm have been observed. Compared to the size of aircraft power distribution systems, an arc this large might damage or be attached to close equipment. In the present study, the depth of the arc is not acquired, however it is planned for further investigation to use another high speed camera to get 3D information of the plasma dimensions.

To conclude on the outputs of the statistical study of the arc behavior, the data to extract in order to design power distribution system is the maximum dimension of the area that can be potentially occupied by an arc and the time the arc takes to reach these dimensions. This gives the available duration to extinguish the arc.

This leads to the second part of the study corresponding to the detection of the arc events thanks to their electrical signatures, in order to design active protection.

B. Arc electrical signatures

In the first part of the study the arc events were observed thanks to a high speed camera. In this section, the electrical measurements made in parallel are studied to identify the signatures linked to the arc behavior. From the statistical investigation carried out, several event were distinguished:

-Arc ignition -Arc expansion -Arc contained -Arc moving freely -Arc extinction While the first and last event correspond to the beginning and the end of the arc, the three others correspond to different arc behaviors.

Three of these events can be observed in the Fig. 7. The event in the picture (1) corresponding to the arc ignition is characterized for the voltage and current by a non-zero value. The electric arc is barely visible at this state. The current amplitude is here limited (about 185 A to 220 A) by the power supply and not by the limiting resistance. In this test, the FOD > REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < was completely ejected leading to a large expansion of the arc. This behavior can be observed in the pictures ( 2) and (3). From the voltage point of view, it is traduced by a monotonous increase in magnitude and for the current by a decrease in magnitude. The current decrease is likely to be due to the power limitation of the power supply. The current stays pretty stable until the voltage reach a certain amplitude and then starts to rapidly drop until the arc extinction that is observed in the picture [START_REF] Naef | Proposed Transient Recovery Voltage Ratings for Power Circuit Breakers[END_REF]. In this test it is hard to state weather the arc is extinguished due to the voltage limitation or the power limitation. In the voltage case, the extinction mechanism involved would be the arc expansion until the source wouldn't be able to provide enough voltage for the arc to grow. In the power limitation case, the mechanism could be the decrease in current due to that limit, leading to the arc extinction when the current is not sufficient to sustain the arc.

The next waveform is the one corresponding to the arc contained, meaning that the arc length remains stable over time. The typical waveform associated can be observed in Fig. 8. The voltage wasn't null before the arc because the FOD in this test was ejected by a previous electric arc, the arc in this figure is the one ignited when the FOD fell on the bars. The arc ignition correspond to the picture [START_REF] Roboam | New trends and challenges of electrical networks embedded in 'more electrical aircraft[END_REF]. Only one luminous spot can be observed due to the fact that the object is soldered on one electrode, so only one arc occurred between the object and the other electrode. For the whole arc duration the voltage remain pretty stable (here around 15 -25 V), this correspond to the image (3), except of the Transient Recovery Voltage (TRV) as in [START_REF] Teulet | One-dimensional hydrokinetic modelling of the decaying arc in air-PA66-copper mixtures: II. Study of the interruption ability[END_REF]. The TRV is a fast rise of the voltage reaching high voltage values. This voltage is not correlated to the arc length as indicated by the image (2) where the plasma area remains small. In this case, the TRV is always greater than the dielectric strength of the medium leading to a reignition at each polarity change [START_REF] Naef | Proposed Transient Recovery Voltage Ratings for Power Circuit Breakers[END_REF]. The last waveform observed was the one associated to a free burning arc. This behavior occurred when the arc was no longer attached to the FOD, but moved freely on the electrodes. It generally generates longer arcs (in duration), the arc is initiated at the FOD location, then moves toward the end of the electrodes due to the Laplace force. In order to avoid damages induced by the arc, heat shrink sleeve was placed near the end of the electrodes. Generally, free burning arcs were stopped and then extinguished by the heat shrink sleeve. From the voltage measurement in Fig. 9 it is possible to observe the classical arc ignition and extinction signatures in (1) and ( 4). On the other hand, the typical waveform of the free burning arc is also visible on the voltage. It is characterized by large variations of the amplitude but, at the contrary of the expansion behavior, the changes in magnitude are not monotonous. This behavior is present before (2) with a steady increase of the voltage and then between (2) and (3) with a sharp decrease of the voltage. This corresponds to a fast motion of the arc that can also be observed in the pictures (2) and [START_REF] André | Arc tracking energy balance for copper and aluminum aeronautic cables[END_REF], where the arc right spot "jumped" from the back to the front of the image.

C. Arc fault detection

The three examples of electrical measurements presented in the previous section prove that typical signatures exist, allowing to identify the arc behavior from the voltage and the current waveforms. If these information bring more tools to understand the electric arc, they cannot help directly to develop arc detection method.

Two difficulties rise when trying to detect arc fault. The variety of waveform generated by the large number of loads present in the electrical network makes complex the detection of arc. Some loads might have an electrical signature very close from the arc one. This leads to the second difficulty, the necessity to avoid false positive. In an aircraft, false positive must be avoided. As can be imagined, the loss of a critical equipment due to a false alarm can be unsafe for the flight.

Current protection devices against arc fault (AFCB and SSPC) are generally based on current monitoring [5]- [START_REF] Leprettre | Detection of electrical series arcs using time-frequency analysis and mathematical morphology[END_REF]. A way to avoid false positives and to increase the chances of arc fault detection is likely to be the tracking of multiple arc markers [START_REF] Andrea | Principle of Arc Fault Detection for Solid State Power Controller[END_REF]. Some authors also studied the impact of the arc on the voltage waveform for arc fault detection purposes [START_REF] Momoh | Design and analysis of aerospace DC arcing faults using fast fourier transformation and artificial neural network[END_REF].

In order to investigate the impact of the arc on the voltage shape, the frequency content of the dV/dt was extracted before, during and after the electric arc. The spectrograms, plotted in Fig. 10, were obtained thanks to Fast Fourier Transform (FFT) applied on windows with 2.5 ms duration, corresponding to one period of the signal. The arc signature is clear from the three spectrograms, before and after the arc, the frequency content of the signal doesn't span above 2.6 kHz, whereas during the arc it spans above 100 kHz. In the present case, using the frequency analysis of the dV/dt would allow differentiating the electric arc from the zero voltage, open circuit and short circuit.

To validate the reliability of the dV/dt frequency content as one of the relevant markers of the electric arc, it would necessitate to carry out arc tests initiated by FOD with different kind of loads connected in parallel and/or in series and also to use electrical sources more representative of aircraft power generation.

IV. CONCLUSION

The electric arc initiated by a FOD in aeronautical conditions has been studied. From the statistical study of the FOD behavior it was determined that a transition occurred between the object being ejected and welded on the electrodes when increasing the mass of the object. The hypothesis that the object behavior would impact the arc dimensions was invalidated thanks to the characterization of the arc length and width from the high speed camera acquisitions. As the FOD mass and nature didn't seemed to impact the arc damaging potential, it appears to be difficult to give guidelines for the design of a mechanical protection for aircraft power distribution systems. From this consideration, the need for active protection was confirmed.

Three main conclusions can be drawn from the present study. In order to determine the time in which a protection device should be able to detect and secure a system, the analysis of the arc dimensions over time, compared to the dimensions of an aircraft power distribution system can provide very interesting data. If detected, the arc behavior can be identified from the voltage waveform. This information can be useful as the arc damaging potential differ from a contained arc to a free burning arc. The best strategy for the detection of the arc with the higher success ratio and the less false positive detection is likely to be achievable thanks to the monitoring and comparison of several arc electrical markers or even other type of markers. Even if generally the studies focus on the current, the voltage contains many information on the arc and the dV/dt frequency content seems to be a relevant marker of the arc.

Further steps of this study include carrying out this experiment in presence of other system parts close to the electrodes to determine the behavior of the arc and using two high speed cameras to get volume information on the arc.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Predominant factor impacting the arc ignition, behavior and extinction in aeronautical environment.From the industrial angle, two strategies are being developed, passive:-Segregation distances between bars -Insulating material -Mechanical protection to prevent FOD from entering in contact with electrodes and active:-Arc detection devices (Arc Fault Circuit Breaker and Solid State Power Controller)

Fig. 3 .

 3 Fig. 3. Scheme of the experimental setup. Geometrical details are not given due to confidentiality matters.

Fig. 4 .

 4 Fig. 4. Probability of FOD ejection as a function of their mass.

Fig. 5 .

 5 Fig. 5. (a) Arc width and (b) height as a function of the FOD mass. The bars and the point represent the results distribution and the mean value respectively.

Fig. 6 .

 6 Fig.6. (top) processed and (bottom) raw images of the electric arc, FOD (left) aluminum 3.12 mm diameter, (center) aluminum 1.92 mm diameter and (right) aluminum 0.7 mm diameter. The processed images correspond to successive integrations of the area occupied by the plasma. The color scale correspond to the time information. The image on the center presents two very distinct color areas, this traduces the occurrence of two arcs, a first one small, leading to the partial ejection of the FOD, followed by a second larger when the FOD fell on the electrodes. On the picture on the right, straight lines can be observed on the top, this is not plasma but the path of ejected molten metal.

Fig. 7 .Fig. 8 .

 78 Fig. 7. (Top) Images of the electric arc. (Middle) Voltage measurement. (Bottom) Current measurement. The FOD was aluminum solid wire with 0.7 mm in diameter.

Fig. 9 .

 9 Fig. 9. (Top) Images of the electric arc. (Middle) Voltage measurement. (Bottom) Current measurement. The FOD was aluminum solid wire with 0.4 mm in diameter.

  . > REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

TABLE I Summary

 I

		of the FOD tested
	Type	Metal nature	Diameter (mm) / Gauge
	Stranded	Aluminum	0,7 / 22
	Stranded	Ni-platted copper	1 / 20 1,2 / 18 1,4 / 16	1,7 / 14 2,7 / 10 3,7 / 08
			0,3	1,6
			0,4	1,8
	Solid	Aluminum	0,7	2
			1	3
			1,5	4
	Solid	Stainless steel	0,7 1,5	2
	Solid	Copper	0,7 1,5	2

III. EXPERIMENTAL RESULTS AND DISCUSSION


François Clautiaux

Ruslan Sadykov

François Vanderbeck

Quentin Viaud

Pattern based diving heuristics

Keywords: Cutting and Packing, Dynamic Programming, Diving Heuristic

come

Introduction

The two-dimensional guillotine cutting-stock problem is a classical problem that occurs in industry when one has to cut rectangular pieces from identical large objects (plates) in such a way that the number of plates used is minimum. In wood, glass, steel, or paper industries, most cutting devices impose so-called ˚Corresponding author. Tel.: +33 5 40 00 21 37 Email addresses: francois.clautiaux@u-bordeaux.fr (François Clautiaux), ruslan.sadykov@inria.fr (Ruslan Sadykov), francois.vanderbeck@u-bordeaux.fr (François Vanderbeck), quentin.viaud@u-bordeaux.fr (Quentin Viaud) guillotine cuts. These cuts are made in a straight line from one border of the plate to the other. Motivated by an industrial application, we restrict our attention to cutting patterns that emanate from a 4-stage guillotine-cut process: i.e. each cut piece can be obtained trough the process of applying at most 4 guillotine cuts.

In our problem, we are given an item set I. Each item i P I is a rectangle with dimensions w i ˆhi , called width and height, and has a demand (or number of copies) equal to d i . Each item i can be rotated, in which case, its dimensions become h i ˆwi . To cut these items, an unlimited number of identical plates (bins) of dimension W ˆH are available. In the remainder of the paper, we assume that all input data are integer.

In practice, set I is partitioned into distinct batches tI 1 , . . . , I B u, which have to be processed independently. This saves operational costs, but leads to a larger amount of leftovers (also called waste). To reduce the loss of material, it is allowed to use the leftover from the last plate of the previous batch to initialize the current batch. A leftover is the result of a first stage cut, i.e. it is a rectangular piece of height H and width W 1 ă W . Therefore, for a given batch, a solution consists of a number of cutting patterns for standard plates, plus a first-plate cutting pattern that has a smaller width, and a last-plate cutting pattern for which the objective is to minimize the width. We call this problem two-dimensional guillotine cutting-stock problems with leftover (2DGCSPL). A solution is depicted in Figure 1.

The problem is NP-hard, since it generalizes cutting-stock. It is also combinatorially complex and the practical instances that we need to solve are large: 10 to 15 batches each having about 150 different items, for which the total demand can be as large as 400 copies, and the plates have a large size p6000ˆ3000q relative to the item size. Hence, we consider only so-called restricted cuts, i.e. cuts of length equal to an item width w i or an item height h i . Special requirement related to restricted cuts is that one of the two sub-plates is immediately initialized with an item of width or height equal to the length of the performed cut. Our initial experimental observations as well as those reported by Furini et al. [START_REF] Furini | Modeling two-dimensional guillotine cutting problems via integer programming[END_REF] for a related problem confirm that using only restricted cuts does not deteriorate the solutions in most cases. The paper is organised as follows. We first give an overview of the literature on related problems in Section 2. Then in Section 3, we recall how it can be modelled by an exponentially large ILP model in which each variable represents a cutting pattern, and how the pricing subproblem can be solved by dynamic programming. Constructive and evolutionary heuristics inspired from the literature for the 2DGCSPL are discussed in Section 4. In Section 5 we present our pattern based diving heuristics for the 2DGCSPL. To obtain better solutions by these diving heuristics, we also propose a partial enumeration technique that is embedded in the dynamic program for solving the pricing subproblem. In Section 6, we report results of computational experiments in which we compare different heuristic algorithms on classical data sets and real-life instances. Conclusions are drawn in Section 7.

I b I b I b I b

Literature review

Cutting-stock problems have drawn a large attention because of their significance for the industry, and their theoretical and practical difficulty. When all demands are unitary (d i " 1, @i P I), the problem is also referred in the literature as the bin-packing problem. From a theoretical point of view, one can reformulate a bin-packing problem as a cutting-stock problem, therefore through the rest of this paper, we will use both names.

The first study on two-dimensional packing problems appeared in Gilmore and Gomory [START_REF] Gilmore | Multistage cutting stock problems of two and more dimensions[END_REF]. Therein, the approach is to formulate the two-dimensional bin-packing problem (2BP) using Linear Programming (LP) and solve it using column generation. Indeed, the Dantzig-Wolfe decomposition of the 2BP gives a nice reformulation with a covering problem as the master and a two-dimensional knapsack as pricing problem. In Gilmore and Gomory [START_REF] Gilmore | Multistage cutting stock problems of two and more dimensions[END_REF] the two-dimensional knapsack pricing problem is tackled by dynamic programming. For the threestage version of the problem, Vanderbeck [START_REF] Vanderbeck | A nested decomposition approach to a three-stage, two-dimensional cutting-stock problem[END_REF] used a nested decomposition, solving the pricing problem using Dantzig-Wolfe reformulation. When the pricing problem is too hard to solve, a level approach can be used instead as outlined by Puchinger and Raidl [START_REF] Puchinger | Models and algorithms for three-stage two-dimensional bin packing[END_REF]. It combines different methods such as heuristics, meta-heuristics or ILP models. Methods are called one after another and when an acceptable solution is found, it is returned. Another way to solve 2BP was proposed by Silva et al. [START_REF] Silva | An integer programming model for two-and three-stage two-dimensional cutting stock problems[END_REF]. It uses an enumeration of the possible patterns occuring during the cutting process. These patterns are fed into a direct Integer Linear Programming (ILP) formulation, which is solved with a general purpose ILP solver. Another approach based on a large ILP model has been developed by Macedo et al. [START_REF] Macedo | Arc-flow model for the two-dimensional guillotine cutting stock problem[END_REF]. Authors used the original arc-flow model of Valério de Carvalho [START_REF] Valério De Carvalho | Exact solution of bin-packing problems using column generation and branch-and-bound[END_REF] for the one-dimensional bin-packing problem (1BP) and extended it to the 2BP.

Most exact methods for the 2BP are not able to scale up to larger instances. This comes from the pseudo-polynomial nature of the models used (either in an ILP formulation, or in the subproblem). To overcome this difficulty, heuristics are often preferred to exact methods. New advanced heuristics such as the agentbased heuristic of Polyakovsky and M'Hallah [START_REF] Polyakovsky | An agent-based approach to the twodimensional guillotine bin packing problem[END_REF] or the repairing heuristic of Fleszar [START_REF] Fleszar | Three insertion heuristics and a justification improvement heuristic for two-dimensional bin packing with guillotine cuts[END_REF] have been shown to be effective.

Column generation or pattern based heuristics have also been proposed. Alvarez-Valdes et al. [START_REF] Alvarez-Valdes | A computational study of lp-based heuristic algorithms for two-dimensional guillotine cutting stock problems[END_REF] used a simple and a more elaborated ad-hoc rounding heuristics. The same authors also used a truncated restricted master heuristic. Furini et al. [START_REF] Furini | A column generation heuristic for the two-dimensional two-staged guillotine cutting stock problem with multiple stock size[END_REF] employed pure diving, restricted master, and diving with sub-MIPing heuristics, if one uses the terminology defined by Sadykov et al. [START_REF] Sadykov | Primal heuristics for branch-and-price: the assets of diving methods[END_REF].

Since the subproblem is time consuming, the previous authors used heuristics to solve it most of the time. Finally, Cintra et al. [START_REF] Cintra | Algorithms for two-dimensional cutting stock and strip packing problems using dynamic programming and column generation[END_REF] worked with column generation approach with non-proper columns. They solved the residual problem obtained after initial rounding using an ad-hoc heuristic.

The main bottleneck of the Dantzig-Wolfe reformulation of the 2BP is the two-dimensional knapsack problem (2KP) pricing subproblem. When the pricing problem is unbounded (i.e. when there are no upper bounds on item productions), dynamic programming (DP) is a natural choice to handle pricing. Recurrence relations were initially introduced by Beasley [START_REF] Beasley | Algorithms for unconstrained two-dimensional guillotine cutting[END_REF]. Russo et al. [START_REF] Russo | An exact dynamic programming algorithm for large-scale unconstrained two-dimensional guillotine cutting problems[END_REF] extended them to tackle huge pricing problems. However it is not straightforward to extend DP to the bounded case, or to a larger number of stages. In Dolatabadi et al. [START_REF] Dolatabadi | Exact algorithms for the twodimensional guillotine knapsack[END_REF], the authors proposed a DP method for the case where the number of cutting stages is unlimited. They combine dynamic programming and implicit enumeration of patterns: assuming a threshold value as a valid lower bound, they enumerate all feasible cutting patterns providing a value better than the threshold. If there are no such patterns, they correct the threshold value and run the method again. For the bounded version of the four-stage version of the problem, Clautiaux et al. [START_REF] Clautiaux | Combining dynamic programming with filtering to solve a four-stage two-dimensional guillotine-cut bounded knapsack problem[END_REF] proposed a hypergraph based label setting algorithm. An initial dynamic program, represented by a hypergraph, is created to solve the unbounded problem. Then a filtering procedure based on Lagrangian relaxation is used to remove non-promising hyperarcs from the hypergraph. Finally, a label setting algorithm is applied to solve the 2KP by handling incrementally item upper bound constraints one by one. Note that when the number of cutting stages is limited to two or three, ILP models perform well as shown by Lodi and Monaci [START_REF] Lodi | Integer linear programming models for 2-staged two-dimensional knapsack problems[END_REF] and Puchinger and Raidl [START_REF] Puchinger | Models and algorithms for three-stage two-dimensional bin packing[END_REF].

Leftovers complicate significantly the problem, mainly because the objective function has to account for the total length of the last cutting pattern. The rounded up lower bound obtained by column generation is seldom equal to the optimal solution value. An application of the cutting-stock problem with leftovers to glass cutting was treated by Puchinger et al. [START_REF] Puchinger | Solving a real-world glass cutting problem[END_REF]. The problem they considered is to minimize the number of plates and the length of the used part of the last plate. Extra constraints are imposed on the order in which items are cut due to storage considerations. The authors have designed heuristics, meta-heuristics and heuristic branch-and-bound methods to solve the problem. Recently, similar approaches have been proposed by Dusberger and Raidl [START_REF] Dusberger | A variable neighborhood search using very large neighborhood structures for the 3-staged 2-dimensional cutting stock problem[END_REF]. They developed a Variable Neighbourhood Search based on a Ruin-and-Recreate principle. A cutting pattern for a plate is represented by a cutting tree. Heuristics aim to remove some parts of the cutting tree for some plates (ruin step), then a new solution is obtained from items removed by the ruin step (recreate step). Recreate step is based on a constructive heuristic; the authors provide an ILP model to handle this part. Dusberger and Raidl [START_REF] Dusberger | Solving the 3-staged 2-dimensional cutting stock problem by dynamic programming and variable neighborhood search[END_REF] then extended their own work by adding dynamic programming in their recreate step. Instead of using heuristics, Andrade et al. [START_REF] Andrade | Two-stage two-dimensional guillotine cutting stock problems with usable leftover[END_REF] proposed direct MILP models for the 2BP with leftover inspired from the works of Lodi and Monaci [START_REF] Lodi | Integer linear programming models for 2-staged two-dimensional knapsack problems[END_REF] and with variables and constraints to handle leftover.

We are not aware of existing works in the literature studying our variant of the cutting-stock problem with four stage cutting, possible rotation, and demands. Moreover, we want to handle leftovers and solve real-life size instances with large plates (typically 6000 ˆ3000) and many items (typically more than 100).

Extended formulations for 2DGCSPL

A classical approach solving for hard packing problems is to use extended formulations resulting from a Dantzig-Wolfe decomposition. In this section, we first describe how this decomposition can be applied to the 2DGCSPL, yielding an exponential size reformulation. The latter is to be solved using column generation and we explain how the pricing subproblem can be solved. Finally, we describe a pseudo-polynomial size extended formulation.

Recall that for 2DGCSPL, one needs to pack set I b of items to three type of plates: a leftover plate from batch b ´1 of given dimension W 1 ˆH (type 1), an unlimited number of standard plates of dimension W ˆH (type 2), and a single potential leftover plate of dimension W ˆH (type 3); and only a part of the width of the plate of type 3 is used (it is equals to the width of the cutting pattern assigned to it). Using these three types, the objective is to minimize the total used plate width and assuring that exactly one plate of type 1 and at most one plate of type 3 are used respectively. For the first batch, the plate of type 1 does not exist; the model can easily be adapted to this special case.

An extended formulation obtained by Dantzig-Wolfe decomposition

To formulate 2DGCSPL, let P t be the set of all valid cutting patterns for a plate of type t P t1, 2, 3u. Let a p i and w p be the number of items i P I in cutting pattern p P P t and its width. Let also λ p be an integer variable which is equal to the number of times pattern p P P t is used in the solution. Then the 2DGCSPL can be formulated as

min ÿ pPP1 W 1 λ p `ÿ pPP2 W λ p `ÿ pPP3 w p λ p (1)
ÿ pPP1YP2YP3 a p i λ p " d i , @i P I, (2)
ÿ pPPt λ p " 1, @t P t1, 3u, (3)
λ p P Z `, @p P P 2 , (4)
λ p P t0, 1u, @p P P t , t P t1, 3u.

Objective function (1) minimizes the total used plate width. Note that the first term is a constant, since the first plate has to be used. Constraints [START_REF] Andrade | Two-stage two-dimensional guillotine cutting stock problems with usable leftover[END_REF] guarantee that the item demands are satisfied. Constraint (3) forces one to use exactly one plate of the first and the third type (the latter can have the width equal to zero, so implicitly we enforce an upper bound of one on type 3 plates). The remaining constraints are added to ensure integrality of the variables. Note that the number of variables λ is exponential. A standard way to tackle formulation (1)-( 5) is to solve its linear relaxation, called Master Problem (MP) using a column generation procedure: one alternates between solving the Restricted Master Problem (RMP), which is the MP with a restricted number of variables λ, and the pricing problem. The latter searches for negative reduced cost variables λ to be added to the RMP. Consider a dual solution pπ, µq of the RMP, where vector π is associated to constraints (2) and vector µ " pµ 1 , µ 3 q is associated to constraints [START_REF] Beasley | Algorithms for unconstrained two-dimensional guillotine cutting[END_REF]. The pricing problem is to find a cutting pattern p P P 1 Y P 2 Y P 3 that minimizes the following reduced cost cp " ´ÿ iPI a p i π i `$ & % W 1 ´µ1 , p P P 1 , W, p P P 2 , w p ´µ3 , p P P 3 .

(6)

Obviously, the pricing problem decomposes into three subproblems, one for each plate type.

Solving the column-generation pricing problem

Each pricing subproblem is a bounded four-stage restricted guillotine-cut two-dimensional knapsack problem. A dynamic programming algorithm (DP) for the unbounded case of this problem (based on the one of Beasley [START_REF] Beasley | Algorithms for unconstrained two-dimensional guillotine cutting[END_REF]) was presented by us in Clautiaux et al. [START_REF] Clautiaux | Combining dynamic programming with filtering to solve a four-stage two-dimensional guillotine-cut bounded knapsack problem[END_REF]. To be self-contained, we recall the dynamic programming recursion here. Then, we restate the complete formulation of the bounded case as an ILP model.

When item rotations are allowed, set I of items is duplicated in the subproblem to include rotated item copies. The two copies of an item i share the same dual value π i from the RMP, and the same upper bound d i .

Since we use restricted cuts, we consider two types of dynamic programming states: those related to a restricted cut (used to initiate a strip), and those that are used to complete a strip.

Let pw, h, sq, w P t1, . . . , W u, h P t1, . . . , Hu, s P t1, . . . , 3u, be the state corresponding to the situation in which a rectangular part of dimension w ˆh of the plate is to be separated from the current plate using a guillotine cut of stage s. Let pw, h, sq, w P t1, . . . , W u, h P t1, . . . , Hu, s P t2, . . . , 4u, be the state corresponding to the same situation with the additional restriction that the next cut should obtain a single item copy. For a given state pw, h, sq, let U pw, h, sq be the maximum value of a configuration obtained from this state. Note that this value only depends on the dual values π associated to item copies, and an additional term for type 3 plate only, which is equal to the total width used. It is sufficient to account for this restricted width in stage 1 cuts only. Then U pW, H, 1q equals to the pricing subproblem optimum value (ignoring constant values in [START_REF] Clautiaux | Combining dynamic programming with filtering to solve a four-stage two-dimensional guillotine-cut bounded knapsack problem[END_REF]).

Let Wpw, hq and Hpw, hq be the set of all possible widths and heights of items which fit into rectangle w ˆh. For the type 1 plate generation, the recursions are identical, only initialization needs to be adapted to the available width. For type 3 plates, the recursive formula for U pw, h, 1q, becomes U pw, h, 1q " max " 0, max

w 1 PWpw,hq
U pw 1 , h, 2q ´w1 `U pw ´w1 , h, 1q ( * where term ´w1 represents the penalty for the width consumed. Using the formalism of Martin et al. [START_REF] Martin | Polyhedral characterization of discrete dynamic programming[END_REF], our dynamic program may be reformulated as a max cost flow problem in a directed acyclic hypergraph. In Clautiaux et al. [START_REF] Clautiaux | Combining dynamic programming with filtering to solve a four-stage two-dimensional guillotine-cut bounded knapsack problem[END_REF], this transformation as well as pre-processing procedures for this hypergraph are presented. Let G " pV, Aq denote this directed acyclic hypergraph. Vertex set V is composed of all states from the dynamic program and also of special states related to single items or waste. Set S of these special states are sources of the hypergraph. Its sink u is the vertex associated to state pW, H, 1q. Each hyperarc a has a head hpaq, which contains a unique vertex, and a tail multiset T paq, which contains one or more vertices representing subplate cutting patterns, each of which can occur more than once. A hyperarc a represents the cutting decision that turns state hpaq into states in T paq. A vertex v has a hyperarc predecessor (resp. successor) set Γ ´pvq (resp. Γ `pvq).

Let Apiq be the set of hyperarcs whose tail sets include the source node representing item i P I. Let x a , for a P A, be the integer variable representing the flow value going through hyperarc a. The resulting flow model associated to the pricing problem can therefore be written as an Integer Linear Program:

max f p " ÿ iPI π i ÿ aPApiq x a (7)
s.t.

ÿ aPΓ ´pvq x a ´ÿ a 1 PΓ `pvq x a 1 " 0, @v P Vztu Y I Y Hu (8) ÿ aPΓ ´puq x a " 1 (9)
ÿ aPApiq

x a ď d i , @i P I (10)

x a P N, @a P A (11)
Objective function [START_REF] Dolatabadi | Exact algorithms for the twodimensional guillotine knapsack[END_REF] aims to maximise the reduced cost of items. This objective yields the non-constant term of the reduced cost f p " ř i a p i π i in [START_REF] Clautiaux | Combining dynamic programming with filtering to solve a four-stage two-dimensional guillotine-cut bounded knapsack problem[END_REF]. Constraints ( 8)-( 9) are classical flow conservation constraints. They ensure that a valid pattern is built according to cutting properties. Constraint set [START_REF] Fleszar | Three insertion heuristics and a justification improvement heuristic for two-dimensional bin packing with guillotine cuts[END_REF] avoids overproduction of an item. When p P P 3 , value w p in ( 6) is computed as: w p " ř aP Ā w a x a , with Ā the set of hyperarcs related to a cut of first stage. Objective function [START_REF] Dolatabadi | Exact algorithms for the twodimensional guillotine knapsack[END_REF] therefore becomes max pPP3 tf p ´wp u.

Tackling the pricing problem by applying an ILP solver to formulation ( 7)-( 11) leads to large computing times. We have shown in Clautiaux et al. [START_REF] Clautiaux | Combining dynamic programming with filtering to solve a four-stage two-dimensional guillotine-cut bounded knapsack problem[END_REF] that an iterative method based on state-space relaxations allows one to solve the pricing problem significantly faster. However, even with the latter approach remains too slow to be used inside column generation. As explained in Section 5, our method will rather rely on solving the unbounded version of the pricing problem.

Deriving a pseudo-polynomial flow formulation for 2DGCSPL

Building on our pricing problem network flow reformulation, one can formulate our cutting-stock problem as an arc-flow model yielding an ILP formulation of pseudo-polynomial size (as done by Valério de Carvalho [START_REF] Valério De Carvalho | Exact solution of bin-packing problems using column generation and branch-and-bound[END_REF] for the 1BP).

Let G t " pV t , A t q be the hypergraph related to a plate of type t P t1, 2, 3u. Let z t , t P t1, 2, 3u be the total number of patterns of type t that are used (z 1 and z 3 are set to be equal to one, whereas z 2 is a variable). Then, the 2DGCSPL model can take the form:

min W 1 z 1 `W z 2 `z3 ÿ aP Ā3 w a x a (12)
s.t.

ÿ aPΓ ´pvq x a ´ÿ a 1 PΓ `pvq x a 1 " 0, @v P V t ztu t Y I Y Hu, @t P t1, 2, 3u, (13)
ÿ aPΓ ´putq x a " z t , @t P t1, 2, 3u (14)
ÿ tPt1,2,3u ÿ aPAtpiq x a " d i , @i P I (15) z 1 " z 3 " 1, z 2 P N (16)
x a P N, @a P A t , @t P t1, 2, 3u

Objective ( 12) is equivalent to [START_REF] Alvarez-Valdes | A computational study of lp-based heuristic algorithms for two-dimensional guillotine cutting stock problems[END_REF]. Constraints ( 13)-( 14) are flow conservation constraints among the different plate types. Constraints [START_REF] Lodi | Recent advances on two-dimensional bin packing problems[END_REF] ensure that the demand is fulfilled. Variable definitions [START_REF] Lodi | Integer linear programming models for 2-staged two-dimensional knapsack problems[END_REF] ensure to select only one pattern from types 1 and 3. Finally, each hyperarc variable is integer [START_REF] Macedo | Arc-flow model for the two-dimensional guillotine cutting stock problem[END_REF].

Constructive and evolutionary heuristics

Due to the large size of instances we aim to solve, a natural approach is to use heuristics, which are usualy used in the literature. In this section, we first propose heuristics to solve the bounded case of the pricing problem. Such a choice is motivated by the fact that solving this problem to optimality is computationally expensive in the column generation context. Second, we also derive constructive heuristics for 2DGCSPL. They will be used in computational comparisons with diving heuristics presented below.

We designed two heuristics for the pricing problem. The first is based on a hypergraph exploration. The second is an evolutionary algorithm.

Hypergraph based heuristics for the pricing problem with bounded production

As described in Section 3.2, the dynamic program seeks a max-cost flow to the sink in the acyclic hypergraph, where vertices correspond to sub-plate cutting patterns and hyperarcs correspond to combinations of these patterns. Thus a cutting pattern can be represented by a set of hyperarcs (how sub-plates were combined) and a set of vertices (intermediate states). As a first step of the heuristic, we run the dynamic program for the unbounded problem. It allows us to compute the best reduced cost cv associated with each vertex v P V, i.e. the dynamic programming value associated to the corresponding state. The reduced cost associated with an hyperarc ca is the sum of the reduced costs of its tails : ca " ř vPT paq cv . We also define as ypA 1 q P Z |I| `the partial solution corresponding to the source vertices in the tail sets of hyperarcs in multiset A 1 :

y i pA 1 q "| A 1 X Apiq |, @i P I.
The constructive stage of the first heuristic is run for a given hyperarc a 1 P A and a given partial solution y 1 P Z |I| `, y 1 ď d. We start by adding ypta 1 uq to y 1 . At any time, we keep the current set L of open vertices, initialized with all non-source vertices in the tail set of hyperarc a 1 . In every iteration of the algorithm, a vertex v 1 is selected in L and removed from it. Then, a hyperarc a incoming to v 1 with the smallest reduced cost such that y 1 `yptauq ď d is chosen, all non-source vertices in the tail set of a are added to L, and we pass to the next iteration. The algorithm stops when set L becomes empty. It returns the obtained solution y 1 and the corresponding multiset of hyperarcs A 1 . The pseudocode of this heuristic is presented as Function HG-Constr-Heur.

Function HG-Constr-Heur(a 1 , y 1 )

L Ð H; A 1 Ð H; if y 1 `ypta 1 uq ď d then y 1 Ð y 1 `ypta 1 uq; A 1 Ð ta 1 u; L Ð T pa 1 qzS; while L ‰ H do Pick v P L; L Ð Lztvu ; A 2 Ð ta P Γ ´pvq : y 1 `yptauq ď du; if A 2 ‰ H then a 2 Ð argmax aPA 2 tc a u; y 1 Ð y 1 `ypta 2 uq; A 1 Ð A 1 Y ta 2 u; L Ð L Y T pa 2 qzS ; return py 1 , A 1 q
The constructive heuristic can then be embedded in a local search method. First the constructive heuristic finds the best feasible solution among all hyperarc a 1 P Γ ´puq incoming to the sink node u, and best solution py ˚, A ˚q is obtained. Then the following local search algorithm is applied to this solution. At every iteration, each hyperarc â P A 1 is replaced by another hyperarc ã P Γ ´pHpâqq, ã ‰ â incoming to the same vertex as â, and solution y 1 is modified accordingly. If an improved solution is found, next iteration starts. The process stops when no improvement occurred or the iteration number limit is reached (minimum between |I| and 50). The formal presentation of this heuristic is given as Algorithm 1.

Evolutionary heuristic for the pricing problem with bounded production

The second heuristic is an evolutionary algorithm inspired from Hadjiconstantinou and Iori [START_REF] Hadjiconstantinou | A hybrid genetic algorithm for the two-dimensional single large object placement problem[END_REF]. Let Ī be an unary representation of item set i P I, i.e. d i copies of each item i are created. We use an indirect encoding: each individual (or genome) g is represented by a sequence of integers g 1 , . . . , g npgq , each of them refers to the g j -th item to cut in Ī. Size npgq of a genome g is usually smaller than | Ī|, as, generally, not all item copies fit into a single plate. To obtain the solution and its value from a genome, we run the following so-called First-Fit heuristic. First-Fit heuristic cuts items in the order given by genome g. Let R be the stack of available sub-plates initialized with initial plate pW, H, 1q at stage 1. First-Fit heuristic takes the first available plate in R, removes it from R, and then tries to cut the first available item l j in the order given by the genome. Since guillotine cuts are considered here, this implies that this cut always divides the current sub-plate r into two smaller sub-plates r 1 and r 2 which are added to stack R. Note that if no remaining item in the genome fits in a given plate r, we use other items which fit into it if possible, thus increasing the solution quality. If no items fit into a sub-plate, it is discarded. The process stops when R " H.

The purpose of the evolutionary algorithm is to create a new population containing improving solutions from the previous one. We initialize our evolutionary algorithm with a population of |I| individuals corresponding to different random permutations of multiset Ī. To ensure good quality solutions, the initial population is then reduced using an elitist strategy. This one aims to keep the pool F of p size best individuals from the initial population. Starting from this pool, a new population is created in two phases: new individual generation and crossover operation. The first phase is the same as the initial population generation. This enriches the pool by p size new individuals. In the second phase, the offspring set O of individuals is produced using a two-point crossover where only 25% of offspring individuals are randomly kept.

Two-point crossover uses as input two parent individuals represented by their genomes g 1 and g 2 and create a new one g 3 by selecting parts of parent genomes.

Two positions in the genome of the first parent p 1 and p 2 are randomly picked. The offspring is then created by merging genome of parent g 1 between positions 0 and p 1 , genome of parent g 2 between positions p 1 and p 2 and then genome of parent g 1 between positions p 2 and npg 1 q.

After a crossover phase, all remaining individuals solutions are moved in pool F . Then only p size best individuals in F are kept for the next population generation. To have enough solutions from one population generation to another, the population size p size is set to 20. The total number of population generations is set to |I|{2.

Heuristics for the 2DGCSPL

A direct way to obtain a heuristic solution for the 2DGCSPL is to use the evolutionary algorithm described in Section 4.2 iteratively for each plate. We call this heuristic peaq.

Alternatively one can use standard bin-packing list heuristics such as Next-Fit (NF), Best-Fit (BF), First-Fit (FF) or Bottom-Left (BL). These heuristics run in polynomial time and produce a feasible packing for our problem. The main difference between them is the packing strategy for a given item. Assume a set of items to cut I and a list of empty sub-plates R. In Next-Fit heuristic, a given item i P I is packed in the next sub-plate r P R in which it fits. If there are previous sub-plates in which the item does not fit, they are simply removed from R. In First-Fit heuristic, a given item i P I is packed in the first sub-plate r P R in which it fits. In Best-Fit heuristic, a given item i P I is packed in the best sub-plate r P R in which it fits. The best sub-plate is the one such that after packing an item its remaining horizontal or vertical unused space is minimised. In Bottom-Left heuristic, a given item i P I is packed in the subplate r P R in which it fits such that position of the item is the closest to the bottom left corner of the bin. At any step of these heuristics, if no sub-plates can accommodate the item, a new one is created, which triggers the addition of a new bin. More details about these heuristics can be found in Lodi et al. [START_REF] Lodi | Recent advances on two-dimensional bin packing problems[END_REF].

We now present heuristic algorithm piubq which is more time consuming than peaq. It combines several list heuristics as well the evolutionary algorithm. We consider 10 ˚|I| random permutations of set Ī. Four list heuristics mentioned above are applied for each of these permutations. These heuristics can be used in several ways. The first option is to fill one plate at a time, i.e. we switch to the next plate only when there are no items in the current list that fit into the current plate. The second option is to open all bins at the same time.

The third option is to implement the following two-phase process. A plate of infinite width is created and items are packed in it using the list heuristics above. According to the guillotine cut property, the obtained solution can be represented by a certain amount of vertical strips. We collect the obtained strips and assign them to original plates using BF, FF and NF list heuristics as in 1BP.

At the same time, we also generate an extra bin-packing solution by solving iteratively 2KPs with our evolutionary algorithm. For the latter we generate 10 times an initial random population and improve it during |I|{10 iterations.

Using the three described options as well as the evolutionary algorithm, we obtain a set of feasible solutions for each random permutation of set Ī. We first record the best solution providing us a valid upper bound for our cutting problem. Then, among all cutting-stock solutions we select a plate of smallest waste, and we add it to the partial solution. We then reiterate the whole process on the residual problem. The algorithm terminates when there are no more items to cut.

A diving heuristic for the 2DGCSPL

In this section we present our diving heuristic to solve 2DGCSPL. We first introduce standard column generation based diving algorithms, following Sadykov et al. [START_REF] Sadykov | Primal heuristics for branch-and-price: the assets of diving methods[END_REF]. Then we show how these algorithms can be adapted when production bound constraints are not enforced in the pricing subproblem. Several additional ingredients are presented to improve the quality of the solutions obtained.

Diving heuristics

A pure column generation based diving heuristic, as presented in Sadykov et al. [START_REF] Sadykov | Primal heuristics for branch-and-price: the assets of diving methods[END_REF], is a depth-first search heuristic in a branch-and-price enumeration tree obtained when branching on the master problem variables λ. It can be directly applied to our problem in the following way. At each iteration one (re)optimizes the master linear problem; as a result, one obtains a fractional solution λ. One chooses a component λ p of the solution vector such that its current value λp is the closest to its nearest non-zero integer r λp u. Then one sets λ p " r λp u in the partial solution, and updates the demand of items: d Ð d ´ap r λp u. The algorithm proceeds to the next iteration with the updated (or residual) master problem. The process terminates when the partial solution becomes complete (i.e. when all master constraints are satisfied). Note that Furini et al. [START_REF] Furini | A column generation heuristic for the two-dimensional two-staged guillotine cutting stock problem with multiple stock size[END_REF] used this pure diving heuristic for two-stage 2D guillotine cutting-stock problem.

As shown in Sadykov et al. [START_REF] Sadykov | Primal heuristics for branch-and-price: the assets of diving methods[END_REF], effectiveness of the diving heuristic significantly increases when it is coupled with Limited Discrepancy Search. This variant allows one to do a limited number of backtracks in the search tree. During a backtrack, the decision to add a column to the partial solution is reverted. Once this happens, this column is added to the tabu list. Columns in the tabu list cannot be added to the partial solution.

To improve the convergence of the column generation procedure on which the diving heuristic relies, we use the automatic dual price smoothing stabilization proposed in Pessoa et al. [START_REF] Pessoa | Automation and combination of linear-programming based stabilization techniques in column generation[END_REF].

Diving with non proper columns

To achieve primal feasibility, the basic diving heuristic requires to use socalled proper columns, i.e. variables that could take a non-zero value in an integer solution of the residual master problem. In our context, a variable λ p , p P P 1 Y P 2 Y P 3 , is proper if a p i ď d i , @i P I. Therefore, in each pricing subproblem we should impose upper bounds on the number of copies of items in the cutting pattern. However when solving the pricing problem [START_REF] Clautiaux | Combining dynamic programming with filtering to solve a four-stage two-dimensional guillotine-cut bounded knapsack problem[END_REF] in Section 3.2, we consider only the unbounded version. In the presence of these upper bounds, the pricing problem becomes significantly harder to solve to optimality. A possible solution to this issue is to solve the pricing problem only heuristically using algorithms from Sections 4.1 and 4.2. However, our preliminary experiments showed that this approach deteriorates significantly the quality of solutions obtained by the diving heuristic.

In this work, we propose a variant of the diving heuristic which uses nonproper columns. A previous study of Cintra et al. [START_REF] Cintra | Algorithms for two-dimensional cutting stock and strip packing problems using dynamic programming and column generation[END_REF] and our preliminary experiments showed that the lower bound obtained by solving the MP with non-proper columns is close to the one obtained when using exclusively proper columns. As the quality of diving heuristics depend mainly on the strength of the MP bound, we may then expect that a "non-proper" diving heuristic will be efficient for our problem.

Our "non-proper" diving heuristic proceeds as follows. Remember that at each iteration, the residual master problem is solved by column generation. Both proper and non-proper columns may be generated. However, the partial solution can only be augmented with proper columns. Therefore, given a fractional solution λ, we choose a proper variable λ p with a value closest to its nearest non-zero integer. If such variable exists, we proceed the same way as in the basic diving heuristic. If there is no such "proper " variable, we choose a column λ p with the smallest reduced cost (with respect to the optimal dual solution of the master problem) among • all proper columns contained in the current RMP;

• and proper columns generated by solving the bounded pricing problem by heuristic algorithms presented in Sections 4.1 and 4.2.

As usual, this column is then added to the partial solution with the value equal to the nearest non-zero integer of λp . In particular, if λp " 0 then λ p " 1 is included in the partial solution.

Our preliminary experiments showed that fixing a type 3 columns λ p , for p P P 3 , early in the search has a negative impact on the quality of solutions obtained. Therefore, we adopt the following modification. Cutting patterns p P P 3 are never added to the partial solution before patterns p P P 1 Y P 2 . As there is exactly one pattern p P P 3 in any feasible solution, once it is added to a partial solution, the latter should become complete. Therefore, each time the partial solution is augmented with a cutting pattern of type 1 or 2, we verify heuristically whether the remaining item copies can be cut into one plate of dimension W ˆH. If it is possible, we produce a cutting pattern p P P 3 including all remaining item copies and minimizing heuristically its width w p . Then this pattern is used to complete the solution, and the diving heuristic terminates.

One can develop further this idea and, each time the partial solution is augmented, to formulate the residual 2DGCSPL and solve it with heuristics described in Sections 4.1 and 4.2. Calls to these heuristics are done iteratively for each plate until the residual problem instance is closed. This modification can be seen as a combination of the diving and pricing heuristics. This combined heuristic is presented in Algorithm 2. A boolean parameter lastP lateOnly is used to set how to evaluate the residual 2DGCSPL instance: building only a solution for the last plate or a complete solution to the residual 2DGCSPL instance.

Algorithm 2: The combined pricing heuristics and diving heuristic with non-proper columns To introduce a Limited Discrepancy Search (LDS) in our "non-proper" diving heuristic, we need to ensure that at least one non-tabu column is produced by our pricing problem heuristic. This can be achieved by generating a sufficient number of different cutting patterns, i.e. the size of set P heur in Algorithm 2 is strictly larger that the current size of the tabu list.

Note that in addition to accelerating column generation, stabilization has an another positive effect. When it is used, the MP solution is generally more fractional and thus contains a larger number of columns. Therefore, there is a larger probability to have a proper column in the fractional solution.

Solving the pricing problem with partial enumeration

Although we have explained above that our diving heuristic can be adapted to handle non-proper columns, still the quality of the solutions may be decreased in comparison with the "proper" case. In this section, we propose a modification to our dynamic program that partially takes into account upper bounds on the number of item copies in order to favour the generation of proper cutting patterns. As it will be seen from our computational results, this approach allows one to improve the quality of the lower bound obtained by solving the MP as well as the quality of solutions produced by the diving heuristic. This comes at the cost of a slightly larger dynamic program. Nevertheless we show below that this can be controlled by a suitable configuration.

On one hand, complete enumeration of the possible patterns would take into account the bound constraints, but the computational cost would be huge. On the other hand, the dynamic program has a reasonable computational cost, but it does not take into account the bound constraints. Our idea is to mix both approaches, by replacing some parts of DP by a partial enumeration.

This idea is implemented using so-called meta-items, each one representing a partial vertical or horizontal stack of item copies satisfying production upper bounds. When restricted states U pw, h, kq, k " 2, 3, 4 are considered, instead of choosing one item to initiate the stripe, we choose a meta-item (or equivalently the items that it represents). There is potentially an exponential number of possible meta-items to initiate the stripe. Therefore, we introduce an additional parameter δ, which restricts the possible meta-items produced by only considering items whose width/height is close enough to the size of the stripe.

Formally, let ām i be the number of copies of item i P I included into metaitem m. Given three values 0 ă w ď W , 0 ă h ď H, and 0 ă δ ď min iPI w i , we define the following set M x pw, h, δq of vertical meta-items. Each meta-item m P M x pw, h, δq forms in the cutting pattern a partial vertical stack of width w containing copies of items i P I such that w ´δ ă w i ď w and h i ď h. Items that do not belong to m may only be cut in other vertical stacks or in the same stack above the item copies in m. In addition, copies of items i P I such that ām i ą 0 may only be cut in other vertical stacks. Formally:

m P M x pw, h, δq ô $ ' ' & ' ' % Di P I, w i " w : ām i ą 0, ām i ą 0 ñ w ´δ ă w i ď w and h i ď h, @i P I, ām i ď d i , @i P I, ř iPI ām i h i ď H.
The first condition ensures that one item has width w (and thus a restricted pattern is built). The second condition ensures that the size of the items in the meta-item satisfies the requested limitations. The third condition ensures that the meta-item satisfies the production bound constraints. The fourth condition ensures that the meta-item height does not exceed the plate height.

Analogously, given values 0 ă w ď W , 0 ă h ď H, and 0 ă δ ď min iPI h i , we define the following set M y ph, w, δq of horizontal meta-items:

m P M y ph, w, δq ô $ ' ' & ' ' % Di P I, h i " h : ām i ą 0, ām i ą 0 ñ h ´δ ă h i ď h and w i ď w, @i P I, ām i ď d i , @i P I, ř iPI ām i w i ď W.
For each meta-item m P M, we define its total value πm " ř iPI π i ām i . For each vertical meta-item m P M x , we define its total height hm " ř iPI ām i h i , and for each horizontal meta-item m P M y , we define its total width wm " ř iPI ām i w i . Note that, by definition, M x pw, h, δq " H if w R WpW, Hq, and M y ph, w, δq " H if h R HpW, Hq. Suppose now that for each w P WpW, Hq a value δ w , 0 ď δ w ď min iPI w i , is fixed, and for each h P HpW, Hq a value δ h , 0 ď δ h ď min iPI h i , is fixed. Then the recursive formulae for states pw, h, sq can be rewritten in the following way without loss of any proper patterns from the set of feasible solutions.

U pw, h, 2q " $ & % max mPM x pw,h,δwq: hmďh πm `U pw, h ´h m , 2q ( , if δ w ą 0, max iPI: wi"w,hiďh tπ i `U pw, h ´hi , 2qu , if δ w " 0, U pw, h,
π i `U pw, h ´hi , 4q ( * if δ w " 0,
If all values δ are fixed to zero, the modified dynamic program reduces to the original one presented in Section 3.2. The larger the values δ are, the fewer non-proper cutting patterns are generated and, at the same time, the larger is the number of meta-items. So, there is a trade-off between the complexity (or the size) of the dynamic program and the strength of the approximation of the space of proper cutting patterns by the space of feasible solutions of the dynamic program. We parametrize this trade-off by defining thresholds ∆ size ě 0 on the size of the sets of meta-items and ∆ diff ě 0 on values δ, respectively for dimension w and h. Given these thresholds, values δ are determined the following way. For each w P WpW, Hq, we set δ w to the largest value δ ď min ∆ diff w , min iPI w i ( such that | M x pw, H, δq |ď ∆ size w . As W x pw, H, 0q " H for any w, such value δ always exists. Analogously, for each h P HpW, Hq, we set δ h to the largest value δ ď min ∆ diff h , min iPI h i ( such that | M y ph, W, δq |ď ∆ size h . Note again that W y ph, W, 0q " H for any h. The sets of meta-items are computed by enumeration.

To illustrate such generation of meta-items, let us consider an instance with a raw plate of size p4, 3q and three items 1ˆi 1 " p4, 1q, 2ˆi 2 " p3, 1q, 1ˆi 3 " p2, 1q. We compute set M x p4, 3, δq with δ " 2. Initially meta-items composed of only items i 1 are created (see Figure 2-(b)). Secondly extra items are added to metaitems derived from item i 1 . Since δ " 2, only item i 2 is added (see Figure 2-(c) and Figure 2-(d)). From each created meta-item, a valid cutting pattern can be obtained (see Figure 3). At this point, item i 3 is not cut because δ " 2. If the cutting process continues with horizontal cuts it will be possible to cut item i 3 . A representation of the resulting patterns is given in Figure 4. Note that all meta-items containing item i 2 are created, therefore it is impossible to add it again to patterns of Figure 3-(a) and Figure 3-(b).

Computational experiments

In this section, we report the results of experiments that we conducted on instances from the literature and on real-world data. The goal of our experiments is threefold: piq to evaluate the impact of the pricing problem with partial enumeration on lower and upper bounds produced by our algorithm; piiq to evaluate the quality of the solutions produced by different variants of the diving heuristic and compare them to heuristics described in Section 4.3; piiiq to measure the impact of re-using leftovers on real-world instances with batches.

All experiments are run using a 2.5 Ghz Haswell Intel Xeon E5-2680 with 128Go of RAM. We used CPLEX 12.6 to solve linear programs. The time limit to solve one instance using one algorithm is set to one hour.

Instances and objective functions

Two objective functions are considered. The first one is the classical cuttingstock objective (minimizing the number of plates). The second one is the total used plate width (takes into account leftovers).

We use two groups of instances. The first group has been proposed in the literature for the classical non-guillotine bin-packing problem. These instances are divided in ten classes. We denote them as χ " U rα, βs if value χ is uniformly generated in interval rα, βs. Six instance classes were proposed by Berkey and Wang [START_REF] Berkey | Two-dimensional finite bin-packing algorithms[END_REF] and are built as follows:

Class 1: W " H " 10, w i , h i " U r1, 10s Class 2: W " H " 30, w i , h i " U r1, 10s Class 3: W " H " 40, w i , h i " U r1, 35s Class 4: W " H " 100, w i , h i " U r1, 35s Class 5: W " H " 100, w i , h i " U r1, 100s Class 6: W " H " 300, w i , h i " U r1, 100s

The four remaining classes with W " H " 100 were created by Martello and Vigo [START_REF] Martello | Exact solution of the two-dimensional finite bin packing problem[END_REF]. Authors divided items in types:

Type 1: w i " U r2{3W, W s, h i " U r1, 1{2Hs Type 2: w i " U r1, 1{2W s, h i " U r2{3H, Hs Type 3: w i " U r1{2W, W s, h i " U r1{2H, Hs Type 4: w i " U r1, 1{2W s, h i " U r1, 1{2Hs
The four classes are: Class 7: Type 1 with probability 70%, Type 2, 3, 4 with probability 10% each Class 8: Type 2 with probability 70%, Type 1, 3, 4 with probability 10% each Class 9: Type 3 with probability 70%, Type 1, 2, 4 with probability 10% each Class 10: Type 4 with probability 70%, Type 1, 2, 3 with probability 10% each In these instances, the total number of item copies ( ř iPI d i ) in each instance belongs to t20, 40, 60, 80, 100u. 10 instances for each class are generated. Thus, the first group contains 500 instances. They are available at http://or.dei. unibo.it/library/two-dimensional-bin-packing-problem. We denote as C ´In the set of all instances in the first group (in classes 1 ´10) having n items.

The second set of industrial instances is retrieved from the real-world data. In our instances H ˆW P t50 ˆ100, 500 ˆ1000, 3000 ˆ6000u, |I| P t25, 50, 100u. The average demand of items is between 2 and 3. 15 instances were generated for each pair plate size/number of items. Thus, the second group contains 135 instances. We denote as R ´In the set of all instances in the second group having n items.

The industrial instances are available on the ESICUP web-site: https:// paginas.fe.up.pt/~esicup/datasets.

Impact of the partial enumeration

Remember that the configuration of the partial enumeration in the pricing problem is characterized by four values: ∆ size w , ∆ size h , ∆ diff w , and ∆ diff h . We consider here three settings. The first one corresponds to the standard dynamic program with no enumeration: ∆ 1 " p∆ size w " 0, ∆ size h " 0, ∆ diff w " 0, ∆ diff h " 0q. The second one corresponds to the enumeration of items with the same h i for odd cutting stages and with the same w i for even cutting stages: ∆ 2 " p∆ size w " 1000, ∆ size h " 1000, ∆ diff w " 1, ∆ diff h " 1q. Our third setting corresponds to the enumeration of items with the same h i for odd cutting stages and enumeration of items with different w i at the second cutting stage: ∆ 3 " p∆ size w " 1000, ∆ size h " 1000, ∆ diff w " 8, ∆ diff h " 1q. For each setting, we report in Table 1 the geometric mean of the number (in thousands) of vertices (V), hyperarcs (A) and the total number of hyperarc tails ( ř |T paq|). It can be seen from this table that the partial enumeration impact on the size of the hypergraph is very reasonable (at most 33% of increase). We do not report the construction time since it is negligible (less than 0.2 seconds on average and at most 6 seconds). Hypergraph building is a one time operation. Therefore a stricter comparison to measure hypergraph configuration has to been done when used to solve the pricing problem. When partial enumeration is applied, many non-proper patterns are excluded from the solution space of the pricing problem. Therefore, the lower bound obtained by column generation may be improved. In the next set of experiments we computationally estimate this improvement. We have implemented four variants of column generation. The first one, denoted as cg mip , uses MILP flow formulation ( 7)- [START_REF] Furini | A column generation heuristic for the two-dimensional two-staged guillotine cutting stock problem with multiple stock size[END_REF] to solve the pricing problem and generates only proper columns. Thus, the best possible "proper" lower bound is obtained at the expense of a large solution time. The three other variants use dynamic programming with partial enumeration to solve the pricing problem. We denote them as cg ∆1 , cg ∆2 , and cg ∆3 depending on the partial enumeration configuration.

Instances |V| |A| ř |T paq| group ∆ 1 ∆ 2 ∆ 3 ∆ 1 ∆ 2 ∆ 3 ∆ 1 ∆ 2 ∆ 3 C 2.
The results are reported in Tables 2 and3. In both tables, the first column reports the instance class and its total number of instances between brackets. The column #pp shows the number of instances solved to optimality using the fast heuristic peaq : the obtained solution value equals to the value of the trivial "surface" lower bound. The next three columns present results for the variant cg mip : number of instances not solved by peaq for which the column generation converged within one hour, average time (t) in seconds, and the average primaldual gap (gap) in percentage from the best known solution. Next three columns give the average gap (gap #opt ) for other three variants of column generation. In order to have a correct comparison, averages in columns gap #opt are calculated only for instances for which the variant cg mip converged. Note that the column generation variants with dynamic programming converged within the time limit for all instances. Thus in columns gap all and t all , we give the average gap and the average solution time among all instances not solved optimally by heuristic peaq for all hypergraph configurations. Reported times are in seconds. In One can see from the results in Tables 2 and3 that the dynamic program is orders of magnitude faster than MIP for solving the pricing problem. The partial enumeration increases the running time of column generation only marginally. Moreover, this technique allows one to obtain a lower bound which is close to the "proper" lower bound, at least for the easiest instances that can be tackled by cg mip . Unfortunately, we were not able to determine how close is the lower bound obtained by cg dp∆3 to the "proper" bound for larger and harder instances, as the latter bound is very time consuming to obtain. Anyway, application of the partial enumeration technique significantly increases the quality of lower bounds obtained by column generation at virtually no cost. Therefore, it offers a good trade-off between the quality of lower bounds and the column generation running time.

It can also be seen that column generation is slower for real-life instances. This is expected as the running time of the dynamic program depends on the plate size, which is larger for the instances in the second group.

In the next experiment, we estimate the impact of the partial enumeration on the pseudo-polynomial formulation ( 12)- [START_REF] Macedo | Arc-flow model for the two-dimensional guillotine cutting stock problem[END_REF] described in Section 3.3. By the result of Martin et al. [START_REF] Martin | Polyhedral characterization of discrete dynamic programming[END_REF], in the absence of constraints [START_REF] Lodi | Recent advances on two-dimensional bin packing problems[END_REF], this formu-lation has the integrality property. Thus the value of its linear programming relaxation is equal to the lower bound obtained by column generation with "nonproper" columns. As we have just seen, the latter may be increased using partial enumeration. Thus, the strength of the linear programming relaxation of formulation ( 12)-( 17) based on partly enumerated hypergraph may be improved too. Moreover, variables x may have coefficients greater than one in constraints [START_REF] Lodi | Recent advances on two-dimensional bin packing problems[END_REF] if there are hyperarcs which correspond to cutting several copies of one item, as in the case of cutting a meta-item. Thus, partial enumeration enables the possibility of adding knapsack cutting planes when solving the formulation by a MIP solver.

The results for direct solution by the CPLEX MIP solver of formulation ( 12)-( 17) with different settings of the partial enumeration of hypergraph are shown in Tables 4 and5. Here we first report the number #pp of instances solved by heuristic peaq. Then for each way to build the hypergraph, we report the number of remaining instances solved in one hour (#opt) with MIP formulation. In the last part of the table, the average time (t) in seconds required to solve all instances is reported. If an instance is not solved within the time limit, we use the time limit value in the calculation of the average time. Reported times are in seconds. According to the results, partial enumeration indeed increases the efficiency of the formulation, as more instances are solved within the time limit and the average solution time is decreased. However, the improvement is not radical. One can also notice that the instances with the total used plate width objective as well as real-life instances are significantly harder to solve.

#opt t, sec. Instances (#) #pp M IP ∆ 1 M IP ∆ 2 M IP ∆ 3 M IP ∆ 1 M IP ∆ 2 M IP ∆ 3 C-I20 (
(#) M IP ∆ 1 M IP ∆ 2 M IP ∆ 3 M IP ∆ 1 M IP ∆ 2 M IP ∆ 3 C-I20 (

Comparison of heuristics

In this experiments, we compare five heuristics for the 2DGCSPL: (i) evolutionary algorithm peaq;

(ii) algorithm piubq, which combines evolutionary algorithm with list heuristics;

(iii) a diving heuristic denoted pdiv H q without partial enumeration in the pricing problem and with simple evaluation of the residual problem in the diving (parameter lastP lateOnly " true);

(iv) a diving heuristic denoted pdivq with partial enumeration ∆ 3 in the pricing problem and with simple evaluation of the residual problem in the diving (parameter lastP lateOnly " true);

(v) a combination of the diving heuristic and the evolutionary algorithm with complete evaluation of the residual problem in the diving (parameter lastP lateOnly " f alse), which we denote as pedivq.

The first two heuristics are presented in Section 4.3. The other three heuristics are described in Section 5. We also consider variants with Limited Discrepancy Search for the last two heuristics and denote them as pdiv 32 q and pediv 32 q. In these variants, the backtrack is allowed up to the depth of 2 of the search tree and the maximum size of the tabu list is 3. This LDS parametrisation results in at most 10 dives in the search tree. Note that diving heuristics are always initialized with the solution produced by heuristic peaq.

In Tables 6 and7, we report the average gap in percentage from the best known solution and the average time in seconds for the five heuristics and two more variants. As observed from Table 6, diving heuristics with LDS led to the best results for the total number of used plates objective: they produced solutions with the best average quality within a few seconds. For group of instances R, a good quality solution is almost always obtained by fast heuristic peaq.

Nevertheless, instances with the total used plate width objective are more interesting for us, as they represent batches in the real problem to solve. From Table 7 one can see that diving algorithms clearly outperform the first two heuristics. Heuristic peaq is the fastest but also produces solutions of the worst quality. Heuristic piubq improves the solution quality at the expense of much larger running time. However, it struggles with real-life instances, as the solu-tion improvement over peaq for this instances is very small. Diving algorithms pdiv H q and pdivq significantly outperform heuristic piubq both in terms of running time and solution quality. The partial enumeration technique increases the effectiveness of the diving heuristic for most instances at a very small cost. This technique is especially useful for large instances with the total used plate width objective. The combination pedivq of the diving heuristic and the evolutionary algorithm further improves the quality of the obtained solutions at a cost of a reasonable increase of running time. The best solutions on average are obtained by diving heuristics with Limited Discrepancy Search. However, the running time of these heuristics is quite long especially for large real-life instances. Thus, to our opinion, the heuristics pdivq and pedivq offer the best tradeoff between solution quality and running time.

In the industrial context, the 2DGCSPL is solved iteratively for different batches of the global problem. To validate our methodology, we consider now real-world instances with batches and plate dimension H ˆW " p3000 ˆ6000q. Each instance is composed of 10 or 15 batches, each batch b is composed of |I b | P t100, 150u items. As previously, average demand of items is between 2 and 3. For each class, we retrieved 25 instances. The class named LmIn contains instances with m batches, each of which having n items. All these instances are are available on the ESICUP web-site: https://paginas.fe.up. pt/~esicup/datasets.

In Table 8, we report the results obtained by the same heuristics tested before. Here we do not report results only for heuristic pdiv H q as it was shown to be dominated by pdivq. We also do not report results obtained by div 32 and ediv 32 since required computation time is very long and obtained solutions are only slightly better than the solutions obtained after only one dive. For each heuristic, we show the solution value obtained which equals to the average number of plates needed to cut items from all batches, as well as the average solution time in minutes. We also give the average lower bound plbq on the optimal value. For each instance, this value is obtained by iteratively computing the rounded-up column generation lower bound for each batch and determining the length of the leftover plate for the next bin based on this bound. To measure the impact of using leftover plates, in column "w/o lo" we give the average of the best solution values obtained by our algorithms for the problem variant in which the batch leftovers cannot be reused in the next batch. As observed in Table 8, keeping a leftover plate for the next batch allows one to save up to 10.6 plates or 3.7% of plates on average. Again, heuristic peaq is the fastest. The more expensive heuristic piubq improves on peaq only marginally. Much better results are obtained by diving heuristics. The standard diving heuristic pdivq saves up to 4.1 plates or 1.4% of plates on average. The extended diving heuristic pedivq saves up to 4.3 plates or 1.5% of plates on average. Moreover the gap with the lower bound is at most 1.9 plates or 0.7% of plates on average using pdivq. For the pedivq this drops to 1.7 plates or 0.6% on average. In our opinion heuristic pdivq offers the best solution quality -running time trade-off. Even if its running time reaches 2 hours for the largest instances, its application in practice is still realistic. These instances correspond to a one day planning horizon. Therefore, spending two hours to obtain a solution seems to be reasonable.

Conclusions

In this work we studied the two-dimensional guillotine cutting-stock problem with leftovers (2DGCSPL) which is solved consecutively when the set of items is partitioned in batches. Given the difficulty and large size of real-life instances, we considered the special case with restricted cuts.

To solve the 2DGCSPL, we have proposed column generation based diving heuristics. An originality of our work is that these heuristics work with nonproper cutting patterns. This variant simplifies the pricing problem but makes it more difficult to obtain feasible solutions for the 2DGCSPL. We proposed several ways to overcome this difficulty, including combination with an evolutionary algorithm and a partial enumeration technique. The latter reduces the number of generated non-proper cutting patterns, tightens the column generation lower bounds and improves the quality of solutions obtained by the diving heuristics.

The computational experiments on the literature and real-life instances revealed that the proposed diving heuristics outperformed significantly the constructive and evolutionary heuristics. The largest improvements are achieved on large instances, as well as on real-life instances and for the instances with the total width objective. The experiments on real-life production plant instances with batches showed that our heuristics run in a reasonable time and allow the decision maker to save up to 1.5% of raw material on average.

Our diving heuristic with non-proper columns is generic and could be applied to other two-dimensional cutting-stock problems. It would be interesting to see how the heuristic running time evolves when considering instances with a different number of stages or with non-restricted cuts. Generalisation to the case with plate defects is especially useful for practical purposes.

We also hope that the publication of real-life instances will inspire more research on this family of problems and will allow the researchers to perform a fair computational comparison between future algorithmic approaches.

Figure 1 :

 1 Figure 1: Representation of a solution for a batch I b . The first plate has a smaller width, since it is the leftover from the previous batch. The right-hand part of the last plate will be re-used in the next batch.

2 :Figure 3 :Figure 4 :

 234 Figure 3: Vertical cutting patterns for a bin p4, 3q using meta-items (b),(c) and (d) from Figure 2

  We now give the recursion for computing values U pw, h, sq and U pw, h, sq for type 2 plates.

			"		*
	U pw, h, 1q " max	0, max w 1 PWpw,hq	U pw 1 , h, 2q `U pw ´w1 , h, 1q	(
			"		*
	U pw, h, 2q " max	0, max h 1 PHpw,hq	U pw, h 1 , 3q `U pw, h ´h1 , 2q (
			"		*
	U pw, h, 3q " max	0, max w 1 PWpw,hq	U pw 1 , h, 4q `U pw ´w1 , h, 3q	(
	U pw, h, 2q "	max iPI: wi"w,hiďh	tπ i `U pw, h ´hi , 2qu
	U pw, h, 3q "	max iPI: hi"h,wiďw	tπ i `U pw ´wi , h, 3qu
			"	*
	U pw, h, 4q " max	0,	max iPI: wi"w,hiďh	π i `U pw, h ´hi , 4q (
	In the recursive formulae for U pw, h, sq and U pw, h, 4q, the alternative with zero
	value corresponds to turning the remaining rectangular part of the plate into
	waste.
			ď		ď
	Wpw, hq "				tw i u, Hpw, hq "	th i u.
	iPI: wiďw,hiďh	iPI: wiďw,hiďh

Algorithm 1 :

 1 Hypergraph based heuristic py ˚, A ˚q " p0, Hq; for a 1 P Γ ´puq do py, Aq Ð HG-Constr-Heurpa 1 , 0q; if πy ą πy ˚then py ˚, A ˚q " py, Aq ; , A 1 q " py ˚, A ˚q; for â P A 1 do Let  be the part of solution (flow) A 1 coming to Hpâq;

	k Ð 0;
	repeat
	improve Ð f alse; k Ð k `1;
	py 1

for ã P Γ ´pHpâqq, ã ‰ â do pỹ, Ãq Ð HG-Constr-Heur `ã, y 1 ´yp Âq ˘; if π ¨py 1 ´yp Âq `ỹq ą πy ˚then y ˚Ð y 1 ´yp Âq `ỹ; A ˚Ð A 1 z  Y Ã; improve Ð true; until improve or k " mint|I|, 50u; Store solution py ˚, A ˚q;

  ˚Ð solution of the 2DGCSPL with demands d by the evolutionary algorithm 2 d 1 Ð d, P part Ð H 3 repeat

	4	Solve the MP with production bounds d 1 by column generation and
		obtain solution	λ

1 P 5 P prop Ð tp P P 1 Y P 2 : λp ą 0, a p ď d 1 u (set of proper patterns in the solution) 6 if P prop ‰ H then 7 p 1 Ð argmin pPP prop t| λp ´rλ p u |u 8 else 9 P RMP Ð set of patterns of type 1 and 2 in the RMP 10 P heur Ð set of heuristic solutions to the pricing problem of type 1 and 2 11 p 1 Ð argmin pPP RMP YP heur tc p u 12 P part Ð P part Y tp 1 u, d 1 Ð d 1 ´ap 1 13 if lastP lateOnly " true then 14 p 3 Ð heuristic solution to the pricing problem of type 3 15 if d 1 ´ap 3 " 0 then 16 d 1 Ð 0 17 if cost `Ppart Y tp 3 u ˘ă costpP ˚q then P ˚Ð P part Y tp 3 u 18 else 19 P ev Ð heuristic solution to the residual 2DGCSPL with demand d 1 20 if cost `Ppart Y P ev ˘ă costpP ˚q then P ˚Ð P part Y P ev 21 until d 1 " 0 22 return P

Table 1 :

 1 Size of the hypergraph (geo. mean., in thousands of vertices, arcs, and tails) with different configuration of the partial enumeration

		6 3.5	3.5	25.2	29.5	29.6	50.0	57.7	57.9
	R	8.3 10.9 10.4 123.9 139.7 162.0 254.0 256.8 340.1

Table 3

 3 column (#pp) is deliberately omitted since no instances are solved by heuristic peaq.% cg ∆1 cg ∆2 cg ∆3 cg ∆1 cg ∆2 cg ∆3 cg ∆1 cg ∆2 cg ∆3

				cg mip			gap #opt ,%		gap all , %	t all
	Instances #pp #opt t gap, C-I20 (100) 63 37 85 6.94 8.63 8.38 8.29 8.63 8.38 8.29 0.1 0.1 0.2
	C-I40 (100)	37	54 583 2.48 2.97 2.80 2.79 6.03 5.82 5.81 0.6 0.6 0.6
	C-I60 (100)	35	33 921 1.26 1.63 1.48 1.48 3.30 3.14 3.14 1.2 1.3 1.3
	C-I80 (100)	29	28 1079 1.01 1.23 1.09 1.09 2.56 2.42 2.42 2.5 2.6 2.6
	C-I100 (100) 31	18 731 0.53 0.66 0.61 0.61 2.83 2.72 2.72 4.5 4.6 4.6
	R-I25 (45)	37	2 1271 10.90 10.96 10.91 10.91 22.22 22.18 22.17 0.4 0.4 0.4
	R-I50 (45)	34	0	-	-	-	-	-	13.85 13.82 13.82 1.9 1.9 2.0
	R-I100 (45)	33	0	-	-	-	-	-	6.53 6.53 6.53 10.7 10.6 10.9

Table 2 :

 2 Comparison of different column generation variants for the total number of used plates objective % cg ∆1 cg ∆2 cg ∆3 cg ∆1 cg ∆2 cg dp∆3 cg ∆1 cg ∆2 cg ∆3

			cg mip			gap #opt ,%		gap all , %	t all
	Instances gap, C-I20 (100) #opt t 62 249.0 1.99 3.82 3.14 3.04 3.79 3.07 3.00 0.3 0.3 0.3
	C-I40 (100)	55 1358.1 1.15 1.86 1.53 1.52 2.73 2.39 2.38 1.4 1.4 1.4
	C-I60 (100)	23 1047.0 0.68 1.00 0.79 0.79 1.82 1.57 1.57 3.0 3.1 3.1
	C-I80 (100)	16 506.4 0.60 0.88 0.65 0.65 1.48 1.29 1.29 5.3 5.5 5.5
	C-I100 (100) 14 503.3 0.59 0.69 0.64 0.64 1.27 1.14 1.14 9.0 9.4 9.3
	R-I25 (45)	20 1486.0 0.66 0.96 0.74 0.72 1.98 1.69 1.65 1.0 1.0 1.0
	R-I50 (45)	4 1879.5 0.60 0.60 0.60 0.60 1.30 1.21 1.20 4.8 4.7 4.9
	R-I100 (45)	0	-	-	-	-	-	0.63 0.59 0.58 26.6 25.2 27.8

Table 3 :

 3 Comparison of different column generation variants for the total used plate width objective

Table 4 :

 4 Comparison of hypergraph-based MIP formulations for the 2DGCSPL with different partial enumeration for the total number of used plates objective

	#opt	t, sec.
	Instances

Table 5 :

 5 Comparison of hypergraph-based MIP formulations for the 2DGCSPL with different partial enumeration for the total used plate width objective

Table 6 :

 6 Comparison of heuristics for the 2DGCSPL for the total number of used plates objective

		ea	iub	div H	div	div 32	ediv	ediv 32
	Instances	gap t gap t gap t gap t gap t gap t gap	t
	C-I20 (100) 3.47 1 2.11 1 2.09 1 2.26 1 1.73 1 1.87 1 1.64	1
	C-I40 (100) 3.64 1 1.79 4 1.59 2 1.42 2 0.55 4 0.86 3 0.37	8
	C-I60 (100) 3.68 1 1.80 13 1.13 3 1.09 4 0.39 9 0.55 8 0.21 27
	C-I80 (100) 4.14 2 1.94 32 0.99 6 0.78 6 0.19 17 0.36 19 0.07 79
	C-I100 (100) 4.47 4 1.84 65 0.75 10 0.56 10 0.15 28 0.32 36 0.02 170
	average	3.88 2 1.90 23 1.31 5 1.22 5 0.60 12 0.79 14 0.46 57
	R-I25 (45)	2.66 1 2.56 2 1.73 2 1.40 2 0.91 5 1.17 2 0.75	7
	R-I50 (45)	1.93 2 1.83 16 1.02 8 0.54 10 0.17 53 0.33 19 0.01 112
	R-I100 (45) 1.51 8 1.37 175 0.61 48 0.30 55 0.08 369 0.18 155 0.00 1224
	average	2.03 4 1.92 65 1.12 20 0.75 23 0.38 143 0.56 59 0.26 448

Table 7 :

 7 Comparison of heuristics for the 2DGCSPL for the total used plate width objective

Table 8 :

 8 Comparison of heuristics on the real-world instances with batches

Acknowledgment

Experiments presented in this paper were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux, CNRS and ANR in accordance to the programme d'investissements d'Avenir (see https://www. plafrim.fr/).