{'fraction_non_alphanumeric': 0.045941098188971824, 'fraction_numerical': 0.014754382919632008, 'mean_word_length': 4.108475578216212, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The fragility-constrained vehicle routing problem with time windows

We study a new variant of the well-studied Vehicle Routing Problem with Time Windows (VRPTW), called the fragility-constrained VRPTW, which assumes that 1) the capacity of a vehicle is organized in multiple identical stacks; 2) all items picked up at a customer are either"fragile"or not; 3) no non-fragile items can be put on top of a fragile item (the fragility constraint) ; and 4) no en-route load rearrangement is possible. We first characterize the feasibility of a route with respect to this fragility constraint. Then, to solve this new problem, we develop an exact branch-price-and-cut (BPC) algorithm that includes a labeling algorithm exploiting this feasibility characterization to efficiently generate feasible routes. This algorithm is benchmarked against another BPC algorithm that deals with the fragility constraint in the column generation master problem through infeasible path cuts. Our computational results show that the former BPC algorithm clearly outperforms the latter in terms of computational time and that the fragility constraint has a greater impact on the optimal solution cost (compared to that of the VRPTW) when vehicle capacity decreases, stack height increases and for a more balance mix of customers with fragile and non-fragile items.


Introduction
This paper introduces a new variant of the classical Vehicle Routing Problem with Time Windows (VRPTW, see Desaulniers et al. 2014) where a particular fragility loading constraint is enforced. This variant is called the Fragility-constrained VRPTW (F-VRPTW). The classical VRPTW considers a set of customers, each specifying a demand volume to be picked up and a pickup time window. Pickups are performed by a fleet of homogeneous vehicles with given capacity, initially located at a depot. Travel times and costs among customers and between each customer and the depot are assumed to be known. The VRPTW consists in computing least-cost vehicle routes such that each customer is visited exactly once within the indicated time window, while the vehicle capacity is not exceed for each route. It has been extensively studied in the literature (see Desaulniers et al. 2014) and the current state-of-the-art exact methodology for solving it is branchprice-and-cut (BPC, see Costa et al. 2019). In particular, the sophisticated BPC algorithms of Pecin et al. (2017a) and Sadykov et al. (2020) are able to solve to proven optimality most tested instances with up to 200 customers.
For the F-VRPTW, the vehicle capacity is organized in multiple identical stacks of a given maximal height. Customers might have items labeled as "fragile". To avoid potential damages to these items, the F-VRPTW imposes the fragility constraint, which implies the following: 1) non-fragile items are prevented to be stacked on top of fragile ones, however 2) fragile items are allowed to be stacked on top of other fragile items. Furthermore, no en-route load rearrangement is allowed. This problem is cast as a pickup one but the developed methodology can easily be adapted for a delivery problem, simply by inverting the role of the fragile and non-fragile items.
Our motivation in studying the F-VRPTW is its relevance for naval shipping companies delivering freight in northern Quebec and, to the best of our knowledge, this problem has never been studied before. However, although in a different loading framework (three-dimensional loading), the fragility constraint, in the exact sense used in the present work, was already introduced by Gendreau et al. (2006). Furthermore, the F-VRPTW is related to a number of studies addressing routing problems combined with several types of loading constraints. For an extensive literature review, the reader is referred to Pollaris et al. (2015) and references therein.
A first stream of work related to the F-VRPTW focuses on routing problems where the vehicle capacity is organized in multiple stacks (see Côté et al. 2012;Carrabs et al. 2013;Cherkesly et al. 2016;Veenstra et al. 2017, to cite a few). Handling multiple stacks responds to practical concerns in a variety of contexts, such as grocery, pharmaceutical and naval applications where either multicompartments vehicles are required, or freight need to be organized into pallets or containers which are then stacked for shipping. From the operational viewpoint, structuring the load in stacks limits the access to the only items located on the top of the stacks, unless time-consuming rearrangements of the freight are operated. In real-world operations, this limitation results in the so-called LIFO policy. Typically, the LIFO policy is a concern when the underlying routing decisions belong to the family of the Pickup and Delivery Problems (PDPs). In fact, the pickup order will influence which items must be delivered first, thus deeply impacting routing possibilities. In the VRPTW literature, however, the LIFO policy is generally not a concern because all the freight either originates from the depot (in case of deliveries) where it can be well positioned in the vehicle before departure or is destined to the depot (in case of pickups) where it can be typically unloaded in any order. Due to the structural differences between the PDPs and the VRPTW, providing an extensive literature review on the PDPs is out of the scope, and we again refer the interested reader to Pollaris et al. (2015). However, we notice that, given the additional fragility constraint in the F-VRPTW, the LIFO rule may render infeasible some customer sequences, and this needs to be taken into account when solving the F-VRPTW.
Another related stream of literature addresses the load-bearing strength constraint, which can be seen, to some extent, as a generalization of the fragility constraint introduced here. This constraint assures that the total weight stacked on top of a given item does not exceed a maximal value, thus preventing item damages. Junqueira et al. (2013) proposes an integer linear programming model for the vehicle routing problem with three-dimensional loading constraints (3L-CVRP), including the load-bearing strength. By adopting a commercial mixed-integer programming solver, the authors obtained the solution for some small problem instances. More recently, a single-stack variant of this problem was addressed by Chabot et al. (2017) in the context of warehouse order picking activities.
The authors propose several heuristics as well as exact branch-and-cut algorithms. As previously mentioned, Gendreau et al. (2006) address the 3L-CVRP where the fragility constraint is handled as in our setting. For this problem, the authors developed an efficient Tabu Search (TS) algorithm. Different metaheuristics have been developed for the same problem by several other authors, see Fuellerer et al. (2010); Tarantilis et al. (2009); Tao and Wang (2015), to cite a few.
In this paper, we propose a set partitioning formulation of the F-VRPTW, which we solve via a BPC algorithm. The most challenging aspect is in the solution of the column generation pricing problem, which turns out to be a variant of the Elementary Shortest Path Problem with Resource Constraints (ESPPRC). The main contribution of this paper is in the way we handle the fragility constraint. We solve the pricing problem by a labeling algorithm and provide a formal characterization of partial routes allowing us to easily check their feasibility with respect to the fragility constraint. In particular, this characterization allows us to avoid the intuitive but naive approach of duplicating labels to account for alternative freight configurations. To improve the performance of the labeling algorithm, we develop an efficient dominance rule with the purpose of discarding proven dominated labels.
To analyze the performance of our algorithm, we perform an extensive computational campaign.
We first introduce F-VRPTW instances that are obtained from the Solomon's VRPTW benchmark instances by varying the vehicle capacity and the maximal stack height. Our computational results show that the fragility constraint is more binding when vehicle capacity is smaller (i.e., for shorter routes) and when stacks are higher. Furthermore, to assess the performance of the proposed BPC algorithm, we have developed an alternative solution method, which is a classic BPC algorithm for the VRPTW (possibly generating routes violating the fragility constraint), augmented with infeasible path cuts to eliminate integer solutions violating the fragility constraint. Our test results show that the proposed BPC algorithm consistently outperforms the alternative one.
The rest of the paper is organized as follows. In Section 2, we formally state the F-VRPTW and present a mathematical formulation. In Section 3, we prove the main theoretical result of the paper, that gives a formal characterization of a route feasible with respect to the fragility constraint. This result is then exploited in Section 4 to develop an efficient BPC algorithm. We report and analyze computational results in Section 5, before providing concluding remarks in Section 6.

