Improved Flow-based Formulations for the Skiving Stock Problem

Thanks to the rapidly advancing development of (commercial) MILP software and hardware components, pseudo-polynomial formulations have been established as a powerful tool for solving cutting and packing problems in recent years. In this paper, we focus on the one-dimensional skiving stock problem (SSP), where a given inventory of small items has to be recomposed to obtain a maximal number of larger objects, each satisfying a minimum threshold length. In the literature, diﬀerent modeling approaches for the SSP have been proposed, and the standard ﬂow-based formulation has turned out to lead to the best trade-oﬀ between eﬃciency and solution time. However, especially for instances of practically meaningful sizes, the resulting models involve very large numbers of variables and constraints, so that appropriate reduction techniques are required to decrease the numerical eﬀorts. For that reason, this paper introduces two improved ﬂow-based formulations for the skiving stock problem that are able to cope with much larger problem sizes. By means of extensive experiments, these new models are shown to possess signiﬁcantly fewer variables as well as an average better computational performance compared to the standard arcﬂow formulation.


Introduction
Modeling frameworks based on graph theory have been investigated and successfully applied in different fields of cutting and packing [8,12,23,26,29,35].Probably, one of the earliest references is given by [39] where theoretical foundations on the relationship between certain integer linear programs (ILPs) and corresponding flow formulations have been proposed.However, at that time, computers were much slower than they are today, and addressing the exact solution of discrete optimization problems was not an easy task.More precisely, the computational times required to enumerate all the variables and constraints and to solve then the complete models exactly were too large.Hence, for a long time, research on these alternative formulations was rather theoretically motivated, so that, for instance, the numerical behavior and benefits of flow-based approaches could not be properly discussed in [39]; nevertheless, some properties suggesting computational advantages of such formulations are already presented therein.
Mainly after the publication of the famous book by Nemhauser and Wolsey [31] in 1988, also computational aspects (e.g., the strength of the considered models) became a central concern in ILP modeling.In recent years, the rapidly advancing development of (commercial) Mixed ILP (MILP) software and the steady progress in terms of powerful hardware components have successively contributed to the scientific importance of pseudo-polynomial formulations for the solution of cutting and packing problems [12,30,37].In particular, a significant body of today's research is dealing with reduction techniques, as well as theoretical and algorithmic improvements [5,9,11] for the existing approaches.
In this paper, we want to promote this general idea by presenting two new and improved flowbased formulations for the one-dimensional skiving stock problem (or SSP for short).In its classical formulation, a threshold length L ∈ N and m ∈ N (different) item types that are specified by their length l i and frequency (of occurence) b i , i ∈ I := {1, . . ., m}, are considered.The SSP requires to recompose the given items in order to obtain a maximal number of objects, each having a length at least L, see also Fig. 1.
given items some possible item combinations L = 10 0 L = 10 0 L = 10 0 In [4], this problem was introduced as the dual bin packing problem (DBPP) whenever b i = 1 holds for all i ∈ I.The term skiving stock problem (SSP) has been proposed in the literature [40] for a particular practical application in the field of paper recycling [18].Nowadays, both names are usually separated according to the typology presented in [38], meaning that • the DBPP refers to highly heterogeneous instances with b i , i ∈ I, close to (or equal to) one, see [22,32], and • the SSP describes instances where i∈I b i is large compared to the parameter m, see [26,40].
The skiving stock problem plays an important role whenever an efficient and sustainable use of limited resources is intended.By way of example, such computations can be found in industrial production processes [18,40], manufacturing and inventory plannings [2,3], wireless communications [25,34], or multiprocessor scheduling problems [1].Further scientific relevance of the SSP is given by the fact that it may appear in holistic cutting-and-skiving scenarios [7,18].Observe that, although its connection to the extensively studied cutting stock problem (CSP) [12,15,19,37] is obvious, both problems are not dual formulations in the sense of mathematical optimization.
Originally, the skiving stock problem was modeled by a pattern-based approach introduced in [40].Although this model is known to possess a very tight LP relaxation [20,27,28], similiar to the CSP context [6,21], the number of variables is exponential with respect to the number of items.Consequently, pseudo-polynomial models for the SSP (most notably the arcflow model and the one-stick model ) have been established in the literature and were shown to always exhibit an equally good relaxation, as well as a reasonable computational performance for most of the considered instances [26].However, especially for instances of practically meaningful sizes, the resulting models involve (very) large numbers of variables and constraints, so that appropriate reduction techniques are required to decrease the numerical efforts for solving the corresponding ILPs.For that reason, we introduce two improved graph-theoretic ILP formulations for the SSP that are based on the ideas of reversed loss arcs and reflected arcs [11], respectively.Both of these new approaches are tailored for addressing specific drawbacks of the existing flow-based modeling framework.More precisely, the first idea mainly eliminates superfluous sinks (and possibly also further redundant vertices) of the graph, whereas the main novelty of the second formulation is to only consider half of the bin capacity so that significantly reduced sets of vertices and arcs are obtained.As a consequence of these reductions, both approaches will be shown to lead to a better numerical behavior (on average) compared to the arcflow model from [26].Moreover, dominance relations between the LP bounds obtained by the various models will be discussed.Altogether, it should be emphasized that our computational study also serves as a first systematic approach to collect and describe benchmarking instances (and their performance) for the skiving stock problem, which can form an experimental basis for future research articles.A preliminary version of this research, containing just one new model and very limited computational results, was presented at the International Conference on Operations Research 2018 (OR 2018, Brussels) as [24].
The paper is organized as follows: in the next section, we review the standard flow-based formulation for the SSP and discuss its most important reduction techniques.Afterwards, in Sect.3, two improved modeling frameworks and related theoretical properties are presented.Then, the computational performance of the three models is compared in Sect. 4. Finally, we summarize the main ideas of this paper and give an outlook on future research.

