On solving the cross-dock door assignment problem

A class of strong lower bounds on the solution value of a Linearised Integer Programming reformulation is introduced for the binary quadratic optimisation model to assign origin and destination nodes to strip and stack doors, resp., in a cross-dock infrastructure. The goal is to minimise the transportation cost of the commodities to be handled at the cross-dock. Strip- and stack-related decomposition submodels are developed by taking benefit of the Integer Linearisation Property that appears in the new model. A linear search heuristic is also provided for obtaining feasible solutions by exploiting the special structure of the problem. We present an extensive computational study on a testbed of 55 instances to show that the proposed joint scheme for lower bounding and feasible solution providing is very efficient. The results obtained by our proposal are similar to those obtained by CPLEX in the 50 instances taken from the literature, but usually the time is much smaller for the large instances. The comparison with CPLEX is particularly good for the new five largest instances, which make it very promising for its application in real-life cases.


Introduction and motivation
Given a network with commodities to deliver from origin nodes to destination ones, the cross-docking problem basically consists of using a entity to serve as a consolidation point for helping on supply chain distribution.The origin nodes can deliver the commodities in the cross-dock, so that, after being classified by type and destination, they can be transported (usually, in smaller quantities) to the destination nodes.A cross-dock has a set of inbound infrastructures (say, doors, so-named strip ones, vehicles, usually, trucks, etc.), and a set of outbound ones (say, doors, so-named stack doors, vehicles, etc.), each of them with a capacity for commodity handling during a given time period.The goal consists of assigning the origin and destination nodes to the inbound and outbound infrastructures, resp., so that the commodity handling and transportation cost in the cross-docking is minimised.
As a matter of fact, industry detected the importance of cross-docking problem earlier than academia.It seems that Wal-Mart and some others in the retailer industry, see Stalk, Evans, and Shulman (1992), were pioneering its applicability in the late 80s, jointly with UPS, in particular, for its parcel delivering system, CONTACT Laureano F. Escudero laureano.escudero@urjc.esStatistics and Operations Research area, Universidad Rey Juan Carlos, URJC, 28933 Móstoles (Madrid), Spain Supplemental data for this article can be accessed online at https://doi.org/10.1080/00207543.2023.2180307.see Forger (1995).Office Depot is another company that took benefit of cross-docking in the earlier years, see Ross (1997).Cross-docking implementations have started to be frequent in fast-moving consumer-goods, manufacturing, automation and, in general, commodity distribution, say, from big companies to small ones; for specific implementations, see for example, Kinnear (1997), Witt (1998), andNapolitano (2011).Serrano, Delorme, and Dolgui (2021) present an application of cross-docking dynamic operations in a Renault's International Logistic Network.See also this type of cross-docking application in Coindreau et al. (2021) for another large car manufacturer.Both works present different mixed integer linear programming (MILP) models for minimising the commodity inventory penalisation when the truck origin node flow is not synchronised, in the same time period, with the container destination nodes.For comprehensive overviews on cross-docking applications, see Ladier andAlpan (2006), VanBelle, Valckenaers, andCattrysse (2012), Fathollahi-Fard et al. (2019), Gelareh et al. (2020), Kiani Mavi et al. (2020), Shahmardan and Sajadieh (2020), Castellucci, Costa, and Toledo (2021), Dulebenets (2021), and Akkerman et al. (2022).
The cross-docking problem solving by scientific means is a difficult task.It was not until the first decade of the current century when a strong advancement in combinatorial theory and algorithms has allowed to develop tools for problem solving.Birim (2106), Enderer, Contardo, and Contreras (2017), Nassief (2017), Nasiri et al. (2018), and Gunawan et al. (2021) integrate the cross-docking problem with the vehicle scheduling one from origin to destination such that the routes always begin and end at the doors of the cross-dock infrastructure.Simulated annealing and reduced variable neighbourhood search algorithms, among others, are considered for problem solving.Yu, Jewpanya, and Kachitvichyanukul (2016) integrate cross-docks in a multi-period product distribution network operations planning.Ye et al. (2018) present an approach for sequencing vehicles -in fact, trucks-to doors in order to minimise waiting time makespan; a particle swarm algorithm is considered jointly with a repair heuristic.Wisittipanich, Irohara, and Hengmeechai (2019), Shahmardan and Sajadieh (2020), and Castellucci, Costa, and Toledo (2021) propose MILP models for the truck scheduling in a multiple cross-docking facility network multi-period setting to minimise commodity storage cost and customer delivery time via a makespan minimisation in order to synchronising the inbound and outbound trucks in the network.The doors so-named the broader term 'boxes' in the cross-docks, are taken into account as a cluster for the related trucks assignment in Wisittipanich, Irohara, and Hengmeechai (2019).On the other hand, the inbound trucks can be partially unloaded and are also used as outbound ones in Shahmardan and Sajadieh (2020), where a hybrid heuristic-simulated annealing is introduced.Castellucci, Costa, and Toledo (2021) present a logic-based Benders Decomposition approach for problem solving.Gaudioso, Monaco, and Sammarra (2021) also deal with the same type of sequencing and scheduling cross-docking problem; a Lagrangean decomposition scheme is considered jointly with two repair heuristics embedded in the Lagrangean process.Dulebenets (2018aDulebenets ( , 2018bDulebenets ( , 2021) ) present MILP models where a single cross-dock facility is considered for inbound and outbound truck scheduling, and the trucks are individually assigned to doors.The aim is to minimise the overall cost of weighted truck handling, waiting time for inbound trucks and delaying departure cost for outbound trucks.Some types of evolutionary metaheuristics are introduced.Amini and Tavakkoli-Moghaddam (2016), Fathollahi-Fard et al. (2019), Keshtzari, Naderi, andMehdizadeh (2016), andLiao, Egbelu, andChang (2012), among others, present approaches for inbound and outbound truck scheduling in a single cross-dock facility, such that the doors are taken into account in some of those works as a cluster where no individual truck-door assignment is performed.The goal is to minimise the makespan.Some metaheuristic approaches for problem solving are considered (as evolutionary, social engineering optimiser, particle swarm algorithms, among others); Amini and Tavakkoli-Moghaddam (2016), additionally, consider trucks breakdown which uncertainty is modelled by considering the Normal distribution as well as commodity storage and outbound due date constraints.Tadumadze et al. (2019) propose a MILP model for integrating the inbound truck scheduling and the workforce to manage the commodity unloading time that is required in a multi-period environment, where the doors are individually considered for inbound truck assignment.The goal is to minimise the commodity unloading time.Theophilus et al. (2021) propose a MILP model for also a single cross-dock facility in a cold-chain assignment application for inbound and outbound truck scheduling, where the doors are individually considered for the related trucks assignment.Temperature-controlled storage areas are considered for perishable commodities.An evolutionary metaheuristic algorithm is introduced.Recently, Faghih-Mohammadi, Nasiri, and Konur (2022) present an approach for disaster relief operations along a 72 h period, as well as a comprehensive overview on the topic, including uncertainty on commodity inbound and demand.The MILP model that is introduced considers individual node-door assignment.A broad computational experiment is performed for validating the matheuristic approaches as fix-and-relax, local search and others.
Besides the works in Amini and Tavakkoli-Moghaddam (2016) and Faghih-Mohammadi, Nasiri, and Konur (2022) mentioned above, not many other works in the literature present schemes dealing with the uncertainty, probably due to the difficulty to be added to the problem.See in Soanpet (2012), Mousavi et al. (2014), andFaghih-Mohammadi, Nasiri, andKonur (2022) interesting overviews on the subject.
It can be observed that the optimisation literature in cross-docking has mostly focused on providing feasible solutions to the inbound and outbound trucks scheduling problem in one or several facilities, being single-period or multi-period oriented applications, usually without detailing the trucks assignment to the doors.So, the general aim is to address the problem solving of supply and distribution network management.Given the difficulty of the problem, several heuristic resolution strategies are considered.It remains as an open problem in many of those cross-docking models to estimate the solutions quality by considering approaches to either providing the optimal solution or obtaining strong lower bounds, at least, to validate the heuristic solutions.On the other hand, the literature for problem solving is scarce on the operations in a cross-dock facility, where the inbound and outbound commodity is known at a given time period.In this case, the aim is the assignment of origin nodes to strip doors and destination nodes to stack doors.It is clear that this proper operation has to be performed anyway, perhaps in a daily basis, once the big planning in single or multi terminals, single or multi-period crossdocking facility network is carried out.It is known as Cross-dock Door Assignment Problem (CDAP), where the aim is to perform such assignment at the minimum cost or, at least, providing a feasible assignment which quality is assessed by strong cost lower bounds in reasonable computing time, given the type of application that is involved.It is a very difficult binary quadratic problem, in fact, it is a NP-hard combinatorial optimisation (CO) problem.There are good works already published on CDAP in the open literature, see Guignard et al. (2012), Nassief, Contreras, andAs'ad (2016), Gelareh et al. (2020), andGuignard (2020), among others.See also comprehensive reviews in Agustina, Lee, andPipiani (2010), VanBelle, Valckenaers, andCattrysse (2012), Buijs, Vis, and Carlo (2014), Nassief (2017), and Goodarzi et al. (2020).Indeed, CDAP is the subject of the current work, since we believe that the solution quality can still be improved requiring small computing effort.It seems that the promising CO solution methodologies for CDAP could be of further interest in some respects.Some of them are related to providing insights and increasing a general understanding on the behaviour of larger cross-docking systems, besides CDAP being a benchmark object where to experiment new CO schemes.
The organisation of the rest of the work is as follows: Section 2 formally presents CDAP as a pure binary quadratic combinatorial optimisation one.It also presents the main contributions of this work.Section 3 deals with the models LIP and LIP-e, and introduces the Lagrangean Decomposition of the split variable reformulation of the second one, the strip-and stack-related submodels and the further tightening ones on one hand and the decoupling scheme on the other one.Section 4 presents the Local Search Heuristic (LSH) for obtaining feasible solutions.Section 5 reports the results of a broad computational experiment that has been conducted for performing a benchmark among the schemes presented in the work by considering large-scale instances.And Section 6 draws the main conclusions and outlines a future research agenda, mainly for dealing with CDAP under uncertainty.As a supplementary file, Appendix A details the proof of equivalence of the new formulations presented in Section 3; Appendix B details the results of the computational experimentation over a testbed of small and middle size instances taken from Guignard et al. (2012); Appendix C illustrates the procedure LSH by using a toy instance; and Appendix D gives the above equivalent formulations of CDAP by using that toy instance.

CDAP description and binary quadratic modelling
CDAP, as one of the key structures of cross-docking, can formally be stated as follows.Let a set of origin nodes M, a set of destination nodes N, a set of strip doors I, and a set of stack doors J.Each strip/stack door may serve more than one origin/destination node satisfying the doors' capacity constraints.On the other hand, each origin/destination must be allocated to one strip/stack door.The classical operation process of a cross-dock is depicted in Figure 1.
commodity handling and transportation unit cost (so-named distance, in the literature) from strip door i to stack door j, for i ∈ I, j ∈ J. w mn , commodity volume from origin node m to be delivered at destination node n, for m ∈ M, n ∈ N, so that d ij w mn is the commodity handling and transportation total cost from node i to node j using the strip door i and stack door j. s m and r n , commodity total volume arriving from origin node m, for m ∈ M, and delivering to node n, for n ∈ N, resp., such that s m = n∈N w mn and r n = m∈M w mn .S i and R j , commodity handling capacity of strip door i, for i ∈ I, and stack door j, for j ∈ J, resp.
The binary variables are as follows: x mi = 1, its value 1 means that origin m is assigned to strip door i; otherwise, 0.
y nj = 1, its value 1 means that destination n is assigned to stack door j; otherwise, 0.
The assignment of origin and destination nodes to strip and stack doors, resp., can be represented by Generalised Assignment Problems (GAPs), to be expressed The objective function that links the GAPs is a binary quadratic one on x mi y nj .So, model CDAP can be expressed The model belongs to the category so-named BQP (Binary Quadratic Problem) as a special mixed integer nonlinear problem.There are several alternatives for BQP solving, some of them have already been implemented in state-of-the-art solvers as CPLEX (the basic solver used in this work); it considers the so-named Fortet inequalities, see Fortet (1960).Other alternatives use reformulations as the so-named convexification of the quadratic objective function, see Billionnet, Elloumi, and Plateau (2009).A recent alternative is based on the Dantzig-Wolfe Reformulation, where the constraint subsets in the model are replaced with the convex hull of its extreme points, see Cesselli, Létocart, and Traversi (2022).
As mentioned above, CDAP (3) is a NP-hard and, then, difficult, combinatorial optimisation problem.Its treatment for problem solving is a young discipline, as a matter of fact, most of the literature has been published during the last 15 years, see Section 1. Cohen and Keren (2009) present one of the first models and a heuristic for cross-docking.Zhu et al. (2009) introduce the binary quadratic model CDAP (3), whose notation and modelling have become a reference in the crossdocking literature.To the best of our knowledge, Guignard et al. (2012) present the first MILP mathematically equivalent model to CDAP (3); for that purpose, the Reformulation Linearisation Technique (RLT1) introduced in Adams and Sherali (1986) has been used.The work outlines a LSH that we have taken benefit from, and a convex hull-based metaheuristic designed explicitly for quadratic 0-1 problems with linear constraints.The testbed of instances that have been generated has become a reference for validating proposals in the crossdocking literature.Nassief, Contreras, and As'ad (2016) introduce a modification of model CDAP by converting it in a mathematically equivalent binary linear one.For that purpose, a new structure is considered, where k is one of the, say, |K| commodities to be sent from origin, say m(k), to destination, say, n(k).A new binary variable z kij is introduced for the triplet (k, i, j), whose value 1 means commodity k enters the cross-dock through strip door i and exiting it through stack door j and otherwise, 0. The objective function (3a) and the constraints (1b) and (2b) are replaced in the new model with z-based ones, and redundant constraints based on (1a) and (2a) are appended for the triplet (z kij , x m(k),i , y n(k),j ).A Lagrangean Relaxation of those constraints results in a model to be decomposed in three types of independent submodels, one of them is further decomposed in |K| trivially solved semi-assignment problems, and the other two are further decomposed in |I| + |J| binary knapsack problems.The solution value at each Lagrangean iteration is a lower bound of the solution value of model CDAP.A lazy primal algorithm obtains a feasible solution from the Lagrangean one to be improved by a local search.The results are very good on the testbed in Guignard et al. (2012) as well as on new generated instances with larger dimensions.Nassief, Contreras, and Jaumard (2018) present two MILP models which are mathematically equivalent to model CDAP.One of such requires a column generation algorithm, so that new columns are appended to the associated so-named master problem relaxation.Guignard (2020) report some computational experiments on RTL1-based cross-docking.A very broad computational experiment have been reported in Gelareh et al. (2020) by considering alternatives to model CDAP, where different RLT1 models are compared which main differences lie on the appending of a variety of tight redundant cuts.

Main contributions of the work and outline
The solving approach for CDAP considered in this work takes benefit of the GAP structures (x, S)-( 1a

Decomposition based lower bounds to the solution value of an enlarged equivalent linearised integer reformulation
The binary quadratic model CDAP (3a) can be transformed through the introduction of new variables into the mathematically equivalent one so-named Linearised mixed Integer Programming Problem (LIP).For that purpose, let us replace the quadratic term x mi y nj with the continuous linear variable v minj .The RTL1 model LIP can be expressed Based on the computational results reported in Gelareh et al. (2020), model LIP (4a) is the most efficient reformulation among the alternatives that have been experimented with by a straightforward use of a state-of-the-art solver in a broad testbed.However, we will also show that the redundant constraints (5b) and (5c) can tighten model (4a), so that the resulting enlarged reformulation, so-named LIP-e (5), could be reformulated itself to consider a Lagrangean Decomposition (LD) scheme.

Lagrangean decomposition (LD)
Let us dualise the split-variable constraint where v minj is a copy of variable v minj in the split-variable reformulation of (5a), for m ∈ M, i ∈ I, n ∈ N, j ∈ J, and δ minj ∀m, i, n, j denote the unsigned Lagrangean multipliers vector of those constraints, to be updated by considering the time-honoured subgradient method, see Geoffrion (1974) and Guignard (2003).The resulting Lagrangean objective function can be expressed Thus, by exploiting the x-related constraints on one side and the y-related constraints on the other one in the splitvariable reformulation of model (5a), it results that the dual model can be decomposed in the submodels (8a) and (9a) to be expressed as follows: Strip-related submodel with additional stack-related constraints Stack-related submodel with additional strip-related constraints Notice that the solution value z LD (δ k+1 ) (7) of the δparametric model at a given Lagrangean iteration, say k (computed as the sum z strip LD (δ k+1 ) + z stack LD (δ k+1 )), is a lower bound of the solution value z * of the original model LIP (4a).
Additionally, both submodels (8a) and (9a) can be tightened by appending the constraints (10a) and (10b), resp.Those new constraints are obtained by taking into account the information implied by constraints (4c) and (4b).The rationale behind it is as follows: It can be seen in constraints (4c) that the sum in the left-hand side (lhs) have the same value for all the indices m, m ∈ M : m = m , independently of the 0-1 value of variable y nj , so, constraint (10a) is satisfied.In the same way, it can be seen in constraints (4b) that the lhs takes the same 0-1 value, independently of the 0-1 value of variable x mi .
Constraint system (10a) is redundant in case it is appended to the split-variable reformulation of (5a), without having any add-value.Moreover, its appending to the strip-and stack-related submodels (8a) and (9a) can have a tightening effect at the price of requiring more computing effort, see next.

Strip-and stack-related tightened models. New equivalent formulations of model LIP (4)
Let us consider the strip-related tightened submodel, for zero Lagrange multipliers, to be expressed as Notice that the solution value of the split-variable reformulation of (5a) is also given by zero-Lagrangean solution value z * LD (0).

Using the ILP to decompose submodels (8) and (9)
The ILP has proved to be very useful in the decoupling of subproblems into smaller ones, see Geoffrion (1974) and Guignard (2003).Note that the strip-related submodel ( 8a) is in the v-and x-binary variables, so that for any pair (m, i) it happens that if x mi = 0 then, constraints (8b) and (8d) force v minj = 0 ∀n ∈ N, j ∈ J. On the other hand, if x mi = 1, then those constraints define a feasible region merely dependent upon pair (m, i).So, the value of the objective function related to pair (m, i) takes the value 0 for x mi = 0 and a value, say β mi , when x mi = 1.
Then, an x mi -parametric subproblem for each pair (m, i) based on submodel (8a) can be expressed as β mi = min n∈N j∈J ( 1 2 d ij w mn + δ minj ) v minj , subject to the corresponding constraints (8b) and (8d), which solution value is 0 for x mi = 0. Thus, the optimal value of submodel (8a) can be found by solving the following subproblem z strip LD (δ) = min m∈M i∈I β mi x mi , subject to the GAP constraints (x, S)-(1a).Note that, after computing all the β mi values, one per each pair (m, i), submodel (8a) reduces to a GAP over the x-variables, which obviously is smaller in size and simpler in structure than the original one.Additionally, the subproblems to solve for obtaining the β-values are also GAPs, thus, they are as small and simple as the other one.A similar decomposition scheme can be applied to solve the stack-related submodel (9a).

On GAP-related sequential decoupling of the strip-and stack-submodels-based lower bound
By considering the decoupling scheme of submodels (8a) and (9a) as outlined above, this section presents the scheme for solving the related GAPs, at a given Lagrangean iteration.
Decoupling strip-related submodel (8a) Step 1 The solution value β mi of the N × J-binary variables in the (m, i)-GAP can be expressed Step 2 Given the set of optimal values {β mi ∀m ∈ M, i ∈ I} retrieved from solving the models (13a) in Step 1, the solution of submodel (8a) can be equivalently obtained by the M × I-binary variables-related GAP to be expressed as Thus, the solution of submodel (8a) requires to solve M × I + 1 GAPs.

Decoupling stack-related submodel (9a)
A similar procedure can be carried out as for decoupling the strip-related submodel (8a).
Step 1 The solution value α nj of the M × I-binary variables in the (n, j)-GAP can be expressed Step 2 Given the set of optimal values {α nj ∀n ∈ N, j ∈ J} retrieved from solving the models (15a), the solution of submodel (9a) can be equivalently obtained by the N × J-binary variables-related GAP to be expressed as Thus, the solution of submodel (9a) requires to solve N × J + 1 GAPs.As summary, and after executing the two-step scheme, a bound, say z, of the optimal solution value z * of model CDAP (3a), obtained at a given Lagrangean iteration, can be expressed as z = z strip LD (δ) + z stack LD (δ).

The LSH for obtaining feasible solutions
Given the difficulty of BQP, a metaheuristic algorithm is considered in this section for obtaining feasible solutions for CDAP also taking benefit of the special structure of its two GAP constraint systems (x, S)-( 1a) and (y, R)-( 2a).The algorithm LSH is an improvement of the algorithm versions, so-named LS1 and LS2, presented in Guignard et al. (2012).It is proposed to split the original CDAP model (3a) into two simpler submodels with linear objective functions, smaller dimensions and the same constraints and variables.
In particular, the quadratic objective function (3a) can be viewed either as a function of the y-variables or as a function of the x-ones, to be expressed as Thus, for a given solution, say, ŷ ≡ (ŷ nj ∀n ∈ N, j ∈ J) of variables vector y, the solution value z(x, ŷ) can be retrieved from the model that can be expressed In a similar way, for a given solution, say, x ≡ (x mi ∀m ∈ M, i ∈ J) of variables vector x, the solution value z(x, y) can be retrieved from the model that can be expressed In case that the vectors x and ŷ are feasible for their respective GAP models ( 18) and ( 19), then a feasible solution is obtained for CDAP model (3a); its objective function value is min{z(x, ŷ), z(x, y)}.Note that in an iterative procedure the value ŷ of variables vector y can be generated at random at the first trial, but the x-solution of model ( 18) is still a feasible one, so that at the next trial the y-solution as retrieved from model ( 19) is also feasible.
Continuing the procedure on that way, a local optimal solution of model (3a) can be obtained through several trials until the solution value z(x, ŷ) cannot be improved any more, based on the initial solution ŷ.At each iteration, say k, a new initial solution ŷ is generated at random.The procedure may continue until a given bound, say k, on the number of iterations is reached.Let us notate the incumbent values of the variables vectors x and y of the procedure as x and y, resp., such that the related cost can be expressed Note 1: The incumbent cost z can be useful as an upper bound in the B&C phase of the straightforward use of the state-of-the-art MILP optimiser for solving some models.Note 2: A similar procedure can be considered by taking x as the initial x-vector, being generated at random.Note: The estimation of t z can be expressed as #g • max{t a , t b }, where t a and t b are the highest computing time required by any iteration of a ≡Steps 0-4 and b ≡Steps 5-9 of LHS, #th denotes de number of available threads,and #g denotes the number of groups (where, sequentially, the parallel execution of #th iterations of algorithm LSH is to be performed, by executing either Steps 0-4 or Steps 5-9) such that #g processors ES-2670 (20 threads), 2.50 Ghz and 128Gb RAM.The CPLEX v12.8 Concert Technology library has been used, embedded in a C/C++ experimental code, for solving the decomposition submodels as well as straightforward solving the full models.
For testing purposes a broad computational experiment has been performed by considering 55 instances in two testbeds.Testbed 1 is composed of a set of 35 benchmark instances in the literature, in fact proposed by Guignard et al. (2012).Testbed 2 is composed of a set of 20 much larger and difficult instances, where the first 15 ones are also taken from Guignard et al. (2012), and the other instances are generated for the experiment to validate the proposal in larger instances.In all of them |M| = |N| and |I| = |J|.The results obtained for the smaller ones, i.e. testbed 1, are presented in Appendix B. The rest of the section shows the results that have been obtained for testbed 2.

LSH performance
As it is shown in Appendix B, see Table B.1, heuristic LSH provides the optimal solution for testbed 1 and nondifficult instances in testbed 2, improving on the results obtained in Guignard et al. (2012).
On the other hand, Table 1 shows the LSH's results for the difficult instances in testbed 2. The headings are as follows: xxxyySzz, instance code, where xx is the number of origin nodes and the number of destination nodes, yy is the number of strip doors and number of stack doors, and zz is the capacity slackness (in percentage); z, incumbent solution value; k, number of LSH iterations that are considered; #reruns, the number of reruns of the algorithm; z, average of the incumbent solution values obtained in the reruns.t z and tt z , estimation of the total computing time required by the reruns in case of parallelisation and total computing time without parallelisation, resp.; GR LSH , goodness ratio of the LSH incumbent solution versus the one provided by straightforward CPLEX, expressed as z z *

CPX
, where z * CPX is the value provided by CPLEX, see also Table 3; and LS1, LS2, CH, incumbent solution values taken from Guignard et al. (2012) as obtained by the related heuristics.
It can be observed in the table that the ratio GR LSH is always slightly smaller than 1, which means that the solution quality is practically the same as the CPLEX one, even when considering the 2h time limit that has been allowed for the latter, see also Subsection 5.2.
Moreover, heuristic LSH provides a slighty smaller incumbent solution value z than any of the heuristics LS1, LS2 or CH for all the difficult instances in Guignard et al. (2012), see the values in bold.
Table 2 shows the main results for assessing the LSH's solution quality by considering the lower bound z retrieved from solving the strip-and stack-related submodels (8a) and (9a), resp., at the Lagrangean iteration where the multiplier updatings have been stabilised.Hence, the new headings are as follows: z LD (0), lower bound provided by the strip-and stack-related submodels (8a) and (9a), with Lagrangean vector δ = 0, by decoupling them into the submodels (13a) and (14a) for the former, and (15a) and (16a) for the latter; t z LD (0) and tt z LD (0) , estimation of the required time (s), for the first Lagrangean iteration in case of parallelisation, and total computing time (s) without parallelisation, resp., for obtaining lower bound z LD (0); z, lower bound provided by the summation of the solution values of the submodels (14a) and (16a) (recall that it is tighter than z LD (0)); t z and tt z , estimation of the required time (s), in case of parallelisation, and total computing time (s) without parallelisation, resp., for obtaining z; and GAP LSH %, optimality gap of the LSH incumbent solution z, expressed as 100 • z−z z .It is worth to point out that z is obtained as the highest summation of the optimal solution values of the decoupled submodel (14a) (that uses as data the solution of the submodels (13a)) and the decoupled  submodel (16a) (that uses as data the solution of the submodels (15a)), among the iterations.
Observe that even the total computing time (s) without parallelisation for obtaining the LSH incumbent solution value z and the lower bound z LD (0), say, tt z + tt z LD (0) is much smaller than the CPLEX time t * CPX−(4) , for the 10 largest instances, where tt z and tt z LD (0) are reported in Tables 1 and 2, resp.
Additionally, we can also notice that the LD approach based on the split-variable reformulation of model (5a) strengths the zero-Lagrangean-based bound (i.e.z is higher than z LD (0)).However, it is not worth the effort (represented in the alternative computing times t z and tt z versus t z LD (0) and tt z LD (0) ), given the tightness of z LD (0) with respect to the solution value z * of the original model.

Performance of straightforward CPLEX
For the original binary quadratic model CDAP (3), the dimensions in the experiment by considering testbed 2 are from 180 up to 500 0-1 variables, from 180 to 500 constraints, and from 8100 up to 62500 quadratic terms.
Table 3  .Note that the higher GR CPX > 1 and the higher TR CPX > 1, the worse performance of CPLEX versus LSH.
Observe in the table that CPLEX obtains the optimality in 50% of the instances (i.e.GAP CPX % = 0).On the other hand, given GR CPX = 1, notice that LSH provides the same solution value with a smaller time (i.e.TR > 1).However, for the largest instances, GAP CPX is reasonable (as GAP LSH is, see Table 2), and GR CPX is close to 1, so, both approaches have practically the same solution quality.Moreover, TR is much higher than 1, i.e. the time required by CPLEX is much higher than the one that does it by LSH.

Conclusions and future research agenda
Several schemes for obtaining lower bounds of the optimal solution value of CDAP, a difficult binary quadratic combinatorial problem, have been introduced.Those schemes take benefit of the GAP structures that are embedded in the model, in particular, in its seemingly most efficient linearised integer programming reformulation LIP (4).
Additionally, a heuristic algorithm has been provided for obtaining good feasible solutions in the two testbeds composed of 55 instances to consider for validating the proposal made in this work.It has been computationally proved that the overall scheme for providing the heuristic incumbent solution and strong lower bounds is very efficient.Those bounds are obtained from an efficient decomposition of model LIP-e (5), an enlarging of model LIP (4) that exploits its GAP structures.The proposed scheme obtains and proves the optimal solution in testbed 1, see Appendix B, composed of 35 instances in the experiment; those instances are a reference taken from the literature.Straightforward CPLEX also does the same for those instances, and requiring smaller computing time.The proposal obtains and proves an optimal solution in 10 out of the 20 instances included in testbed 2, the 5 largest ones are generated for testing it in larger instances.CPLEX also obtains and proves the optimal solution for those instances but requiring much higher computing time in most of them.
Neither the proposal, nor CPLEX can prove the optimality of the solution in the 10 largest instances in testbed 2; anyway, the optimality gaps are reasonable in both systems, although the proposal provides slightly better incumbent solution values, and it requires much less computing time than CPLEX.
Based on the computational experiment performed in this work, a provisional conclusion would be that the subgradient method-related LD lower bound of the splitvariable reformulation of the LIP-e (5a) of the original model CDAP is stronger than the zero-Lagrangean one.However, it is not worth the effort, given the computing time required by the former and the tightness of the latter.
Finally, two new equivalent MILP models of CDAP are introduced.Both of them are obtained by tightening the corresponding strip-and stack-related submodels.It is also proved that any of these new formulations allow to obtain the solution value of the original model CDAP in reasonable computing time for the instances in testbed 1.

Main items in the future research agenda
• The Simplicial Decomposition and Proximal Level methodology types would be investigated looking for even more efficient decomposition approaches for solving the original model CDAP.• A further tightening of the decomposed submodels ( 8) and ( 9).• Cross-dock Door Design Planning under uncertainty, where a two-stage distributionally robust optimisation binary mixed quadratic approach can be considered, where the number and capacity of the strip and stack doors are first stage decisions.The uncertainty could lie on the sets of origin and destination nodes in a network and the door capacity disruption; it would be represented in a finite set of scenarios for the ambiguity set.

Disclosure statement
No potential conflict of interest was reported by the author(s).
Lagrangean decomposition, by dualising the split-variable constraints in the split-variable reformulation of model LIP-e, is introduced in order to obtain lower bounds in reasonable computing time.That model is decomposed in two submodels, such that one is related to the strip doors, to be named strip submodel, and the other one is related to the stack doors, to be named stack submodel.
Guignard et al. (2012)ose many tested ones, in the testbed of instances inGuignard et al. (2012).An enlarged redundant model is introduced in our work, let us name it model LIP-e.(ii)A min Proof: It is trivial by swapping the roles of M for N, I for J and x mi for y nj in the proof of Proposition 3.1.Corollary 3.3: An alternative formulation of model LIP (4a) is the one that results as the decomposition in the submodels (8a) plus (10a) and (9a) plus (10b), with zero Lagrange multipliers, resp.In this case, the solution value coincides with z * LD (0).Proof: It is followed from Propositions 3.1 and 3.2.

Table 1 .
LSH incumbent solution value for difficult instances in testbed 2.

Table 2 .
LSH solution quality for the instances in testbed 2.The estimation of t z LD (0) and t z can be expressed as #it • #g strip • t mod−(13) + t mod−(14) + #g stack • t mod−(15) + t mod−(16) , where #it is the number of Lagrangean iterations, #th denotes the number of threads in the available HW, #g strip and #g stack denote the number of groups (where, sequentially, the parallel execution of #th submodels (13a) and (14a) can be performed, resp., for each group) such that #g strip = |M|•|I|

Table 3 .
CPLEX performance for models LIP for the instances in testbed 2.
shows the main results when solving model LIP (4) by straightforward CPLEX for this set of instances.The headings are as follows: z CPX , lower bound obtained by CPLEX; z * CPX , and t * CPX−(4) , incumbent solution value and required computing time (s), resp.; z LP−(4) and t LP−(4) , related LP lower bound and required computing time (s); and symbol means that the LP solution value as well as the related time have the same value for the set of instances with the same parameters, but the flow S-slackness percentage; GAP CPX %, optimality gap of the incumbent solution values z * CPX , expressed as 100 •