Problem statement and mathematical formulation
The classical VRPTW can be described as follows. Consider a directed graph G = (N , A) where N = {0, 1, . . . , n, n + 1} is the node set and A = {(i, j) | i, j ∈ N, i = j, (i, j) = (n + 1, 0)} is the arc set. We denote N c = {1, . . . , n} the subset of nodes representing the n customers, while nodes 0 and n + 1 represent the same depot at the beginning and the end of a route, respectively. A fleet of homogeneous vehicles with capacity Q is initially located at the depot. A demand volume q i (an integer number of items of the same size) to be collected and a time window [w i , w i ] are associated with each customer i ∈ N c . A vehicle is allowed to arrive at customer i earlier than w i , but must leave in the specified time window. A cost c ij and a travel time t ij (possibly including a service time at i) are associated with each arc (i, j) ∈ A. We assume that these costs and travel times satisfy the triangle inequality. Observe that set A can be reduced by eliminating all arcs (i, j) with i, j ∈ N c that are load or time infeasible, i.e., such that q i + q j > Q or w i + t ij > w j .
The VRPTW calls for finding a minimum-cost set of routes starting and ending at the depot, such that each customer is visited exactly once, and the time windows and vehicle capacity constraints are fulfilled.
The F-VRPTW differs from the VRPTW in that the vehicle capacity is organized in a set of m identical stacks M = {1, . . . , m}, each with a set of positions K = {1, . . . , k} (see Figure   1). Consequently, k is the height of a stack and the vehicle capacity is set as Q = km. Items are assumed to be put in the vehicle from the top and each of them occupies one position in one stack. Furthermore, the customer set N c is partitioned in two subsets L and H, identifying customers with fragile (light) and non-fragile (heavy) items, respectively. The fragility constraint requires that no item from customers in H can be stacked on the top of items from customers in L.
However, fragile items can be stacked on top of other fragile ones. For example, the positioning of the items in stack 2 in Figure 1 is infeasible, while the positioning in the other stacks is legitimate.
No en-route rearrangement of items is permitted. It should be noticed that the case of customers having both fragile and non-fragile items can be easily reduced to the present setting by adding suitable dummy nodes. We formulate the F-VRPTW as a set-partitioning problem. To this scope, consider a route r as a sequence of nodes r = (v 0 , v 1 , . . . , v p , v p+1 ), where v 0 = 0 and v p+1 = n + 1 represent the depot, and let Ω be set of all feasible routes. For all i ∈ N c , let a ir be a parameter with value 1 if route r visits customer i and 0 otherwise. The cost c r of route r ∈ Ω is computed as By introducing a binary variable x r for each route r ∈ Ω that takes value 1 if route r is chosen and 0 otherwise, the F-VRPTW can be formulated as follows: where the objective function (2) minimizes the total cost, constraints (3) ensure that each customer is visited exactly once, and constraints (4) impose that the variables are binary.