The Arcflow Model
Throughout this paper we assume the item lengths to be ordered non-increasingly, i.e., If required, we can easily obtain this order by a sorting algorithm, such as merge sort, with O(m • log m) operations.
Let E = (m, l, L, b) be an instance of the SSP.Any combination of items whose lengths sum up to at least L is called (packing) pattern of E.More formally, a pattern can be referred to as an integer vector a ∈ Z m + where the ith component states how many items of type i ∈ I are involved.Then, the pattern set is given by P (E) = a ∈ Z m + l a ≥ L .As this set (at least theoretically) contains an infinite number of elements, normally the restriction to minimal patterns (collected in the set P (E)) is sufficient.Any minimal pattern possesses the property that none of its items can be removed without underrunning the threshold L.
Based on these preliminaries, we can define the quantity v := max l a a ∈ P (E) , which represents the largest total length of a minimal pattern.
Remark 1.It is not necessary to know all the (exponentially many) minimal patterns to calculate v. Instead, as we will see later in Algorithm 1, this value is "automatically" found while generating the arcflow graph.
As a starting point, we consider the directed graph G = (V, E) given by the set of vertices V = {0, 1, . . ., v} and the set of arcs In this setting, an arc (p, q) ∈ E represents the positioning1 of an item of length q − p = l j ∈ {l 1 , . . ., l m } with its left point at vertex p ∈ V, see Fig. 2 for an example.Then, any path s = (v 0 , v 1 , . . ., v k ) in G from v 0 = 0 to some vertex v k ≥ L corresponds to a pattern a ∈ P (E).However, this basic approach still has some major drawbacks:  • There are paths not corresponding to minimal patterns, such as the sequence 0 → 5 → 7 → 9 → 12.
• There are redundant vertices that cannot be contained in paths starting at v 0 = 0 (e.g., v = 1).
Consequently, reduction techniques are required to tackle these problems: i) Obviously, only nonnegative integer linear combinations of the given item lengths have to be included in the set of vertices, since all other nodes cannot occur in a path s = (v 0 , v 1 , . . ., v k ) from v 0 = 0 to some v k ≥ L. Thus, the set V can be replaced by a reduced one (whose elements are mostly called raster points) This technique originates from the consideration of so-called normal patterns in early publications on two-dimensional cutting problems [10,17], and refers to the observation [36,Criterion 3].Even nowadays the theory of potential allocation points is a quite active field of research in cutting and packing [9,14].
ii) In many cases, the previous reduction is not sufficient for vertices v ≥ L, because the minimality of the corresponding item combinations is not ensured.By way of example, the instance E = (2, (7, 2), 10, (7, 7)) certainly leads to 12 ∈ V , but there is no minimal pattern a ∈ P (E) with l a = 12.For that reason, let be the set of all possible total lengths of minimal patterns.Then, we can replace V by These observations also lead to a reduced set of arcs If we want to assign a path s in G = (V , E ) to a minimal pattern a ∈ P (E), this identification is not necessarily unique, because a ∈ Z m + does not carry any information about the particular item order.For that reason, the consideration of monotonically decreasing paths (i.e., paths whose corresponding item lengths are sorted in non-increasing order) is proposed in the literature, see [33,Subsection 4.8.4.] or [35] for an exemplary application to the related CSP.Formally, this can be done by defining an index µ(p) ∈ I ∪ {0} for all p ∈ V \ L by Remark 2. Practically, this additional constraint can be incorporated by processing the item types in the order specified by assumption (1) within the graph generation, see also Algorithm 1.
This leads to a new set of arcs and a reduced graph G = (V , E ), see Fig. 3 for an example.Remark 3. The introduction of µ may not solve the problem of equivalent paths within G entirely.Consider, by way of example, the instance E 0 where G (illustrated in Fig. 3) contains the two paths 0 → 3 → 5 → 10 and 0 → 5 → 8 → 10 which both refer to the pattern a = (1, 1, 1) .Moreover, the first path does not respect the monotonicity requirement.However, this problem cannot be avoided if, for instance, a (large) item length can be obtained by summing up two (or more) smaller item lengths.
Similar to [9, Algorithm 1], where a possible construction of the vertex set for the CSP scenario is described, the reduced arcflow graph can be generated in O(mL) by the procedure given in Algorithm 1.
for j = 1 to b i do end if end for 18: end for 19: return V , E In order to formulate the arcflow model in a preferably convenient way, we introduce the auxiliary index sets for every q ∈ V and the set for every i ∈ I.
Remark 4. Note that A − (q) and A + (q) model the predecessors and successors of vertex q ∈ V , respectively, whereas E(i) collects all arcs referring to an item of length l i .Moreover, the cases A − (q) = ∅ and A + (q) = ∅ can occur for some q ∈ V .
Let x pq ∈ Z + denote the flow along arc (p, q) ∈ E , then we can state the Arcflow Model of the Skiving Stock Problem x 0q → max s.t.
Constraints (3) can be interpreted as a flow conservation: at every interior vertex, which is not a possible endpoint, the units of incoming flow (i.e., the number of items "terminating" at this position) have to equal the units of emanating flow (that is the number of items "starting" at this position).Thus, the objective function maximizes the total flow within G , that is, the total number of objects obtained.Conditions (4) ensure that the given item supply is not exceeded.In this arcflow model, the numbers of variables and constraints are O(mL) and O(m + L), respectively.
Although this model has turned out to be competitive in solving randomly generated instances [26], there are still some drawbacks connected to different aspects of this formulation: (a) The size of the vertex set V (and thus the size of the graph G ) is strongly determined by the quantity v, which is not known until the whole arcflow graph has been generated.Moreover, G normally consists of multiple sinks (i.e., many different endpoints of patterns) whenever v > L holds.
(b) As indicated by Remark 3, in some cases, the graph still contains some symmetries or non-monotonous paths.
(c) The reduction caused by raster points is particularly successful at the beginning of the graph (i.e., in the first half of V) whereas it tends to be nearly irrelevant in the second part.
Whereas the second problem may not be addressed in the general context (see Remark 3), we aim at tackling the remaining unfavorable properties by two new flow-based approaches in the following section.