Characterization of fragility-feasible routes
The present section is devoted to introducing a formal characterization of a route that is feasible with respect to the fragility constraint (hereafter called a fragility-feasible route). As previously mentioned, the proposed solution method heavily stands on this characterization. In particular, as it will be clear in Section 4, this characterization allows us to drastically reduce the number of labels generated in the labeling algorithm developed to solve the column generation pricing problem.
The main observation behind the proposed characterization is that a route can be infeasible with respect to the fragility constraint only if, in the collecting process, we must form an incomplete stack containing fragile items, and then we keep collecting a large enough number of non-fragile items so that some of them must be placed on the incomplete stack with fragile items, hence violating the fragility constraint. As a consequence, the number of non-fragile items that can be collected after a given customer is bounded by a function depending on the number of fragile and non-fragile items collected so far and the order in which they were collected. However, when it is possible to complete all the stacks containing fragile items at a given customer with fragile items, the number of non-fragile items that can still be collected is only bounded by the vehicle capacity constraint.
Given a route r = (v 0 , v 1 , . . . , v p , v p+1 ), let us consider the induced sequence of individual collected items S r = (b 1 , b 2 , . . . , b s ). Given that customers belong either to set L or H, the order used to embed the items of a given customer in sequence S r is irrelevant. We denote by a L (i) the number of fragile items in S r until and including item b i , by a H (i) the corresponding number of non-fragile items, and by l(i) = a L (i) + a H (i) the total number of collected items. Even if l(i) = i here, notation l(i) is used because the items will not always be considered individually in the rest of this paper. Let us also define a modified modulo function as: Finally, by a slight abuse of notation, we will write b i ∈ L or b i ∈ H to express that b i is a fragile or non-fragile item, respectively.
where a H (s) denotes the total number of non-fragile items and is an upper bound on the number of non-fragile items that can be collected after item b i .
Before proving this proposition, let us present the intuition behind conditions (c1) and (c2), and how they can be interpreted. Each fragile item b i ∈ L in S r may induce an upper bound on the number of non-fragile items that can be collected after it. This upper bound corresponds to the optimal value of the problem of finding a feasible loading configuration for the first i items in S r that maximizes the number of positions available to load additional non-fragile items. Note that this upper bound does not depend on the items to be collected after item b i and the corresponding configuration is not necessarily a sub-configuration of a feasible configuration for the whole sequence S r . Consequently, this configuration is called, hereafter, a locally-optimal configuration for item b i i a H (i) + F (a L (i)) U (i) mk − l(i) max H (i)   1  1  9  11  9  3  3  -9  8  4  4  -8  8  5  2  6  7  6  6  3  -6  6  8  5  -4  4  9 3 -3 3 Table 1: Values for items in L for an example with m = 4 and k = 3 (even if b i ∈ H). When condition (c1) is not met for b i , i.e., a H (i)+F (a L (i)) < k, then there remain in any locally-optimal configuration for b i exactly k − l(i) mod k empty positions above a fragile item to which no non-fragile item can be assigned. Therefore, the maximum number of non-fragile items that can be loaded after item is the total number of empty positions. On the other hand, when condition (c1) holds for b i , we can show that either 1) item b i has no impact on the number of additional non-fragile items that can be loaded, or 2) the corresponding upper bound is equal to mk − l(i) and, thus, redundant with the vehicle capacity constraint. In both cases, there is no need to compute the upper bound.
To illustrate these different cases, let us consider the following example. Let m = 4, k = 3, where the upper index indicated if the item belongs to L or H. Table 1 displays relevant values for each item b L i : the left-hand side of condition (c1), i.e., a H (i) + F (a L (i)); the upper bound U (i) of condition (c2) whenever condition (c1) does not hold; the upper bound mk − l(i) induced by the vehicle capacity; and the exact maximum number of non-fragile items that can be loaded after item b L i (max H (i)). The latter is computed as: where B i is the subset of fragile item indices j ∈ {1, . . . , i} such that a H (j) + F (a L (j)) ≥ k. This maximum is obtained by considering the upper bounds induced by the vehicle capacity and by the fragile items b L j , j < i, for which condition (c1) does not hold. In the latter case, the corresponding upper bound U (j) is adjusted by subtracting the number of non-fragile items loaded between b L j and b L i .
Let us discuss the information reported in Table 1 for some specific items b L i . For i = 1, the values are obvious because item b L 1 must be assigned to a stack, leaving the other three stacks (i.e., 9 positions) to potentially load non-fragile items. For i = 3, condition (c1) is satisfied (a H (3) + F (a L (3)) = 3 ≥ 3). In this case, b L 3 can be positioned on top of a fragile item (b L 1 ) in the locally-optimal configuration for b H 2 (see Figure 2a) and, therefore, does not impact the number of non-fragile items that can still be loaded. Note that max H (3) = 8 stems from the upper bound U (1) = 9 from which item b H 2 has been subtracted. For i = 5, condition (c1) is not satisfied (a H (5) + F (a L (5)) = 2 < 3) and the locally-optimal configuration (see Figure 2b) Figure 2: Locally-optimal configurations for some items of the example as at least one (= k − l(5) mod k) of the seven (= mk − l(5)) unoccupied positions cannot be assigned to a non-fragile item (e.g., the empty position in stack 2 in Figure 2b). Note that item b L 5 induces a drop of two (from max H (4) = 8 to max H (5) = 6) in the exact maximum number of non-fragile items that can still be loaded. Finally, for i = 9, condition (c1) holds. In this case, the locally-optimal configuration for b L 9 (see Figure 2c) provides an upper bound that is equal to mk − l(9) = 3 and is, thus, redundant with the vehicle capacity bound.
To prove Proposition 3.1, we proceed by induction on the number of stacks m. To make explicit the dependence of the proposition on the number of stacks, we will use the notation U (i, m) instead of U (i). The main line of the proof consists first in showing that the basis of the induction P (m) with m = 1 is true. In particular we show that the feasibility and characterization conditions are equivalent: F eas(1) ⇐⇒ Char(1). We then prove the induction step P (m) =⇒ P (m + 1).
To prove that P (1) is true, we need to show that F eas(1) ⇐⇒ Char(1). Observe that the case with no fragile item is trivial because all routes are feasible with respect to the fragility constraint and, in this case, the characterization requires no condition to be verified. We can then concentrate on routes including fragile items (a L (s) > 0). We prove the sufficiency part first and the necessity second.
• Feas(1) =⇒ Char(1) Feasibility implies that non-fragile items are not stacked on top of fragile ones, i.e., a H (s) − a H (i) = 0 for all b i ∈ L. For all b i ∈ L such that l(i) < k, we have that U (i, 1) = 0, and condition (c2) is satisfied. Furthermore, if there exists b i ∈ L such that l(i) = k, then F (a L (i)) = a L (i) and a H (i) + F (a L (i)) = a H (i) + a L (i) = l(i) ≥ k and condition (c1) is verified.
• Feas(1) ⇐= Char(1) We prove the following counterpositive version of this statement: If sequence S r is infeasible with respect to the fragility constraint, then there exists an item b i ∈ L such that If S r is infeasible, then a non-fragile item is collected after a fragile one and there exists i ∈ {1, . . . , s} such that b i ∈ L and a H (s) − a H (i) > 0. Furthermore, l(i) < k because an item is loaded after b i . Consequently, a H (i) + F (a L (i)) = a H (i) + a L (i) = l(i) < k, thus condition (c3) is satisfied. Also, U (i, 1) = 0, and condition (c4) follows.
Here we prove that P (m + 1) is true when assuming that P (m) is true, and in particular we need to prove that F eas(m + 1) ⇐⇒ Char(m + 1). We do this by proving the necessity condition first and sufficiency second.
Feas(m + 1) ⇐= Char(m + 1) We start by considering a route r collecting at most (m+1)k items and complying with Char(m+ 1). We build a suitable subsequenceS r of at most mk items and we show thatS r complies with Char(m). By the induction hypothesis, this implies thatS r is a feasible sequence. The feasibility of S r will be trivially obtained by construction.
We distinguish two cases: 1. a L (s) ≥ k. In this case, route r collects enough fragile items to be able to complete a stack made of fragile items only. Consider the sequenceS r of at most mk items obtained from S r by eliminating the first k fragile items. Formally,S r = (b i | b i ∈ S r and either b i ∈ H or b i ∈ L and a L (i) > k). We use bars to indicate that a given entity refers toS r ; for example,ā H (i) is the number of non-fragile items collected inS r until and including item b i .
Given that the number of non-fragile items collected in S r andS r is the same and that exactly k fragile items have been deleted from S r , for all i ∈ {1, . . . , s} with b i ∈S r , we havē To prove thatS r is feasible, we show that it satisfies Char(m). To this scope, let us consider a generic fragile item b i inS r . Because we assume that Char(m + 1) holds for S r , b i satisfies condition (c1) or condition (c2) for S r , yielding the following two cases.
As a consequence of the induction hypothesis, the sequenceS r is feasible when considering m stacks. Then, for m + 1 stacks, S r is also feasible because it is always possible to stack the first k fragile items (that were put aside when buildingS r ) in a dedicated stack.
2. a L (s) < k. We use a similar technique as for the previous case. Because it is not possible to fill a stack with fragile items only, the strategy consists in building a dedicated stack with a number of non-fragile items collected before the first fragile item, and with all fragile items.
In case the number of non-fragile items collected before the first fragile item plus the total number of fragile items is less than k, the dedicated stack will be incomplete.
Let f be the index in S r of the first fragile item. Consider the subsequenceS r obtained from S r by eliminating the first min{a H (f ), k − a L (s)} non-fragile items and a L (s) fragile item. BecauseS r is only made of non-fragile items, it is certainly feasible with respect to the fragility constraint. However, we need to ensure the capacity constraint ofS, i.e., that l(s) ≤ mk is fulfilled. We distinguish two cases: • k − a L (s) ≤ a H (f ). In this case, the construction eliminates exactly k items andl(s) = l(s) − k ≤ mk.
• k − a L (s) > a H (f ). In this case, the construction eliminates a H (f ) + a L (s) items. Hence, ≤ U (f, m + 1) Equation (16) holds because l(s) = a H (s)+a L (s). Inequality (17) is due to the hypothesis that S r complies with Char(m + 1) and, in particular, with condition (c2) for b f because Consequently,S r is feasible and S r is also feasible, because it is always possible to build a dedicated stack with the eliminated items. This completes the proof of the necessity part of the statement.
Feas(m + 1) =⇒ Char(m + 1) We distinguish two cases: 1. a L (s) ≥ k. Here, we use some of the constructions developed in the above necessity proof.
We start with a sequence S r of at most (m + 1)k items, which we assume to be feasible. Then, we construct a suitable sequenceS r of at most mk items and show that it is feasible. By using the induction hypothesis (F eas(m) ⇐⇒ Char(m)), we then prove that S r complies with Char(m + 1).
As in the necessity proof, letS r be the subsequence obtained by eliminating the first k fragile items from S r . Because we assume that S r is feasible and it is always possible to put the first k fragile items in a dedicated stack without affecting the feasibility of S r , the subsequenceS r is also feasible and, by the induction hypothesis,S r complies with Char(m). We also observe that relations (7)-(9) hold by construction.
To prove that S r complies with Char(m + 1), we first consider the fragile items b i belonging to both S r andS r and observe that, by the same reasoning done for the necessity proof, they all satisfy at least one of the conditions of Char(m + 1). In fact, for a given i ∈ {1, . . . , s} such that b i ∈ L ∩S r , we have two possible cases: Then equations (7) and (8) imply Hence, condition (c2) holds for S r .
Now, let us verify that one of the conditions of Char(m + 1) also holds for the first k frag- where the inequality arises from the fact that a L (s) ≥ k and the equality holds because 2. a L (s) < k. In this case we can prove the statement directly, without employing the induction To prove that this condition is met, we start by observing We consider two cases: Therefore, condition (c2) holds.
ii) a H (f ) ≥ k − a L (s). In this case, a H (s) ≤ mk + k − a L (s) and we find where the last inequality comes from the assumption that This completes the proof of the sufficiency part of the statement.
The characterization provided in Proposition 3.1 requires that conditions (c1) and (c2) are satisfied for all fragile items b i in S r . The following proposition ensures that, if the characterization conditions are satisfied for a suitable subset of fragile items, they are satisfied for all fragile items.
In fact, it turns out that it is sufficient to check the conditions of the characterization only for b i ∈ L such that i = s or b i+1 ∈ L. In particular, when route r visits a customer i ∈ L with demand q i , there is no need to consider these items individually. Moreover, if r visits several customers in L consecutively, it is sufficient that the last item of the last of these customers verifies the conditions.