An Arcflow Model with Reversed Loss Arcs
As a starting point, we want to deal with problem (a) mentioned above.To this end, we assume that L ∈ V holds without loss of generality.If this is not the case, then the given instance E = (m, l, L, b) can be replaced by an equivalent instance E = (m, l, L , b) with L ≥ L + 1.Now, instead of allowing multiple endpoints within the graph, we force any path to end at vertex v = L.This can be obtained by the following three steps: (1) Any arc (p, q) ∈ E with q = L + r for some r ≥ 1 is shifted back to (p − r, L).After having completely executed this step, we let v denote the lowest-indexed vertex that is the tail of a shifted arc.
(2) Add reversed loss arcs (e, d) with v ≤ d < e < L and e = succ V (d), where succ M (m) describes the successor of element m ∈ M within the canonically ordered set M ⊆ N.
(3) Now, any vertex that is not part anymore of at least one arc is deleted.
Let us refer to the resulting graph as G L := (V L , E L ), then principally the same optimization model as in Sect. 2 can be applied (now, of course, based on the updated sets V L and E L ), if the artificial loss arcs are also respected in the flow conservation conditions.We will briefly refer to this model as the loss arcflow model (or simply larcflow ), and define by z L its objective function.
Example 1.The modified graph for the instance E 0 = (3, (5, 3, 2), 10, (3,4,4)) is depicted in Fig. 4.  Here, the formerly existing arcs (9,12), (8,11), (9,11) ∈ E have been shifted back to the arcs (7, 10), (7,10), (8,10) ∈ E L , and two additional reversed loss arcs have been introduced, namely (8, 7) and (9,8).In this new setting, the vertices 11 and 12 do not appear anymore; and, by way of example, the path 0 → 3 → 6 → 9 → 12 is now represented by 0 → 3 → 6 → 9 Besides having eliminated the additional sinks of the graph (so that any path now terminates at v = L), there are some further advantageous, neutral and disadvantageous properties (marked by (+), (±) and (−)) that we want to mention here as a concluding summary of this subsection: (+) In our example, see Fig. 4, the new arcflow graph G L contains the same number of arcs as before since the number of additional loss arcs is equal to the savings that are obtained by shifting arcs back (onto some already existing arcs).However, for practically meaningful problem sizes, this effect is not likely to happen so that, normally, a smaller number of variables can be expected (see also Sect.4).
(±) Typically, the loss arcflow model has (roughly) the same number of constraints as the standard arcflow formulation because the number of interior vertices is not likely to change considerably.However, in some (rather exceptional) cases, it may happen that interior vertices of V can be deleted after having shifted the arcs according to Step (1) from above.On the other hand, for instances with sparse vertex set, it is also possible that the shifting of arcs can activate additional interior vertices.As we will see in our computational part (see Sect. 4), both of these phenomena do not play a very important role for larger instances, in general.
(−) Possibly, the shifting of some arcs can produce additional symmetries or can destroy some monotonicity properties of the graph.For instance, the graph G L depicted in Fig. 4 now contains a further possibility to represent the pattern a = (1, 1, 1) .Moreover, at least from a theoretical point of view, arbitrarily long paths (patterns) can be obtained in G L because the graph does not have to be acyclic anymore, see also Fig. 4.

The Reflect Arcflow Model
In this subsection, the problems (a) and (c) from the list presented in Sect. 2 shall be dealt with collectively.For that reason, we consider an approach basically introduced in [11] (for the CSP scenario), whose key idea is to reduce the set of vertices by decomposing any path into two subpaths (providing roughly half of the desired total length) that are linked by a reflected arc.
Example 2. Let us again consider the instance E 0 = (3, (5, 3, 2), 10, (3,4,4)) and the corresponding graph G depicted in Fig. 3.By way of example, the description of the path 0 → 3 → 6 → 9 → 12 actively uses five vertices and four different arcs.In a new setting, that only takes the vertices of the first half (of this path) into account, we could also describe this path as illustrated in Fig. 5.
Figure 5: A preliminary schema to represent paths on only one half of the vertex set.The third index (s or r) of an arc is used to distinguish between standard arcs and reflected arcs, where in the latter case also the coordinate of the reflection point is mentioned as an auxiliary information.
Here, the two (identical) subpaths (0, 3, s) → (3, 6, s) (each referring to two items of length l 2 = 3) are connected by a reflected arc (6, 6, r = 6) (where the reflection point r = 6 is given by half of the previous total path length of l a = 12).Hence, we obtain an equivalent description of the pattern a = (0, 4, 0) by only using three vertices and three different arcs.
As skiving stock patterns can exhibit several lengths L t ∈ L := l a a ∈ P (E) (which are identical to the sinks of G ), we usually would have to cope with multiple reflection points r t = L t /2 if the idea presented in [11] and the previous example is simply overtaken.Moreover, the elements of L are not known in advance and, hence, cannot be used to efficiently implement a reflect graph.Due to these difficulties, we propose another strategy to represent a pattern in the SSP context.
In order to avoid multiple reflection points, we again try to force any path (or pattern) to have exactly the length L by introducing appropriate loss arcs (but now in forward direction).Then, Algorithm 2 can be applied to build the reflect graph G R := (U, A) (for the SSP) with only one reflection point2 r := L/2.
More precisely, we have that any arc (u, v) ∈ E (appearing in Sec. 2) After having executed these steps, let v denote the lowest-indexed vertex that is the head of a reflected arc.Then, we add the loss arcs (d, e, l) with d ∈ U, v ≤ d < L/2 and e = succ U (d) acting as the direct successor of d.Finally, we have to introduce a special reflected arc (not referring to any item length) (L/2, L/2, r) to enable the combination of two subpaths providing exactly half of the bin capacity.
Example 3. The reflect arcflow graph G R for the instance E 0 = (3, (5, 3, 2), 10, (3,4,4)) is depicted in Fig. 6.For the sake of a better comprehension, the colors of the arcs are now referring to the different arc types (standard, reflected and loss).Here, the path 0 → 3 → 6 → 9 → 12 from  for k = L/2 − 1 down to 0 do 10: if else if the arc is reflected end if

22:
end for

23:
end for 24: end for (4, 4, r) (5, 5, r)  , so that (at the moment) also non-minimal patterns can be discovered in the graph.However, this possibility will later be excluded by adapted flow conservation constraints.
Observe that, even for this quite small example, the numbers of vertices and arcs could be reduced by roughly 50 percent compared to the standard arcflow model (from Fig. 3) or the model with reversed loss arcs (from Fig. 4).As mentioned in Sect.2, these significant savings are obtained because the effect of normal patterns (or raster points) is noticed almost only at the beginning of the graph, whereas it tends to be irrelevant for the second half of the vertices (which is not modeled explicitly anymore).

Let us define
(0,e,s)∈A (d,e,l)∈A The objective function maximizes the total flow on the reflected arcs.As (on balance) any unit of flow on a reflected arc is responsible for linking two subpaths to one single pattern, this quantity actually displays the number of obtained objects (with length ≥ L).For this reason, we further have to demand that the total flow in the system is twice the objective value, see constraint (7).Moreover, conditions (6) require each item to be used at most b i times, whereas constraints (8) model the (adapted) flow conservation at each interior vertex: note that (in a pattern sense) each flow entering a node through a standard arc has to be continued by a flow leaving the node through a standard or a reflected arc (classical flow conservation), or by a flow entering the node through a reflected arc (connection between two subpaths), including the additional flow brought and removed by the loss arcs.
Example 4. The conditions explained so far would entail the following properties: i) Fortunately, the non-minimal pattern a = (0, 4, 1) from the previous example cannot be represented by a feasible path (consisting of two copies of 0 → 3 reflect → 4 linked by (4, 4, r) ∈ A) in G R anymore.To this end, let us consider the flow propagation for the arcs and vertices depicted in Fig. 7.For the first interior node e = 3, we have two units of flow on the incoming standard arc (0, 3, s) and two units of flow on the leaving reflected arc (3, 4, r), so that condition (8) is satisfied.However, when looking at e = 4 there is no flow corresponding to the lefthand side of (8), whereas on the right-hand side we have two units of flow for the incoming reflected arc (3, 4, r) and an additional nonnegative flow on the reflected arc (4, 4, r).Hence, the flow conservation does not hold for the vertex e = 4, so that the non-minimal pattern cannot be part of a feasible flow (in G R ) anymore.
ii) Contrary to the first point, at the moment, it would still be possible to link the two identical subpaths 0 → 2 → 4 loss → 5 by means of the reflected arc (5, 5, r) without violating constraints (8) for any of the involved interior vertices.As this would refer to a vector a = (0, 0, 4) (which does not represent a pattern) some further conditions are needed to limit the use of loss arcs.
Based on the second observation given in the previous example, loss arcs may only be used for special purposes.More precisely, they either have to directly follow another loss arc or a reflected arc, see (9).Then, by way of example, the vector a = (0, 0, 4) cannot be obtained anymore as a feasible connection of two subpaths.
One particularity of this model is that we have to allow negative flows on the (artificial) reflected arc (L/2, L/2, r), see (11), to also enable two reflected paths to be merged together without violating the flow conservation and without miscounting the number of obtained patterns in the objective function.
Example 5.As long as a nonnegative flow on the artificial reflected arc (L/2, L/2, r) ∈ A is demanded, the flow conservation conditions for e = 5 in the scenario depicted in Fig. 8 would lead to the same contradiction as in the previous example.Moreover, as already two units of flow are required on the reflected arc (3, 4, r) ∈ A, this pattern would not be counted correctly (as one single pattern) in the objective function.To solve these two problems, the flow variable corresponding to (L/2, L/2, r) ∈ A is allowed to possess negative values.Then, with ξ (0,3,s) = ξ (3,4,r) = ξ (4,5,l) = 2 and ξ (5,5,r) = −1 all constraints are satisfied and the total flow on all the involved reflected arcs is equal to one, so that exactly one pattern is added to the objective value.Example 6.Another typical case where linking two reflected paths can sometimes be meaningful occurs if very large items may appear together in a single pattern, see Fig. 9 for an illustration of the exemplary instance E 1 = (3, (18,16,8), 20, (10, 10, 10)).Observe that, also in this setting, ξ (L/2,L/2,r) will become negative if, for instance, the two identical paths 0   We also note that the reduction criteria proposed in [11, Proposition 5] stating that any feasible pattern can be represented by a pair of colliding paths whose reflected arc (d, e, r) has d < e cannot be used for the SSP.Indeed, let us consider the pattern a = (0, 0, 3) from the previous example: it can only be obtained by the two colliding paths 0 → 8 reflect → 4 loss → 8 and 0 → 8, in which the reflected arc (8, 4, r) has to be used.
In summary, modeling and generating the reflect arcflow formulation for the SSP is more challenging than in the cutting scenario [11] and requires (partly) independent techniques as well as considerable modifications of the corresponding ILP.
Altogether, this model possesses much fewer variables than the previous flow-based formulations (in Sect. 2 and Subsect.3.1).As regards the number of constraints, no clear prediction can be made.On the one hand, our new approach allows for using much fewer flow conservation constraints (in most cases only roughly 50% of the original number) of type (8), whereas, on the other hand, an additional constraint (see ( 9)) has to be included for any vertex.Depending on the particular instance, these two effects may lead to slight differences (in both directions) if compared to the number of constraints of the standard arcflow model, see also Sect. 4.