Proof
Let us denote the last two conditions (c5) and (c6). Observe that a H (i + 1) = a H (i), a L (i + 1) = a L (i) + 1 and l(i + 1) = l(i) + 1. Furthermore, it is easy to prove that U Let us now consider two cases: where the second equality is valid because F (a L (i)) < k − 1. Thus, condition (c5) is also satisfied.

If a H
which is not possible given the vehicle capacity. Consequently, case 2 is not possible and the proof is complete.
The labeling algorithm described in Section 4.1.2 largely stands on Propositions 3.1 and 3.2.
To conclude this section, we prove two other propositions concerning specific cases. Even if the stated results are not exploited in the proposed BPC algorithm, they are useful to understand when it is more difficult to satisfy the fragility constraint. Furthermore, they served to design our test instances and to analyze the obtained computational results.
From Propositions 3.1 and 3.2, we can prove that r is fragility-feasible by showing that, for the last item b i j of every customer v j ∈ L, condition (c2) is where the first inequality comes from the fact that some customers v j ′ , j ′ ∈ {j + 1, . . . , p}, might belong to set L and the second inequality from the assumption. Therefore, condition (c2) holds. Proposition 3.3 indicates that routes must be sufficiently filled to be fragility-infeasible. For example, if we consider m = 12 stacks of height k = 4, then only routes with a total load of 45 or more (i.e., with a filling rate exceeding 91.7%) can violate the fragility constraint. Next, we discuss the case k = 2.
Proposition 3.4 Let r = (v 0 , v 1 , . . . , v p , v p+1 ) be a route that is feasible with respect to the capacity constraint and let k = 2. Route r does not satisfy the fragility constraint if and only if p ≥ 2, Proof ( =⇒ ) Assume that r does not satisfy the fragility constraint. From Propositions 3.1 and 3.2, it means that there exists an index j * ∈ {1, . . . , p} such that v j * ∈ L and both conditions (c1) and (c2) are not met for the last item , if a H (s) + a L (i j * ) = Q or, equivalently, p j=1 q v j = Q and v j ∈ H for all j ∈ {j * + 1, . . . , p}.
( ⇐= ) Assume that p ≥ 2, p j=1 q v j = Q and there exists an index j * ∈ {1, . . . , p − 1} such that v j ∈ L for all j ∈ {1, . . . , j * }, v j ∈ H for all j ∈ {j * + 1, . . . , p} and j * j=1 q v j is odd. Denote by b i j * the last item at customer v j * . Let us show that both conditions (c1) and (c2) of Proposition 3.1 do not hold for b i j * , indicating that route r is fragility-infeasible. From the assumption, we can easily deduce that a H (i j * ) = 0 and F (a L (i j * )) = 1. Thus, (c1) is violated.
Proposition 3.4 describes explicitly the routes that do not satisfy the fragility constraint when k = 2. In particular, it shows that, if a customer in H is visited before a customer in L, then the route is necessary fragility-feasible. Indeed, it gives the opportunity to start a stack with a nonfragile item and fill it with a fragile item if needed. This principle, which is expressed by condition (c1), is also valid for larger stack heights k. From it, we can deduce that routes visiting customers in H with sufficiently large demands before customers in L have high chances to be fragility-feasible.
This may not be true for a large stack height k but, in practice, k is not too large.

The branch-price-and-cut framework
BPC is a variant of the branch-and-bound algorithm where the lower bounds at the nodes of the search tree are computed via column generation. Lower bounds are then strengthened by the dynamic generation of valid inequalities. In Section 4.1, we present the proposed column generation algorithm for solving the linear relaxation of the F-VRPTW formulation (2)-(4). In Section 4.2, we describe some acceleration strategies. Cutting planes and branching strategies are then discussed in Sections 4.3 and 4.4, respectively. Finally, an alternative BPC algorithm which does not exploit the fragility-feasible route characterization is introduced in Section 4.5. This algorithm will be used as a benchmark against the proposed BPC algorithm.

Column generation
At each node of the search tree, the BPC algorithm must compute a lower bound by solving the linear relaxation of (2)-(4) possibly augmented by the constraints implied by the branching decisions or valid inequalities previously added. However, due to the extremely large number of variables in (2)-(4), this linear problem cannot be solved directly and, consequently, an iterative column generation procedure needs to be invoked (see, e.g., Lübbecke and Desrosiers 2005). A column generation iteration first solves a Restricted Master Problem (RMP) where only a small subset Ω ′ ⊂ Ω of feasible routes is considered. Letx be the solution of the RMP. The iteration then proceeds by solving a pricing problem to verify the optimality ofx for the whole linear relaxation.
When the optimality check fails, a new set of routes (identified by the pricing problem) is added to the RMP, and the process iterates. Otherwise, the column generation procedure stops as the current RMP solution yields a valid lower bound.
In Section 4.1.1, we define the pricing problem, which turns out to be a variant of the ESPPRC.
This problem is solved using the labeling algorithm described in Section 4.1.2.

Pricing problem formulation
To verify the optimality of the current RMP solution in the column generation procedure, the pricing problem searches for routes with a negative reduced cost. Let π i for all i ∈ N c be the dual variables associated with constraints (3). The reduced cost of a route can be expressed as follows: Thus, the pricing problem minimizes (21) over the set of all feasible routes Ω.
Any route r ∈ Ω can be represented as a path in network G. By combining expressions (1) and (21), it is possible to express the reduced costc r of a route r = (v 0 , v 1 , . . . , v p , v p+1 ) as the sum of the contributions of its arcs: This property can be suitably used to express the column generation pricing problem as an arc-flow formulation which is presented in Appendix A. This formulation is, however, a non-linear integer program that is hard to solve via state-of-the-art solvers. On the other hand, we observe that the pricing problem corresponds to an ESPPRC, where the resources ensure route feasibility with respect to the time windows, the vehicle capacity, and the fragility constraint.