Relations between the LP bounds
Whenever there are different modeling approaches for one and the same optimization problem, an important question deals with the strength of the bounds obtained by the respective LP relaxations.It is well known that the tightness of a relaxation is a crucial factor in the size of branch-and-bound trees, because it significantly influences the efforts to solve the ILP.To this end, we will investigate the relationships between the optimal objective values of the three LP relaxations z LP A := z LP A (E), z LP L := z LP L (E), and z LP R := z LP R (E) and discuss possible dominance relations.
A first important observation is given by: Proof.For the given instance, we have z LP R = 1.5, z LP A = 1.6, and z LP L = 5/3, obtained by the patterns (paths): The main reason for the occuring differences is given by the fact that some patterns (paths) can be feasible in one graph, while they cannot be composed in another graph.In the example presented in the previous theorem, the patterns a = (1, 0, 0, 2) and a = (0, 1, 0, 3) cannot appear in the reflect graph because it is impossible to start a subpath with more than one copy of the item length l 4 = 2. On the other hand, the pattern a = (0, 2, 0, 0) can for instance appear in the reflect framework, whereas it cannot be found in the standard arcflow formulation.
Given the information of Theorem 5, in the following result we clarify that there is no dominance relation between the reflect model on the one hand and either the standard arcflow model or the loss arcflow model on the other hand.Proof.Indeed, for the given instance we have z LP R = 1.5 and z LP A = z LP L = 1.4,obtained by the patterns (paths): × (0, 5) .
Again, the reason for the differing optimal values is based on different sets of feasible patterns (paths).For instance, the pattern a = (2, 0) used in the solution of the reflect model does not appear in the other two graphs due to the construction presented in Algorithm 1, where the for-loop for i = 1 (with 1 ≤ j ≤ b 1 implying j = 1) only generates one single arc (in the whole graph) corresponding to this item.Of course, also the reflect graph does only contain this arc once, but here two identical copies of this subpath (0 → 5) can be linked to obtain a feasible pattern (that can then be used 0.5 times).
Based on the two instances given in the previous theorems, only the relationship between the optimal values z LP A and z LP L is not fully discussed yet.
To answer this open question, we present the following concluding result.Proof.Consider a given instance E = (m, l, L, b) and a corresponding feasible solution x = (x pq ) (p,q)∈E of the LP relaxation of the standard arcflow model.Since the directed graph G of that instance is acyclic, any flow (i.e., any feasible solution) decomposes into feasible paths γ 1 , . . ., γ s from v 0 = 0 to some v ≥ L with, in general, fractional flow values f 1 , . . ., f s ≥ 0 (being constant along the arcs of a fixed path).Now note that, by construction, any path γ in arcflow has a corresponding path γ in loss arcflow, where all the arcs except the last one(s) are the same: • If the last arc (p, q) of γ ends after L, then loss arcs have to be added from the penultimate node p of γ to the node L − (q − p) in loss arcflow.These loss arcs exist since loss arcs are created (between all active nodes) from L − 1 to L − max{l i : i ∈ I}.The path γ is then completed by an arc (L − (q − p), L) which exists as any arc ending after L in arcflow for item i ∈ I is transposed to (L − l i , L) in loss arcflow.
• If the last arc of γ ends in L, then the paths γ and γ are exactly the same.
It is straightforward to show that these transformed paths γ 1 , . . ., γ s (with the same flow values f 1 , . . ., f s ≥ 0 as before) also satisfy the constraints of loss arcflow: • As we did not touch arcs starting from v 0 = 0, the objective function (and also the objective value) remains the same.
• Any of the transformed paths γ 1 , . . ., γ s satisfies the flow conservation constraints by construction; so the sum of all these paths will also satisfy it.
• The flow on the possibly new arcs of type (L − l i , L) in loss arcflow equals the flow on the arcs (d, d + l i ) with d + l i ≥ L in arcflow.Consequently, also the item limitations are obeyed.
• All flow variables are still nonnegative.
Altogether, we obtain a feasible solution x of the LP relaxation of loss arcflow with the same objective value.
Summarizing the main observations of this subsection, we can finally point out that arcflow dominates loss arcflow (in terms of the LP bound), whereas there are no further dominance relations among the three optimal objective values z LP A , z LP L , and z LP R .

Computational Experiments
In addition to the theoretical presentation that we provided for the new mathematical models of the SSP, their computational performance and average efficiency shall be discussed based on numerical experiments.For that purpose, three general categories of instances are considered: (A) datasets designed for the SSP (or dual bin packing problem) from the literature, (B) randomly generated instances, (C) benchmarking instances originally introduced for the CSP (or for the bin packing problem) 3 .
Based on the fact that the SSP is a relatively young and rather unknown optimization problem, only a few families of instances of the first type can be found in the relevant literature, whereas much more test sets are available for the CSP.As both problems share the same input data (but with a partly different meaning), we decided to also include these CSP benchmarks into our experimental study in order to provide a wider variety of data sets.It is not guaranteed that these typically challenging instances for the CSP will be difficult to solve in the SSP context as well; but at least some of them (as, for instance, the AI/ANI sets introduced in [12, Subsection 7.1.3])will definitly preserve their hardness so that the consideration of CSP benchmarks is justified also from this point of view.Altogether, besides purely providing information on the practical behavior of our different models, this section also serves as a first systematic approach to collect and describe benchmarking instances for the SSP that can form an experimental basis for future research articles.