Labeling algorithm
In the vehicle routing literature, the ESPPRCs are frequently solved by dynamic programming, which is usually implemented by means of a labeling algorithm (Irnich and Desaulniers 2005). In our specific problem setting, given that the fragility constraint is dealt with in the pricing problem, the labeling algorithm is potentially very challenging. In fact, given a route visiting a set of customers in a specific order, there are many possible ways to position the freight in the vehicle.
Although some of these configurations may be infeasible, there could still exist a large number of feasible loading configurations. A potential labeling algorithm would create one new label for each feasible configuration or, if clever modeling is used to avoid symmetry between the identical stacks as in Cherkesly et al. (2016), one new label for each feasible stack-anonymous configuration.
However, both intuition and preliminary computational tests for the simplest case with stacks of height k = 2 suggest that this method is not efficient due to the large number of generated labels (Altman 2017).
We then directed our research efforts towards the theoretical results of Propositions 3.1 and 3.2.
In particular, the characterization of a fragility-feasible route enables us to carry information about the existence of a feasible loading configuration for a given partial route. In practice, as long as at least one feasible configuration exists, our labeling algorithm extends partial routes. However, no attempt is made to build a specific configuration, which can be trivially retrieved a posteriori.
A labeling algorithm (see Irnich and Desaulniers 2005) uses labels to represent partial paths (routes) in a network. Paths are enumerated by extending recursively an initial label E 0 from the source node towards the destination. Labels are extended according to extension functions. A dominance rule is applied to eliminate partial routes that cannot yield complete optimal routes. These algorithmic components are described next.
with the following n + 6 components: • one component C r accounting for the reduced cost; • n binary components V l r , l ∈ N c , indicating whether or not customer node l is unreachable, i.e., it cannot be visited in any feasible extension of route r because it has already been visited or because its time window or the vehicle capacity cannot be satisfied; • one component T r indicating the (earliest) service starting time at node r; • one component L r accounting for total collected load; • two components A H r and A L r providing the total number of non-fragile and fragile items collected, respectively; • one component X r indicating the maximum number of non-fragile items that can still be collected according to Propositions 3.1 and 3.2.
The initial label at node 0 is Label extension functions. Given a partial route r = (0, . . . , i), i ∈ N \ {n + 1}, with label , it can be extended along an arc (i, j) ∈ A using the following label extension functions to yield a new partial route r ′ = (0, . . . , i, j) represented by the label where Extension function (29) relies on Proposition 3.1 and the vehicle capacity constraint to compute the maximum number of non-fragile items that can still be loaded after visiting a customer. In the first case, given that A H r ′ + F (A L r ′ ) ≥ k, the loading of the q j fragile items at node j does not affect the maximum number of non-fragile items that can still be loaded in virtue of the fragility constraint. However, if the capacity constraint becomes binding, the term Q − L r ′ in the minimum function ensures that X r ′ does not exceed the residual available space (i.e., X r ′ ≤ Q − L r ′ ). The second case is similar but, when A H r ′ + F (A L r ′ ) < k, Proposition 3.1 imposes a maximum number of non-fragile items U r ′ . In this case, the capacity limit is implicit because Q − L r ′ > U r ′ . Finally, in the third case, because j ∈ H or j = n + 1, the extension function simply decreases the number of non-fragile items by q j .
The obtained label E r ′ is declared feasible if j = n + 1 or if V l r ′ ≤ 1 for all l ∈ N c , T r ′ ≤ w j , L r ′ ≤ Q, and X r ′ ≥ 0. Otherwise, it is discarded. Note that all these conditions except the last one can be verified before performing the label extension to a node j ∈ N c by checking the value of V j r as follows. If V j r = 1, then at least one of these conditions will not be met and there is no need to compute E r ′ . Otherwise, the extension must be performed and the condition X r ′ ≥ 0 must be checked to determine the feasibility of label E r ′ .
Dominance rule. Let r d , d = 1, 2, be two partial routes, both ending at the same node i ∈ N and represented by the labels (c7) any feasible (single-or multiple-arc) extension e of r 2 is also feasible for r 1 , and (c8) for any such extension e, the inequality C r 1 ⊕e ≤ C r 2 ⊕e holds, where symbol ⊕ denotes the concatenation operator.
The dominated label E r 2 can then be discarded. However, when multiple labels dominate each other, one of them must be kept.
The dominance conditions (c7) and (c8) cannot be checked easily in practice. In general, sufficient conditions defining a so-called dominance rule are used instead. In particular, when all label extension functions are monotone (either non-decreasing or non-increasing), it is easy to provide such a dominance rule (see, e.g., Irnich and Desaulniers 2005). In our case, monotonicity does not hold for extension function (29). Consequently, the dominance rule that we introduce in the following proposition is more complex and will be proven.
Proposition 4.1 Route r 1 dominates route r 2 if the following relations hold: In the following, we use R(σ 1 , σ 2 ) as a shorthand notation for the set of conditions (30)-(34) when r 1 = σ 1 and r 2 = σ 2 . Before proving Proposition 4.1, we state and prove the following preliminary lemma.
Lemma 4.1 If relations R(r 1 , r 2 ) are true, then relations R(r 1 ⊕ e, r 2 ⊕ e) also hold for any extension e of r 2 .

Proof.
We begin by showing that the statement holds true for any single-arc extension e made of a generic arc (i, j). To highlight that routes r d , d = 1, 2, end at node i, we denote their labels Similarly, the labels for routes r d ⊕ e are written Given that the extension functions (23)- (25) and (28) are non-decreasing, it is trivial to show that the statement holds for relations (30)-(33). We then concentrate on relation (34) for which we want to prove that: We distinguish three cases related to the extension of r 1 , one for each case in (29).
• If j ∈ H or j = n + 1, then where the equality ensues from X 1j = X 1i +q j and L 1j −A H 1j = L 1i +q j −(A H 1i +q j ) = L 1i −A H 1i , and the second inequality from the assumption.
where the equality is derived from X 1j = min{X 1i , Q−L 1i −q j } and L 1j −A H 1j = L 1i −A H 1i +q j , the second inequality from the assumption, and the last from L 1i ≤ L 2i .
where the equality arises from X 1j = min{X 1i , U j } and 1i − 1) − q j (as discussed next), and the third inequality from the assumption. Showing that where the second equality stems from A L 1j ≥ 1 (as j ∈ L), yielding A L 1j = F (A L 1j ), the third equality from A H 1j + F (A L 1j ) < k, and the inequality from This proves that the implication (35) holds and that the statement is true for single-arc extensions.
By induction on the number of single-arc extensions, the lemma is also true for any extension e with an arbitrary number of arcs.

Proof of Proposition 4.1.
To prove this proposition, we show that, if relations R(r 1 , r 2 ) are satisfied, then conditions (c7) and (c8) are also satisfied. Let e be a feasible extension of r 2 that ends at a node j ∈ N . Denote by Assuming that relations R(r 1 , r 2 ) hold, we deduce by Lemma 4.1 that From inequality (40), we can easily deduce that because Q − L 2j ≥ X 2j by definition of X 2j . Therefore, inequalities (37)- (39) and (41) indicate that r 1 ⊕ e is feasible, i.e., that condition (c7) holds, whereas inequality (36) implies condition (c8).

Acceleration strategies
Acceleration strategies play a key role for the development of efficient BPC algorithms. In this section we describe the adopted procedures, all of them geared towards improving the performance of the labeling algorithm in the column generation step.

Decremental state space relaxation
An efficient technique to generate negative reduced cost elementary paths is the decremental state space relaxation (DSSR) introduced independently by Boland et al. (2006) and Righini and Salani (2008). The pricing problem is initially solved without the label components V 1 , . . . , V n , thus completely relaxing the elementarity requirements. If no negative reduced cost paths are found in this way, some of the corresponding label components are dynamically added, and the process is iterated until elementary paths with a negative reduced cost are found or the reduced cost of the shortest path is non-negative. As proposed in Desaulniers et al. (2008), at each iteration of the column generation, instead of initializing the algorithm with an empty set of visit-label components, we use the components generated in the previous iteration.