Data Sets
In what follows we aim at describing in more detail the three instance categories mentioned above.
(A) As regards datasets particularly designed for the SSP, the literature only offers a very few of them.Moreover, given the year (and the computational limits at that time) when they were published, some of these instances (e.g., those presented in [40, Table 3] with L = 100 and m ≤ 10) nowadays have to be considered as too easy to notice any substantial differences between the various approaches.To the best of our knowledge, the only considerable systematic description of (larger) SSP benchmarking instances is contained in [32], where also some test sets previously appearing in [22] are proposed.Contrary to our definitions, these constructions mostly fix the total number n = i∈I b i (instead of the number m of different item types) as an input parameter of a test class.Theoretically, given the fact that the numbers of variables and constraints (and, hence, the size of the ILP model) rather depend on m and L, this could lead to quite heterogeneous instances within one and the same test class.However, in our experiments we always have b i , i ∈ I, equal to (or at least very close to) one so that n ≈ m holds, and this leads to instance sets possessing comparably complex ILP models for fixed input parameters n and L.
Unfortunately, the original instances from [32] are not available anymore, so we had to build our own instances based on the constructions described in that paper.More precisely, this first category consists of two test sets with the following input parameters: (A1) These rather small instances were proposed in [22] and are described by l max = 99 and the following integers: n ∈ {10, 20, 40, 60, 80, 100}, L ∈ {100, 120, 150, 200, 300, 400, 500}, For any of these 126 combinations, we randomly generated 10 instances, each with uniformly distributed item lengths l i ∈ [l min , l max ], i ∈ I.Note that the values of b i do not have to be specified because any of the n items is drawn individually.
(B) Randomly generated instances for the SSP have previously been addressed in [26].As all these problems were solved in reasonably short computation times, we intend to consider similar, but much larger instances in this second category.More precisely, any class is parametrized by a pair (m, L) in the following way: m ∈ {200, 300, 400, 500}, L ∈ {10000, 20000, 30000, 50000}.
In accordance with [26], for any of the 16 combinations we randomly generated 10 instances, each with uniformly distributed item lengths (C) The third category of test sets consists of 1895 challenging instances that were originally invented for the CSP (or for the bin packing problem) over the past years.All these instances are well known benchmarks whose specifications are described in the relevant literature.More precisely, as similarly pinpointed in [11,Section 7.1], this category contains the three subclasses: (C1) Classical: A set of 1615 instances presented in various articles in the last decades and having highly different characteristics.A complete description of these data sets can be found in [12].
(C2) GI: A total number of 80 instances (divided into four classes of 20 instances each) proposed in [16] partly involving extremely large threshold lengths L ≤ 1500000.
(C3) AI/ANI: Two sets containing 100 challenging instances each that have been proposed in [12].Here, the word "challenging" either means that already finding reasonable candidates for an optimal solution may be difficult (as in the AI case), or that (as regards the ANI instances) proving the optimality of a given solution becomes really hard since the additive integrality gap (measured with respect to the LP bound) is at least one.
Instances in set C have been gathered together in the BPPLIB [13] and can be found on the Internet by http://or.dei.unibo.it/library/bpplib.Instances from sets A and B are available from the authors upon request.

Methodology
All our experiments were executed on an AMD A10-5800K CPU with 16 GB RAM, using Gurobi 8.0.1 (imposed to run on a single thread) as MILP solver with a time limit of one hour.The model generation itself was performed in Python 3.6.4and the Gurobi Python module.For any instance and for any formulation, we collected the following data: • opt: an indicator stating whether the instance could be tackled successfully (opt = 1) or not (opt = 0) within the given time limit, • z , z c : optimal value of the ILP and the LP relaxation, respectively (if solvable), • t, t LP : total time (in seconds) needed to solve the ILP and the LP relaxation, respectively, • n var , n con : number of variables and constraints appearing in the respective ILP, • n nz : number of nonzero entries in the constraint matrix of the respective ILP, • LB, U B: best lower (upper) bound obtained for the optimal objective value.
Note that, given the large number of instances attempted, we will only refer to aggregate (in terms of opt) or averaged values for each considered class of instances instead of providing the individual results for any single example.Whenever an instance could not be solved to optimality, the time limit (one hour) will be used as the computation time t.
In addition to the comparison of the three (pure) ILP formulations, we also measure the effects of passing an initial feasible solution x heu to the Gurobi solver.This feasible solution was always obtained by a heuristic proposed in [32] whose iterations can be described as follows: Choose the largest item type (say k) that is still available (i.e., with residual b k > 0) and allocate a kj = min{ (L − 1)/l k , b k } items of type k to the current bin j.
(I) If a kj < b k holds, then we choose the smallest possible item type i ∈ I with b i > 0 so that increasing a ij by one unit leads to a feasible pattern.(This step can always be performed by choosing i = k.) (II) If a kj = b k holds, then we continue with item type k + 1 and try to assign as many items of that type to bin j without reaching the threshold L. Depending on the situation, then we have to apply step (I) or (II) with respect to the item index k + 1.
This general procedure is repeated (while successively updating the values b i , i ∈ I), until all items have been assigned or no feasible bin can be found anymore.
Based on the collected data, we will evaluate and compare the computational behavior of six solution strategies (three ILP formulations with or without the additional application of the heuristic) in this article.Note that, in [26], the original arcflow model turned out to provide the best computational performance for practically meaningful problem sizes (among the tested formulations) although it possesses slightly more variables and constraints compared to the onestick formulation, see [26,Table 6].Consequently, the conventional arcflow formulation (introduced in Sect.2) can be considered as the state-of-the-art solution strategy for skiving stock problems.Hence, our comparison also involves the most promising existing ILP formulation, so that a reasonable investigation of the various solution approaches is conducted.

Results and Discussion
In this subsection, we aim at presenting some of the results obtained by the computations described above.Note that, for the sake of simplicity, the arcflow model of Sect.2, the loss arcflow model of Subsect.3.1 and the reflect model of Subsect.3.2 will mostly be referred to by their respective abbreviations "arcflow", "larcflow" and "reflect".Moreover, for the tables, we always use the short form "heu" instead of heuristic.
At first, we consider the instances of category (A) and state that the subset (A1) could be solved to optimality by any of the six variants within really short times, see Tab. 1.Hence, these instances are rather small to observe significant time differences among the considered formulations.However, having a look at especially the numbers of variables and nonzeros, we can notice a considerably less complex structure of the ILP model for reflect (having roughly a third of the variables and half of the nonzero matrix elements compared to arcflow) and larcflow.Moreover, slight differences for the various LP relaxation values could be observed underlining the (theoretical) fact that the considered formulations are not equivalent if continuous variables are considered.To obtain a more precise picture of the six solution variants, we move forward to the harder instances contained in category (A2).At first, we summarize the averaged (or aggregated, in terms of opt) data in the same way as for the subset (A1), see Tab. 2. Here, we can see that both alternative formulations can solve more instances (out of a total number of 1050) than arcflow while also being faster on average (despite having a slightly worse LP bound).Again, this behavior is mainly due to significant savings related to the ILP model's complexity.Moreover, a very interesting observation is given by the usefulness of the heuristic for any of the three different approaches: both the number of instances solved to optimality and the average time needed to solve them become better if a starting point is passed to the ILP solver.In Tab. 2 we are averaging over a wide variety of very heterogeneous instances.To highlight possible computational insights, we decided to perform more detailed analyses.To this end, Tab.3-5 report on the values opt and t for one input parameter being fixed.Obviously, in almost all of the experiments and for almost any modeling approach, the application of the heuristic leads to better results with respect to the considered data.We also observe that the difficulty of the problem increases as: (i) the size of the smallest item decreases, (ii) the threshold length increases, and (iii) the number of items increases.At this position, we refer the interested reader to the appendix (Tab.A.15-A.17) to have a look at an even more detailed presentation of the results for the A2-instances solved by any of the six formulations.In any of the experiments related to the subset (A2), the heuristic could be performed in much less than one second.Consequently, in addition to the advantageous effects observed earlier, the heuristic should always be applied to obtain a feasible starting point for the ILP solver.By way of example, Fig. 10 compares the solution times for the arcflow model with and without applying the heuristic.It can clearly be seen that there are many instances possessing roughly the same solution time in both cases, but there are also many instances whose solution time decreases drastically (partly from the time limit, indicating that there was no solution found previously) if the heuristic is applied prior to the optimization.Moreover, note that the opposite effect (i.e., a significant increase of the computation time) is also possible, but happens very rarely.
As a summary, we would like to further underline the positive effects of the heuristic by the data collected in Tab.6-7.For any combination (F 1 , F 2 ) of the six solution variants, we are counting the number of A2-instances: • that are solved faster by formulation F 1 than by F 2 , see Tab. 6. (Here, the word "faster" means that a time difference of at least 0.5 seconds could be observed.) • that are solved by formulation F 1 , but not by F 2 , see Tab. 7.  In particular, if considering a fixed model type, than the entry at position (F 1 , F 2 ) is always better than that at (F 2 , F 1 ) if F 1 applies the heuristic and F 2 does not.Moreover, we can again state that the reflect formulation is superior to the other two approaches.This observation can also be verified by the data contained in Tab. 8. There, for any of the six formulations we counted the number of A2-instances where the respective approach possessed the fastest solution time.Again, being faster than another formulation is connected with a time difference of at least 0.5 seconds.Altogether, given the reasons presented (based on extensive computations) so far, for the remaining two categories of instances, (B) and (C), we will later only consider the three different models with the application of the heuristic.
One last thing we would like to report for the A2-instances deals with the strength of the respective LP relaxations.As to be clearly seen in Tab. 9, apart from the theoretically proved statement that arcflow dominates the loss arcflow model, any relation is possible among the other pairs of optimal values.This also emphasizes the fact that differences between the LP relaxation values do not only appear for the very simple instances particularly constructed for the proofs within the previous section, but also for larger instances of practically meaningful complexity.Now, let us consider the instances from subset (B), see Tab. 10.Note that, although they are also randomly generated, they are completely different from the previous set (A) because the values b i are usually much larger than one, so that typical SSP instances are built.Given the extensive discussion for the set (A) and the fact that randomly generated instances have already been partly adressed in the relevant literature [26], we only consider the average solution times and the numbers of solved instances per feasible combination of input data.Contrary to the A-instances, the model with loss arcs is usually not better than the standard arcflow model, meaning that, in many cases, they possess a comparable performance, while in a few scenarios arcflow is also performing strictly better.A first explanation for this observation is given by the high degree of randomization within the Gurobi solution procedure.Hence, even if the same solver is called in both situations, different sub-routines or heuristics can appear in different orders, and this may influence the total solution efforts.On the other hand, even though the loss arcflow formulation basically has fewer variables (and roughly the same number of constraints), from our point of view, it does not counter balance the fact that its feasible integer solutions have a more complex structure that may be harder to deal with, compared to the standard flow propagation of arcflow, for a general ILP solver.Moreover, note that the B-instances are typically more complex and harder (on average) than those from the subset (A), and this especially leads to larger branching trees with much more nodes to visit.This is also due to the fact that the starting point provided by the heuristic leads to a lower bound whose absolute difference to the true optimal value4 is much larger than in the previous settings.Altogether, under these conditions, finding better integer solutions becomes more and more important, so that the arcflow model can show advantages compared to the loss arcflow formulation.Similar argumentations can be applied to further explain the very few cases where the arcflow model is beating the reflect formulation: also here, a good integer solution found quickly can lead to significant improvements in dealing with the branch-and-bound tree.However, it should be stated that also for this kind of instances the reflect formulation is (on average) the best one, with significantly lower computation times and a larger total number of instances solved to optimality.
For the sake of completeness, as a last experiment we are also considering benchmark instances that have been presented for the related bin packing problem (or for the CSP) over the last couple of decades.Note that an instance which is hard in the bin packing case is not guaranteed to preserve its difficulty if interpreted as a dual bin packing instance, and vice versa.However, at least for some of the hard instances from the C-part, it is clear that they will remain hard in the SSP scenario.By way of example, let us particularly mention the ANI-instances, whose integrality gaps will also be at least one if interpreted as an SSP instance.The precise results can be found in Tab.11.Note that, for some of the very large instances proposed by [16], we observed memory issues (i.e., the branching trees became too large), so that no data could be obtained.These situations are indicated by "-" in the table.Also here, the impression we got from the previous instance sets is again emphasized.For any subcategory (C1), (C2), and (C3), the reflect formulation is able to cope with the most instances, and possesses the lowest solution time.Again, the behavior between arcflow and loss arcflow is changing depending on the particular problem set: for instances where the main difficulty is the model size (e.g., Scholl 3, GI AA125 and GI BA125), loss arcflow seems better, while for instances instances where the As a summary of all instance types, we would like to report about the average relations between the numbers of variables, constraints and nonzeros for the three models in Tab.12.Both alternative formulations always possess much fewer variables and nonzeros compared to the standard arcflow model5 .As regards the number of constraints, they are always roughly similar in any case, but we also observe that tiny differences (in both directions) compared to the standard arcflow model are possible.Altogether, note that in any of the three categories we observed instances that cannot be solved (in one hour and with the given computational environment) even by the most promising of the considered formulations.Consequently, our experiments do not only cover a significant set of instances, but also clearly point out the limits of what can currently be solved within a reasonable amount of time.We better highlight this fact in Tab. 13, where we show the behavior of the reflect formulation for the very hard instances of the categories (C2) and (C3).
Finally, we would like to refer the interested reader to Tab.A.14 in our appendix, where the quality of the heuristic has been evaluated also for the sets (B) and (C).