The ng-path relaxation
The technique described in this section consists in partially relaxing the elementarity requirements of partial routes. The adopted path relaxation, called ng-path, was introduced by Baldacci et al.
(2011). With each node i ∈ N c , the approach associates a subset of nodes N i ⊆ N c such that i ∈ N i and |N i | ≤ ∆ 0 , where ∆ 0 is a predefined integer parameter. Typically, N i contains the closest customers of i. Given a partial route r = (v 0 , . . . , v p ), the subsets N i allow us to define a new subset Π(r) whose elements are prevented to be direct extension candidates for r. The subset Π(r) is defined as Π The value of ∆ 0 can influence the degree of elementarity of an ng-path. A small value of ∆ 0 allows the ng-paths to contain many cycles, while ∆ 0 = |N c | imposes elementarity. Allowing cycles generally makes the corresponding ESPPRC easier to solve, but the quality of the provided lower bound deteriorates.
Hence the right choice of ∆ 0 is key for the overall computational efficiency of the BPC algorithm.
Instead of choosing neighborhoods of the same size for all nodes, Contardo et al. (2015) propose to adjust them dynamically and individually in a DSSR fashion as follows. First, the size of each neighborhood is initially set to a relatively small value ∆ 0 (10 in our tests). Column generation is then applied to solve a linear relaxation. If the computed linear relaxation solution contains routes with a cycle, neighbors are added to the nodes in these cycles in such a way that they become forbidden. The corresponding columns are then removed from the RMP and column generation is started over again. The process is repeated until the linear relaxation solution is free of cycles or when it is not possible to forbid the remaining cycles because a maximum of ∆ + additional neighbors has already been added to a node (∆ + = 10 for our tests).
To implement ng-path, the extension function of the label components V 1 , . . . , V n needs to be suitably modified (see Desaulniers et al. 2014;Costa et al. 2019, for a detailed description of the technique).

Heuristic dynamic programming
To rapidly obtain routes with negative reduced costs, we adopt graphs with reduced size. In particular, at each iteration, the labeling algorithm is first executed on a simplified graph G ′ containing only a subset A ′ of the arcs in A. In case no route with negative reduced cost is found, the labeling algorithm is run again, but on the complete graph. As it is detailed in Desaulniers et al. (2008), the choice of the initial set of arcs A ′ depends on the current modified costs of the arcs (22). Thus, the set A ′ changes at every iteration, according to the dual values corresponding to RMP solution. More precisely, for every node i ∈ N , all incoming arcs and all outgoing arcs are sorted according to their modified cost in two separated lists I i and O i . If an arc (i, j) is such that its position in both I j and O i is larger than a given threshold τ , then arc (i, j) is eliminated. In this study, we set τ = 10.

Cutting planes
At a given node of the branch-and-bound tree, the column generation process provides a linear relaxation solution that may be fractional, as well as a valid lower bound on the value of the node.
When this solution is fractional, it is possible to strengthen the obtained lower bound by looking for potential violated valid inequalities.
We consider two families of inequalities. The first was introduced in Jepsen et al. (2008) and is known as subset-row inequalities, which are Chvátal-Gomory rank-one cuts defined over subsets of the set partitioning constraints (3). In the general case, subset-row inequalities are expressed as: Following Jepsen et al. (2008), we only consider cuts obtained with |W | = 3 and k = 2, resulting where Ω W is the subset of paths visiting at least two customers in W . In a column generation method, the addition of subset-row inequalities to the RMP requires several adjustments to the pricing problem (see, e.g., Desaulniers et al. 2011), as well as careful management when the number of added cuts increases. The implementation proposed in this study, follows the one described in Desaulniers et al. (2008). Note that, given the length of the routes generated for our tests, it was not worthwhile to apply a more sophisticated variant of the subset-row cuts such as the limitednode-memory subset-row cuts (Pecin et al. 2017b).
The second family of inequalities we consider is the 2-path cuts, introduced by Kohl et al. (1999) for the VRPTW. Let W ⊆ N c be a subset of nodes such that it is not possible (e.g., because of the time windows or the vehicle capacity) to serve all customers in W by a single route. The corresponding 2-path inequality is given by where µ W r is the number of times route r "enters" the set W , i.e., µ W r = i∈N \W i∈W η ijr , where non negative parameter η ijr is equal to the number of times route r traverses arc (i, j) ∈ A.
We separate these cuts by a heuristic procedure similar to that of Desaulniers et al. (2008). In particular, we enumerate all subsets W ⊆ N c such that 1) |W | is less than or equal to a given value (15 in our studies), 2) the flow entering W is strictly less than two, and 3) the nodes in W are connected in the support graph of the current linear relaxation solution. If i∈W q i > Q, then a violated inequality is found. Otherwise, we use dynamic programming to solve a minimum Hamiltonian path problem with time windows over the set W ∪ {o, d}. If the problem is infeasible, a violated inequality is found. The new violated inequalities are added to the RMP. Furthermore, the dual variable associated with the 2-path cut for subset W must be subtracted from the modified arc costc ij for all arcs (i, j) ∈ A with i ∈ N \ W and j ∈ W .

Branching strategies
To enforce integrality, we apply the following branching scheme whenever required. First, we branch on the total number of vehicles used. If this number is integer, we branch on the arc-flow variables.
In this case, we select the arc (i, j) with flow closest to 0.5. To fix the flow on this arc to 0, we simply remove (i, j) from set A. To fix it to 1, we remove from A all arcs (i, ℓ), ℓ = j if i = 0 and all arcs (ℓ, j), ℓ = i if j = n + 1. The columns in the RMP containing any removed arc are deleted.
Finally, the branch-and-bound tree is explored using a best-first strategy.

Alternative algorithm
To assess the effectiveness of the proposed BPC algorithm (denoted BPC-FC-PP because it handles the fragility constraint in the pricing problem), we implemented another BPC algorithm that does not exploit the fragility-feasible route characterization developed in Section 3. In fact, one may suspect that fragility-infeasible routes are relatively rare and should not often be encountered during the solution process. Therefore, it might not be worthwhile to consider the fragility constraint while solving the pricing problem and infeasible path cuts should rather be added to the linear relaxation whenever a linear relaxation solution contains a fragility-infeasible route. Infeasible path cuts have been used in the vehicle routing literature to handle difficult-to-model constraints pertaining to individual routes (see, e.g., Cordeau 2006;Cherkesly et al. 2016). Given an infeasible route r * composed of a subset of arcs A r * , the corresponding infeasible path cut writes as: where Ω ′ is the set of routes that are feasible with respect to the time windows and the vehicle capacity and, for each arc (i, j) ∈ A, non negative parameter η ijr is again equal to the number of times route r ∈ Ω ′ traverses arc (i, j).
In summary, the alternative BPC algorithm is identical to the BPC algorithm described above except that 1) the labeling algorithm does not consider the A H , A L and X label components, nor the condition (34) in the dominance rule; and 2) infeasible path cuts of the form (44) are generated as needed whenever a linear relaxation solution contains a fragility-infeasible route. To highlight that the fragility constraint is treated in the master problem, this BPC algorithm is denoted BPC-FC-MP.

Computational experiments
In this section, we present the results of the computational experiments that we conducted to assess the effectiveness of the proposed BPC algorithm and to compare the solutions of the F-VRPTW obtained when varying the number of stacks and against the VRPTW solutions. In Section 5.1, we present the instances used for these tests. In Section 5.2, we describe the tests performed and which algorithms were involved in these tests. In Section 5.3, we report the computational results obtained and discuss computational performance. Finally, a comparison between the costs of the F-VRPTW and VRPTW solutions is realized in Section 5.4.