Conclusions
In this paper, we investigated an existing and two new arcflow formulations for the skiving stock problem.The main idea of the two new approaches is given by avoiding superfluous nodes that are larger than the capacity L or by only considering the first half of the vertex set and to model a pattern as the sum of two subpaths that are connected by a reflected arc.For these three formulations, the quality of the continuous relaxation has been investigated from a theoretical point of view.Moreover, based on extensive numerical tests with a wide variety of differently characterized instances, we have shown that especially the new reflect model possesses significantly fewer variables and much better solution times (thus also resulting to a larger number of instances solved to optimality) compared to the original formulation.Consequently, the reflect arcflow model may be seen as a powerful tool for solving (large) instances of the SSP in reasonably short time.As future research, we envisage the use of these enhanced arcflow models to other relevant cutting and packing problems.

Theorem 7 .
For any instance E = (m, l, L, b) we have z LP A ≤ z LP L .

Figure 10 :
Figure 10: Comparison of the ILP solution times for the A2-instances and the arcflow model: for each instance, a is drawn at position (x, y), where the x-and y-values contain the times with and without heuristic, respectively.The red line displays the function f (x) = x.

Table 1 :
Computational results for the 1260 A1-instances: opt presents an aggregated number whereas all other data are average values

Table 2 :
Computational results for the 1050 A2-instances: opt presents an aggregated number whereas all other data are average values

Table 3 :
Averaged computation times and total numbers (in parentheses) of solved A2-instances for different values of l min

Table 4 :
Averaged computation times and total numbers (in parentheses) of solved A2-instances for different

Table 6 :
Comparison between the solution times for the A2-instances: an entry x at position (F 1 , F 2 ) gives the number of instances where formulation F 1 was faster than formulation F 2 (i.e., where a time difference of at least 0.5 seconds can be observed)

Table 7 :
Comparison between the A2-instances solved to optimality: an entry x at position (F 1 , F 2 ) states the number of instances which were solved by formulation F 1 but not by F 2

Table 8 :
Number of A2-instances where the respective model was the fastest of the six solution approaches

Table 9 :
Comparison of the LP relaxation values of the three formulations: an entry x at position (F 1 , F 2 ) gives the number of A2-instances with z c (F 1 ) ≤ z c (F 2 ) (measured with a tolerance of ε = 10 −8 )

Table 10 :
Results for the 160 B-instances (average solution times for the ILP and numbers of solved instances) always with application of heuristic

Table 11 :
Results for 1855 (out of 1895) C-instances (average solution times for the ILP and numbers of solved instances)

Table 12 :
Relative numbers of variables, constraints and nonzeros averaged over all instances of the respective category (where the value of the arcflow model is normalized to 1.00)

Table 13 :
Some further results for very hard CSP-bechmark instances also showing the limits of the reflect formulation

Table A .
14:Relative gap ∆r = (z − z heu )/z between heuristic value and optimal value for all instance categories.If optimality was not proved, z ≈ U B with the upper bound U B (from the reflect model) found by Gurobi was taken instead.

Table A .
15: Arcflow Results for A2-instances (average solution times for the ILP and number of solved instances)

Table A .
16: Loss-Arcflow Results for A2-instances (average solution times for the ILP and number of solved instances)