Test instances
To perform our tests, we first selected the 12 instances in the class R1 of the well-known Solomon's VRPTW benchmark instances (that can be found at http://w.cba.neu.edu/∼msolomon/r101.htm).
In these instances involving 100 customers each, the customers are randomly located in a square grid and the time windows are relatively narrow with respect to the traveling times. These instances differ only by their customer time windows. Because the fragility constraint has more chances to be binding if vehicle capacity is relatively tight, i.e., if the filling rate of the vehicles is close to 100%, and the routes do not contain too many customers in H and L intertwined (see the end of Section 3), we have considered different tight vehicle capacities, namely, Q = 48, 60, 72. To study the impact of the number of stacks on the computational times and solution costs, we have also considered different stack heights relevant in practice, namely, k = 2, 3, 4, 6, yielding varying numbers of stacks m = Q/k for each tested vehicle capacity Q.
Finally, for each instance, we assumed that a proportion ρ of the customers belong to H (with non-fragile items) and, thus, (1 − ρ) belong to set L (with fragile items). We considered three different ρ values, namely, 25%, 50%, and 75%. To determine to which set each customer in N c belongs to, we used the following procedure. Let us assume that the customers are numbered from 1 to n according to the order provided by Solomon. Then, if ρ = 50%, the customers with a number equal to 2i, i = 1, 2, 3, . . . belongs to H; if ρ = 25%, the customers with a number equal to 4i − 3, i = 1, 2, 3, . . . belongs to H; and, if ρ = 75%, the customers with a number equal to 4i − 1, i = 1, 2, 3, . . . belongs to L.
Given the difficulty to solve some of these 100 customers instances to optimality, we created smaller instances by extracting the first 50 customers or the first 75 customers of each 100-customer instances. Thus, for ρ = 50%, we ran tests on a total of 12 × 3 × 4 × 3 = 432 instances. For ρ = 25% and ρ = 75%, we also performed tests on instances with 75 customers, Q = 60 and k = 2, 3, 4, 6.

Computational campaign and environment
Our computational experiments proceeded as follows. First, to assess the efficiency of algorithm BPC-FC-PP, we solved the 432 instances with ρ = 50% using both BPC algorithms, i.e., BPC-FC-PP and BPC-FC-MP. Then, to evaluate the impact of considering the fragility constraint on the computational times and solution costs, we also solved the corresponding VRPTW solutions using a classic BPC algorithm, i.e., algorithm BPC-FC-MP without generating any infeasible path cuts.
Note that all these algorithms do not include all state-of-the-art strategies, such as route enumeration, variable fixing, or limited-memory Chvátal-Gomory rank-one inequalities (see Costa et al. 2019, for details). Therefore, the computational times achieved for the VRPTW instances are not competitive with those obtained by the state-of-the-art algorithms of Pecin et al. (2017a) and Sadykov et al. (2020).
The results collected from these experiments allow to perform sensitivity analyses with respect to the number of customers, the vehicle capacity, and the stack height, but not with respect to the proportion of customers in sets L and H. For the latter, we also solved the F-VRPTW instances with ρ = 25% and ρ = 75% (for 75 customers, Q = 60 and k = 2, 3, 4, 6) using algorithm BPC-FC-PP only.
All algorithms were coded in C/C++ using the Gencol 4.5 library and the primal simplex algorithm of Cplex 12.6.3 for solving the RMPs. Our tests were performed on a PC equipped with eight Intel Core i7-4770 processors clocked at 3.40 GHz, and 16 Gb of RAM. For all runs, a single core was used and a 2-hour time limit was imposed for solving each instance.

Computational performance comparison
We start by presenting the results of our first experiment where, as mentioned above, we solved 432 instances using three algorithms. The detailed results for each individual instance can be found in Appendix B. We report in Table 2  In fact, the former succeeded to solve 115 instances with 100 customers, compared to only 50 instances for the latter. Thus, we can say that BPC-FC-PP clearly outperforms BPC-FC-MP. This is obviously explained by how the fragility constraint is treated in each algorithm because this is the only difference between these two algorithms. Handling it in the pricing problem increases the difficulty of solving this pricing problem but, as the results show, it yields substantially smaller average integrality gaps (note that the gap differences would probably be much larger for the instances with 75 and 100 customers if the same number of instances were solved to optimality by both algorithms) and, therefore, less branch-and-bound nodes. Also, the algorithm BPC-FC-MP   in BPC-FC-PP increases the difficulty of solving the pricing problem.
Next, we study separately the impact of the vehicle capacity and the stack height on the computational performance of algorithm BPC-FC-PP. Table 3 reports the average results obtained by BPC-FC-PP by group of 144 instances with the same vehicle capacity (still with ρ = 50%). For each instance group, the same statistics as in Table 2 are provided. These results (in particular, the total number of instances solved to optimality) indicate that the difficulty of solving the F-VRPTW increases with the vehicle capacity. Indeed, increasing vehicle capacity allows longer routes to be generated and makes the pricing problem harder to solve. Furthermore, as shown by the results, the average integrality gaps tend to increase with vehicle capacity, requiring more effort to derive an optimal integer solution. However, this is not corroborated by the average number of nodes explored that decreases when Q increases. This may be explained by the fact that, compared to the case Q = 48, less instances were solved to optimality and more subset-row cuts were generated for the cases Q = 60 and 72. Note that, from these results, it is not easy to see if the impact of handling the fragility constraint varies with the vehicle capacity. Nonetheless, given that the labeling algorithm used for the F-VRPTW requires an additional dominance condition (34) compared to the algorithm used for the VRPTW, we can speculate that handling the fragility constraint contributes to the increased difficulty of solving the pricing problem.
Using the same format as in Table 3, Table 4 displays the average BPC-FC-PP results for groups of 108 instances with the same stack height. Here, the variation between the results obtained for the different stack heights is much less pronounced. In fact, the number of solved instances for each stack height is the same except for k = 6 for which three additional instances were solved. On the other hand, a slight increase of the average computational time is observed when k increases.
Given that the average integrality gaps are relatively constant with respect to stack height, this small computational time increase might be due to slightly harder-to-solve pricing problems. In fact, one can observe that, for a given route r represented by label E r , if the number of collected non-fragile items A H r is greater than or equal to k − 1, then condition (c2) of Proposition 3.1 is always fulfilled for any subsequently visited customer in L and X r is often equal to Q − L r . In this case, dominance condition (32) becomes redundant with condition (34) and the chances that E r dominates another label are increased. Thus, according to this hypothesis, a larger value of k may reduce the number of dominated labels.
We conducted a second series of experiments to see if the proportion of customers in set H influences the computational time. In this series, we solved all 75-customer instances with Q = 60 k #Opt T (s) Gap (%) #Nodes   and three different ρ values, namely, 25%, 50%, and 75%, using the BPC-FC-PP algorithm. Using the same statistics as above, the average results (over 48 instances per row) are reported in Table 5.
These results indicate that the instances with ρ = 50% seem easier to solve than those with ρ = 25% and ρ = 75%: compared to the latter, the former exhibit average time reductions of 40% and 27%, respectively. It should be noticed that the high average time for the instances with ρ = 25% is mostly due to one outlier (instance R104 with k = 4) which required more than 6700 seconds compared to less than 1800 seconds for all others. Removing this outlier would bring down the average time to around 290 seconds for the instances with ρ = 25%, that is, still a 12% difference.
Consequently, a balanced mix of customers with fragile and non-fragile items tends to ease the solution process, but not substantially. A possible reason for explaining this tendency is that the fragility constraint is more binding when ρ is close to 50% and, thus, reduces more the solution space as discussed at the end of the next section.

Solution cost comparison
In this section, we compare the solutions obtained for the VRPTW and the F-VRPTW, namely, we compute the number of times that the optimal solution computed for the VRPTW is not feasible for the F-VRPTW and the induced average cost increase. As in Section 5.3, we analyze the results with respect to the number of customers, the vehicle capacity, the stack height, and the proportion of customers in H. The results are presented in Tables 6 to 9. In all these tables, we report the following information by group of instances (the group definition depends on the table): the total number of instances for which the F-VRPTW was solved to optimality by either the BPC-FC-PP or the BPC-FC-MP algorithm (#Opt); the total number of VRPTW solutions computed for these instances that are not feasible for the F-VRPTW (#ModS) and, in parentheses, the percentage of instances with a modified solution, i.e., 100 #M odS #Opt ; and the average cost increase (in percentage) induced by the modified solutions (∆z IP ).   Table 7: Average cost increase by capacity (ρ = 50%) Table 6 provides the results with respect to the number of customers, i.e., each row summarizes the results for all instances with the same number of customers. From these results, we can clearly observe that the number of instances impacted by the fragility constraint increases with the number of customers, going from 44% when n = 50 to 89% when n = 100. This can be explained by the fact that the average number of routes in an optimal solution increases with the number of customers yielding more chances that at least one route in the solution is fragility-infeasible. The results also show that the average cost increase is relatively low (between 0.35% to 0.61%) and tends to decrease with the number of customers. Here, it is likely that more customers offer more options to find alternative fragility-feasible routes.
In Table 7, we present the results obtained when grouping the instances by vehicle capacity.
From these results, we observe that the number of instances with a modified solution tends to decrease with vehicle capacity. This is to be expected because, as routes contain more customers, it becomes easier to find a feasible loading configuration. However, given that the proportion of instances with a modified solution seems to stall between Q = 60 (59%) and Q = 72 (58%), we ran additional tests with Q = 84. For this capacity, the proportion decreases to 38%, confirming the observed trend. On the other hand, the results also indicate that the fragility constraint impacts less the solution cost as the vehicle capacity increases. This may be explained by the fact that, when capacity increases, the solutions contain less routes and, therefore, less of them are fragilityinfeasible and need to be changed.
Next, we study the impact that handling the fragility constraint has on the cost of the solutions with respect to the stack height. The results for instances grouped by stack height are reported in Table 8. Without surprise, they indicate that the number of instances with a modified solution and the average cost increase are influenced by stack height. Over all instances with the same height, we can observe an increase from 55% to 78% of the #ModS statistic when k goes from 2 to 6, and an increase from 0.30% to 0.55% of ∆z IP . These increases are simply due to a reduction of the solution space that occurs when k increases.   Finally, Table 9 provides results for instances grouped by proportion of customers in H. They indicate that the fragility constraint is much more binding when there is a balance mix of customers with fragile and non-fragile items as 71% of the optimal VRPTW solutions are fragility-infeasible when ρ = 50%, whereas only 23% and 33% of the optimal solutions are infeasible for ρ = 25% and ρ = 75%, respectively. Furthermore, the cost increase is also larger for ρ = 50% as more routes in the VRPTW solutions become infeasible for the F-VRPTW. We explain this behavior as follows.
As ρ tends toward 0%, the optimal solution routes tend to contain only customers in L. They are, thus, fragility-feasible. Symmetrically, when ρ tends toward 100%, the routes tend to be composed of customers in H only and are, thus, also fragility-feasible. Therefore, the maximum number of routes that include both customer types and that are at risk of being fragility-infeasible occur when ρ = 50%.

Conclusions
In this paper, we studied the F-VRPTW, a new variant of the VRPTW that considers a fragility constraint restricting the positioning of the items in the stacks of the vehicles. We first established necessary and sufficient conditions that characterize the feasibility of a route with respect to the fragility constraint. Then, we developed a BPC algorithm, denoted BPC-FC-PP, which exploits this characterization to handle efficiently the fragility constraint in the pricing problem. To evaluate this algorithm, we performed computational experiments on instances derived from some well-known VRPTW benchmark instances and that have been designed to yield a relatively tight fragility constraint. We compared the computational performance of BPC-FC-PP with that of another BPC algorithm, denoted BPC-FC-MP, that handles the fragility constraint in the master problem through infeasible path cuts. This comparison shows that handling the fragility constraint in the pricing problem is much more efficient than in the master problem: out of 432 instances, BPC-FC-PP solved 391 instances to optimality within a 2-hour time limit, whereas BPC-FC-MP could only solve 301 instances. Finally, we analyzed the impact that the fragility constraint has on the optimal solution cost with respect to the VRPTW solutions. Our results indicate that, for the tested instances, the computed VRPTW optimal solution is infeasible for the F-VRPTW in 67% of the cases. This proportion increases with the number of customers and the stack height, but decreases with vehicle capacity. We also observed that the fragility constraint is at its maximal tightness when there is a balance mix of customers with fragile items and with non-fragile items.
For the instances requiring a modified solution, an average cost increase of around 36% is incurred.
As possible future works, we can think about tackling more general problem variants. In particular, stacks of various heights could be considered as well as fragile items and non-fragile items that do not take the same space in a stack.

A An arc-flow formulation for the pricing problem
In this appendix, we propose an arc-flow formulation for the column generation pricing problem.
For each (i, j) ∈ A, we introduce a binary variable y ij that takes value 1 if the route traverses arc (i, j), and 0 otherwise. For each node i ∈ N , let T i be a continuous variable indicating the service starting time. Furthermore, for each (γ, β) ∈ M × K potential position of an item in the vehicle and each customer i ∈ N c , we introduce a binary variable z iγβ which is equal to 1 if an item of customer i is placed in position (γ, β), and 0 otherwise.
The pricing problem can then be formulated as follows: whereT is an upper bound on the service starting time at a customer and T 0 = 0.
Objective function (45) minimizes the total reduced cost of the route, where parametersc ij are computed according to (22). Constraints (46)-(48) ensure that the y ij variables describe a path in G between the origin depot node 0 and the destination depot node n + 1. Constraints (49) impose lower bounds on the service starting times, while (50) enforce time window restrictions. Relations (51)-(55) regulate the loading rules. In particular, inequalities (51) ensure that no more than one item is put in a given position, equalities (52) impose that the demand of each visited customer is fully satisfied, inequalities (53) imply that items are stacked on the top of others, while constraints (54) ensure coherence between the order of the customer visits and the relative positions of the items in each stack. Furthermore, relations (55) impose that no fragile item will be put under a non-fragile one. Finally, (56)-(58) define the variable domains.
Note that the capacity constraint is not expressed explicitly but rather implicitly by the number of available positions in M × K.

B Detailed computational results
This appendix presents the detailed results obtained from our computational experiments. They are reported in Tables 10-20, where each table is  ∆z IP : The difference in percentage between the optimal values of the F-VRPTW and the VRPTW, whenever z IP is an optimal value. Note that all VRPTW optimal values have been computed even if some of them required to exceed the time limit.
Note that Tables 19 and 20 do not report any results for algorithm BPC-FC-MP.       Table 17: Results for n = 100, Q = 60 and ρ = 50%